Script UQ basique

diff --git a/ReacteurNumerique_UQ/CasJouet_cath_OT.py b/ReacteurNumerique_UQ/CasJouet_cath_OT.py
new file mode 100644 (file)
index 0000000..48007cc
--- /dev/null
+++ b/ReacteurNumerique_UQ/CasJouet_cath_OT.py
@@ -0,0 +1,97 @@
+#module load openturns/1.14.0 (à charger avant)
+import sys, os
+import openturns as ot
+#import matplotlib.pyplot as plt
+
+sys.path.insert(0, '/data/projets/projets.002/rn_gt_incertitudes.281/D - Divers/cathWrap2')
+from cathWrap import Cathare_Wrapper
+
+orginal_dir = os.getcwd()
+
+def _exec(P_pa,T_C):
+
+    inputs = {'Pstart':P_pa,
+              'Pstart_I':P_pa,
+              'Tstart':T_C}
+    
+    c = Cathare_Wrapper(jdd = '/data/projets/projets.002/rn_gt_incertitudes.281/D - Divers/C3_CANONFORMAT-132bar.dat',
+                        in_dict = inputs,
+                        main_dir = '/home/i46979/Documents/Persalys/Test_Cathare',
+                        #name='mod_jdd',
+                        cath_link='/home/i46979/Documents/Cathare_jdd/INSTALLATION_C3/init_c3.sh')
+
+    c.main()
+
+    time, alpha = c.res['TOPA']
+
+    #i = next(i for i,v in enumerate(alpha) if v > 0.8) #commenté au 12/04/2022
+    #t_s = time[i]
+    #affichage graphique à résoudre?
+    #plt.figure()
+    #plt.title('titre?')
+    
+    t_s = time[-1]
+    #alpha_end = alpha[-1]
+
+    return t_s
+
+if __name__ == '__main__':
+    
+    #Définition du vecteur d'entrée à partir des arguments en ligne de commande
+    #P = 13 200 000 Pa & T = 300°C (rappel)
+    #E = [float(sys.argv[1]),float(sys.argv[2])]
+
+    #Evaluation avec openturns
+    #Définition de la loi pour la pression
+    m_P = float(sys.argv[1])
+    loi_P = ot.Uniform(0.99*m_P,1.01*m_P)
+    #Définition de la loi pour la température
+    m_T = float(sys.argv[2])
+    loi_T = ot.Uniform(0.99*m_T,1.01*m_T)
+    #Définition de la distribution d'entrée
+    dependance = ot.IndependentCopula(2)
+    loi_E = ot.ComposedDistribution([loi_P,loi_T],dependance)
+    #Création du vecteur d'entrée aléatoire
+    random_loi_E = ot.RandomVector(loi_E)
+
+    #Conversion de la fonction _exec vers une fonction openturns
+    def fct_conv(X):
+        return [_exec(X[0],X[1])]
+    #Définition de la fonction python à partir de la fonction convertie
+    fct_cath = ot.PythonFunction(2, 1, fct_conv)
+
+    #Définition du gradient de fct_cath
+    pasGradient = [10, 0.05]
+    def fct_cath_grad(X):
+        y = ot.NonCenterdFinitedDifferenceGradient(pasGradient,fct_cath.getEvaluation())
+        y_grad = y.gradient
+        return y_grad([X[0],X[1]])
+
+    #Définition de la fonction python définitive à manipuler
+    fct_ot_def = ot.PythonFunction(2, 1, fct_cath, gradient = fct_cath_grad)
+
+    #Définition du hessien de la fonction python définitive (optionnelle)
+    pasHessien = [1, 0.005]
+    h = ot.CenteredFiniteDifferenceHessian(pasHessien,fct_ot_def.getEvaluation())
+    fct_ot_def.setHessian(h)
+
+    #Calcul de la sortie aléatoire
+    fct_ot_def = ot.MemoizeFunction(fct_ot_def)#facultatif, pour consultation historique
+    random_S = ot.CompositeRandomVector(fct_ot_def,random_loi_E)
+
+    #Obtention d'une sortie échantillonnée
+    #ech_S = random_S.getSample(10)
+    #mc_1_S = ech_S.computeMean()
+    #mc_2_S = ech_S.computeStandardDeviation()
+    
+    #Obtention d'une sortie de type Taylor/Cumul quadratique
+    taylor_S = ot.TaylorExpansionMoments(random_S)
+    taylor_1_S = taylor_S.getMeanFirstOrder()
+    taylor_2_S = taylor_S.getCovariance()
+    
+    #Obtention d'une sortie dépassement de seuil avec Monte-Carlo
+    #non retenue dans les fonctionnalités prioritaires
+
+    #print("tirages aléatoires de S = ", ech_S)
+    print("moyenne = ",mc_1_S," ou", taylor_1_S)
+    print("variance = ",mc_2_S," ou", taylor_2_S)