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 with heating source term (phi) in one dimension on regular domain [a,b]\r
   9 Riemann problemn with ghost cell boundary condition\r
  10 Left : Inlet boundary condition (velocity and temperature imposed)\r
  11 Right : Outlet boundary condition (pressure imposed)\r
  12 Roe scheme\r
  13 Regular square mesh\r
  14 \r
  15 State law Stiffened gaz : p = (gamma - 1) * rho * (e - q) - gamma * p0\r
  16 4 choices of parameters gamma and p0 are available : \r
  17   - Lagrange interpolation (with q=0)\r
  18   - Hermite interpolation with reference point at 575K (with q=0)\r
  19   - Hermite interpolation with reference point at 590K (with q=0)\r
  20   - Hermite interpolation with reference point at 617.94K (saturation at 155 bar)  with q=0\r
  21   \r
  22 Linearized enthalpy : h = h_sat + cp * (T - T_sat)\r
  23 Values for cp and T_sat parameters are taken at the reference point chosen for the state law\r
  24 \r
  25 To do correct the computation of the time step : lambda_max (maximum eigenvalue) should be computed first)\r
  26 """\r
  27 \r
  28 import cdmath\r
  29 import numpy as np\r
  30 import matplotlib\r
  31 \r
  32 matplotlib.use("Agg")\r
  33 import matplotlib.pyplot as plt\r
  34 import matplotlib.animation as manimation\r
  35 import sys\r
  36 from math import sqrt, atan, pi\r
  37 from numpy import sign\r
  38 \r
  39 \r
  40 #### Initial and boundary condition (T in K, v in m/s, p in Pa)\r
  41 T_inlet  = 565.\r
  42 v_inlet  = 5.\r
  43 p_outlet = 155.0 * 10**5\r
  44 \r
  45 #initial parameters are determined from boundary conditions\r
  46 p_0   = p_outlet       #initial pressure\r
  47 v_0   = v_inlet        #initial velocity\r
  48 T_0   = T_inlet        #initial temperature\r
  49 ### Heating source term\r
  50 phi=1.e8\r
  51 \r
  52 ## Numerical parameter\r
  53 precision = 1e-6\r
  54 \r
  55 #state law parameter : can be 'Lagrange', 'Hermite590K', 'Hermite617K', or 'FLICA'\r
  56 state_law = "Hermite575K"\r
  57 \r
  58 def state_law_parameters(state_law):\r
  59         #state law Stiffened Gaz : p = (gamma - 1) * rho * e - gamma * p0\r
  60         global gamma\r
  61         global p0\r
  62         global q\r
  63         global c0\r
  64         global cp\r
  65         global h_sat\r
  66         global T_sat\r
  67         \r
  68         if state_law == "Lagrange":\r
  69                 # reference values for Lagrange interpolation\r
  70                 p_ref = 155. * 10**5     #Reference pressure in a REP 900 nuclear power plant     \r
  71                 p1    = 153. * 10**5     # value of pressure at inlet of a 900 MWe PWR vessel\r
  72                 rho_ref = 594.38        #density of water at saturation temperature of 617.94K and 155 bars\r
  73                 rho1 = 742.36           # value of density at inlet of a 900 MWe PWR vessel (T1 = 565K)\r
  74                 e_ref = 1603.8 * 10**3  #internal energy of water at saturation temperature of 617.94K and 155 bars\r
  75                 e1 = 1273.6 * 10**3     # value of internal energy at inlet of a 900 MWe PWR vessel\r
  76                 \r
  77                 gamma = (p1 - p_ref) / (rho1 * e1 - rho_ref *e_ref) + 1.\r
  78                 p0 = - 1. / gamma * ( - (gamma - 1) * rho_ref * e_ref + p_ref)\r
  79                 q=0.\r
  80                 c0 = sqrt((gamma - 1) * (e_ref + p_ref / rho_ref))\r
  81                 \r
  82                 cp = 8950.\r
  83                 h_sat = 1.63 * 10 ** 6     # saturation enthalpy of water at 155 bars\r
  84                 T_sat = 617.94 \r
  85 \r
  86         elif state_law == "Hermite617K":\r
  87                 # reference values for Hermite interpolation\r
  88                 p_ref = 155. * 10**5     #Reference pressure in a REP 900 nuclear power plant\r
  89                 T_ref = 617.94          #Reference temperature for interpolation at 617.94K\r
  90                 rho_ref = 594.38        #density of water at saturation temperature of 617.94K and 155 bars\r
  91                 e_ref = 1603.8 * 10**3  #internal energy of water at saturation temperature of 617.94K and 155 bars\r
  92                 h_ref   = e_ref + p_ref / rho_ref\r
  93                 c_ref = 621.43          #sound speed for water at 155 bars and 617.94K\r
  94 \r
  95                 gamma = 1. + c_ref * c_ref / (e_ref + p_ref / rho_ref)           # From the sound speed formula\r
  96                 p0 = 1. / gamma * ( (gamma - 1) * rho_ref * e_ref - p_ref)       \r
  97                 q=0.\r
  98                 c0 = sqrt((gamma - 1) * (e_ref + p_ref / rho_ref))\r
  99                 \r
 100                 cp = 8950.                  # value at 155 bar and 617.94K\r
 101                 h_sat = 1.63 * 10 ** 6     # saturation enthalpy of water at 155 bars\r
 102                 T_sat = 617.94  \r
 103         \r
 104         elif state_law == 'Hermite590K':\r
 105                 # reference values for Hermite interpolation\r
 106                 p_ref = 155. * 10**5     #Reference pressure  in a REP 900 nuclear power plant\r
 107                 T_ref = 590.             #Reference temperature for interpolation at 590K\r
 108                 rho_ref = 688.3         #density of water at 590K and 155 bars\r
 109                 e_ref = 1411.4 * 10**3  #internal energy of water at 590K and 155 bars\r
 110                 h_ref   = e_ref + p_ref / rho_ref\r
 111                 c_ref = 866.29          #sound speed for water at 155 bars and 590K\r
 112                 \r
 113                 gamma = 1. + c_ref * c_ref / (e_ref + p_ref / rho_ref)           # From the sound speed formula\r
 114                 p0 = 1. / gamma * ( (gamma - 1) * rho_ref * e_ref - p_ref)       \r
 115                 q=0.\r
 116                 c0 = sqrt((gamma - 1) * (e_ref + p_ref / rho_ref))\r
 117                 \r
 118                 cp = 5996.8                  # value at 155 bar and 590K\r
 119                 h_sat = 1433.9 * 10 ** 3     # saturation enthalpy of water at 155 bars\r
 120                 T_sat = 590.  \r
 121 \r
 122         elif state_law == 'Hermite575K':\r
 123                 # reference values for Hermite interpolation\r
 124                 p_ref = 155 * 10**5     #Reference pressure  in a REP 900 nuclear power plant\r
 125                 T_ref = 575             #Reference temperature at inlet in a REP 900 nuclear power plant\r
 126                 #Remaining values determined using iapws python package\r
 127                 rho_ref = 722.66        #density of water at 575K and 155 bars\r
 128                 e_ref = 1326552.66  #internal energy of water at 575K and 155 bars\r
 129                 h_ref   = e_ref + p_ref / rho_ref\r
 130                 c_ref = 959.28          #sound speed for water at 155 bars and 575K\r
 131                 \r
 132                 gamma = 1 + c_ref * c_ref / (e_ref + p_ref / rho_ref)           # From the sound speed formula\r
 133                 p0 = 1 / gamma * ( (gamma - 1) * rho_ref * e_ref - p_ref)       \r
 134                 q=0.\r
 135                 c0 = sqrt((gamma - 1) * (e_ref + p_ref / rho_ref))#This is actually c_ref\r
 136                 \r
 137                 cp = 5504.05                 # value at 155 bar and 590K\r
 138                 h_sat = h_ref     # saturation enthalpy of water at 155 bars\r
 139                 T_sat = T_ref\r
 140         else:\r
 141                 raise ValueError("Incorrect value for parameter state_law")\r
 142                 \r
 143 def initial_conditions_Riemann_problem(a, b, nx):\r
 144         print("Initial data Riemann problem")\r
 145         dx = (b - a) / nx  # space step\r
 146         x = [a + 0.5 * dx + i * dx for i in range(nx)]  # array of cell center (1D mesh)\r
 147 \r
 148         p_initial = np.array([ p_0 for xi in x])\r
 149         v_initial = np.array([ v_0 for xi in x])\r
 150         T_initial = np.array([ T_0 for xi in x])\r
 151 \r
 152         rho_initial = p_to_rho_StiffenedGaz(p_initial, T_initial)\r
 153         q_initial = rho_initial * v_initial\r
 154         rho_E_initial = T_to_rhoE_StiffenedGaz(T_initial, rho_initial, q_initial)\r
 155 \r
 156         return rho_initial, q_initial, rho_E_initial, p_initial, v_initial, T_initial\r
 157 \r
 158 def p_to_rho_StiffenedGaz(p_field, T_field):\r
 159         rho_field = (p_field + p0) * gamma / (gamma - 1) * 1. / (h_sat + cp * (T_field - T_sat))\r
 160         return rho_field\r
 161         \r
 162 def T_to_rhoE_StiffenedGaz(T_field, rho_field, q_field):\r
 163         rho_E_field = 1. / 2. * (q_field) ** 2 / rho_field + p0 + rho_field / gamma * (h_sat + cp * (T_field- T_sat))\r
 164         return rho_E_field\r
 165 \r
 166 def rhoE_to_T_StiffenedGaz(rho_field, q_field, rho_E_field):\r
 167         T_field = T_sat + 1 / cp * (gamma * (rho_E_field / rho_field - 1 / 2 * (q_field / rho_field) ** 2) - gamma * p0 / rho_field - h_sat)\r
 168         return T_field\r
 169 \r
 170 def rho_to_p_StiffenedGaz(rho_field, q_field, rho_E_field):\r
 171         p_field = (gamma - 1) * (rho_E_field - 1. / 2 * q_field ** 2 / rho_field) - gamma * p0\r
 172         return p_field\r
 173 \r
 174 def T_to_E_StiffenedGaz(p_field, T_field, v_field):\r
 175         rho_field = p_to_rho_StiffenedGaz(p_field, T_field)\r
 176         E_field = (p_field + gamma * p0) / ((gamma-1) * rho_field) + 0.5 * v_field **2\r
 177         return E_field\r
 178                 \r
 179 def dp_drho_e_const_StiffenedGaz( e ):\r
 180         return (gamma-1)*(e-q)\r
 181 \r
 182 def dp_de_rho_const_StiffenedGaz( rho ):\r
 183         return (gamma-1)*rho\r
 184 \r
 185 def sound_speed_StiffenedGaz( h ):\r
 186         return np.sqrt((gamma-1)*(h-q))\r
 187 \r
 188 def rho_h_to_p_StiffenedGaz( rho, h ):\r
 189         return (gamma - 1) * rho * ( h - q ) / gamma - p0\r
 190 \r
 191 def MatRoe_StiffenedGaz( rho_l, q_l, rho_E_l, rho_r, q_r, rho_E_r):\r
 192         RoeMat = cdmath.Matrix(3, 3)\r
 193 \r
 194         u_l = q_l / rho_l\r
 195         u_r = q_r / rho_r\r
 196         p_l = rho_to_p_StiffenedGaz(rho_l, q_l, rho_E_l)\r
 197         p_r = rho_to_p_StiffenedGaz(rho_r, q_r, rho_E_r)\r
 198         H_l = rho_E_l / rho_l + p_l / rho_l\r
 199         H_r = rho_E_r / rho_r + p_r / rho_r\r
 200 \r
 201         # Roe averages\r
 202         rho = np.sqrt(rho_l * rho_r)\r
 203         u   = (u_l * sqrt(rho_l) + u_r * sqrt(rho_r)) / (sqrt(rho_l) + sqrt(rho_r))\r
 204         H   = (H_l * sqrt(rho_l) + H_r * sqrt(rho_r)) / (sqrt(rho_l) + sqrt(rho_r))\r
 205 \r
 206         p = rho_h_to_p_StiffenedGaz( rho, H - u**2/2. )\r
 207         e = H - p / rho - 1./2 * u**2\r
 208         dp_drho = dp_drho_e_const_StiffenedGaz( e )\r
 209         dp_de   = dp_de_rho_const_StiffenedGaz( rho )\r
 210 \r
 211         RoeMat[0, 0] = 0\r
 212         RoeMat[0, 1] = 1\r
 213         RoeMat[0, 2] = 0\r
 214         RoeMat[1, 0] = dp_drho - u ** 2 + dp_de / rho * (u**2/2 - e)\r
 215         RoeMat[1, 1] = 2 * u - u * dp_de / rho\r
 216         RoeMat[1, 2] = dp_de / rho\r
 217         RoeMat[2, 0] = -u * ( -dp_drho + H - dp_de / rho * (u**2/2 - e) )\r
 218         RoeMat[2, 1] = H - dp_de / rho * u ** 2\r
 219         RoeMat[2, 2] = (dp_de / rho +1) * u\r
 220         \r
 221         return(RoeMat)\r
 222 \r
 223         \r
 224 def Droe_StiffenedGaz( rho_l, q_l, rho_E_l, rho_r, q_r, rho_E_r):\r
 225     Droe   = cdmath.Matrix(3, 3)\r
 226 \r
 227     u_l = q_l / rho_l\r
 228     u_r = q_r / rho_r\r
 229     p_l = rho_to_p_StiffenedGaz(rho_l, q_l, rho_E_l)\r
 230     p_r = rho_to_p_StiffenedGaz(rho_r, q_r, rho_E_r)\r
 231     H_l = rho_E_l / rho_l + p_l / rho_l\r
 232     H_r = rho_E_r / rho_r + p_r / rho_r\r
 233 \r
 234     # Roe averages\r
 235     rho = np.sqrt(rho_l * rho_r)\r
 236     u = (u_l * sqrt(rho_l) + u_r * sqrt(rho_r)) / (sqrt(rho_l) + sqrt(rho_r))\r
 237     H = (H_l * sqrt(rho_l) + H_r * sqrt(rho_r)) / (sqrt(rho_l) + sqrt(rho_r))\r
 238 \r
 239     c = sound_speed_StiffenedGaz( H - u**2/2. )\r
 240     \r
 241     lamb  = cdmath.Vector(3)\r
 242     lamb[0] = u-c\r
 243     lamb[1] = u\r
 244     lamb[2] = u+c    \r
 245 \r
 246     r   = cdmath.Matrix(3, 3)\r
 247     r[0,0] = 1.\r
 248     r[1,0] = u-c\r
 249     r[2,0] = H-u*c    \r
 250     r[0,1] = 1.\r
 251     r[1,1] = u   \r
 252     r[2,1] = H-c**2/(gamma-1)    \r
 253     r[0,2] = 1.\r
 254     r[1,2] = u+c\r
 255     r[2,2] = H+u*c         \r
 256 \r
 257     l   = cdmath.Matrix(3, 3)\r
 258     l[0,0] = (1./(2*c**2))*(0.5*(gamma-1)*u**2+u*c)\r
 259     l[1,0] = (1./(2*c**2))*(-u*(gamma-1)-c)\r
 260     l[2,0] = (1./(2*c**2))*(gamma-1)\r
 261     l[0,1] = ((gamma-1)/c**2)*(H-u**2)\r
 262     l[1,1] = ((gamma-1)/c**2)*u   \r
 263     l[2,1] = -((gamma-1)/c**2)    \r
 264     l[0,2] = (1./(2*c**2))*(0.5*(gamma-1)*u**2-u*c)\r
 265     l[1,2] = (1./(2*c**2))*(c-u*(gamma-1))\r
 266     l[2,2] = (1./(2*c**2))*(gamma-1)\r
 267 \r
 268     M1 = cdmath.Matrix(3, 3) #abs(lamb[0])*r[:,0].tensProduct(l[:,0])\r
 269     M2 = cdmath.Matrix(3, 3) #abs(lamb[1])*r[:,1].tensProduct(l[:,1])   \r
 270     M3 = cdmath.Matrix(3, 3) #abs(lamb[2])*r[:,2].tensProduct(l[:,2])\r
 271     for i in range(3):\r
 272         for j in range(3):\r
 273             M1[i,j] = abs(lamb[0])*r[i,0]*l[j,0]\r
 274             M2[i,j] = abs(lamb[1])*r[i,1]*l[j,1]            \r
 275             M3[i,j] = abs(lamb[2])*r[i,2]*l[j,2]\r
 276             \r
 277     Droe = M1+M2+M3 \r
 278     \r
 279     return(Droe)    \r
 280 \r
 281 \r
 282 def jacobianMatricesm(coeff, rho_l, q_l, rho_E_l, rho_r, q_r, rho_E_r):\r
 283 \r
 284         if rho_l < 0 or rho_r < 0:\r
 285                 print("rho_l=", rho_l, " rho_r= ", rho_r)\r
 286                 raise ValueError("Negative density")\r
 287         if rho_E_l < 0 or rho_E_r < 0:\r
 288                 print("rho_E_l=", rho_E_l, " rho_E_r= ", rho_E_r)\r
 289                 raise ValueError("Negative total energy")\r
 290 \r
 291         RoeMat = MatRoe_StiffenedGaz( rho_l, q_l, rho_E_l, rho_r, q_r, rho_E_r)\r
 292         \r
 293         Droe = Droe_StiffenedGaz( rho_l, q_l, rho_E_l, rho_r, q_r, rho_E_r)    \r
 294         return (RoeMat - Droe) * coeff * 0.5\r
 295 \r
 296 \r
 297 def jacobianMatricesp(coeff, rho_l, q_l, rho_E_l, rho_r, q_r, rho_E_r):\r
 298         if rho_l < 0 or rho_r < 0:\r
 299                 print("rho_l=", rho_l, " rho_r= ", rho_r)\r
 300                 raise ValueError("Negative density")\r
 301         if rho_E_l < 0 or rho_E_r < 0:\r
 302                 print("rho_E_l=", rho_E_l, " rho_E_r= ", rho_E_r)\r
 303                 raise ValueError("Negative total energy")\r
 304 \r
 305         RoeMat = MatRoe_StiffenedGaz( rho_l, q_l, rho_E_l, rho_r, q_r, rho_E_r)\r
 306         Droe = Droe_StiffenedGaz( rho_l, q_l, rho_E_l, rho_r, q_r, rho_E_r)    \r
 307 \r
 308         return (RoeMat + Droe) * coeff * 0.5\r
 309 \r
 310 \r
 311 def FillEdges(j, Uk, nbComp, divMat, Rhs, Un, dt, dx, isImplicit):\r
 312         dUi1 = cdmath.Vector(3)\r
 313         dUi2 = cdmath.Vector(3)\r
 314         temp1 = cdmath.Vector(3)\r
 315         temp2 = cdmath.Vector(3)\r
 316 \r
 317         if (j == 0):\r
 318                 rho_l   = Uk[j * nbComp + 0]\r
 319                 q_l     = Uk[j * nbComp + 1]\r
 320                 rho_E_l = Uk[j * nbComp + 2]\r
 321                 rho_r   = Uk[(j + 1) * nbComp + 0]\r
 322                 q_r     = Uk[(j + 1) * nbComp + 1]\r
 323                 rho_E_r = Uk[(j + 1) * nbComp + 2]      \r
 324 \r
 325                 Am = jacobianMatricesm(dt / dx, rho_l, q_l, rho_E_l, rho_r, q_r, rho_E_r)\r
 326                 divMat.addValue(j * nbComp, (j + 1) * nbComp, Am)\r
 327                 divMat.addValue(j * nbComp,       j * nbComp, Am * (-1.))\r
 328 \r
 329                 p_inlet = rho_to_p_StiffenedGaz(Uk[j * nbComp + 0], Uk[j * nbComp + 1], Uk[j * nbComp + 2])# We take p from inside the domain\r
 330                 rho_l=p_to_rho_StiffenedGaz(p_inlet, T_inlet) # rho is computed from the temperature BC and the inner pressure\r
 331                 #rho_l   = Uk[j * nbComp + 0]                            # We take rho from inside the domain\r
 332                 q_l     = rho_l * v_inlet                               # q is imposed by the boundary condition v_inlet\r
 333                 rho_E_l = T_to_rhoE_StiffenedGaz(T_inlet, rho_l, q_l)   #rhoE is obtained  using the two boundary conditions v_inlet and e_inlet\r
 334                 rho_r   = Uk[j * nbComp + 0]\r
 335                 q_r     = Uk[j * nbComp + 1]\r
 336                 rho_E_r = Uk[j * nbComp + 2]\r
 337 \r
 338                 Ap = jacobianMatricesp(dt / dx, rho_l, q_l, rho_E_l, rho_r, q_r, rho_E_r)\r
 339                 divMat.addValue(j * nbComp, j * nbComp, Ap)\r
 340         \r
 341                 if(isImplicit):\r
 342                         dUi1[0] = Uk[(j + 1) * nbComp + 0] - Uk[j * nbComp + 0]\r
 343                         dUi1[1] = Uk[(j + 1) * nbComp + 1] - Uk[j * nbComp + 1]\r
 344                         dUi1[2] = Uk[(j + 1) * nbComp + 2] - Uk[j * nbComp + 2]\r
 345                         temp1 = Am * dUi1\r
 346                 \r
 347                         dUi2[0] = rho_l   -  Uk[(j ) * nbComp + 0]\r
 348                         dUi2[1] = q_l     -  Uk[(j ) * nbComp + 1]\r
 349                         dUi2[2] = rho_E_l -  Uk[(j ) * nbComp + 2]\r
 350                         temp2 = Ap * dUi2\r
 351                 else:\r
 352                         dUi2[0] = rho_l   \r
 353                         dUi2[1] = q_l     \r
 354                         dUi2[2] = rho_E_l \r
 355                         temp2 = Ap * dUi2\r
 356 \r
 357         elif (j == nx - 1):\r
 358                 rho_l   = Uk[(j - 1) * nbComp + 0]\r
 359                 q_l     = Uk[(j - 1) * nbComp + 1]\r
 360                 rho_E_l = Uk[(j - 1) * nbComp + 2]\r
 361                 rho_r   = Uk[j * nbComp + 0]\r
 362                 q_r     = Uk[j * nbComp + 1]\r
 363                 rho_E_r = Uk[j * nbComp + 2]\r
 364 \r
 365                 Ap = jacobianMatricesp(dt / dx, rho_l, q_l, rho_E_l, rho_r, q_r, rho_E_r)\r
 366                 divMat.addValue(j * nbComp, j * nbComp, Ap)\r
 367                 divMat.addValue(j * nbComp, (j - 1) * nbComp, Ap * (-1.))\r
 368 \r
 369                 rho_l   = Uk[j * nbComp + 0]\r
 370                 q_l     = Uk[j * nbComp + 1]\r
 371                 rho_E_l = Uk[j * nbComp + 2]\r
 372                 rho_r   = Uk[(j ) * nbComp + 0]                               # We take rho inside the domain\r
 373                 q_r     = Uk[(j ) * nbComp + 1]                               # We take q from inside the domain\r
 374                 rho_E_r = (p_outlet+gamma*p0)/(gamma-1) + 0.5*q_r**2/rho_r    # rhoE is obtained using the boundary condition p_outlet\r
 375 \r
 376                 Am = jacobianMatricesm(dt / dx, rho_l, q_l, rho_E_l, rho_r, q_r, rho_E_r)       \r
 377                 divMat.addValue(j * nbComp, j * nbComp, Am * (-1.))\r
 378 \r
 379                 if(isImplicit):\r
 380                         dUi1[0] = rho_r   - Uk[j * nbComp + 0]\r
 381                         dUi1[1] = q_r     - Uk[j * nbComp + 1]\r
 382                         dUi1[2] = rho_E_r - Uk[j * nbComp + 2]\r
 383                         temp1 = Am * dUi1\r
 384 \r
 385                         dUi2[0] = Uk[(j - 1) * nbComp + 0] - Uk[j * nbComp + 0]\r
 386                         dUi2[1] = Uk[(j - 1) * nbComp + 1] - Uk[j * nbComp + 1]\r
 387                         dUi2[2] = Uk[(j - 1) * nbComp + 2] - Uk[j * nbComp + 2]\r
 388                         temp2 = Ap * dUi2\r
 389                 else:\r
 390                         dUi1[0] = rho_r   \r
 391                         dUi1[1] = q_r     \r
 392                         dUi1[2] = rho_E_r \r
 393                         temp1 = Am * dUi1\r
 394         \r
 395         if(isImplicit):#implicit scheme, contribution from the Newton scheme residual\r
 396                 Rhs[j * nbComp + 0] = -temp1[0] + temp2[0] - (Uk[j * nbComp + 0] - Un[j * nbComp + 0])\r
 397                 Rhs[j * nbComp + 1] = -temp1[1] + temp2[1] - (Uk[j * nbComp + 1] - Un[j * nbComp + 1])\r
 398                 Rhs[j * nbComp + 2] = -temp1[2] + temp2[2] - (Uk[j * nbComp + 2] - Un[j * nbComp + 2])\r
 399         else:#explicit scheme, contribution from the boundary data the right hand side\r
 400                 Rhs[j * nbComp + 0] = -temp1[0] + temp2[0] \r
 401                 Rhs[j * nbComp + 1] = -temp1[1] + temp2[1] \r
 402                 Rhs[j * nbComp + 2] = -temp1[2] + temp2[2] \r
 403 \r
 404 def FillInnerCell(j, Uk, nbComp, divMat, Rhs, Un, dt, dx, isImplicit):\r
 405 \r
 406         rho_l   = Uk[(j - 1) * nbComp + 0]\r
 407         q_l     = Uk[(j - 1) * nbComp + 1]\r
 408         rho_E_l = Uk[(j - 1) * nbComp + 2]\r
 409         rho_r   = Uk[j * nbComp + 0]\r
 410         q_r     = Uk[j * nbComp + 1]\r
 411         rho_E_r = Uk[j * nbComp + 2]\r
 412         Ap = jacobianMatricesp(dt / dx, rho_l, q_l, rho_E_l, rho_r, q_r, rho_E_r)\r
 413         \r
 414         rho_l   = Uk[j * nbComp + 0]\r
 415         q_l     = Uk[j * nbComp + 1]\r
 416         rho_E_l = Uk[j * nbComp + 2]\r
 417         rho_r   = Uk[(j + 1) * nbComp + 0]\r
 418         q_r     = Uk[(j + 1) * nbComp + 1]\r
 419         rho_E_r = Uk[(j + 1) * nbComp + 2]\r
 420         Am = jacobianMatricesm(dt / dx, rho_l, q_l, rho_E_l, rho_r, q_r, rho_E_r)\r
 421 \r
 422         divMat.addValue(j * nbComp, (j + 1) * nbComp, Am)\r
 423         divMat.addValue(j * nbComp, j * nbComp, Am * (-1.))\r
 424         divMat.addValue(j * nbComp, j * nbComp, Ap)\r
 425         divMat.addValue(j * nbComp, (j - 1) * nbComp, Ap * (-1.))\r
 426 \r
 427         if(isImplicit):\r
 428                 dUi1 = cdmath.Vector(3)\r
 429                 dUi2 = cdmath.Vector(3)\r
 430                 \r
 431                 dUi1[0] = Uk[(j + 1) * nbComp + 0] - Uk[j * nbComp + 0]\r
 432                 dUi1[1] = Uk[(j + 1) * nbComp + 1] - Uk[j * nbComp + 1]\r
 433                 dUi1[2] = Uk[(j + 1) * nbComp + 2] - Uk[j * nbComp + 2]\r
 434         \r
 435                 dUi2[0] = Uk[(j - 1) * nbComp + 0] - Uk[j * nbComp + 0]\r
 436                 dUi2[1] = Uk[(j - 1) * nbComp + 1] - Uk[j * nbComp + 1]\r
 437                 dUi2[2] = Uk[(j - 1) * nbComp + 2] - Uk[j * nbComp + 2]\r
 438                 \r
 439                 temp1 = Am * dUi1\r
 440                 temp2 = Ap * dUi2\r
 441 \r
 442                 Rhs[j * nbComp + 0] = -temp1[0] + temp2[0] - (Uk[j * nbComp + 0] - Un[j * nbComp + 0])\r
 443                 Rhs[j * nbComp + 1] = -temp1[1] + temp2[1] - (Uk[j * nbComp + 1] - Un[j * nbComp + 1])\r
 444                 Rhs[j * nbComp + 2] = -temp1[2] + temp2[2] - (Uk[j * nbComp + 2] - Un[j * nbComp + 2])\r
 445         else:\r
 446                 Rhs[j * nbComp + 0] = 0\r
 447                 Rhs[j * nbComp + 1] = 0\r
 448                 Rhs[j * nbComp + 2] = 0\r
 449 \r
 450 def SetPicture(rho_field_Roe, q_field_Roe, h_field_Roe, p_field_Roe, v_field_Roe, T_field_Roe, dx):\r
 451         max_initial_rho = max(rho_field_Roe)\r
 452         min_initial_rho = min(rho_field_Roe)\r
 453         max_initial_q = max(q_field_Roe)\r
 454         min_initial_q = min(q_field_Roe)\r
 455         min_initial_h = min(h_field_Roe)\r
 456         max_initial_h = max(h_field_Roe)\r
 457         max_initial_p = max(p_field_Roe)\r
 458         min_initial_p = min(p_field_Roe)\r
 459         max_initial_v = max(v_field_Roe)\r
 460         min_initial_v = min(v_field_Roe)\r
 461         max_initial_T = max(T_field_Roe)\r
 462         min_initial_T = min(T_field_Roe)\r
 463 \r
 464         fig, ([axDensity, axMomentum, axh],[axPressure, axVitesse, axTemperature]) = plt.subplots(2, 3,sharex=True, figsize=(14,10))\r
 465         plt.gcf().subplots_adjust(wspace = 0.5)\r
 466 \r
 467         lineDensity_Roe, = axDensity.plot([a+0.5*dx + i*dx for i in range(nx)], rho_field_Roe, label='Roe')\r
 468         axDensity.set(xlabel='x (m)', ylabel='Densité (kg/m3)')\r
 469         axDensity.set_xlim(a,b)\r
 470         axDensity.set_ylim(680, 800)\r
 471         axDensity.legend()\r
 472 \r
 473         lineMomentum_Roe, = axMomentum.plot([a+0.5*dx + i*dx for i in range(nx)], q_field_Roe, label='Roe')\r
 474         axMomentum.set(xlabel='x (m)', ylabel='Momentum (kg/(m2.s))')\r
 475         axMomentum.set_xlim(a,b)\r
 476         axMomentum.set_ylim(3500,       4000)\r
 477         axMomentum.legend()\r
 478 \r
 479         lineh_Roe, = axh.plot([a+0.5*dx + i*dx for i in range(nx)], h_field_Roe, label='Roe')\r
 480         axh.set(xlabel='x (m)', ylabel='h (J/m3)')\r
 481         axh.set_xlim(a,b)\r
 482         axh.set_ylim(1.2 * 10**6, 1.5*10**6)\r
 483         axh.legend()\r
 484         \r
 485         linePressure_Roe, = axPressure.plot([a+0.5*dx + i*dx for i in range(nx)], p_field_Roe, label='Roe')\r
 486         axPressure.set(xlabel='x (m)', ylabel='Pression (bar)')\r
 487         axPressure.set_xlim(a,b)\r
 488         axPressure.set_ylim(155, 155.5)\r
 489         axPressure.ticklabel_format(axis='y', style='sci', scilimits=(0,0))\r
 490         axPressure.legend()\r
 491 \r
 492         lineVitesse_Roe, = axVitesse.plot([a+0.5*dx + i*dx for i in range(nx)], v_field_Roe, label='Roe')\r
 493         axVitesse.set(xlabel='x (m)', ylabel='Vitesse (m/s)')\r
 494         axVitesse.set_xlim(a,b)\r
 495         axVitesse.set_ylim(v_0-1, v_0+1)\r
 496         axVitesse.ticklabel_format(axis='y', style='sci', scilimits=(0,0))\r
 497         axVitesse.legend()\r
 498 \r
 499         lineTemperature_Roe, = axTemperature.plot([a+0.5*dx + i*dx for i in range(nx)], T_field_Roe, label='Roe')\r
 500         axTemperature.set(xlabel='x (m)', ylabel='Température (K)')\r
 501         axTemperature.set_xlim(a,b)\r
 502         axTemperature.set_ylim(T_0-10, T_0+30)\r
 503         axTemperature.ticklabel_format(axis='y', style='sci', scilimits=(0,0))\r
 504         axTemperature.legend()\r
 505         \r
 506         return(fig, lineDensity_Roe, lineMomentum_Roe, lineh_Roe, linePressure_Roe, lineVitesse_Roe, lineTemperature_Roe)\r
 507 \r
 508 \r
 509 def EulerSystemRoe(ntmax, tmax, cfl, a, b, nbCells, output_freq, meshName, state_law, isImplicit):\r
 510         state_law_parameters(state_law)\r
 511         dim = 1\r
 512         nbComp = 3\r
 513         dt = 0.\r
 514         time = 0.\r
 515         it = 0\r
 516         isStationary = False\r
 517         dx = (b - a) / nx\r
 518         dt = cfl * dx / c0\r
 519         #dt = 5*10**(-6)\r
 520         nbVoisinsMax = 2\r
 521 \r
 522         # iteration vectors\r
 523         Un_Roe  = cdmath.Vector(nbCells * (nbComp))\r
 524         dUn_Roe = cdmath.Vector(nbCells * (nbComp))\r
 525         dUk_Roe = cdmath.Vector(nbCells * (nbComp))\r
 526         Rhs_Roe = cdmath.Vector(nbCells * (nbComp))\r
 527 \r
 528         # Initial conditions\r
 529         print("Construction of the initial condition …")\r
 530 \r
 531         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
 532         h_field_Roe = (rho_E_field_Roe + p_field_Roe) / rho_field_Roe - 0.5 * (q_field_Roe / rho_field_Roe) **2\r
 533         p_field_Roe = p_field_Roe * 10 ** (-5)\r
 534         \r
 535 \r
 536         for k in range(nbCells):\r
 537                 Un_Roe[k * nbComp + 0] = rho_field_Roe[k]\r
 538                 Un_Roe[k * nbComp + 1] = q_field_Roe[k]\r
 539                 Un_Roe[k * nbComp + 2] = rho_E_field_Roe[k]\r
 540 \r
 541         divMat_Roe = cdmath.SparseMatrixPetsc(nbCells * nbComp, nbCells * nbComp, (nbVoisinsMax + 1) * nbComp)\r
 542         \r
 543         # Picture settings\r
 544         fig, 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
 545 \r
 546         plt.savefig("EulerComplet_HeatedChannel_" + str(dim) + "D_Roe" + meshName + "0" + ".png")\r
 547         iterGMRESMax = 50\r
 548         newton_max = 100\r
 549 \r
 550         print("Starting computation of the non linear Euler non isentropic system with Roe scheme …")\r
 551         # STARTING TIME LOOP\r
 552         while (it < ntmax and time <= tmax and not isStationary):\r
 553                 dUn_Roe = Un_Roe.deepCopy()\r
 554                 Uk_Roe  = Un_Roe.deepCopy()\r
 555                 residu_Roe = 1.\r
 556                 \r
 557                 k_Roe = 0\r
 558                 while (k_Roe < newton_max and residu_Roe > precision):\r
 559                         # STARTING NEWTON LOOP\r
 560                         divMat_Roe.zeroEntries()  #sets the matrix coefficients to zero\r
 561                         for j in range(nbCells):\r
 562                                 \r
 563                                 # traitements des bords\r
 564                                 if (j == 0):\r
 565                                         FillEdges(j, Uk_Roe, nbComp, divMat_Roe, Rhs_Roe, Un_Roe, dt, dx, isImplicit)\r
 566                                 elif (j == nbCells - 1):\r
 567                                         FillEdges(j, Uk_Roe, nbComp, divMat_Roe, Rhs_Roe, Un_Roe, dt, dx, isImplicit)\r
 568 \r
 569                                 # traitement des cellules internes\r
 570                                 else:\r
 571                                         FillInnerCell(j, Uk_Roe, nbComp, divMat_Roe, Rhs_Roe, Un_Roe, dt, dx, isImplicit)\r
 572                                         \r
 573                                 Rhs_Roe[j * nbComp + 2] += phi*dt\r
 574                         \r
 575                         if(isImplicit):\r
 576                                 #solving the linear system divMat * dUk = Rhs\r
 577                                 divMat_Roe.diagonalShift(1.)\r
 578                                 LS_Roe = cdmath.LinearSolver(divMat_Roe, Rhs_Roe, iterGMRESMax, precision, "GMRES", "LU")\r
 579                                 dUk_Roe = LS_Roe.solve()\r
 580                                 vector_residu_Roe = dUk_Roe.maxVector(nbComp)\r
 581                                 residu_Roe = max(abs(vector_residu_Roe[0])/rho0, abs(vector_residu_Roe[1])/(rho0*v_0), abs(vector_residu_Roe[2])/rhoE0 )\r
 582                         else:\r
 583                                 dUk_Roe=Rhs_Roe - divMat_Roe*Un_Roe\r
 584                                 residu_Roe = 0.#Convergence schéma Newton\r
 585                         \r
 586                         if (isImplicit and (it == 1 or it % output_freq == 0 or it >= ntmax or isStationary or time >= tmax)):\r
 587                                 print("Residu Newton at iteration ",k_Roe, " :   ", residu_Roe)\r
 588                                 print("Linear system converged in ", LS_Roe.getNumberOfIter(), " GMRES iterations")\r
 589 \r
 590                         #updates for Newton loop\r
 591                         Uk_Roe += dUk_Roe\r
 592                         k_Roe = k_Roe + 1\r
 593                         if (isImplicit and not LS_Roe.getStatus()):\r
 594                                 print("Linear system did not converge ", LS_Roe.getNumberOfIter(), " GMRES iterations")\r
 595                                 raise ValueError("No convergence of the linear system")\r
 596                         \r
 597                         if k_Roe == newton_max:\r
 598                                 raise ValueError("No convergence of Newton Roe Scheme")\r
 599 \r
 600                 #updating fields\r
 601                 Un_Roe = Uk_Roe.deepCopy()\r
 602                 dUn_Roe -= Un_Roe\r
 603 \r
 604                 #Testing stationarity\r
 605                 residu_stat = dUn_Roe.maxVector(nbComp)#On prend le max de chaque composante\r
 606                 if (it % output_freq == 0 ):\r
 607                         print("Test de stationarité : Un+1-Un= ", max(abs(residu_stat[0])/rho0, abs(residu_stat[1])/(rho0*v_0), abs(residu_stat[2])/rhoE0 ))\r
 608 \r
 609                 if ( it>1 and abs(residu_stat[0])/rho0<precision  and abs(residu_stat[1])/(rho0*v_0)<precision and abs(residu_stat[2])/rhoE0<precision):\r
 610                                 isStationary = True\r
 611                 \r
 612                 for k in range(nbCells):\r
 613                         rho_field_Roe[k]   = Un_Roe[k * nbComp + 0]\r
 614                         q_field_Roe[k]     = Un_Roe[k * nbComp + 1]\r
 615                         rho_E_field_Roe[k] = Un_Roe[k * nbComp + 2]\r
 616 \r
 617                 v_field_Roe = q_field_Roe / rho_field_Roe\r
 618                 p_field_Roe = rho_to_p_StiffenedGaz(rho_field_Roe, q_field_Roe, rho_E_field_Roe)\r
 619                 T_field_Roe = rhoE_to_T_StiffenedGaz(rho_field_Roe, q_field_Roe, rho_E_field_Roe)\r
 620                 h_field_Roe = (rho_E_field_Roe + p_field_Roe) / rho_field_Roe - 0.5 * (q_field_Roe / rho_field_Roe) **2\r
 621                 p_field_Roe = p_field_Roe * 10 ** (-5)\r
 622                 \r
 623                 if( min(p_field_Roe)<0) :\r
 624                         raise ValueError("Negative pressure, stopping calculation")\r
 625 \r
 626                 #picture and video updates\r
 627                 lineDensity_Roe.set_ydata(rho_field_Roe)\r
 628                 lineMomentum_Roe.set_ydata(q_field_Roe)\r
 629                 lineRhoE_Roe.set_ydata(h_field_Roe)\r
 630                 linePressure_Roe.set_ydata(p_field_Roe)\r
 631                 lineVitesse_Roe.set_ydata(v_field_Roe)\r
 632                 lineTemperature_Roe.set_ydata(T_field_Roe)\r
 633                 \r
 634                 time = time + dt\r
 635                 it = it + 1\r
 636 \r
 637                 # Savings\r
 638                 if (it == 1 or it % output_freq == 0 or it >= ntmax or isStationary or time >= tmax):\r
 639         \r
 640                         print("-- Time step : " + str(it) + ", Time: " + str(time) + ", dt: " + str(dt))\r
 641 \r
 642                         print("Temperature gain between inlet and outlet is ", T_field_Roe[nbCells-1]-T_field_Roe[0],"\n")\r
 643 \r
 644                         plt.savefig("EulerComplet_HeatedChannel_" + str(dim) + "D_Roe" + meshName + "Implicit"+str(isImplicit)+ str(it) + '_time' + str(time) + ".png")\r
 645 \r
 646         print("-- Iter: " + str(it) + ", Time: " + str(time) + ", dt: " + str(dt)+"\n")\r
 647         \r
 648         if (it >= ntmax):\r
 649                 print("Maximum number of time steps ntmax= ", ntmax, " reached")\r
 650                 return\r
 651 \r
 652         elif (isStationary):\r
 653                 print("Stationary regime reached at time step ", it, ", t= ", time)\r
 654                 print("------------------------------------------------------------------------------------")\r
 655                 np.savetxt("EulerComplet_HeatedChannel_" + str(dim) + "D_Roe" + meshName + "_rho_Stat.txt", rho_field_Roe, delimiter="\n")\r
 656                 np.savetxt("EulerComplet_HeatedChannel_" + str(dim) + "D_Roe" + meshName + "_q_Stat.txt", q_field_Roe, delimiter="\n")\r
 657                 np.savetxt("EulerComplet_HeatedChannel_" + str(dim) + "D_Roe" + meshName + "_rhoE_Stat.txt", rho_E_field_Roe, delimiter="\n")\r
 658                 np.savetxt("EulerComplet_HeatedChannel_" + str(dim) + "D_Roe" + meshName + "_p_Stat.txt", p_field_Roe, delimiter="\n")\r
 659                 plt.savefig("EulerComplet_HeatedChannel_" + str(dim) + "D_Roe" + meshName + "Implicit"+str(isImplicit)+"_Stat.png")\r
 660                 return\r
 661         else:\r
 662                 print("Maximum time Tmax= ", tmax, " reached")\r
 663                 return\r
 664 \r
 665 \r
 666 def solve(a, b, nx, meshName, meshType, cfl, state_law, isImplicit):\r
 667         print("Simulation of a heated channel in dimension 1 on " + str(nx) + " cells")\r
 668         print("State Law Stiffened Gaz, " + state_law)\r
 669         print("Initial data : ", "constant fields")\r
 670         print("Boundary conditions : ", "Inlet (Left), Outlet (Right)")\r
 671         print("Mesh name : ", meshName, ", ", nx, " cells")\r
 672         # Problem data\r
 673         tmax = 10.\r
 674         ntmax = 100000\r
 675         output_freq = 1000\r
 676         EulerSystemRoe(ntmax, tmax, cfl, a, b, nx, output_freq, meshName, state_law, isImplicit)\r
 677         return\r
 678 \r
 679 def FillMatrixFromEdges(j, Uk, nbComp, divMat, dt, dx):\r
 680 \r
 681         if (j == 0):\r
 682                 rho_l   = Uk[j * nbComp + 0]\r
 683                 q_l     = Uk[j * nbComp + 1]\r
 684                 rho_E_l = Uk[j * nbComp + 2]\r
 685                 rho_r   = Uk[(j + 1) * nbComp + 0]\r
 686                 q_r     = Uk[(j + 1) * nbComp + 1]\r
 687                 rho_E_r = Uk[(j + 1) * nbComp + 2]      \r
 688 \r
 689                 Am = jacobianMatricesm(dt / dx, rho_l, q_l, rho_E_l, rho_r, q_r, rho_E_r)\r
 690                 divMat.addValue(j * nbComp, (j + 1) * nbComp, Am)\r
 691                 divMat.addValue(j * nbComp,       j * nbComp, Am * (-1.))\r
 692         \r
 693                 p_inlet = rho_to_p_StiffenedGaz(Uk[j * nbComp + 0], Uk[j * nbComp + 1], Uk[j * nbComp + 2])# We take p from inside the domain\r
 694                 rho_l=p_to_rho_StiffenedGaz(p_inlet, T_inlet) # rho is computed from the temperature BC and the inner pressure\r
 695                 #rho_l   = Uk[j * nbComp + 0]                            # We take rho from inside the domain\r
 696                 q_l     = rho_l * v_inlet                               # q is imposed by the boundary condition v_inlet\r
 697                 rho_E_l = T_to_rhoE_StiffenedGaz(T_inlet, rho_l, q_l)   #rhoE is obtained  using the two boundary conditions v_inlet and e_inlet\r
 698                 rho_r   = Uk[j * nbComp + 0]\r
 699                 q_r     = Uk[j * nbComp + 1]\r
 700                 rho_E_r = Uk[j * nbComp + 2]\r
 701 \r
 702                 Ap = jacobianMatricesp(dt / dx, rho_l, q_l, rho_E_l, rho_r, q_r, rho_E_r)\r
 703                 divMat.addValue(j * nbComp, j * nbComp, Ap)\r
 704 \r
 705         elif (j == nx - 1):\r
 706                 rho_l   = Uk[j * nbComp + 0]\r
 707                 q_l     = Uk[j * nbComp + 1]\r
 708                 rho_E_l = Uk[j * nbComp + 2]\r
 709                 rho_r   = Uk[(j ) * nbComp + 0]                               # We take rho inside the domain\r
 710                 q_r     = Uk[(j ) * nbComp + 1]                               # We take q from inside the domain\r
 711                 rho_E_r = (p_outlet+gamma*p0)/(gamma-1) + 0.5*q_r**2/rho_r    # rhoE is obtained using the boundary condition p_outlet\r
 712 \r
 713                 Am = jacobianMatricesm(dt / dx, rho_l, q_l, rho_E_l, rho_r, q_r, rho_E_r)       \r
 714                 divMat.addValue(j * nbComp, j * nbComp, Am * (-1.))\r
 715         \r
 716                 rho_l   = Uk[(j - 1) * nbComp + 0]\r
 717                 q_l     = Uk[(j - 1) * nbComp + 1]\r
 718                 rho_E_l = Uk[(j - 1) * nbComp + 2]\r
 719                 rho_r   = Uk[j * nbComp + 0]\r
 720                 q_r     = Uk[j * nbComp + 1]\r
 721                 rho_E_r = Uk[j * nbComp + 2]\r
 722 \r
 723                 Ap = jacobianMatricesp(dt / dx, rho_l, q_l, rho_E_l, rho_r, q_r, rho_E_r)\r
 724                 divMat.addValue(j * nbComp, j * nbComp, Ap)\r
 725                 divMat.addValue(j * nbComp, (j - 1) * nbComp, Ap * (-1.))\r
 726 \r
 727 \r
 728 def FillMatrixFromInnerCell(j, Uk, nbComp, divMat, dt, dx):\r
 729 \r
 730         rho_l   = Uk[(j - 1) * nbComp + 0]\r
 731         q_l     = Uk[(j - 1) * nbComp + 1]\r
 732         rho_E_l = Uk[(j - 1) * nbComp + 2]\r
 733         rho_r   = Uk[j * nbComp + 0]\r
 734         q_r     = Uk[j * nbComp + 1]\r
 735         rho_E_r = Uk[j * nbComp + 2]\r
 736         Ap = jacobianMatricesp(dt / dx, rho_l, q_l, rho_E_l, rho_r, q_r, rho_E_r)\r
 737         \r
 738         rho_l   = Uk[j * nbComp + 0]\r
 739         q_l     = Uk[j * nbComp + 1]\r
 740         rho_E_l = Uk[j * nbComp + 2]\r
 741         rho_r   = Uk[(j + 1) * nbComp + 0]\r
 742         q_r     = Uk[(j + 1) * nbComp + 1]\r
 743         rho_E_r = Uk[(j + 1) * nbComp + 2]\r
 744         Am = jacobianMatricesm(dt / dx, rho_l, q_l, rho_E_l, rho_r, q_r, rho_E_r)\r
 745 \r
 746         divMat.addValue(j * nbComp, (j + 1) * nbComp, Am)\r
 747         divMat.addValue(j * nbComp, j * nbComp, Am * (-1.))\r
 748         divMat.addValue(j * nbComp, j * nbComp, Ap)\r
 749         divMat.addValue(j * nbComp, (j - 1) * nbComp, Ap * (-1.))\r
 750                         \r
 751 def FillRHSFromEdges(j, Uk, nbComp, Rhs, Un, dt, dx, isImplicit):\r
 752         dUi1 = cdmath.Vector(3)\r
 753         dUi2 = cdmath.Vector(3)\r
 754         temp1 = cdmath.Vector(3)\r
 755         temp2 = cdmath.Vector(3)\r
 756 \r
 757         if (j == 0):\r
 758                 rho_l   = Uk[j * nbComp + 0]\r
 759                 q_l     = Uk[j * nbComp + 1]\r
 760                 rho_E_l = Uk[j * nbComp + 2]\r
 761                 rho_r   = Uk[(j + 1) * nbComp + 0]\r
 762                 q_r     = Uk[(j + 1) * nbComp + 1]\r
 763                 rho_E_r = Uk[(j + 1) * nbComp + 2]      \r
 764 \r
 765                 Am = jacobianMatricesm(dt / dx, rho_l, q_l, rho_E_l, rho_r, q_r, rho_E_r)\r
 766 \r
 767                 p_inlet = rho_to_p_StiffenedGaz(Uk[j * nbComp + 0], Uk[j * nbComp + 1], Uk[j * nbComp + 2])# We take p from inside the domain\r
 768                 rho_l=p_to_rho_StiffenedGaz(p_inlet, T_inlet) # rho is computed from the temperature BC and the inner pressure\r
 769                 #rho_l   = Uk[j * nbComp + 0]                            # We take rho from inside the domain\r
 770                 q_l     = rho_l * v_inlet                               # q is imposed by the boundary condition v_inlet\r
 771                 rho_E_l = T_to_rhoE_StiffenedGaz(T_inlet, rho_l, q_l)   #rhoE is obtained  using the two boundary conditions v_inlet and e_inlet\r
 772                 rho_r   = Uk[j * nbComp + 0]\r
 773                 q_r     = Uk[j * nbComp + 1]\r
 774                 rho_E_r = Uk[j * nbComp + 2]\r
 775 \r
 776                 Ap = jacobianMatricesp(dt / dx, rho_l, q_l, rho_E_l, rho_r, q_r, rho_E_r)\r
 777         \r
 778                 if(isImplicit):\r
 779                         dUi1[0] = Uk[(j + 1) * nbComp + 0] - Uk[j * nbComp + 0]\r
 780                         dUi1[1] = Uk[(j + 1) * nbComp + 1] - Uk[j * nbComp + 1]\r
 781                         dUi1[2] = Uk[(j + 1) * nbComp + 2] - Uk[j * nbComp + 2]\r
 782                         temp1 = Am * dUi1\r
 783                 \r
 784                         dUi2[0] = rho_l   -  Uk[(j ) * nbComp + 0]\r
 785                         dUi2[1] = q_l     -  Uk[(j ) * nbComp + 1]\r
 786                         dUi2[2] = rho_E_l -  Uk[(j ) * nbComp + 2]\r
 787                         temp2 = Ap * dUi2\r
 788                 else:\r
 789                         dUi2[0] = rho_l   \r
 790                         dUi2[1] = q_l     \r
 791                         dUi2[2] = rho_E_l \r
 792                         temp2 = Ap * dUi2\r
 793 \r
 794         elif (j == nx - 1):\r
 795                 rho_l   = Uk[(j - 1) * nbComp + 0]\r
 796                 q_l     = Uk[(j - 1) * nbComp + 1]\r
 797                 rho_E_l = Uk[(j - 1) * nbComp + 2]\r
 798                 rho_r   = Uk[j * nbComp + 0]\r
 799                 q_r     = Uk[j * nbComp + 1]\r
 800                 rho_E_r = Uk[j * nbComp + 2]\r
 801 \r
 802                 Ap = jacobianMatricesp(dt / dx, rho_l, q_l, rho_E_l, rho_r, q_r, rho_E_r)\r
 803 \r
 804                 rho_l   = Uk[j * nbComp + 0]\r
 805                 q_l     = Uk[j * nbComp + 1]\r
 806                 rho_E_l = Uk[j * nbComp + 2]\r
 807                 rho_r   = Uk[(j ) * nbComp + 0]                               # We take rho inside the domain\r
 808                 q_r     = Uk[(j ) * nbComp + 1]                               # We take q from inside the domain\r
 809                 rho_E_r = (p_outlet+gamma*p0)/(gamma-1) + 0.5*q_r**2/rho_r    # rhoE is obtained using the boundary condition p_outlet\r
 810 \r
 811                 Am = jacobianMatricesm(dt / dx, rho_l, q_l, rho_E_l, rho_r, q_r, rho_E_r)       \r
 812 \r
 813                 if(isImplicit):\r
 814                         dUi1[0] = rho_r   - Uk[j * nbComp + 0]\r
 815                         dUi1[1] = q_r     - Uk[j * nbComp + 1]\r
 816                         dUi1[2] = rho_E_r - Uk[j * nbComp + 2]\r
 817                         temp1 = Am * dUi1\r
 818         \r
 819                         dUi2[0] = Uk[(j - 1) * nbComp + 0] - Uk[j * nbComp + 0]\r
 820                         dUi2[1] = Uk[(j - 1) * nbComp + 1] - Uk[j * nbComp + 1]\r
 821                         dUi2[2] = Uk[(j - 1) * nbComp + 2] - Uk[j * nbComp + 2]\r
 822                         temp2 = Ap * dUi2\r
 823                 else:\r
 824                         dUi1[0] = rho_r   \r
 825                         dUi1[1] = q_r     \r
 826                         dUi1[2] = rho_E_r \r
 827                         temp1 = Am * dUi1\r
 828         \r
 829         if(isImplicit):\r
 830                 Rhs[j * nbComp + 0] = -temp1[0] + temp2[0] - (Uk[j * nbComp + 0] - Un[j * nbComp + 0])\r
 831                 Rhs[j * nbComp + 1] = -temp1[1] + temp2[1] - (Uk[j * nbComp + 1] - Un[j * nbComp + 1])\r
 832                 Rhs[j * nbComp + 2] = -temp1[2] + temp2[2] - (Uk[j * nbComp + 2] - Un[j * nbComp + 2])\r
 833         else:#explicit scheme, contribution from the boundary data the right hand side\r
 834                 Rhs[j * nbComp + 0] = -temp1[0] + temp2[0] \r
 835                 Rhs[j * nbComp + 1] = -temp1[1] + temp2[1] \r
 836                 Rhs[j * nbComp + 2] = -temp1[2] + temp2[2] \r
 837 \r
 838 def FillRHSFromInnerCell(j, Uk, nbComp, Rhs, Un, dt, dx, isImplicit):\r
 839 \r
 840         rho_l   = Uk[(j - 1) * nbComp + 0]\r
 841         q_l     = Uk[(j - 1) * nbComp + 1]\r
 842         rho_E_l = Uk[(j - 1) * nbComp + 2]\r
 843         rho_r   = Uk[j * nbComp + 0]\r
 844         q_r     = Uk[j * nbComp + 1]\r
 845         rho_E_r = Uk[j * nbComp + 2]\r
 846         Ap = jacobianMatricesp(dt / dx, rho_l, q_l, rho_E_l, rho_r, q_r, rho_E_r)\r
 847         \r
 848         rho_l   = Uk[j * nbComp + 0]\r
 849         q_l     = Uk[j * nbComp + 1]\r
 850         rho_E_l = Uk[j * nbComp + 2]\r
 851         rho_r   = Uk[(j + 1) * nbComp + 0]\r
 852         q_r     = Uk[(j + 1) * nbComp + 1]\r
 853         rho_E_r = Uk[(j + 1) * nbComp + 2]\r
 854         Am = jacobianMatricesm(dt / dx, rho_l, q_l, rho_E_l, rho_r, q_r, rho_E_r)\r
 855 \r
 856         if(isImplicit):#Contribution to the right hand side if te scheme is implicit\r
 857                 dUi1 = cdmath.Vector(3)\r
 858                 dUi2 = cdmath.Vector(3)\r
 859                 \r
 860                 dUi1[0] = Uk[(j + 1) * nbComp + 0] - Uk[j * nbComp + 0]\r
 861                 dUi1[1] = Uk[(j + 1) * nbComp + 1] - Uk[j * nbComp + 1]\r
 862                 dUi1[2] = Uk[(j + 1) * nbComp + 2] - Uk[j * nbComp + 2]\r
 863         \r
 864                 dUi2[0] = Uk[(j - 1) * nbComp + 0] - Uk[j * nbComp + 0]\r
 865                 dUi2[1] = Uk[(j - 1) * nbComp + 1] - Uk[j * nbComp + 1]\r
 866                 dUi2[2] = Uk[(j - 1) * nbComp + 2] - Uk[j * nbComp + 2]\r
 867 \r
 868                 temp1 = Am * dUi1\r
 869                 temp2 = Ap * dUi2\r
 870 \r
 871                 Rhs[j * nbComp + 0] = -temp1[0] + temp2[0] - (Uk[j * nbComp + 0] - Un[j * nbComp + 0])\r
 872                 Rhs[j * nbComp + 1] = -temp1[1] + temp2[1] - (Uk[j * nbComp + 1] - Un[j * nbComp + 1])\r
 873                 Rhs[j * nbComp + 2] = -temp1[2] + temp2[2] - (Uk[j * nbComp + 2] - Un[j * nbComp + 2])\r
 874         else:\r
 875                 Rhs[j * nbComp + 0] = 0\r
 876                 Rhs[j * nbComp + 1] = 0\r
 877                 Rhs[j * nbComp + 2] = 0\r
 878 \r
 879 \r
 880 def computeSystemMatrix(a,b,nx, cfl, Uk, isImplicit):\r
 881         dim = 1\r
 882         nbComp = 3\r
 883         dx = (b - a) / nx\r
 884         dt = cfl * dx / c0\r
 885         nbVoisinsMax = 2\r
 886 \r
 887         nbCells = nx\r
 888         divMat = cdmath.SparseMatrixPetsc(nbCells * nbComp, nbCells * nbComp, (nbVoisinsMax + 1) * nbComp)\r
 889 \r
 890         divMat.zeroEntries()  #sets the matrix coefficients to zero\r
 891         for j in range(nbCells):\r
 892                 \r
 893                 # traitements des bords\r
 894                 if (j == 0):\r
 895                         FillMatrixFromEdges(j, Uk, nbComp, divMat, dt, dx)\r
 896                 elif (j == nbCells - 1):\r
 897                         FillMatrixFromEdges(j, Uk, nbComp, divMat, dt, dx)\r
 898                 # traitement des cellules internes\r
 899                 else:\r
 900                         FillMatrixFromInnerCell(j, Uk, nbComp, divMat, dt, dx)\r
 901         \r
 902         if(isImplicit):         \r
 903                 divMat.diagonalShift(1.)  # add one on the diagonal\r
 904         else:\r
 905                 divMat*=-1.\r
 906                 divMat.diagonalShift(1.)  # add one on the diagonal\r
 907 \r
 908         return divMat\r
 909 \r
 910 def computeRHSVector(a,b,nx, cfl, Uk, Un, isImplicit):\r
 911         dim = 1\r
 912         nbComp = 3\r
 913         dx = (b - a) / nx\r
 914         dt = cfl * dx / c0\r
 915         nbVoisinsMax = 2\r
 916 \r
 917         nbCells = nx\r
 918         Rhs = cdmath.Vector(nbCells * (nbComp))\r
 919 \r
 920         for j in range(nbCells):\r
 921                 \r
 922                 # traitements des bords\r
 923                 if (j == 0):\r
 924                         FillRHSFromEdges(j, Uk, nbComp, Rhs, Un, dt, dx, isImplicit)\r
 925                 elif (j == nbCells - 1):\r
 926                         FillRHSFromEdges(j, Uk, nbComp, Rhs, Un, dt, dx, isImplicit)\r
 927                 # traitement des cellules internes\r
 928                 else:\r
 929                         FillRHSFromInnerCell(j, Uk, nbComp, Rhs, Un, dt, dx, isImplicit)\r
 930                         \r
 931         return Rhs\r
 932 \r
 933 \r
 934 if __name__ == """__main__""":\r
 935         nbComp=3 # number of equations \r
 936         a = 0.# domain is interval [a,b]\r
 937         b = 4.2# domain is interval [a,b]\r
 938         nx = 10# number of cells\r
 939         dx = (b - a) / nx  # space step\r
 940         x = [a + 0.5 * dx + i * dx for i in range(nx)]  # array of cell center (1D mesh)\r
 941         state_law = "Hermite575K"\r
 942         state_law_parameters(state_law)\r
 943         rho0=p_to_rho_StiffenedGaz(p_0, T_0)\r
 944         rhoE0=T_to_rhoE_StiffenedGaz(T_0, rho0, rho0*v_0)\r
 945 \r
 946 \r
 947 #### initial condition (T in K, v in m/s, p in Pa)\r
 948         p_initial   = np.array([ p_outlet      for xi in x])\r
 949         v_initial   = np.array([ v_inlet       for xi in x])\r
 950         T_initial   = np.array([ T_inlet       for xi in x])\r
 951         \r
 952         rho_field = p_to_rho_StiffenedGaz(p_initial, T_initial)\r
 953         q_field = rho_field * v_initial\r
 954         rho_E_field = rho_field * T_to_E_StiffenedGaz(p_initial, T_initial, v_initial)\r
 955 \r
 956         U = cdmath.Vector(nx * (nbComp))#Inutile à terme mais nécessaire pour le moment\r
 957 \r
 958         for k in range(nx):\r
 959                 U[k * nbComp + 0] = rho_field[k]\r
 960                 U[k * nbComp + 1] = q_field[k]\r
 961                 U[k * nbComp + 2] = rho_E_field[k]\r
 962         print("\n Testing function computeSystemMatrix \n")\r
 963         cfl = 0.5\r
 964         computeSystemMatrix(a, b, nx, cfl, U,True)  #Implicit matrix\r
 965         computeSystemMatrix(a, b, nx, cfl, U,False) #Explicit matrix\r
 966 \r
 967         print("\n Testing function computeRHSVector \n")\r
 968         cfl = 0.5\r
 969         computeRHSVector(a, b, nx, cfl, U, U,True)  #Implicit RHS\r
 970         computeRHSVector(a, b, nx, cfl, U, U,False) #Explicit RHS\r
 971 \r
 972         print("\n Testing function solve (Implicit scheme) \n")\r
 973         isImplicit=True\r
 974         cfl = 1000.\r
 975         solve(a, b, nx, "RegularSquares", "", cfl, state_law, isImplicit)\r
 976         \r
 977         print("\n Testing function solve (Explicit scheme) \n")\r
 978         isImplicit=False\r
 979         cfl = 0.5\r
 980         solve(a, b, nx, "RegularSquares", "", cfl, state_law, isImplicit)\r
 981 \r