SALOME platform Git repositories - tools/solverlab.git/blob

   1 #!/usr/bin/env python3\r
   2 # -*-coding:utf-8 -*\r
   3 \r
   4 """\r
   5 Created on Mon Aug 30 2021\r
   6 @author: Michael NDJINGA, Katia Ait Ameur, Coraline Mounier\r
   7 \r
   8 Euler system without source term in one dimension on regular domain [a,b]\r
   9 Riemann problemn with ghost cell boundary condition\r
  10 Left and right: Neumann boundary condition \r
  11 Roe scheme\r
  12 Regular square mesh\r
  13 \r
  14 State law Stiffened gaz : p = (gamma - 1) * rho * (e - q) - gamma * p0\r
  15 4 choices of parameters gamma and p0 are available : \r
  16   - Lagrange interpolation (with q=0)\r
  17   - Hermite interpolation with reference point at 575K (with q=0)\r
  18   - Hermite interpolation with reference point at 590K (with q=0)\r
  19   - Hermite interpolation with reference point at 617.94K (saturation at 155 bar)  with q=0\r
  20   \r
  21 Linearized enthalpy : h = h_sat + cp * (T - T_sat)\r
  22 Values for cp and T_sat parameters are taken at the reference point chosen for the state law\r
  23 \r
  24 To do correct the computation of the time step : lambda_max (maximum eigenvalue) should be computed first)\r
  25 """\r
  26 \r
  27 import cdmath\r
  28 import numpy as np\r
  29 import matplotlib\r
  30 \r
  31 matplotlib.use("Agg")\r
  32 import matplotlib.pyplot as plt\r
  33 import matplotlib.animation as manimation\r
  34 import sys\r
  35 from math import sqrt, atan, pi\r
  36 from numpy import sign\r
  37 \r
  38 \r
  39 ## Numerical parameter\r
  40 precision = 1e-5\r
  41 \r
  42 #state law parameter : can be 'Lagrange', 'Hermite590K', 'Hermite617K', or 'FLICA'\r
  43 state_law = "Hermite590K"\r
  44 \r
  45 #indicates with test case is simulated to compare with FLICA5 results\r
  46 #num_test = 0 means there are no FLICA5 results given here\r
  47 num_test = 0\r
  48 \r
  49 #def state_law_parameters(state_law):\r
  50 #state law Stiffened Gaz : p = (gamma - 1) * rho * e - gamma * p0\r
  51 global gamma\r
  52 global p0\r
  53 global q\r
  54 global c0\r
  55 global cp\r
  56 global h_sat\r
  57 global T_sat\r
  58 \r
  59 if state_law == "Lagrange":\r
  60         # reference values for Lagrange interpolation\r
  61         p_ref = 155. * 10**5     #Reference pressure in a REP 900 nuclear power plant     \r
  62         p1    = 153. * 10**5     # value of pressure at inlet of a 900 MWe PWR vessel\r
  63         rho_ref = 594.38        #density of water at saturation temperature of 617.94K and 155 bars\r
  64         rho1 = 742.36           # value of density at inlet of a 900 MWe PWR vessel (T1 = 565K)\r
  65         e_ref = 1603.8 * 10**3  #internal energy of water at saturation temperature of 617.94K and 155 bars\r
  66         e1 = 1273.6 * 10**3     # value of internal energy at inlet of a 900 MWe PWR vessel\r
  67         \r
  68         gamma = (p1 - p_ref) / (rho1 * e1 - rho_ref *e_ref) + 1.\r
  69         p0 = - 1. / gamma * ( - (gamma - 1) * rho_ref * e_ref + p_ref)\r
  70         q=0.\r
  71         c0 = sqrt((gamma - 1) * (e_ref + p_ref / rho_ref))\r
  72         \r
  73         cp = 8950.\r
  74         h_sat = 1.63 * 10 ** 6     # saturation enthalpy of water at 155 bars\r
  75         T_sat = 617.94 \r
  76 \r
  77 elif state_law == "Hermite617K":\r
  78         # reference values for Hermite interpolation\r
  79         p_ref = 155. * 10**5     #Reference pressure in a REP 900 nuclear power plant\r
  80         T_ref = 617.94          #Reference temperature for interpolation at 617.94K\r
  81         rho_ref = 594.38        #density of water at saturation temperature of 617.94K and 155 bars\r
  82         e_ref = 1603.8 * 10**3  #internal energy of water at saturation temperature of 617.94K and 155 bars\r
  83         h_ref   = e_ref + p_ref / rho_ref\r
  84         c_ref = 621.43          #sound speed for water at 155 bars and 617.94K\r
  85 \r
  86         gamma = 1. + c_ref * c_ref / (e_ref + p_ref / rho_ref)           # From the sound speed formula\r
  87         p0 = 1. / gamma * ( (gamma - 1) * rho_ref * e_ref - p_ref)       \r
  88         q=0.\r
  89         c0 = sqrt((gamma - 1) * (e_ref + p_ref / rho_ref))\r
  90         \r
  91         cp = 8950.                  # value at 155 bar and 617.94K\r
  92         h_sat = 1.63 * 10 ** 6     # saturation enthalpy of water at 155 bars\r
  93         T_sat = 617.94  \r
  94 \r
  95 elif state_law == 'Hermite590K':\r
  96         # reference values for Hermite interpolation\r
  97         p_ref = 155. * 10**5     #Reference pressure  in a REP 900 nuclear power plant\r
  98         T_ref = 590.             #Reference temperature for interpolation at 590K\r
  99         rho_ref = 688.3         #density of water at 590K and 155 bars\r
 100         e_ref = 1411.4 * 10**3  #internal energy of water at 590K and 155 bars\r
 101         h_ref   = e_ref + p_ref / rho_ref\r
 102         c_ref = 866.29          #sound speed for water at 155 bars and 590K\r
 103         \r
 104         gamma = 1. + c_ref * c_ref / (e_ref + p_ref / rho_ref)           # From the sound speed formula\r
 105         p0 = 1. / gamma * ( (gamma - 1) * rho_ref * e_ref - p_ref)       \r
 106         q=0.\r
 107         c0 = sqrt((gamma - 1) * (e_ref + p_ref / rho_ref))\r
 108         \r
 109         cp = 5996.8                  # value at 155 bar and 590K\r
 110         h_sat = 1433.9 * 10 ** 3     # saturation enthalpy of water at 155 bars\r
 111         T_sat = 590.  \r
 112 \r
 113 elif state_law == 'Hermite575K':\r
 114         # reference values for Hermite interpolation\r
 115         p_ref = 155 * 10**5     #Reference pressure  in a REP 900 nuclear power plant\r
 116         T_ref = 575             #Reference temperature at inlet in a REP 900 nuclear power plant\r
 117         #Remaining values determined using iapws python package\r
 118         rho_ref = 722.66        #density of water at 575K and 155 bars\r
 119         e_ref = 1326552.66  #internal energy of water at 575K and 155 bars\r
 120         h_ref   = e_ref + p_ref / rho_ref\r
 121         c_ref = 959.28          #sound speed for water at 155 bars and 575K\r
 122         \r
 123         gamma = 1 + c_ref * c_ref / (e_ref + p_ref / rho_ref)           # From the sound speed formula\r
 124         p0 = 1 / gamma * ( (gamma - 1) * rho_ref * e_ref - p_ref)       \r
 125         q=0.\r
 126         c0 = sqrt((gamma - 1) * (e_ref + p_ref / rho_ref))#This is actually c_ref\r
 127         \r
 128         cp = 5504.05                 # value at 155 bar and 590K\r
 129         h_sat = h_ref     # saturation enthalpy of water at 155 bars\r
 130         T_sat = T_ref\r
 131 else:\r
 132         raise ValueError("Incorrect value for parameter state_law")\r
 133 \r
 134 \r
 135 #initial parameters for Riemann problem (p in Pa, v in m/s, T in K)\r
 136 p_L = 155. * 10**5  \r
 137 p_R = 150. * 10**5\r
 138 v_L = 0.\r
 139 v_R = 0.\r
 140 h_L = 1.4963*10**6\r
 141 h_R = 1.4963*10**6\r
 142 \r
 143 T_L = (h_L - h_sat ) / cp + T_sat\r
 144 T_R = (h_R - h_sat ) / cp + T_sat\r
 145 \r
 146 def initial_conditions_Riemann_problem(a, b, nx):\r
 147         print("Initial data Riemann problem")\r
 148         dx = (b - a) / nx  # space step\r
 149         x = [a + 0.5 * dx + i * dx for i in range(nx)]  # array of cell center (1D mesh)        \r
 150 \r
 151         p_initial = np.array([ (xi < (a + b) / 2) * p_L + (xi >= (a + b) / 2) * p_R for xi in x])\r
 152         v_initial = np.array([ (xi < (a + b) / 2) * v_L + (xi >= (a + b) / 2) * v_R for xi in x])\r
 153         T_initial = np.array([ (xi < (a + b) / 2) * T_L + (xi >= (a + b) / 2) * T_R for xi in x])\r
 154 \r
 155         rho_initial = p_to_rho_StiffenedGaz(p_initial, T_initial)\r
 156         q_initial = rho_initial * v_initial\r
 157         rho_E_initial = T_to_rhoE_StiffenedGaz(T_initial, rho_initial, q_initial)\r
 158 \r
 159         return rho_initial, q_initial, rho_E_initial, p_initial, v_initial, T_initial\r
 160 \r
 161 def rho_to_p_StiffenedGaz(rho_field, q_field, rho_E_field):\r
 162         p_field = (gamma - 1) * ( rho_E_field - 1. / 2 * q_field ** 2 / rho_field - rho_field * q) - gamma * p0\r
 163         return p_field\r
 164         \r
 165 \r
 166 def p_to_rho_StiffenedGaz(p_field, T_field):\r
 167         rho_field = (p_field + p0) * gamma / (gamma - 1) * 1 / (h_sat + cp * (T_field - T_sat) - q)\r
 168         return rho_field\r
 169         \r
 170 \r
 171 def T_to_rhoE_StiffenedGaz(T_field, rho_field, q_field):\r
 172         rho_E_field = 1 / 2 * (q_field) ** 2 / rho_field + p0 + rho_field / gamma * (h_sat + cp * (T_field- T_sat) + (gamma - 1) * q)\r
 173         return rho_E_field\r
 174 \r
 175                 \r
 176 def rhoE_to_T_StiffenedGaz(rho_field, q_field, rho_E_field):\r
 177         T_field = T_sat + 1 / cp * (gamma * (rho_E_field / rho_field - 1 / 2 * (q_field / rho_field) ** 2) - gamma * p0 / rho_field -  (gamma - 1) * q - h_sat)\r
 178         return T_field\r
 179 \r
 180 def dp_drho_e_const_StiffenedGaz( e ):\r
 181         return (gamma-1)*(e-q)\r
 182 \r
 183 def dp_de_rho_const_StiffenedGaz( rho ):\r
 184         return (gamma-1)*rho\r
 185 \r
 186 def sound_speed_StiffenedGaz( h ):\r
 187         return np.sqrt((gamma-1)*(h-q))\r
 188 \r
 189 def rho_h_to_p_StiffenedGaz( rho, h ):\r
 190         return (gamma - 1) * rho * ( h - q ) / gamma - p0\r
 191 \r
 192 def MatRoe_StiffenedGaz( rho_l, q_l, rho_E_l, rho_r, q_r, rho_E_r):\r
 193         RoeMat = cdmath.Matrix(3, 3)\r
 194 \r
 195         u_l = q_l / rho_l\r
 196         u_r = q_r / rho_r\r
 197         p_l = rho_to_p_StiffenedGaz(rho_l, q_l, rho_E_l)\r
 198         p_r = rho_to_p_StiffenedGaz(rho_r, q_r, rho_E_r)\r
 199         H_l = rho_E_l / rho_l + p_l / rho_l\r
 200         H_r = rho_E_r / rho_r + p_r / rho_r\r
 201 \r
 202         # Roe averages\r
 203         rho = np.sqrt(rho_l * rho_r)\r
 204         u   = (u_l * sqrt(rho_l) + u_r * sqrt(rho_r)) / (sqrt(rho_l) + sqrt(rho_r))\r
 205         H   = (H_l * sqrt(rho_l) + H_r * sqrt(rho_r)) / (sqrt(rho_l) + sqrt(rho_r))\r
 206 \r
 207         p = rho_h_to_p_StiffenedGaz( rho, H - u**2/2. )\r
 208         e = H - p / rho - 1./2 * u**2\r
 209         dp_drho = dp_drho_e_const_StiffenedGaz( e )\r
 210         dp_de   = dp_de_rho_const_StiffenedGaz( rho )\r
 211 \r
 212         RoeMat[0, 0] = 0\r
 213         RoeMat[0, 1] = 1\r
 214         RoeMat[0, 2] = 0\r
 215         RoeMat[1, 0] = dp_drho - u ** 2 + dp_de / rho * (u**2/2 - e)\r
 216         RoeMat[1, 1] = 2 * u - u * dp_de / rho\r
 217         RoeMat[1, 2] = dp_de / rho\r
 218         RoeMat[2, 0] = -u * ( -dp_drho + H - dp_de / rho * (u**2/2 - e) )\r
 219         RoeMat[2, 1] = H - dp_de / rho * u ** 2\r
 220         RoeMat[2, 2] = (dp_de / rho +1) * u\r
 221         \r
 222         return(RoeMat)\r
 223 \r
 224         \r
 225 def Droe_StiffenedGaz( rho_l, q_l, rho_E_l, rho_r, q_r, rho_E_r):\r
 226     Droe   = cdmath.Matrix(3, 3)\r
 227 \r
 228     u_l = q_l / rho_l\r
 229     u_r = q_r / rho_r\r
 230     p_l = rho_to_p_StiffenedGaz(rho_l, q_l, rho_E_l)\r
 231     p_r = rho_to_p_StiffenedGaz(rho_r, q_r, rho_E_r)\r
 232     H_l = rho_E_l / rho_l + p_l / rho_l\r
 233     H_r = rho_E_r / rho_r + p_r / rho_r\r
 234 \r
 235     # Roe averages\r
 236     rho = np.sqrt(rho_l * rho_r)\r
 237     u = (u_l * sqrt(rho_l) + u_r * sqrt(rho_r)) / (sqrt(rho_l) + sqrt(rho_r))\r
 238     H = (H_l * sqrt(rho_l) + H_r * sqrt(rho_r)) / (sqrt(rho_l) + sqrt(rho_r))\r
 239 \r
 240     c = sound_speed_StiffenedGaz( H - u**2/2. )\r
 241     \r
 242     lamb  = cdmath.Vector(3)\r
 243     lamb[0] = u-c\r
 244     lamb[1] = u\r
 245     lamb[2] = u+c    \r
 246 \r
 247     r   = cdmath.Matrix(3, 3)\r
 248     r[0,0] = 1.\r
 249     r[1,0] = u-c\r
 250     r[2,0] = H-u*c    \r
 251     r[0,1] = 1.\r
 252     r[1,1] = u   \r
 253     r[2,1] = H-c**2/(gamma-1)    \r
 254     r[0,2] = 1.\r
 255     r[1,2] = u+c\r
 256     r[2,2] = H+u*c         \r
 257 \r
 258     l   = cdmath.Matrix(3, 3)\r
 259     l[0,0] = (1./(2*c**2))*(0.5*(gamma-1)*u**2+u*c)\r
 260     l[1,0] = (1./(2*c**2))*(-u*(gamma-1)-c)\r
 261     l[2,0] = (1./(2*c**2))*(gamma-1)\r
 262     l[0,1] = ((gamma-1)/c**2)*(H-u**2)\r
 263     l[1,1] = ((gamma-1)/c**2)*u   \r
 264     l[2,1] = -((gamma-1)/c**2)    \r
 265     l[0,2] = (1./(2*c**2))*(0.5*(gamma-1)*u**2-u*c)\r
 266     l[1,2] = (1./(2*c**2))*(c-u*(gamma-1))\r
 267     l[2,2] = (1./(2*c**2))*(gamma-1)\r
 268 \r
 269     M1 = cdmath.Matrix(3, 3) #abs(lamb[0])*r[:,0].tensProduct(l[:,0])\r
 270     M2 = cdmath.Matrix(3, 3) #abs(lamb[1])*r[:,1].tensProduct(l[:,1])   \r
 271     M3 = cdmath.Matrix(3, 3) #abs(lamb[2])*r[:,2].tensProduct(l[:,2])\r
 272     for i in range(3):\r
 273         for j in range(3):\r
 274             M1[i,j] = abs(lamb[0])*r[i,0]*l[j,0]\r
 275             M2[i,j] = abs(lamb[1])*r[i,1]*l[j,1]            \r
 276             M3[i,j] = abs(lamb[2])*r[i,2]*l[j,2]\r
 277             \r
 278     Droe = M1+M2+M3 \r
 279     \r
 280     return(Droe)    \r
 281 \r
 282 \r
 283 def jacobianMatricesm(coeff, rho_l, q_l, rho_E_l, rho_r, q_r, rho_E_r):\r
 284 \r
 285         if rho_l < 0 or rho_r < 0:\r
 286                 print("rho_l=", rho_l, " rho_r= ", rho_r)\r
 287                 raise ValueError("Negative density")\r
 288         if rho_E_l < 0 or rho_E_r < 0:\r
 289                 print("rho_E_l=", rho_E_l, " rho_E_r= ", rho_E_r)\r
 290                 raise ValueError("Negative total energy")\r
 291 \r
 292         RoeMat = MatRoe_StiffenedGaz( rho_l, q_l, rho_E_l, rho_r, q_r, rho_E_r)\r
 293         \r
 294         Droe = Droe_StiffenedGaz( rho_l, q_l, rho_E_l, rho_r, q_r, rho_E_r)    \r
 295         return (RoeMat - Droe) * coeff * 0.5\r
 296 \r
 297 \r
 298 def jacobianMatricesp(coeff, rho_l, q_l, rho_E_l, rho_r, q_r, rho_E_r):\r
 299         if rho_l < 0 or rho_r < 0:\r
 300                 print("rho_l=", rho_l, " rho_r= ", rho_r)\r
 301                 raise ValueError("Negative density")\r
 302         if rho_E_l < 0 or rho_E_r < 0:\r
 303                 print("rho_E_l=", rho_E_l, " rho_E_r= ", rho_E_r)\r
 304                 raise ValueError("Negative total energy")\r
 305 \r
 306         RoeMat = MatRoe_StiffenedGaz( rho_l, q_l, rho_E_l, rho_r, q_r, rho_E_r)\r
 307         Droe = Droe_StiffenedGaz( rho_l, q_l, rho_E_l, rho_r, q_r, rho_E_r)    \r
 308 \r
 309         return (RoeMat + Droe) * coeff * 0.5\r
 310 \r
 311 \r
 312 def FillEdges(j, Uk, nbComp, divMat, Rhs, Un, dt, dx, isImplicit):\r
 313         dUi = cdmath.Vector(3)\r
 314 \r
 315         if (j == 0):\r
 316                 rho_l   = Uk[j * nbComp + 0]\r
 317                 q_l     = Uk[j * nbComp + 1]\r
 318                 rho_E_l = Uk[j * nbComp + 2]\r
 319                 rho_r   = Uk[(j + 1) * nbComp + 0]\r
 320                 q_r     = Uk[(j + 1) * nbComp + 1]\r
 321                 rho_E_r = Uk[(j + 1) * nbComp + 2]\r
 322                 \r
 323                 Am = jacobianMatricesm(dt / dx, rho_l, q_l, rho_E_l, rho_r, q_r, rho_E_r)\r
 324                 divMat.addValue(j * nbComp, (j + 1) * nbComp, Am)\r
 325                 divMat.addValue(j * nbComp,       j * nbComp, Am * (-1.))\r
 326 \r
 327                 dUi[0] = Uk[(j + 1) * nbComp + 0] - Uk[j * nbComp + 0]\r
 328                 dUi[1] = Uk[(j + 1) * nbComp + 1] - Uk[j * nbComp + 1]\r
 329                 dUi[2] = Uk[(j + 1) * nbComp + 2] - Uk[j * nbComp + 2]\r
 330                 temp = Am * dUi\r
 331 \r
 332                 if(isImplicit):\r
 333                         Rhs[j * nbComp + 0] = -temp[0] - (Uk[j * nbComp + 0] - Un[j * nbComp + 0])\r
 334                         Rhs[j * nbComp + 1] = -temp[1] - (Uk[j * nbComp + 1] - Un[j * nbComp + 1])\r
 335                         Rhs[j * nbComp + 2] = -temp[2] - (Uk[j * nbComp + 2] - Un[j * nbComp + 2])\r
 336 \r
 337         elif (j == nx - 1):\r
 338                 rho_l   = Uk[(j - 1) * nbComp + 0]\r
 339                 q_l     = Uk[(j - 1) * nbComp + 1]\r
 340                 rho_E_l = Uk[(j - 1) * nbComp + 2]\r
 341                 rho_r   = Uk[j * nbComp + 0]\r
 342                 q_r     = Uk[j * nbComp + 1]\r
 343                 rho_E_r = Uk[j * nbComp + 2]\r
 344 \r
 345                 Ap = jacobianMatricesp(dt / dx, rho_l, q_l, rho_E_l, rho_r, q_r, rho_E_r)\r
 346                 divMat.addValue(j * nbComp, j * nbComp, Ap)\r
 347                 divMat.addValue(j * nbComp, (j - 1) * nbComp, Ap * (-1.))\r
 348 \r
 349                 dUi[0] = Uk[(j - 1) * nbComp + 0] - Uk[j * nbComp + 0]\r
 350                 dUi[1] = Uk[(j - 1) * nbComp + 1] - Uk[j * nbComp + 1]\r
 351                 dUi[2] = Uk[(j - 1) * nbComp + 2] - Uk[j * nbComp + 2]\r
 352 \r
 353                 temp = Ap * dUi\r
 354 \r
 355                 if(isImplicit):\r
 356                         Rhs[j * nbComp + 0] = temp[0] - (Uk[j * nbComp + 0] - Un[j * nbComp + 0])\r
 357                         Rhs[j * nbComp + 1] = temp[1] - (Uk[j * nbComp + 1] - Un[j * nbComp + 1])\r
 358                         Rhs[j * nbComp + 2] = temp[2] - (Uk[j * nbComp + 2] - Un[j * nbComp + 2])\r
 359 \r
 360 def FillInnerCell(j, Uk, nbComp, divMat, Rhs, Un, dt, dx, isImplicit):\r
 361 \r
 362         rho_l   = Uk[(j - 1) * nbComp + 0]\r
 363         q_l     = Uk[(j - 1) * nbComp + 1]\r
 364         rho_E_l = Uk[(j - 1) * nbComp + 2]\r
 365         rho_r   = Uk[j * nbComp + 0]\r
 366         q_r     = Uk[j * nbComp + 1]\r
 367         rho_E_r = Uk[j * nbComp + 2]\r
 368         Ap = jacobianMatricesp(dt / dx, rho_l, q_l, rho_E_l, rho_r, q_r, rho_E_r)\r
 369         \r
 370         rho_l   = Uk[j * nbComp + 0]\r
 371         q_l     = Uk[j * nbComp + 1]\r
 372         rho_E_l = Uk[j * nbComp + 2]\r
 373         rho_r   = Uk[(j + 1) * nbComp + 0]\r
 374         q_r     = Uk[(j + 1) * nbComp + 1]\r
 375         rho_E_r = Uk[(j + 1) * nbComp + 2]\r
 376         Am = jacobianMatricesm(dt / dx, rho_l, q_l, rho_E_l, rho_r, q_r, rho_E_r)\r
 377 \r
 378         divMat.addValue(j * nbComp, (j + 1) * nbComp, Am)\r
 379         divMat.addValue(j * nbComp, j * nbComp, Am * (-1.))\r
 380         divMat.addValue(j * nbComp, j * nbComp, Ap)\r
 381         divMat.addValue(j * nbComp, (j - 1) * nbComp, Ap * (-1.))\r
 382         dUi1 = cdmath.Vector(3)\r
 383         dUi2 = cdmath.Vector(3)\r
 384         dUi1[0] = Uk[(j + 1) * nbComp + 0] - Uk[j * nbComp + 0]\r
 385         dUi1[1] = Uk[(j + 1) * nbComp + 1] - Uk[j * nbComp + 1]\r
 386         dUi1[2] = Uk[(j + 1) * nbComp + 2] - Uk[j * nbComp + 2]\r
 387 \r
 388         dUi2[0] = Uk[(j - 1) * nbComp + 0] - Uk[j * nbComp + 0]\r
 389         dUi2[1] = Uk[(j - 1) * nbComp + 1] - Uk[j * nbComp + 1]\r
 390         dUi2[2] = Uk[(j - 1) * nbComp + 2] - Uk[j * nbComp + 2]\r
 391         temp1 = Am * dUi1\r
 392         temp2 = Ap * dUi2\r
 393 \r
 394         if(isImplicit):\r
 395                 Rhs[j * nbComp + 0] = -temp1[0] + temp2[0] - (Uk[j * nbComp + 0] - Un[j * nbComp + 0])\r
 396                 Rhs[j * nbComp + 1] = -temp1[1] + temp2[1] - (Uk[j * nbComp + 1] - Un[j * nbComp + 1])\r
 397                 Rhs[j * nbComp + 2] = -temp1[2] + temp2[2] - (Uk[j * nbComp + 2] - Un[j * nbComp + 2])\r
 398 \r
 399 def SetPicture(rho_field_Roe, q_field_Roe, h_field_Roe, p_field_Roe, v_field_Roe, T_field_Roe, dx):\r
 400         max_initial_rho = max(rho_field_Roe)\r
 401         min_initial_rho = min(rho_field_Roe)\r
 402         max_initial_q = max(q_field_Roe)\r
 403         min_initial_q = min(q_field_Roe)\r
 404         min_initial_h = min(h_field_Roe)\r
 405         max_initial_h = max(h_field_Roe)\r
 406         max_initial_p = max(p_field_Roe)\r
 407         min_initial_p = min(p_field_Roe)\r
 408         max_initial_v = max(v_field_Roe)\r
 409         min_initial_v = min(v_field_Roe)\r
 410         max_initial_T = max(T_field_Roe)\r
 411         min_initial_T = min(T_field_Roe)\r
 412 \r
 413         fig, ([axDensity, axMomentum, axRhoE],[axPressure, axVitesse, axTemperature]) = plt.subplots(2, 3,sharex=True, figsize=(14,10))\r
 414         plt.gcf().subplots_adjust(wspace = 0.5)\r
 415 \r
 416         lineDensity_Roe, = axDensity.plot([a+0.5*dx + i*dx for i in range(nx)], rho_field_Roe, label='Roe')\r
 417         axDensity.set(xlabel='x (m)', ylabel='Densité (kg/m3)')\r
 418         axDensity.set_xlim(a,b)\r
 419         #axDensity.set_ylim(min_initial_rho - 0.1 * (max_initial_rho - min_initial_rho), 700.)\r
 420         axDensity.set_ylim(657, 660.5)\r
 421         axDensity.legend()\r
 422 \r
 423         lineMomentum_Roe, = axMomentum.plot([a+0.5*dx + i*dx for i in range(nx)], q_field_Roe, label='Roe')\r
 424         axMomentum.set(xlabel='x (m)', ylabel='Momentum (kg/(m2.s))')\r
 425         axMomentum.set_xlim(a,b)\r
 426         #axMomentum.set_ylim(min_initial_q - 0.1*(max_initial_q-min_initial_q), max_initial_q * 2.5)\r
 427         axMomentum.set_ylim(-50, 500)\r
 428         axMomentum.legend()\r
 429         \r
 430         lineRhoE_Roe, = axRhoE.plot([a+0.5*dx + i*dx for i in range(nx)], h_field_Roe, label='Roe')\r
 431         axRhoE.set(xlabel='x (m)', ylabel='h (J/m3)')\r
 432         axRhoE.set_xlim(a,b)\r
 433         #axRhoE.set_ylim(min_initial_h - 0.05*abs(min_initial_h), max_initial_h +  0.5*(max_initial_h-min_initial_h))\r
 434         axRhoE.set_ylim(1.495 * 10**6, 1.5*10**6)\r
 435         axRhoE.legend()\r
 436         \r
 437         linePressure_Roe, = axPressure.plot([a+0.5*dx + i*dx for i in range(nx)], p_field_Roe, label='Roe')\r
 438         axPressure.set(xlabel='x (m)', ylabel='Pression (bar)')\r
 439         axPressure.set_xlim(a,b)\r
 440         #axPressure.set_ylim(min_initial_p - 0.05*abs(min_initial_p), max_initial_p +  0.5*(max_initial_p-min_initial_p))\r
 441         axPressure.set_ylim(149, 156)\r
 442         axPressure.ticklabel_format(axis='y', style='sci', scilimits=(0,0))\r
 443         axPressure.legend()\r
 444 \r
 445         lineVitesse_Roe, = axVitesse.plot([a+0.5*dx + i*dx for i in range(nx)], v_field_Roe, label='Roe')\r
 446         axVitesse.set(xlabel='x (m)', ylabel='Vitesse (m/s)')\r
 447         axVitesse.set_xlim(a,b)\r
 448         #axVitesse.set_ylim(min_initial_v - 0.05*abs(min_initial_v), 15)\r
 449         axVitesse.set_ylim(-0.5, 1)\r
 450         axVitesse.ticklabel_format(axis='y', style='sci', scilimits=(0,0))\r
 451         axVitesse.legend()\r
 452 \r
 453         lineTemperature_Roe, = axTemperature.plot([a+0.5*dx + i*dx for i in range(nx)], T_field_Roe, label='Roe')\r
 454         axTemperature.set(xlabel='x (m)', ylabel='Température (K)')\r
 455         axTemperature.set_xlim(a,b)\r
 456         #axTemperature.set_ylim(min_initial_T - 0.005*abs(min_initial_T), 604)\r
 457         axTemperature.set_ylim(600, 601)\r
 458         axTemperature.ticklabel_format(axis='y', style='sci', scilimits=(0,0))\r
 459         axTemperature.legend()\r
 460         \r
 461         return(fig, lineDensity_Roe, lineMomentum_Roe, lineRhoE_Roe, linePressure_Roe, lineVitesse_Roe, lineTemperature_Roe, lineDensity_Roe, lineMomentum_Roe, lineRhoE_Roe, linePressure_Roe, lineVitesse_Roe, lineTemperature_Roe)\r
 462 \r
 463 \r
 464 def EulerSystemRoe(ntmax, tmax, cfl, a, b, nbCells, output_freq, meshName, state_law):\r
 465         #state_law_parameters(state_law)\r
 466         dim = 1\r
 467         nbComp = 3\r
 468         dt = 0.\r
 469         time = 0.\r
 470         it = 0\r
 471         isStationary = False\r
 472         isImplicit = False\r
 473         dx = (b - a) / nx\r
 474         dt = cfl * dx / c0\r
 475         #dt = 5*10**(-6)\r
 476         nbVoisinsMax = 2\r
 477 \r
 478         # iteration vectors\r
 479         Un_Roe  = cdmath.Vector(nbCells * (nbComp))\r
 480         dUn_Roe = cdmath.Vector(nbCells * (nbComp))\r
 481         dUk_Roe = cdmath.Vector(nbCells * (nbComp))\r
 482         Rhs_Roe = cdmath.Vector(nbCells * (nbComp))\r
 483 \r
 484         # Initial conditions\r
 485         print("Construction of the initial condition …")\r
 486 \r
 487         rho_field_Roe, q_field_Roe, rho_E_field_Roe, p_field_Roe, v_field_Roe, T_field_Roe = initial_conditions_Riemann_problem(a, b, nx)\r
 488         h_field_Roe = (rho_E_field_Roe + p_field_Roe) / rho_field_Roe - 0.5 * (q_field_Roe / rho_field_Roe) **2\r
 489         p_field_Roe = p_field_Roe * 10 ** (-5)\r
 490         \r
 491 \r
 492         for k in range(nbCells):\r
 493                 Un_Roe[k * nbComp + 0] = rho_field_Roe[k]\r
 494                 Un_Roe[k * nbComp + 1] = q_field_Roe[k]\r
 495                 Un_Roe[k * nbComp + 2] = rho_E_field_Roe[k]\r
 496 \r
 497         divMat_Roe = cdmath.SparseMatrixPetsc(nbCells * nbComp, nbCells * nbComp, (nbVoisinsMax + 1) * nbComp)\r
 498         \r
 499         # Picture settings\r
 500         fig, lineDensity, lineMomentum, lineRhoE, linePressure, lineVitesse, lineTemperature, lineDensity_Roe, lineMomentum_Roe, lineRhoE_Roe, linePressure_Roe, lineVitesse_Roe, lineTemperature_Roe  = SetPicture( rho_field_Roe, q_field_Roe, h_field_Roe, p_field_Roe, v_field_Roe, T_field_Roe, dx)\r
 501 \r
 502         # Video settings\r
 503         FFMpegWriter = manimation.writers['ffmpeg']\r
 504         metadata = dict(title="Roe scheme for the 1D Euler system", artist="CEA Saclay", comment="Shock tube")\r
 505         writer = FFMpegWriter(fps=10, metadata=metadata, codec='h264')\r
 506         with writer.saving(fig, "1DEuler_complet_RP_Roe" + ".mp4", ntmax):\r
 507                 writer.grab_frame()\r
 508                 plt.savefig("EulerComplet_RP_" + str(dim) + "D_Roe" + meshName + "0" + ".png")\r
 509                 np.savetxt( "EulerComplet_RP_" + str(dim) + "D_Roe" + meshName + "_rho" + "0" + ".txt", rho_field_Roe, delimiter="\n")\r
 510                 np.savetxt( "EulerComplet_RP_" + str(dim) + "D_Roe" + meshName + "_q" + "0" + ".txt", q_field_Roe,  delimiter="\n")\r
 511                 iterGMRESMax = 50\r
 512                 newton_max = 100\r
 513         \r
 514                 print("Starting computation of the non linear Euler non isentropic system with Roe scheme …")\r
 515                 # STARTING TIME LOOP\r
 516                 while (it < ntmax and time <= tmax and not isStationary):\r
 517                         dUn_Roe = Un_Roe.deepCopy()\r
 518                         Uk_Roe  = Un_Roe.deepCopy()\r
 519                         residu_Roe = 1.\r
 520                         \r
 521                         k_Roe = 0\r
 522                         while (k_Roe < newton_max and residu_Roe > precision):\r
 523                                 # STARTING NEWTON LOOP\r
 524                                 divMat_Roe.zeroEntries()  #sets the matrix coefficients to zero\r
 525                                 for j in range(nbCells):\r
 526                                         \r
 527                                         # traitements des bords\r
 528                                         if (j == 0):\r
 529                                                 FillEdges(j, Uk_Roe, nbComp, divMat_Roe, Rhs_Roe, Un_Roe, dt, dx, isImplicit)\r
 530                                         elif (j == nbCells - 1):\r
 531                                                 FillEdges(j, Uk_Roe, nbComp, divMat_Roe, Rhs_Roe, Un_Roe, dt, dx, isImplicit)\r
 532         \r
 533                                         # traitement des cellules internes\r
 534                                         else:\r
 535                                                 FillInnerCell(j, Uk_Roe, nbComp, divMat_Roe, Rhs_Roe, Un_Roe, dt, dx, isImplicit)\r
 536                                                 \r
 537                                 #solving the linear system divMat * dUk = Rhs\r
 538                                 \r
 539                                 if(isImplicit):\r
 540                                         divMat_Roe.diagonalShift(1.)\r
 541                                         LS_Roe = cdmath.LinearSolver(divMat_Roe, Rhs_Roe, iterGMRESMax, precision, "GMRES", "LU")\r
 542                                         dUk_Roe = LS_Roe.solve()\r
 543                                         residu_Roe = dUk_Roe.norm()\r
 544                                 else:\r
 545                                         dUk_Roe=Rhs_Roe - divMat_Roe*Un_Roe\r
 546                                         residu_Roe = 0.#Convergence schéma Newton\r
 547                                 \r
 548                                 if (isImplicit and (it == 1 or it % output_freq == 0 or it >= ntmax or isStationary or time >= tmax)):\r
 549                                         print("Residu Newton :   ", residu_Roe)\r
 550                                         print("Linear system converged in ", LS_Roe.getNumberOfIter(), " GMRES iterations")\r
 551         \r
 552                                 #updates for Newton loop\r
 553                                 Uk_Roe += dUk_Roe\r
 554                                 k_Roe = k_Roe + 1\r
 555                                 if (isImplicit and not LS_Roe.getStatus()):\r
 556                                         print("Linear system did not converge ", LS_Roe.getNumberOfIter(), " GMRES iterations")\r
 557                                         raise ValueError("No convergence of the linear system")\r
 558                                 \r
 559                                 if k_Roe == newton_max:\r
 560                                         raise ValueError("No convergence of Newton Roe Scheme")\r
 561         \r
 562                         #updating fields\r
 563                         Un_Roe = Uk_Roe.deepCopy()\r
 564                         dUn_Roe -= Un_Roe\r
 565                         if (dUn_Roe.norm()<precision):\r
 566                                         isStationary = True\r
 567                         \r
 568                         for k in range(nbCells):\r
 569                                 rho_field_Roe[k]   = Un_Roe[k * nbComp + 0]\r
 570                                 q_field_Roe[k]     = Un_Roe[k * nbComp + 1]\r
 571                                 rho_E_field_Roe[k] = Un_Roe[k * nbComp + 2]\r
 572         \r
 573                         v_field_Roe = q_field_Roe / rho_field_Roe\r
 574                         p_field_Roe = rho_to_p_StiffenedGaz(rho_field_Roe, q_field_Roe, rho_E_field_Roe)\r
 575                         T_field_Roe = rhoE_to_T_StiffenedGaz(rho_field_Roe, q_field_Roe, rho_E_field_Roe)\r
 576                         h_field_Roe = (rho_E_field_Roe + p_field_Roe) / rho_field_Roe - 0.5 * (q_field_Roe / rho_field_Roe) **2\r
 577                         p_field_Roe = p_field_Roe * 10 ** (-5)\r
 578                         \r
 579                         if( min(p_field_Roe)<0) :\r
 580                                 raise ValueError("Negative pressure, stopping calculation")\r
 581         \r
 582                         #picture and video updates\r
 583                         lineDensity_Roe.set_ydata(rho_field_Roe)\r
 584                         lineMomentum_Roe.set_ydata(q_field_Roe)\r
 585                         lineRhoE_Roe.set_ydata(h_field_Roe)\r
 586                         linePressure_Roe.set_ydata(p_field_Roe)\r
 587                         lineVitesse_Roe.set_ydata(v_field_Roe)\r
 588                         lineTemperature_Roe.set_ydata(T_field_Roe)\r
 589                         \r
 590                         writer.grab_frame()\r
 591         \r
 592                         time = time + dt\r
 593                         it = it + 1\r
 594         \r
 595                         # Savings\r
 596                         if (it == 1 or it % output_freq == 0 or it >= ntmax or isStationary or time >= tmax):\r
 597                 \r
 598                                 print("-- Time step : " + str(it) + ", Time: " + str(time) + ", dt: " + str(dt))\r
 599         \r
 600                                 plt.savefig("EulerComplet_RP_" + str(dim) + "D_Roe" + meshName + str(it) + '_time' + str(time) + ".png")\r
 601 \r
 602         print("-- Iter: " + str(it) + ", Time: " + str(time) + ", dt: " + str(dt))\r
 603         if (it >= ntmax):\r
 604                 print("Maximum number of time steps ntmax= ", ntmax, " reached")\r
 605                 return\r
 606 \r
 607         elif (isStationary):\r
 608                 print("Stationary regime reached at time step ", it, ", t= ", time)\r
 609                 print("------------------------------------------------------------------------------------")\r
 610                 np.savetxt("EulerComplet_RP_" + str(dim) + "D_Roe" + meshName + "_rho_Stat.txt", rho_field_Roe, delimiter="\n")\r
 611                 np.savetxt("EulerComplet_RP_" + str(dim) + "D_Roe" + meshName + "_q_Stat.txt", q_field_Roe, delimiter="\n")\r
 612                 np.savetxt("EulerComplet_RP_" + str(dim) + "D_Roe" + meshName + "_rhoE_Stat.txt", rho_E_field_Roe, delimiter="\n")\r
 613                 np.savetxt("EulerComplet_RP_" + str(dim) + "D_Roe" + meshName + "_p_Stat.txt", p_field_Roe, delimiter="\n")\r
 614                 plt.savefig("EulerComplet_RP_" + str(dim) + "D_Roe" + meshName + "_Stat.png")\r
 615                 return\r
 616         else:\r
 617                 print("Maximum time Tmax= ", tmax, " reached")\r
 618                 return\r
 619 \r
 620 \r
 621 def solve(a, b, nx, meshName, meshType, cfl, state_law):\r
 622         print("Resolution of the Euler System in dimension 1 on " + str(nx) + " cells")\r
 623         print("State Law Stiffened Gaz, " + state_law)\r
 624         print("Initial data : ", "Riemann problem")\r
 625         print("Boundary conditions : ", "Neumann")\r
 626         print("Mesh name : ", meshName, ", ", nx, " cells")\r
 627         # Problem data\r
 628         tmax = 10.\r
 629         ntmax = 25\r
 630         output_freq = 1\r
 631         EulerSystemRoe(ntmax, tmax, cfl, a, b, nx, output_freq, meshName, state_law)\r
 632         return\r
 633 \r
 634 \r
 635 if __name__ == """__main__""":\r
 636         a = 0.\r
 637         b = 1.\r
 638         nx = 50\r
 639         cfl = 0.5\r
 640         #state_law = "Hermite590K"\r
 641         #state_law_parameters(state_law)\r
 642         solve(a, b, nx, "RegularSquares", "", cfl, state_law)\r