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

   1 #!/usr/bin/env python3\r
   2 # -*-coding:utf-8 -*\r
   3 \r
   4 #===============================================================================================================================\r
   5 # Name        : Résolution VF des équations d'Euler isotherme 2D sans terme source\r
   6 #                \partial_t rho + \div q = 0\r
   7 #                \partial_t q   + \div q \otimes q/rho   \grad p = 0\r
   8 # Author      : Michaël Ndjinga, Coraline Mounier\r
   9 # Copyright   : CEA Saclay 2020\r
  10 # Description : Propagation d'une onde de choc droite\r
  11 #               Utilisation du schéma upwind explicite ou implicite sur un maillage général\r
  12 #               Initialisation par une discontinuité verticale\r
  13 #               Conditions aux limites de Neumann\r
  14 #                Création et sauvegarde du champ résultant et des figures\r
  15 #================================================================================================================================\r
  16 \r
  17 \r
  18 import cdmath\r
  19 import numpy as np\r
  20 import matplotlib\r
  21 matplotlib.use("Agg")\r
  22 import matplotlib.pyplot as plt\r
  23 import matplotlib.animation as manimation\r
  24 from math import sqrt\r
  25 import sys\r
  26 import PV_routines\r
  27 \r
  28 p0=155.e5 #reference pressure in a pressurised nuclear vessel\r
  29 c0=700.   #reference sound speed for water at 155 bars\r
  30 rho0=p0/(c0*c0)   #reference density\r
  31 precision=1e-5\r
  32 \r
  33 def initial_conditions_Riemann_problem(my_mesh):\r
  34     print( "Initial data : Riemann problem" )\r
  35     dim     = my_mesh.getMeshDimension()\r
  36     nbCells = my_mesh.getNumberOfCells()\r
  37 \r
  38     xcentre = 0.5\r
  39 \r
  40     density_field = cdmath.Field("Density",    cdmath.CELLS, my_mesh, 1)\r
  41     q_x_field     = cdmath.Field("Momentum x", cdmath.CELLS, my_mesh, 1)\r
  42     q_y_field     = cdmath.Field("Momentum y", cdmath.CELLS, my_mesh, 1)\r
  43     #Velocity field with 3 components should be created for streamlines to be drawn\r
  44     velocity_field = cdmath.Field("Velocity", cdmath.CELLS, my_mesh, 3)\r
  45 \r
  46     for i in range(nbCells):\r
  47         x = my_mesh.getCell(i).x()\r
  48 \r
  49         # Initial momentum is zero\r
  50         q_x_field[i] = 0\r
  51         q_y_field[i] = 0\r
  52 \r
  53         # Initial velocity is zero\r
  54         velocity_field[i,0] = 0\r
  55         velocity_field[i,1] = 0\r
  56         velocity_field[i,2] = 0\r
  57 \r
  58         if x < xcentre:\r
  59             density_field[i] = rho0\r
  60             pass\r
  61         else:\r
  62             density_field[i] = rho0/2\r
  63             pass\r
  64         pass\r
  65 \r
  66     return density_field, q_x_field, q_y_field, velocity_field\r
  67 \r
  68 \r
  69 def jacobianMatricesm_x(coeff,rho_l,q_lx,q_ly,rho_r,q_rx,q_ry):\r
  70     RoeMatx = cdmath.Matrix(3,3);\r
  71     Dmacx = cdmath.Matrix(3,3);\r
  72 \r
  73     u_lx=q_lx/rho_l\r
  74     u_rx=q_rx/rho_r\r
  75     u_ly=q_ly/rho_l\r
  76     u_ry=q_ry/rho_r\r
  77     if rho_l<0 or rho_r<0:\r
  78         print("rho_l=",rho_l, " rho_r= ",rho_r)\r
  79         raise ValueError("Negative density")\r
  80     ux = (u_lx*sqrt(rho_l)+u_rx*sqrt(rho_r))/(sqrt(rho_l)+sqrt(rho_r));\r
  81     uy = (u_ly*sqrt(rho_l)+u_ry*sqrt(rho_r))/(sqrt(rho_l)+sqrt(rho_r));\r
  82 \r
  83     RoeMatx[0,0] = 0\r
  84     RoeMatx[0,1] = 1\r
  85     RoeMatx[0,2] = 0\r
  86     RoeMatx[1,0] = c0*c0 - ux*ux\r
  87     RoeMatx[1,1] = 2*ux\r
  88     RoeMatx[1,2] = 0\r
  89     RoeMatx[2,0] = -ux*uy\r
  90     RoeMatx[2,1] = uy\r
  91     RoeMatx[2,2] = ux\r
  92 \r
  93     Dmacx[0,0] = abs(ux)-ux\r
  94     Dmacx[0,1] = 1\r
  95     Dmacx[0,2] = 0\r
  96     Dmacx[1,0] = -c0*c0-ux*ux\r
  97     Dmacx[1,1] = abs(ux)+ux\r
  98     Dmacx[1,2] = 0\r
  99     Dmacx[2,0] = -ux*uy\r
 100     Dmacx[2,1] = uy\r
 101     Dmacx[2,2] = abs(ux)\r
 102 \r
 103     return (RoeMatx-Dmacx)*coeff*0.5\r
 104 \r
 105 def jacobianMatricesp_x(coeff,rho_l,q_lx,q_ly,rho_r,q_rx,q_ry):\r
 106     RoeMatx   = cdmath.Matrix(3,3)\r
 107     Dmacx = cdmath.Matrix(3,3)\r
 108 \r
 109     u_lx=q_lx/rho_l\r
 110     u_rx=q_rx/rho_r\r
 111     u_ly=q_ly/rho_l\r
 112     u_ry=q_ry/rho_r\r
 113     if rho_l<0 or rho_r<0:\r
 114         print("rho_l=",rho_l, " rho_r= ",rho_r)\r
 115         raise ValueError("Negative density")\r
 116     ux = (u_lx*sqrt(rho_l)+u_rx*sqrt(rho_r))/(sqrt(rho_l)+sqrt(rho_r));\r
 117     uy = (u_ly*sqrt(rho_l)+u_ry*sqrt(rho_r))/(sqrt(rho_l)+sqrt(rho_r));\r
 118 \r
 119     RoeMatx[0,0] = 0\r
 120     RoeMatx[0,1] = 1\r
 121     RoeMatx[0,2] = 0\r
 122     RoeMatx[1,0] = c0*c0 - ux*ux\r
 123     RoeMatx[1,1] = 2*ux\r
 124     RoeMatx[1,2] = 0\r
 125     RoeMatx[2,0] = -ux*uy\r
 126     RoeMatx[2,1] = uy\r
 127     RoeMatx[2,2] = ux\r
 128 \r
 129     Dmacx[0,0] = abs(ux)-ux\r
 130     Dmacx[0,1] = 1\r
 131     Dmacx[0,2] = 0\r
 132     Dmacx[1,0] = -c0*c0-ux*ux\r
 133     Dmacx[1,1] = abs(ux)+ux\r
 134     Dmacx[1,2] = 0\r
 135     Dmacx[2,0] = -ux*uy\r
 136     Dmacx[2,1] = uy\r
 137     Dmacx[2,2] = abs(ux)\r
 138 \r
 139     return (RoeMatx+Dmacx)*coeff*0.5\r
 140 \r
 141 def jacobianMatricesm_y(coeff,rho_l,q_lx,q_ly,rho_r,q_rx,q_ry):\r
 142     RoeMaty   = cdmath.Matrix(3,3);\r
 143     Dmacy = cdmath.Matrix(3,3);\r
 144 \r
 145     u_lx=q_lx/rho_l\r
 146     u_rx=q_rx/rho_r\r
 147     u_ly=q_ly/rho_l\r
 148     u_ry=q_ry/rho_r\r
 149     if rho_l<0 or rho_r<0:\r
 150         print("rho_l=",rho_l, " rho_r= ",rho_r)\r
 151         raise ValueError("Negative density")\r
 152     ux = (u_lx*sqrt(rho_l)+u_rx*sqrt(rho_r))/(sqrt(rho_l)+sqrt(rho_r));\r
 153     uy = (u_ly*sqrt(rho_l)+u_ry*sqrt(rho_r))/(sqrt(rho_l)+sqrt(rho_r));\r
 154 \r
 155     RoeMaty[0,0] = 0\r
 156     RoeMaty[0,1] = 0\r
 157     RoeMaty[0,2] = 1\r
 158     RoeMaty[1,0] = -ux*uy\r
 159     RoeMaty[1,1] = uy\r
 160     RoeMaty[1,2] = ux\r
 161     RoeMaty[2,0] = c0*c0-uy*uy\r
 162     RoeMaty[2,1] = 0\r
 163     RoeMaty[2,2] = 2*uy\r
 164 \r
 165     Dmacy[0,0] = abs(uy)-uy\r
 166     Dmacy[0,1] = 0\r
 167     Dmacy[0,2] = 1\r
 168     Dmacy[1,0] = -ux*uy\r
 169     Dmacy[1,1] = abs(uy)\r
 170     Dmacy[1,2] = ux\r
 171     Dmacy[2,0] = -c0*c0-uy*uy\r
 172     Dmacy[2,1] = 0\r
 173     Dmacy[2,2] = abs(uy)+uy\r
 174 \r
 175     return (RoeMaty-Dmacy)*coeff*0.5\r
 176 \r
 177 def jacobianMatricesp_y(coeff,rho_l,q_lx,q_ly,rho_r,q_rx,q_ry):\r
 178     RoeMaty   = cdmath.Matrix(3,3);\r
 179     Dmacy = cdmath.Matrix(3,3);\r
 180 \r
 181     u_lx=q_lx/rho_l\r
 182     u_rx=q_rx/rho_r\r
 183     u_ly=q_ly/rho_l\r
 184     u_ry=q_ry/rho_r\r
 185     if rho_l<0 or rho_r<0:\r
 186         print("rho_l=",rho_l, " rho_r= ",rho_r)\r
 187         raise ValueError("Negative density")\r
 188     ux = (u_lx*sqrt(rho_l)+u_rx*sqrt(rho_r))/(sqrt(rho_l)+sqrt(rho_r));\r
 189     uy = (u_ly*sqrt(rho_l)+u_ry*sqrt(rho_r))/(sqrt(rho_l)+sqrt(rho_r));\r
 190 \r
 191     RoeMaty[0,0] = 0\r
 192     RoeMaty[0,1] = 0\r
 193     RoeMaty[0,2] = 1\r
 194     RoeMaty[1,0] = -ux*uy\r
 195     RoeMaty[1,1] = uy\r
 196     RoeMaty[1,2] = ux\r
 197     RoeMaty[2,0] = c0*c0-uy*uy\r
 198     RoeMaty[2,1] = 0\r
 199     RoeMaty[2,2] = 2*uy\r
 200 \r
 201     Dmacy[0,0] = abs(uy)-uy\r
 202     Dmacy[0,1] = 0\r
 203     Dmacy[0,2] = 1\r
 204     Dmacy[1,0] = -ux*uy\r
 205     Dmacy[1,1] = abs(uy)\r
 206     Dmacy[1,2] = ux\r
 207     Dmacy[2,0] = -c0*c0-uy*uy\r
 208     Dmacy[2,1] = 0\r
 209     Dmacy[2,2] = abs(uy)+uy\r
 210 \r
 211     return (RoeMaty+Dmacy)*coeff*0.5\r
 212 \r
 213 def EulerSystemStaggered(ntmax, tmax, cfl,output_freq, my_mesh, meshName):\r
 214 \r
 215     if not my_mesh.isStructured() :\r
 216         raise ValueError("Euler_ConservativeStaggered2D_RiemannProblem.py expects a structured mesh")\r
 217 \r
 218     dim=my_mesh.getMeshDimension()\r
 219     if dim != 2 :\r
 220         raise ValueError("Euler_ConservativeStaggered2D_RiemannProblem.py expects a 2D mesh")\r
 221 \r
 222     nbComp=dim+1\r
 223     # Mesh parameters\r
 224     nbCells = my_mesh.getNumberOfCells()\r
 225     nbCells_x = my_mesh.getNx()\r
 226     nbCells_y = my_mesh.getNy()\r
 227     dx,dy = my_mesh.getDXYZ()\r
 228     nbVoisinsMax=my_mesh.getMaxNbNeighbours(cdmath.CELLS)\r
 229     dx_min=my_mesh.minRatioVolSurf()\r
 230     \r
 231     # Time evolution parameters\r
 232     time = 0.\r
 233     it=0;\r
 234     isStationary=False\r
 235     iterGMRESMax = 50\r
 236     newton_max   = 50\r
 237 \r
 238     #iteration vectors\r
 239     Un =cdmath.Vector(nbCells*nbComp)\r
 240     dUn=cdmath.Vector(nbCells*nbComp)\r
 241     dUk=cdmath.Vector(nbCells*nbComp)\r
 242     Rhs=cdmath.Vector(nbCells*nbComp)\r
 243 \r
 244     dUi_x=cdmath.Vector(nbComp)\r
 245     dUi_y=cdmath.Vector(nbComp)\r
 246     dUi1_x=cdmath.Vector(nbComp)\r
 247     dUi2_x=cdmath.Vector(nbComp)\r
 248     dUi1_y=cdmath.Vector(nbComp)\r
 249     dUi2_y=cdmath.Vector(nbComp)\r
 250 \r
 251     # Initial conditions #\r
 252     print("Construction of the initial condition …")\r
 253     rho_field, q_x_field, q_y_field, velocity_field= initial_conditions_Riemann_problem(my_mesh)\r
 254 \r
 255     for i in range(nbCells):\r
 256         Un[nbComp*i]   = rho_field[i]\r
 257         Un[nbComp*i+1] = q_x_field[i]\r
 258         Un[nbComp*i+2] = q_y_field[i]\r
 259 \r
 260     #Sauvegarde de la donnée initiale\r
 261     rho_field.setTime(time,it);\r
 262     rho_field.writeVTK("EulerIsothermal_"+str(dim)+"DConservativeStaggered_"+meshName+"_density");\r
 263     q_x_field.setTime(time,it);\r
 264     q_x_field.writeVTK("EulerIsothermal_"+str(dim)+"DConservativeStaggered_"+meshName+"_momentumX");\r
 265     q_y_field.setTime(time,it);\r
 266     q_y_field.writeVTK("EulerIsothermal_"+str(dim)+"DConservativeStaggered_"+meshName+"_momentumY");\r
 267     velocity_field.setTime(time,it);\r
 268     velocity_field.writeVTK("EulerIsothermal_"+str(dim)+"DConservativeStaggered_"+meshName+"_velocity");\r
 269 \r
 270     print("Starting computation of the isothermal Euler system with a conservative staggered scheme …")\r
 271     divMat=cdmath.SparseMatrixPetsc(nbCells*nbComp,nbCells*nbComp,(nbVoisinsMax+1)*nbComp)\r
 272 \r
 273     # Starting time loop\r
 274     while (it<ntmax and time <= tmax and not isStationary):\r
 275         dUn = Un.deepCopy()\r
 276         Uk  = Un.deepCopy()\r
 277         residu = 1e10\r
 278         k=0\r
 279         while (k<newton_max and residu > 1/precision ):\r
 280 \r
 281             dt = cfl * dx_min / c0# This choice should be improved when possible by using the actual eigenvalue abs(u)+c0, that is the time step should be determined avec the computation of the jacobian matrices\r
 282             #DEBUT BOUCLE NEWTON\r
 283             for j in range(nbCells_y):\r
 284                 for i in range(nbCells_x):\r
 285                     #Traitement des coins\r
 286                     if ( j==0 and i==0) :\r
 287                         #Calcul de Am_x\r
 288                         rho_l=Uk[3*nbCells_x*j+3*i]\r
 289                         qx_l =Uk[3*nbCells_x*j+3*i+1]\r
 290                         qy_l =Uk[3*nbCells_x*j+3*i+2]\r
 291                         rho_r=Uk[3*nbCells_x*j+3*(i+1)]\r
 292                         qx_r =Uk[3*nbCells_x*j+3*(i+1)+1]\r
 293                         qy_r =Uk[3*nbCells_x*j+3*(i+1)+2]\r
 294                         Am_x= jacobianMatricesm_x(dt/dx,rho_l,qx_l,qy_l,rho_r,qx_r,qy_r)\r
 295 \r
 296                         #calcul de Am_y\r
 297                         rho_l=Uk[3*nbCells_x*j+3*i]\r
 298                         qx_l =Uk[3*nbCells_x*j+3*i+1]\r
 299                         qy_l =Uk[3*nbCells_x*j+3*i+2]\r
 300                         rho_r=Uk[3*nbCells_x*(j+1)+3*i]\r
 301                         qx_r =Uk[3*nbCells_x*(j+1)+3*i+1]\r
 302                         qy_r =Uk[3*nbCells_x*(j+1)+3*i+2]\r
 303                         Am_y= jacobianMatricesm_y(dt/dx,rho_l,qx_l,qy_l,rho_r,qx_r,qy_r)\r
 304 \r
 305                         divMat.addValue(3*nbCells_x*j+3*i,3*nbCells_x*j+3*(i+1),Am_x)\r
 306                         divMat.addValue(3*nbCells_x*j+3*i,3*nbCells_x*(j+1)+3*i,Am_y)\r
 307                         divMat.addValue(3*nbCells_x*j+3*i,3*nbCells_x*j+3*i,Am_x*(-1.)-Am_y)\r
 308 \r
 309                         dUi_x[0] = Uk[3*nbCells_x*j+3*(i+1)]-Uk[3*nbCells_x*j+3*i]\r
 310                         dUi_x[1] = Uk[3*nbCells_x*j+3*(i+1)+1]-Uk[3*nbCells_x*j+3*i+1]\r
 311                         dUi_x[2] = Uk[3*nbCells_x*j+3*(i+1)+2]-Uk[3*nbCells_x*j+3*i+2]\r
 312                         dUi_y[0] = Uk[3*nbCells_x*(j+1)+3*i]-Uk[3*nbCells_x*j+3*i]\r
 313                         dUi_y[1] = Uk[3*nbCells_x*(j+1)+3*i+1]-Uk[3*nbCells_x*j+3*i+1]\r
 314                         dUi_y[2] = Uk[3*nbCells_x*(j+1)+3*i+2]-Uk[3*nbCells_x*j+3*i+2]\r
 315 \r
 316                         temp_x = Am_x*dUi_x\r
 317                         temp_y = Am_y*dUi_y\r
 318                         #print("Bloc 0 matrix  :   ", Am)\r
 319                         Rhs[3*nbCells_x*j+3*i+0] = -temp_x[0] - temp_y[0]-(Uk[3*nbCells_x*j+3*i+0]-Un[3*nbCells_x*j+3*i+0])\r
 320                         Rhs[3*nbCells_x*j+3*i+1] = -temp_x[1] - temp_y[1]-(Uk[3*nbCells_x*j+3*i+1]-Un[3*nbCells_x*j+3*i+1])\r
 321                         Rhs[3*nbCells_x*j+3*i+2] = -temp_x[2] - temp_y[2]-(Uk[3*nbCells_x*j+3*i+2]-Un[3*nbCells_x*j+3*i+2])\r
 322 \r
 323                     elif(i==0 and j==nbCells_y-1):\r
 324                         #Calcul de Am_x\r
 325                         rho_l=Uk[3*nbCells_x*j+3*i]\r
 326                         qx_l =Uk[3*nbCells_x*j+3*i+1]\r
 327                         qy_l =Uk[3*nbCells_x*j+3*i+2]\r
 328                         rho_r=Uk[3*nbCells_x*j+3*(i+1)]\r
 329                         qx_r =Uk[3*nbCells_x*j+3*(i+1)+1]\r
 330                         qy_r =Uk[3*nbCells_x*j+3*(i+1)+2]\r
 331                         Am_x= jacobianMatricesm_x(dt/dx,rho_l,qx_l,qy_l,rho_r,qx_r,qy_r)\r
 332 \r
 333                         #calcul de Ap_y\r
 334                         rho_l=Uk[3*nbCells_x*(j-1)+3*i]\r
 335                         qx_l =Uk[3*nbCells_x*(j-1)+3*i+1]\r
 336                         qy_l =Uk[3*nbCells_x*(j-1)+3*i+2]\r
 337                         rho_r=Uk[3*nbCells_x*j+3*i]\r
 338                         qx_r =Uk[3*nbCells_x*j+3*i+1]\r
 339                         qy_r =Uk[3*nbCells_x*j+3*i+2]\r
 340                         Ap_y= jacobianMatricesp_y(dt/dx,rho_l,qx_l,qy_l,rho_r,qx_r,qy_r)\r
 341 \r
 342                         divMat.addValue(3*nbCells_x*j+3*i,3*nbCells_x*j+3*(i+1),Am_x)\r
 343                         divMat.addValue(3*nbCells_x*j+3*i,3*nbCells_x*(j-1)+3*i,Ap_y*(-1.))\r
 344                         divMat.addValue(3*nbCells_x*j+3*i,3*nbCells_x*j+3*i,Am_x*(-1.)+Ap_y)\r
 345 \r
 346                         dUi_x[0] = Uk[3*nbCells_x*j+3*(i+1)]-Uk[3*nbCells_x*j+3*i]\r
 347                         dUi_x[1] = Uk[3*nbCells_x*j+3*(i+1)+1]-Uk[3*nbCells_x*j+3*i+1]\r
 348                         dUi_x[2] = Uk[3*nbCells_x*j+3*(i+1)+2]-Uk[3*nbCells_x*j+3*i+2]\r
 349                         dUi_y[0] = Uk[3*nbCells_x*j+3*i]-Uk[3*nbCells_x*(j-1)+3*i]\r
 350                         dUi_y[1] = Uk[3*nbCells_x*j+3*i+1]-Uk[3*nbCells_x*(j-1)+3*i+1]\r
 351                         dUi_y[2] = Uk[3*nbCells_x*j+3*i+2]-Uk[3*nbCells_x*(j-1)+3*i+2]\r
 352 \r
 353                         temp_x = Am_x*dUi_x\r
 354                         temp_y = Ap_y*dUi_y\r
 355                         #print("Bloc 0 matrix  :   ", Am)\r
 356                         Rhs[3*nbCells_x*j+3*i+0] = -temp_x[0]- temp_y[0]-(Uk[3*nbCells_x*j+3*i+0]-Un[3*nbCells_x*j+3*i+0])\r
 357                         Rhs[3*nbCells_x*j+3*i+1] = -temp_x[1]- temp_y[1]-(Uk[3*nbCells_x*j+3*i+1]-Un[3*nbCells_x*j+3*i+1])\r
 358                         Rhs[3*nbCells_x*j+3*i+2] = -temp_x[2]- temp_y[2]-(Uk[3*nbCells_x*j+3*i+2]-Un[3*nbCells_x*j+3*i+2])\r
 359 \r
 360                     elif(i==nbCells_x-1 and j==0):\r
 361 \r
 362                         #Calcul de Ap_x\r
 363                         rho_l=Uk[3*nbCells_x*j+3*(i-1)]\r
 364                         qx_l =Uk[3*nbCells_x*j+3*(i-1)+1]\r
 365                         qy_l =Uk[3*nbCells_x*j+3*(i-1)+2]\r
 366                         rho_r=Uk[3*nbCells_x*j+3*i]\r
 367                         qx_r =Uk[3*nbCells_x*j+3*i+1]\r
 368                         qy_r =Uk[3*nbCells_x*j+3*i+2]\r
 369                         Ap_x= jacobianMatricesp_x(dt/dx,rho_l,qx_l,qy_l,rho_r,qx_r,qy_r)\r
 370 \r
 371                         #calcul de Am_y\r
 372                         rho_l=Uk[3*nbCells_x*j+3*i]\r
 373                         qx_l =Uk[3*nbCells_x*j+3*i+1]\r
 374                         qy_l =Uk[3*nbCells_x*j+3*i+2]\r
 375                         rho_r=Uk[3*nbCells_x*(j+1)+3*i]\r
 376                         qx_r =Uk[3*nbCells_x*(j+1)+3*i+1]\r
 377                         qy_r =Uk[3*nbCells_x*(j+1)+3*i+2]\r
 378                         Am_y= jacobianMatricesm_y(dt/dx,rho_l,qx_l,qy_l,rho_r,qx_r,qy_r)\r
 379 \r
 380                         divMat.addValue(3*nbCells_x*j+3*i,3*nbCells_x*j+3*(i-1),Ap_x*(-1))\r
 381                         divMat.addValue(3*nbCells_x*j+3*i,3*nbCells_x*(j+1)+3*i,Am_y)\r
 382                         divMat.addValue(3*nbCells_x*j+3*i,3*nbCells_x*j+3*i,Ap_x-Am_y)\r
 383 \r
 384                         dUi_x[0] = Uk[3*nbCells_x*j+3*i]-Uk[3*nbCells_x*j+3*(i-1)]\r
 385                         dUi_x[1] = Uk[3*nbCells_x*j+3*i+1]-Uk[3*nbCells_x*j+3*(i-1)+1]\r
 386                         dUi_x[2] = Uk[3*nbCells_x*j+3*i+2]-Uk[3*nbCells_x*j+3*(i-1)+2]\r
 387                         dUi_y[0] = Uk[3*nbCells_x*(j+1)+3*i]-Uk[3*nbCells_x*j+3*i]\r
 388                         dUi_y[1] = Uk[3*nbCells_x*(j+1)+3*i+1]-Uk[3*nbCells_x*j+3*i+1]\r
 389                         dUi_y[2] = Uk[3*nbCells_x*(j+1)+3*i+2]-Uk[3*nbCells_x*j+3*i+2]\r
 390 \r
 391                         temp_x = Ap_x*dUi_x\r
 392                         temp_y= Am_y*dUi_y\r
 393                         #print("Bloc 0 matrix  :   ", Am)\r
 394                         Rhs[3*nbCells_x*j+3*i+0] = -temp_x[0]- temp_y[0]-(Uk[3*nbCells_x*j+3*i+0]-Un[3*nbCells_x*j+3*i+0])\r
 395                         Rhs[3*nbCells_x*j+3*i+1] = -temp_x[1]- temp_y[1]-(Uk[3*nbCells_x*j+3*i+1]-Un[3*nbCells_x*j+3*i+1])\r
 396                         Rhs[3*nbCells_x*j+3*i+2] = -temp_x[2]- temp_y[2]-(Uk[3*nbCells_x*j+3*i+2]-Un[3*nbCells_x*j+3*i+2])\r
 397 \r
 398                     elif ( j==nbCells_y-1 and i==nbCells_x-1) :\r
 399 \r
 400                         #Calcul de Ap_x\r
 401                         rho_l=Uk[3*nbCells_x*j+3*(i-1)]\r
 402                         qx_l =Uk[3*nbCells_x*j+3*(i-1)+1]\r
 403                         qy_l =Uk[3*nbCells_x*j+3*(i-1)+2]\r
 404                         rho_r=Uk[3*nbCells_x*j+3*i]\r
 405                         qx_r =Uk[3*nbCells_x*j+3*i+1]\r
 406                         qy_r =Uk[3*nbCells_x*j+3*i+2]\r
 407                         Ap_x= jacobianMatricesp_x(dt/dx,rho_l,qx_l,qy_l,rho_r,qx_r,qy_r)\r
 408 \r
 409                         #calcul de Ap_y\r
 410                         rho_l=Uk[3*nbCells_x*(j-1)+3*i]\r
 411                         qx_l =Uk[3*nbCells_x*(j-1)+3*i+1]\r
 412                         qy_l =Uk[3*nbCells_x*(j-1)+3*i+2]\r
 413                         rho_r=Uk[3*nbCells_x*j+3*i]\r
 414                         qx_r =Uk[3*nbCells_x*j+3*i+1]\r
 415                         qy_r =Uk[3*nbCells_x*j+3*i+2]\r
 416                         Ap_y= jacobianMatricesp_y(dt/dx,rho_l,qx_l,qy_l,rho_r,qx_r,qy_r)\r
 417 \r
 418                         divMat.addValue(3*nbCells_x*j+3*i,3*nbCells_x*j+3*(i-1),Ap_x*(-1.))\r
 419                         divMat.addValue(3*nbCells_x*j+3*i,3*nbCells_x*(j-1)+3*i,Ap_y*(-1.))\r
 420                         divMat.addValue(3*nbCells_x*j+3*i,3*nbCells_x*j+3*i,Ap_x+Ap_y)\r
 421 \r
 422                         dUi_x[0] = Uk[3*nbCells_x*j+3*i]-Uk[3*nbCells_x*j+3*(i-1)]\r
 423                         dUi_x[1] = Uk[3*nbCells_x*j+3*i+1]-Uk[3*nbCells_x*j+3*(i-1)+1]\r
 424                         dUi_x[2] = Uk[3*nbCells_x*j+3*i+2]-Uk[3*nbCells_x*j+3*(i-1)+2]\r
 425                         dUi_y[0] = Uk[3*nbCells_x*j+3*i]-Uk[3*nbCells_x*(j-1)+3*i]\r
 426                         dUi_y[1] = Uk[3*nbCells_x*j+3*i+1]-Uk[3*nbCells_x*(j-1)+3*i+1]\r
 427                         dUi_y[2] = Uk[3*nbCells_x*j+3*i+2]-Uk[3*nbCells_x*(j-1)+3*i+2]\r
 428 \r
 429                         temp_x = Ap_x*dUi_x\r
 430                         temp_y = Ap_y*dUi_y\r
 431                         Rhs[3*nbCells_x*j+3*i+0] = -temp_x[0]-temp_y[0]-(Uk[3*nbCells_x*j+3*i+0]-Un[3*nbCells_x*j+3*i+0])\r
 432                         Rhs[3*nbCells_x*j+3*i+1] = -temp_x[1]-temp_y[1]-(Uk[3*nbCells_x*j+3*i+1]-Un[3*nbCells_x*j+3*i+1])\r
 433                         Rhs[3*nbCells_x*j+3*i+2] = -temp_x[2]-temp_y[2]-(Uk[3*nbCells_x*j+3*i+2]-Un[3*nbCells_x*j+3*i+2])\r
 434 \r
 435                     #Traitement des bords restants\r
 436                     elif i==0 :\r
 437 \r
 438                         #Calcul de Am_x\r
 439                         rho_l=Uk[3*nbCells_x*j+3*i]\r
 440                         qx_l =Uk[3*nbCells_x*j+3*i+1]\r
 441                         qy_l =Uk[3*nbCells_x*j+3*i+2]\r
 442                         rho_r=Uk[3*nbCells_x*j+3*(i+1)]\r
 443                         qx_r =Uk[3*nbCells_x*j+3*(i+1)+1]\r
 444                         qy_r =Uk[3*nbCells_x*j+3*(i+1)+2]\r
 445                         Am_x= jacobianMatricesm_x(dt/dx,rho_l,qx_l,qy_l,rho_r,qx_r,qy_r)\r
 446 \r
 447                         #calcul de Am_y\r
 448                         rho_l=Uk[3*nbCells_x*j+3*i]\r
 449                         qx_l =Uk[3*nbCells_x*j+3*i+1]\r
 450                         qy_l =Uk[3*nbCells_x*j+3*i+2]\r
 451                         rho_r=Uk[3*nbCells_x*(j+1)+3*i]\r
 452                         qx_r =Uk[3*nbCells_x*(j+1)+3*i+1]\r
 453                         qy_r =Uk[3*nbCells_x*(j+1)+3*i+2]\r
 454                         Am_y= jacobianMatricesm_y(dt/dx,rho_l,qx_l,qy_l,rho_r,qx_r,qy_r)\r
 455 \r
 456                         #calcul de Ap_y\r
 457                         rho_l=Uk[3*nbCells_x*(j-1)+3*i]\r
 458                         qx_l =Uk[3*nbCells_x*(j-1)+3*i+1]\r
 459                         qy_l =Uk[3*nbCells_x*(j-1)+3*i+2]\r
 460                         rho_r=Uk[3*nbCells_x*j+3*i]\r
 461                         qx_r =Uk[3*nbCells_x*j+3*i+1]\r
 462                         qy_r =Uk[3*nbCells_x*j+3*i+2]\r
 463                         Ap_y= jacobianMatricesp_y(dt/dx,rho_l,qx_l,qy_l,rho_r,qx_r,qy_r)\r
 464 \r
 465                         #remplissage de la divMat\r
 466                         divMat.addValue(3*nbCells_x*j+3*i,3*nbCells_x*j+3*(i+1),Am_x)\r
 467                         divMat.addValue(3*nbCells_x*j+3*i,3*nbCells_x*(j+1)+3*i,Am_y)\r
 468                         divMat.addValue(3*nbCells_x*j+3*i,3*nbCells_x*(j-1)+3*i,Ap_y*(-1.))\r
 469                         divMat.addValue(3*nbCells_x*j+3*i,3*nbCells_x*j+3*i,Ap_y-Am_x-Am_y)\r
 470 \r
 471                         #remplissage du membre de droite\r
 472                         dUi_x[0] = Uk[3*nbCells_x*j+3*(i+1)]-Uk[3*nbCells_x*j+3*i]\r
 473                         dUi_x[1] = Uk[3*nbCells_x*j+3*(i+1)+1]-Uk[3*nbCells_x*j+3*i+1]\r
 474                         dUi_x[2] = Uk[3*nbCells_x*j+3*(i+1)+2]-Uk[3*nbCells_x*j+3*i+2]\r
 475                         dUi1_y[0] = Uk[3*nbCells_x*(j+1)+3*i]-Uk[3*nbCells_x*j+3*i]\r
 476                         dUi1_y[1] = Uk[3*nbCells_x*(j+1)+3*i+1]-Uk[3*nbCells_x*j+3*i+1]\r
 477                         dUi1_y[2] = Uk[3*nbCells_x*(j+1)+3*i+2]-Uk[3*nbCells_x*j+3*i+2]\r
 478                         dUi2_y[0] = Uk[3*nbCells_x*j+3*i]-Uk[3*nbCells_x*(j-1)+3*i]\r
 479                         dUi2_y[1] = Uk[3*nbCells_x*j+3*i+1]-Uk[3*nbCells_x*(j-1)+3*i+1]\r
 480                         dUi2_y[2] = Uk[3*nbCells_x*j+3*i+2]-Uk[3*nbCells_x*(j-1)+3*i+2]\r
 481 \r
 482                         temp_x  = Am_x*dUi_x\r
 483                         temp1_y = Am_y*dUi1_y\r
 484                         temp2_y = Ap_y*dUi2_y\r
 485                         Rhs[3*nbCells_x*j+3*i+0] = -temp_x[0]-temp1_y[0]-temp2_y[0]-(Uk[3*nbCells_x*j+3*i+0]-Un[3*nbCells_x*j+3*i+0])\r
 486                         Rhs[3*nbCells_x*j+3*i+1] = -temp_x[1]-temp1_y[1]-temp2_y[1]-(Uk[3*nbCells_x*j+3*i+1]-Un[3*nbCells_x*j+3*i+1])\r
 487                         Rhs[3*nbCells_x*j+3*i+2] = -temp_x[2]-temp1_y[2]-temp2_y[2]-(Uk[3*nbCells_x*j+3*i+2]-Un[3*nbCells_x*j+3*i+2])\r
 488 \r
 489                     elif i==nbCells_x-1:\r
 490 \r
 491                         #Calcul de Ap_x\r
 492                         rho_l=Uk[3*nbCells_x*j+3*(i-1)]\r
 493                         qx_l =Uk[3*nbCells_x*j+3*(i-1)+1]\r
 494                         qy_l =Uk[3*nbCells_x*j+3*(i-1)+2]\r
 495                         rho_r=Uk[3*nbCells_x*j+3*i]\r
 496                         qx_r =Uk[3*nbCells_x*j+3*i+1]\r
 497                         qy_r =Uk[3*nbCells_x*j+3*i+2]\r
 498                         Ap_x= jacobianMatricesp_x(dt/dx,rho_l,qx_l,qy_l,rho_r,qx_r,qy_r)\r
 499 \r
 500                         #calcul de Am_y\r
 501                         rho_l=Uk[3*nbCells_x*j+3*i]\r
 502                         qx_l =Uk[3*nbCells_x*j+3*i+1]\r
 503                         qy_l =Uk[3*nbCells_x*j+3*i+2]\r
 504                         rho_r=Uk[3*nbCells_x*(j+1)+3*i]\r
 505                         qx_r =Uk[3*nbCells_x*(j+1)+3*i+1]\r
 506                         qy_r =Uk[3*nbCells_x*(j+1)+3*i+2]\r
 507                         Am_y= jacobianMatricesm_y(dt/dx,rho_l,qx_l,qy_l,rho_r,qx_r,qy_r)\r
 508 \r
 509                         #calcul de Ap_y\r
 510                         rho_l=Uk[3*nbCells_x*(j-1)+3*i]\r
 511                         qx_l =Uk[3*nbCells_x*(j-1)+3*i+1]\r
 512                         qy_l =Uk[3*nbCells_x*(j-1)+3*i+2]\r
 513                         rho_r=Uk[3*nbCells_x*j+3*i]\r
 514                         qx_r =Uk[3*nbCells_x*j+3*i+1]\r
 515                         qy_r =Uk[3*nbCells_x*j+3*i+2]\r
 516                         Ap_y= jacobianMatricesp_y(dt/dx,rho_l,qx_l,qy_l,rho_r,qx_r,qy_r)\r
 517 \r
 518                         #remplissage de la divMat\r
 519                         divMat.addValue(3*nbCells_x*j+3*i,3*nbCells_x*j+3*(i-1),Ap_x*(-1.))\r
 520                         divMat.addValue(3*nbCells_x*j+3*i,3*nbCells_x*(j+1)+3*i,Am_y)\r
 521                         divMat.addValue(3*nbCells_x*j+3*i,3*nbCells_x*(j-1)+3*i,Ap_y*(-1.))\r
 522                         divMat.addValue(3*nbCells_x*j+3*i,3*nbCells_x*j+3*i,Ap_y+Ap_x-Am_y)\r
 523 \r
 524                         #remplissage du membre de droite\r
 525                         dUi_x[0] = Uk[3*nbCells_x*j+3*i]-Uk[3*nbCells_x*j+3*(i-1)]\r
 526                         dUi_x[1] = Uk[3*nbCells_x*j+3*i+1]-Uk[3*nbCells_x*j+3*(i-1)+1]\r
 527                         dUi_x[2] = Uk[3*nbCells_x*j+3*i+2]-Uk[3*nbCells_x*j+3*(i-1)+2]\r
 528                         dUi1_y[0] = Uk[3*nbCells_x*(j+1)+3*i]-Uk[3*nbCells_x*j+3*i]\r
 529                         dUi1_y[1] = Uk[3*nbCells_x*(j+1)+3*i+1]-Uk[3*nbCells_x*j+3*i+1]\r
 530                         dUi1_y[2] = Uk[3*nbCells_x*(j+1)+3*i+2]-Uk[3*nbCells_x*j+3*i+2]\r
 531                         dUi2_y[0] = Uk[3*nbCells_x*j+3*i]-Uk[3*nbCells_x*(j-1)+3*i]\r
 532                         dUi2_y[1] = Uk[3*nbCells_x*j+3*i+1]-Uk[3*nbCells_x*(j-1)+3*i+1]\r
 533                         dUi2_y[2] = Uk[3*nbCells_x*j+3*i+2]-Uk[3*nbCells_x*(j-1)+3*i+2]\r
 534 \r
 535                         temp_x  = Ap_x*dUi_x\r
 536                         temp1_y = Am_y*dUi1_y\r
 537                         temp2_y = Ap_y*dUi2_y\r
 538                         Rhs[3*nbCells_x*j+3*i+0] = -temp_x[0]-temp1_y[0]-temp2_y[0]-(Uk[3*nbCells_x*j+3*i+0]-Un[3*nbCells_x*j+3*i+0])\r
 539                         Rhs[3*nbCells_x*j+3*i+1] = -temp_x[1]-temp1_y[1]-temp2_y[1]-(Uk[3*nbCells_x*j+3*i+1]-Un[3*nbCells_x*j+3*i+1])\r
 540                         Rhs[3*nbCells_x*j+3*i+2] = -temp_x[2]-temp1_y[2]-temp2_y[2]-(Uk[3*nbCells_x*j+3*i+2]-Un[3*nbCells_x*j+3*i+2])\r
 541 \r
 542                     elif j==0:\r
 543 \r
 544                         #Calcul de Am_x\r
 545                         rho_l=Uk[3*nbCells_x*j+3*i]\r
 546                         qx_l =Uk[3*nbCells_x*j+3*i+1]\r
 547                         qy_l =Uk[3*nbCells_x*j+3*i+2]\r
 548                         rho_r=Uk[3*nbCells_x*j+3*(i+1)]\r
 549                         qx_r =Uk[3*nbCells_x*j+3*(i+1)+1]\r
 550                         qy_r =Uk[3*nbCells_x*j+3*(i+1)+2]\r
 551                         Am_x= jacobianMatricesm_x(dt/dx,rho_l,qx_l,qy_l,rho_r,qx_r,qy_r)\r
 552 \r
 553                         #Calcul de Ap_x\r
 554                         rho_l=Uk[3*nbCells_x*j+3*(i-1)]\r
 555                         qx_l =Uk[3*nbCells_x*j+3*(i-1)+1]\r
 556                         qy_l =Uk[3*nbCells_x*j+3*(i-1)+2]\r
 557                         rho_r=Uk[3*nbCells_x*j+3*i]\r
 558                         qx_r =Uk[3*nbCells_x*j+3*i+1]\r
 559                         qy_r =Uk[3*nbCells_x*j+3*i+2]\r
 560                         Ap_x= jacobianMatricesp_x(dt/dx,rho_l,qx_l,qy_l,rho_r,qx_r,qy_r)\r
 561 \r
 562                         #calcul de Am_y\r
 563                         rho_l=Uk[3*nbCells_x*j+3*i]\r
 564                         qx_l =Uk[3*nbCells_x*j+3*i+1]\r
 565                         qy_l =Uk[3*nbCells_x*j+3*i+2]\r
 566                         rho_r=Uk[3*nbCells_x*(j+1)+3*i]\r
 567                         qx_r =Uk[3*nbCells_x*(j+1)+3*i+1]\r
 568                         qy_r =Uk[3*nbCells_x*(j+1)+3*i+2]\r
 569                         Am_y= jacobianMatricesm_y(dt/dx,rho_l,qx_l,qy_l,rho_r,qx_r,qy_r)\r
 570 \r
 571                         #remplissage de la divMat\r
 572                         divMat.addValue(3*nbCells_x*j+3*i,3*nbCells_x*j+3*(i-1),Ap_x*(-1.))\r
 573                         divMat.addValue(3*nbCells_x*j+3*i,3*nbCells_x*j+3*(i+1),Am_x)\r
 574                         divMat.addValue(3*nbCells_x*j+3*i,3*nbCells_x*(j+1)+3*i,Am_y)\r
 575                         divMat.addValue(3*nbCells_x*j+3*i,3*nbCells_x*j+3*i,Ap_x-Am_x-Am_y)\r
 576 \r
 577                         #remplissage du membre de droite\r
 578                         dUi1_x[0] = Uk[3*nbCells_x*j+3*(i+1)]-Uk[3*nbCells_x*j+3*i]\r
 579                         dUi1_x[1] = Uk[3*nbCells_x*j+3*(i+1)+1]-Uk[3*nbCells_x*j+3*i+1]\r
 580                         dUi1_x[2] = Uk[3*nbCells_x*j+3*(i+1)+2]-Uk[3*nbCells_x*j+3*i+2]\r
 581                         dUi2_x[0] = Uk[3*nbCells_x*j+3*i]-Uk[3*nbCells_x*j+3*(i-1)]\r
 582                         dUi2_x[1] = Uk[3*nbCells_x*j+3*i+1]-Uk[3*nbCells_x*j+3*(i-1)+1]\r
 583                         dUi2_x[2] = Uk[3*nbCells_x*j+3*i+2]-Uk[3*nbCells_x*j+3*(i-1)+2]\r
 584                         dUi_y[0] = Uk[3*nbCells_x*(j+1)+3*i]-Uk[3*nbCells_x*j+3*i]\r
 585                         dUi_y[1] = Uk[3*nbCells_x*(j+1)+3*i+1]-Uk[3*nbCells_x*j+3*i+1]\r
 586                         dUi_y[2] = Uk[3*nbCells_x*(j+1)+3*i+2]-Uk[3*nbCells_x*j+3*i+2]\r
 587 \r
 588                         temp1_x = Am_x*dUi1_x\r
 589                         temp2_x = Ap_x*dUi2_x\r
 590                         temp_y = Am_y*dUi_y\r
 591                         Rhs[3*nbCells_x*j+3*i+0] = -temp1_x[0]-temp2_x[0]-temp_y[0]-(Uk[3*nbCells_x*j+3*i+0]-Un[3*nbCells_x*j+3*i+0])\r
 592                         Rhs[3*nbCells_x*j+3*i+1] = -temp1_x[1]-temp2_x[1]-temp_y[1]-(Uk[3*nbCells_x*j+3*i+1]-Un[3*nbCells_x*j+3*i+1])\r
 593                         Rhs[3*nbCells_x*j+3*i+2] = -temp1_x[2]-temp2_x[2]-temp_y[2]-(Uk[3*nbCells_x*j+3*i+2]-Un[3*nbCells_x*j+3*i+2])\r
 594 \r
 595                     elif j==nbCells_y-1:\r
 596 \r
 597                         #Calcul de Am_x\r
 598                         rho_l=Uk[3*nbCells_x*j+3*i]\r
 599                         qx_l =Uk[3*nbCells_x*j+3*i+1]\r
 600                         qy_l =Uk[3*nbCells_x*j+3*i+2]\r
 601                         rho_r=Uk[3*nbCells_x*j+3*(i+1)]\r
 602                         qx_r =Uk[3*nbCells_x*j+3*(i+1)+1]\r
 603                         qy_r =Uk[3*nbCells_x*j+3*(i+1)+2]\r
 604                         Am_x= jacobianMatricesm_x(dt/dx,rho_l,qx_l,qy_l,rho_r,qx_r,qy_r)\r
 605 \r
 606                         #Calcul de Ap_x\r
 607                         rho_l=Uk[3*nbCells_x*j+3*(i-1)]\r
 608                         qx_l =Uk[3*nbCells_x*j+3*(i-1)+1]\r
 609                         qy_l =Uk[3*nbCells_x*j+3*(i-1)+2]\r
 610                         rho_r=Uk[3*nbCells_x*j+3*i]\r
 611                         qx_r =Uk[3*nbCells_x*j+3*i+1]\r
 612                         qy_r =Uk[3*nbCells_x*j+3*i+2]\r
 613                         Ap_x= jacobianMatricesp_x(dt/dx,rho_l,qx_l,qy_l,rho_r,qx_r,qy_r)\r
 614 \r
 615                         #calcul de Ap_y\r
 616                         rho_l=Uk[3*nbCells_x*(j-1)+3*i]\r
 617                         qx_l =Uk[3*nbCells_x*(j-1)+3*i+1]\r
 618                         qy_l =Uk[3*nbCells_x*(j-1)+3*i+2]\r
 619                         rho_r=Uk[3*nbCells_x*j+3*i]\r
 620                         qx_r =Uk[3*nbCells_x*j+3*i+1]\r
 621                         qy_r =Uk[3*nbCells_x*j+3*i+2]\r
 622                         Ap_y= jacobianMatricesp_y(dt/dx,rho_l,qx_l,qy_l,rho_r,qx_r,qy_r)\r
 623 \r
 624                         #remplissage de la divMat\r
 625                         divMat.addValue(3*nbCells_x*j+3*i,3*nbCells_x*j+3*(i-1),Ap_x*(-1.))\r
 626                         divMat.addValue(3*nbCells_x*j+3*i,3*nbCells_x*j+3*(i+1),Am_x)\r
 627                         divMat.addValue(3*nbCells_x*j+3*i,3*nbCells_x*(j-1)+3*i,Ap_y*(-1.))\r
 628                         divMat.addValue(3*nbCells_x*j+3*i,3*nbCells_x*j+3*i,Ap_x+Ap_y-Am_x)\r
 629 \r
 630                         #remplissage du membre de droite\r
 631                         dUi1_x[0] = Uk[3*nbCells_x*j+3*(i+1)]-Uk[3*nbCells_x*j+3*i]\r
 632                         dUi1_x[1] = Uk[3*nbCells_x*j+3*(i+1)+1]-Uk[3*nbCells_x*j+3*i+1]\r
 633                         dUi1_x[2] = Uk[3*nbCells_x*j+3*(i+1)+2]-Uk[3*nbCells_x*j+3*i+2]\r
 634                         dUi2_x[0] = Uk[3*nbCells_x*j+3*i]-Uk[3*nbCells_x*j+3*(i-1)]\r
 635                         dUi2_x[1] = Uk[3*nbCells_x*j+3*i+1]-Uk[3*nbCells_x*j+3*(i-1)+1]\r
 636                         dUi2_x[2] = Uk[3*nbCells_x*j+3*i+2]-Uk[3*nbCells_x*j+3*(i-1)+2]\r
 637                         dUi_y[0] = Uk[3*nbCells_x*j+3*i]-Uk[3*nbCells_x*(j-1)+3*i]\r
 638                         dUi_y[1] = Uk[3*nbCells_x*j+3*i+1]-Uk[3*nbCells_x*(j-1)+3*i+1]\r
 639                         dUi_y[2] = Uk[3*nbCells_x*j+3*i+2]-Uk[3*nbCells_x*(j-1)+3*i+2]\r
 640 \r
 641                         temp1_x = Am_x*dUi1_x\r
 642                         temp2_x = Ap_x*dUi2_x\r
 643                         temp_y = Ap_y*dUi_y\r
 644                         Rhs[3*nbCells_x*j+3*i+0] = -temp1_x[0]-temp2_x[0]-temp_y[0]-(Uk[3*nbCells_x*j+3*i+0]-Un[3*nbCells_x*j+3*i+0])\r
 645                         Rhs[3*nbCells_x*j+3*i+1] = -temp1_x[1]-temp2_x[1]-temp_y[1]-(Uk[3*nbCells_x*j+3*i+1]-Un[3*nbCells_x*j+3*i+1])\r
 646                         Rhs[3*nbCells_x*j+3*i+2] = -temp1_x[2]-temp2_x[2]-temp_y[2]-(Uk[3*nbCells_x*j+3*i+2]-Un[3*nbCells_x*j+3*i+2])\r
 647 \r
 648                     #Traitement des autres cellules\r
 649                     else:\r
 650 \r
 651                         #Calcul de Am_x\r
 652                         rho_l=Uk[3*nbCells_x*j+3*i]\r
 653                         qx_l =Uk[3*nbCells_x*j+3*i+1]\r
 654                         qy_l =Uk[3*nbCells_x*j+3*i+2]\r
 655                         rho_r=Uk[3*nbCells_x*j+3*(i+1)]\r
 656                         qx_r =Uk[3*nbCells_x*j+3*(i+1)+1]\r
 657                         qy_r =Uk[3*nbCells_x*j+3*(i+1)+2]\r
 658                         Am_x= jacobianMatricesm_x(dt/dx,rho_l,qx_l,qy_l,rho_r,qx_r,qy_r)\r
 659 \r
 660                         #Calcul de Ap_x\r
 661                         rho_l=Uk[3*nbCells_x*j+3*(i-1)]\r
 662                         qx_l =Uk[3*nbCells_x*j+3*(i-1)+1]\r
 663                         qy_l =Uk[3*nbCells_x*j+3*(i-1)+2]\r
 664                         rho_r=Uk[3*nbCells_x*j+3*i]\r
 665                         qx_r =Uk[3*nbCells_x*j+3*i+1]\r
 666                         qy_r =Uk[3*nbCells_x*j+3*i+2]\r
 667                         Ap_x= jacobianMatricesp_x(dt/dx,rho_l,qx_l,qy_l,rho_r,qx_r,qy_r)\r
 668 \r
 669                         #calcul de Am_y\r
 670                         rho_l=Uk[3*nbCells_x*j+3*i]\r
 671                         qx_l =Uk[3*nbCells_x*j+3*i+1]\r
 672                         qy_l =Uk[3*nbCells_x*j+3*i+2]\r
 673                         rho_r=Uk[3*nbCells_x*(j+1)+3*i]\r
 674                         qx_r =Uk[3*nbCells_x*(j+1)+3*i+1]\r
 675                         qy_r =Uk[3*nbCells_x*(j+1)+3*i+2]\r
 676                         Am_y= jacobianMatricesm_y(dt/dx,rho_l,qx_l,qy_l,rho_r,qx_r,qy_r)\r
 677 \r
 678                         #calcul de Ap_y\r
 679                         rho_l=Uk[3*nbCells_x*(j-1)+3*i]\r
 680                         qx_l =Uk[3*nbCells_x*(j-1)+3*i+1]\r
 681                         qy_l =Uk[3*nbCells_x*(j-1)+3*i+2]\r
 682                         rho_r=Uk[3*nbCells_x*j+3*i]\r
 683                         qx_r =Uk[3*nbCells_x*j+3*i+1]\r
 684                         qy_r =Uk[3*nbCells_x*j+3*i+2]\r
 685                         Ap_y= jacobianMatricesp_y(dt/dx,rho_l,qx_l,qy_l,rho_r,qx_r,qy_r)\r
 686 \r
 687                         #remplissage de la divMat\r
 688                         divMat.addValue(3*nbCells_x*j+3*i,3*nbCells_x*j+3*(i-1),Ap_x*(-1.))\r
 689                         divMat.addValue(3*nbCells_x*j+3*i,3*nbCells_x*(j-1)+3*i,Ap_y*(-1.))\r
 690                         divMat.addValue(3*nbCells_x*j+3*i,3*nbCells_x*j+3*(i+1),Am_x)\r
 691                         divMat.addValue(3*nbCells_x*j+3*i,3*nbCells_x*(j+1)+3*i,Am_y)\r
 692                         divMat.addValue(3*nbCells_x*j+3*i,3*nbCells_x*j+3*i,Ap_x+Ap_y-Am_x-Am_y)\r
 693 \r
 694                         #remplissage du membre de droite\r
 695 \r
 696                         dUi1_x[0] = Uk[3*nbCells_x*j+3*(i+1)]-Uk[3*nbCells_x*j+3*i]\r
 697                         dUi1_x[1] = Uk[3*nbCells_x*j+3*(i+1)+1]-Uk[3*nbCells_x*j+3*i+1]\r
 698                         dUi1_x[2] = Uk[3*nbCells_x*j+3*(i+1)+2]-Uk[3*nbCells_x*j+3*i+2]\r
 699 \r
 700                         dUi2_x[0] = Uk[3*nbCells_x*j+3*i]-Uk[3*nbCells_x*j+3*(i-1)]\r
 701                         dUi2_x[1] = Uk[3*nbCells_x*j+3*i+1]-Uk[3*nbCells_x*j+3*(i-1)+1]\r
 702                         dUi2_x[2] = Uk[3*nbCells_x*j+3*i+2]-Uk[3*nbCells_x*j+3*(i-1)+2]\r
 703 \r
 704                         dUi1_y[0] = Uk[3*nbCells_x*(j+1)+3*i]-Uk[3*nbCells_x*j+3*i]\r
 705                         dUi1_y[1] = Uk[3*nbCells_x*(j+1)+3*i+1]-Uk[3*nbCells_x*j+3*i+1]\r
 706                         dUi1_y[2] = Uk[3*nbCells_x*(j+1)+3*i+2]-Uk[3*nbCells_x*j+3*i+2]\r
 707 \r
 708                         dUi2_y[0] = Uk[3*nbCells_x*j+3*i]-Uk[3*nbCells_x*(j-1)+3*i]\r
 709                         dUi2_y[1] = Uk[3*nbCells_x*j+3*i+1]-Uk[3*nbCells_x*(j-1)+3*i+1]\r
 710                         dUi2_y[2] = Uk[3*nbCells_x*j+3*i+2]-Uk[3*nbCells_x*(j-1)+3*i+2]\r
 711 \r
 712                         temp1_x = Am_x*dUi1_x\r
 713                         temp2_x = Ap_x*dUi2_x\r
 714                         temp1_y = Am_y*dUi1_y\r
 715                         temp2_y = Ap_y*dUi2_y\r
 716                         Rhs[3*nbCells_x*j+3*i+0] = -temp1_x[0]-temp2_x[0]-temp1_y[0]-temp2_y[0]-(Uk[3*nbCells_x*j+3*i+0]-Un[3*nbCells_x*j+3*i+0])\r
 717                         Rhs[3*nbCells_x*j+3*i+1] = -temp1_x[1]-temp2_x[1]-temp1_y[1]-temp2_y[1]-(Uk[3*nbCells_x*j+3*i+1]-Un[3*nbCells_x*j+3*i+1])\r
 718                         Rhs[3*nbCells_x*j+3*i+2] = -temp1_x[2]-temp2_x[2]-temp1_y[2]-temp2_y[2]-(Uk[3*nbCells_x*j+3*i+2]-Un[3*nbCells_x*j+3*i+2])\r
 719 \r
 720             divMat.diagonalShift(1)  #only after  filling all coefficients\r
 721             LS=cdmath.LinearSolver(divMat,Rhs,iterGMRESMax, precision, "GMRES","LU")\r
 722             dUk=LS.solve();\r
 723             residu = dUk.norm()\r
 724             Uk+=dUk\r
 725             #print("Newton iteration number ",k, "residu = ",residu)\r
 726             if(not LS.getStatus()):\r
 727                 print("Linear system did not converge ", LS.getNumberOfIter(), " GMRES iterations")\r
 728                 raise ValueError("Pas de convergence du système linéaire");\r
 729             k=k+1\r
 730         Un = Uk.deepCopy()\r
 731         dUn-=Un\r
 732         isStationary = dUn.norm()<precision\r
 733         \r
 734         for i in range(nbCells):\r
 735             rho_field[i] = Un[nbComp*i]\r
 736             q_x_field[i] = Un[nbComp*i+1]\r
 737             q_y_field[i] = Un[nbComp*i+2]\r
 738 \r
 739         time=time+dt;\r
 740         it=it+1;\r
 741         #Sauvegardes\r
 742         if(it==1 or it%output_freq==0 or it>=ntmax or isStationary or time >=tmax):\r
 743             print("-- Iter: " + str(it) + ", Time: " + str(time) + ", dt: " + str(dt) + ", Newton iterations: "+str(k)+", ||dUn|| = "+str(dUn.norm()))\r
 744             print("Linear system converged in ", LS.getNumberOfIter(), " GMRES iterations, résidu = ", residu)\r
 745             for k in range(nbCells):\r
 746                 rho  =rho_field[k]\r
 747                 velocity_field[k,0]=q_x_field[i]/rho\r
 748                 if(dim>1):\r
 749                     velocity_field[k,1]=q_y_field[k]/rho\r
 750             print\r
 751             rho_field.setTime(time,it);\r
 752             rho_field.writeVTK("EulerIsothermal_"+str(dim)+"DConservativeStaggered_"+meshName+"_density",False);\r
 753             q_x_field.setTime(time,it);\r
 754             q_x_field.writeVTK("EulerIsothermal_"+str(dim)+"DConservativeStaggered_"+meshName+"_momentumX",False);\r
 755             q_y_field.setTime(time,it);\r
 756             q_y_field.writeVTK("EulerIsothermal_"+str(dim)+"DConservativeStaggered_"+meshName+"_momentumY",False);\r
 757             velocity_field.setTime(time,it);\r
 758             velocity_field.writeVTK("EulerIsothermal_"+str(dim)+"DConservativeStaggered_"+meshName+"_velocity",False);\r
 759             #Postprocessing : save 2D picture\r
 760             PV_routines.Save_PV_data_to_picture_file("EulerIsothermal_"+str(dim)+"DConservativeStaggered_"+meshName+"_density_"  +str(it)+'.vtu',"Density",   'CELLS',"EulerIsothermal_"+str(dim)+"DConservativeStaggered_"+meshName+"_density"  +str(it))\r
 761             PV_routines.Save_PV_data_to_picture_file("EulerIsothermal_"+str(dim)+"DConservativeStaggered_"+meshName+"_momentumX_"+str(it)+'.vtu',"Momentum x",'CELLS',"EulerIsothermal_"+str(dim)+"DConservativeStaggered_"+meshName+"_momentumX"+str(it))\r
 762             PV_routines.Save_PV_data_to_picture_file("EulerIsothermal_"+str(dim)+"DConservativeStaggered_"+meshName+"_momentumY_"+str(it)+'.vtu',"Momentum y",'CELLS',"EulerIsothermal_"+str(dim)+"DConservativeStaggered_"+meshName+"_momentumY"+str(it))\r
 763     print("-- Iter: " + str(it) + ", Time: " + str(time) + ", dt: " + str(dt))\r
 764     if(it>=ntmax):\r
 765         print("Nombre de pas de temps maximum ntmax= ", ntmax, " atteint")\r
 766         return\r
 767     elif(isStationary):\r
 768         print("Régime stationnaire atteint au pas de temps ", it, ", t= ", time)\r
 769         print("------------------------------------------------------------------------------------")\r
 770         return\r
 771     else:\r
 772         print("Temps maximum Tmax= ", tmax, " atteint")\r
 773         return\r
 774 \r
 775 def solve( my_mesh,filename, resolution):\r
 776     print("Resolution of the Isothermal Euler system in dimension 2 on "+str(my_mesh.getNumberOfCells())+ " cells")\r
 777     print("Initial data : ", "Riemann problem")\r
 778     print("Boundary conditions : ", "Neumann")\r
 779     print("Mesh name : ",filename )\r
 780     # Problem data\r
 781     tmax = 10.\r
 782     ntmax = 10\r
 783     cfl=1\r
 784     output_freq = 1\r
 785     EulerSystemStaggered(ntmax, tmax, cfl, output_freq, my_mesh, filename)\r
 786     return\r
 787 \r
 788 if __name__ == """__main__""":\r
 789     if len(sys.argv) >1 :\r
 790         filename=sys.argv[1]\r
 791         my_mesh = cdmath.Mesh(filename)\r
 792     else :\r
 793         filename = "CartesianGrid"\r
 794         ax=0;bx=1;nx=20;\r
 795         ay=0;by=1;ny=10;        \r
 796         my_mesh = cdmath.Mesh(ax,bx,nx,ay,by,ny)\r
 797 \r
 798     solve(my_mesh,filename,100)\r