--- /dev/null
+#!/usr/bin/env python3
+# -*-coding:utf-8 -*
+
+# Isothermal Euler system
+# d rho/d t + d q/d x =0
+# d q/d t + d (q^2/rho+p)/d x = 0, where q = rho*u and p = c^2*rho
+# UU = (rho,q) : conservative variable
+# Scheme : Conservative stagerred scheme (Ait Ameur et al)
+# Author : Katia Ait Ameur
+# Date : November 2020
+
+import cdmath
+import numpy as np
+import matplotlib
+matplotlib.use("Agg")
+import matplotlib.pyplot as plt
+import matplotlib.animation as manimation
+from math import sqrt
+
+c0=330.#reference sound speed for water at 15 bars
+precision=1e-5
+rho0=1
+
+def initial_conditions_Riemann_problem(a,b,nx):
+ print("Initial data Riemann problem")
+ dx = (b - a) / nx #space step
+ x=[a+0.5*dx + i*dx for i in range(nx)] # array of cell center (1D mesh)
+ rho_initial = [ (xi<(a+b)/2)*rho0 + (xi>=(a+b)/2)*rho0*2 for xi in x]
+ q_initial = [ (xi<(a+b)/2)*rho0*(-300) + (xi>=(a+b)/2)*rho0*2*(-300) for xi in x]
+
+ return rho_initial, q_initial
+
+def matrix_coef(rho_l,q_l,rho_r,q_r):
+ m_coef = -0.5*(q_l-q_r)/rho_r
+ return m_coef
+
+def jacobianMatricesm(coeff,rho_l,q_l,rho_r,q_r):
+ RoeMat = cdmath.Matrix(2,2);
+ Dmac = cdmath.Matrix(2,2);
+
+ u_l=q_l/rho_l
+ u_r=q_r/rho_r
+ Dmac_coef = matrix_coef(rho_l,q_l,rho_r,q_r)
+ if rho_l<0 or rho_r<0 :
+ print( "rho_l=",rho_l, " rho_r= ",rho_r)
+ raise ValueError("Negative density")
+ u = (u_l*sqrt(rho_l)+u_r*sqrt(rho_r))/(sqrt(rho_l)+sqrt(rho_r));
+ RoeMat[0,0] = 0
+ RoeMat[0,1] = 1
+ RoeMat[1,0] = c0*c0 - u*u
+ RoeMat[1,1] = 2*u
+
+ Dmac[0,0]= abs(u)-u;
+ Dmac[0,1]= 1;
+ Dmac[1,0]= -c0*c0-u*u;
+ Dmac[1,1]= abs(u)+u;
+
+ return (RoeMat-Dmac)*coeff*0.5
+
+def jacobianMatricesp(coeff,rho_l,q_l,rho_r,q_r):
+ RoeMat = cdmath.Matrix(2,2);
+ Dmac = cdmath.Matrix(2,2);
+
+ u_l=q_l/rho_l
+ u_r=q_r/rho_r
+ Dmac_coef = matrix_coef(rho_l,q_l,rho_r,q_r)
+ if rho_l<0 or rho_r<0 :
+ print( "rho_l=",rho_l, " rho_r= ",rho_r)
+ raise ValueError("Negative density")
+ u = (u_l*sqrt(rho_l)+u_r*sqrt(rho_r))/(sqrt(rho_l)+sqrt(rho_r));
+ RoeMat[0,0] = 0
+ RoeMat[0,1] = 1
+ RoeMat[1,0] = c0*c0 - u*u
+ RoeMat[1,1] = 2*u
+
+ Dmac[0,0]= abs(u)-u;
+ Dmac[0,1]= 1;
+ Dmac[1,0]= -c0*c0-u*u;
+ Dmac[1,1]= abs(u)+u;
+
+ return (RoeMat+Dmac)*coeff*0.5
+
+def EulerSystemStaggered(ntmax, tmax, cfl, a,b,nx, output_freq, meshName):
+ dim=1
+ nbComp=dim+1
+ nbCells = nx
+ dt = 0.
+ dx=(b-a)/nx
+ dt = cfl * dx / c0
+ nbVoisinsMax=2
+
+ time = 0.
+ it=0;
+ iterGMRESMax=50
+ newton_max = 50
+ isStationary=False
+
+ #iteration vectors
+ Un =cdmath.Vector(nbCells*(dim+1))
+ dUn=cdmath.Vector(nbCells*(dim+1))
+ dUk=cdmath.Vector(nbCells*(dim+1))
+ Rhs=cdmath.Vector(nbCells*(dim+1))
+
+ # Initial conditions #
+ print("Construction of the initial condition …")
+ rho_field, q_field = initial_conditions_Riemann_problem(a,b,nx)
+ v_field = [ q_field[i]/rho_field[i] for i in range(nx)]
+ p_field = [ rho_field[i]*(c0*c0) for i in range(nx)]
+
+ max_initial_rho=max(rho_field)
+ min_initial_rho=min(rho_field)
+ max_initial_q=max(q_field)
+ min_initial_q=min(q_field)
+ max_initial_p=max(p_field)
+ min_initial_p=min(p_field)
+ max_initial_v=max(v_field)
+ min_initial_v=min(v_field)
+
+ for k in range(nbCells):
+ Un[k*nbComp+0] = rho_field[k]
+ Un[k*nbComp+1] = q_field[k]
+
+ print("Starting computation of the non linear Euler system with staggered scheme …")
+ divMat=cdmath.SparseMatrixPetsc(nbCells*nbComp,nbCells*nbComp,(nbVoisinsMax+1)*nbComp)
+
+ # Picture settings
+ fig, ([axDensity, axMomentum],[axVelocity, axPressure]) = plt.subplots(2, 2,sharex=True, figsize=(10,10))
+ fig.suptitle('Conservative staggered scheme')
+ lineDensity, = axDensity.plot([a+0.5*dx + i*dx for i in range(nx)], rho_field, label='Density') #new picture for video # Returns a tuple of line objects, thus the comma
+ axDensity.set(xlabel='x (m)', ylabel='Density')
+ axDensity.set_xlim(a,b)
+ axDensity.set_ylim(min_initial_rho - 0.1*(max_initial_rho-min_initial_rho), max_initial_rho + 0.1*(max_initial_rho-min_initial_rho) )
+ axDensity.legend()
+ lineMomentum, = axMomentum.plot([a+0.5*dx + i*dx for i in range(nx)], q_field, label='Momentum')
+ axMomentum.set(xlabel='x (m)', ylabel='Momentum')
+ axMomentum.set_xlim(a,b)
+ axMomentum.set_ylim(min_initial_q - 0.1*(max_initial_q-min_initial_q), max_initial_q + 0.1*(max_initial_q-min_initial_q) )
+ axMomentum.legend()
+ lineVelocity, = axVelocity.plot([a+0.5*dx + i*dx for i in range(nx)], v_field, label='Velocity')
+ axVelocity.set(xlabel='x (m)', ylabel='Velocity')
+ axVelocity.set_xlim(a,b)
+ axVelocity.set_ylim(min_initial_v - 0.4*abs(min_initial_v), max_initial_v + 0.05*abs(max_initial_v) )
+ axVelocity.legend()
+ linePressure, = axPressure.plot([a+0.5*dx + i*dx for i in range(nx)], p_field, label='Pressure')
+ axPressure.set(xlabel='x (m)', ylabel='Pressure')
+ axPressure.set_xlim(a,b)
+ axPressure.set_ylim(min_initial_p - 0.05*abs(min_initial_p), max_initial_p + 0.05*abs(max_initial_p) )
+ axPressure.ticklabel_format(axis='y', style='sci', scilimits=(0,0))
+ axPressure.legend()
+
+ # Video settings
+ FFMpegWriter = manimation.writers['ffmpeg']
+ metadata = dict(title="Conservative staggered scheme for the 1D isothermal Euler System", artist = "CEA Saclay", comment="Shock tube")
+ writer=FFMpegWriter(fps=10, metadata=metadata, codec='h264')
+ with writer.saving(fig, "1DEuler_System_ConservativeStaggered"+".mp4", ntmax):
+ writer.grab_frame()
+ plt.savefig("EulerSystem"+str(dim)+"D_ConservativeStaggered"+meshName+"_0"+".png")
+
+ # Starting time loop
+ while (it<ntmax and time <= tmax and not isStationary):
+ dUn = Un.deepCopy()
+ Uk = Un.deepCopy()
+ residu = 1.
+ k=0
+ while (k<newton_max and residu > precision ):
+ #DEBUT BOUCLE NEWTON
+ divMat.zeroEntries()#sets the matrix coefficients to zero
+ for j in range(nbCells):#
+ if ( j==0) :
+ rho_r = Uk[j*nbComp+0]
+ q_r = Uk[j*nbComp+1]
+ rho_l = rho_r # Conditions de Neumann
+ q_l = q_r
+ Am= jacobianMatricesm(dt/dx,rho_l,q_l,rho_r,q_r)
+ divMat.addValue(j*nbComp,(j+1)*nbComp,Am)
+ divMat.addValue(j*nbComp, j*nbComp,Am*(-1.))
+ dUi=cdmath.Vector(2)
+ dUi[0] = Uk[(j+1)*nbComp+0]-Uk[j*nbComp+0]
+ dUi[1] = Uk[(j+1)*nbComp+1]-Uk[j*nbComp+1]
+ temp=cdmath.Vector(2)
+ temp = Am*dUi
+ #print("Bloc 0 matrix : ", Am)
+ Rhs[j*nbComp+0] = -temp[0]-(Uk[j*nbComp+0]-Un[j*nbComp+0])
+ Rhs[j*nbComp+1] = -temp[1]-(Uk[j*nbComp+1]-Un[j*nbComp+1])
+ elif ( j==nbCells-1) :
+ rho_l = Uk[j*nbComp+0]
+ q_l = Uk[j*nbComp+1]
+ rho_r = rho_l # Conditions de Neumann
+ q_r = q_l
+ Ap= jacobianMatricesp(dt/dx,rho_l,q_l,rho_r,q_r)
+ divMat.addValue(j*nbComp, j*nbComp,Ap)
+ divMat.addValue(j*nbComp,(j-1)*nbComp,Ap*(-1.))
+ dUi=cdmath.Vector(2)
+ dUi[0] = Uk[(j-1)*nbComp+0]-Uk[j*nbComp+0]
+ dUi[1] = Uk[(j-1)*nbComp+1]-Uk[j*nbComp+1]
+ temp=cdmath.Vector(2)
+ temp = Ap*dUi
+ Rhs[j*nbComp+0] = temp[0]-(Uk[j*nbComp+0]-Un[j*nbComp+0])
+ Rhs[j*nbComp+1] = temp[1]-(Uk[j*nbComp+1]-Un[j*nbComp+1])
+ else :
+ rho_l = Uk[(j-1)*nbComp+0]
+ q_l = Uk[(j-1)*nbComp+1]
+ rho_r = Uk[j*nbComp+0]
+ q_r = Uk[j*nbComp+1]
+ Ap = jacobianMatricesp(dt/dx,rho_l,q_l,rho_r,q_r)
+ ###############################################################
+ rho_l = Uk[j*nbComp+0]
+ q_l = Uk[j*nbComp+1]
+ rho_r = Uk[(j+1)*nbComp+0]
+ q_r = Uk[(j+1)*nbComp+1]
+ Am = jacobianMatricesm(dt/dx,rho_l,q_l,rho_r,q_r)
+ divMat.addValue(j*nbComp,(j+1)*nbComp,Am)
+ divMat.addValue(j*nbComp, j*nbComp,Am*(-1.))
+ divMat.addValue(j*nbComp, j*nbComp,Ap)
+ divMat.addValue(j*nbComp,(j-1)*nbComp,Ap*(-1.))
+ dUi1=cdmath.Vector(2)
+ dUi2=cdmath.Vector(2)
+ dUi1[0] = Uk[(j+1)*nbComp+0]-Uk[j*nbComp+0]
+ dUi1[1] = Uk[(j+1)*nbComp+1]-Uk[j*nbComp+1]
+ dUi2[0] = Uk[(j-1)*nbComp+0]-Uk[j*nbComp+0]
+ dUi2[1] = Uk[(j-1)*nbComp+1]-Uk[j*nbComp+1]
+ temp1 = cdmath.Vector(2)
+ temp2 = cdmath.Vector(2)
+ temp1 = Am*dUi1
+ temp2 = Ap*dUi2
+ Rhs[j*nbComp+0] = -temp1[0] + temp2[0] -(Uk[j*nbComp+0]-Un[j*nbComp+0])
+ Rhs[j*nbComp+1] = -temp1[1] + temp2[1] -(Uk[j*nbComp+1]-Un[j*nbComp+1])
+ divMat.diagonalShift(1)#only after filling all coefficients
+ LS=cdmath.LinearSolver(divMat,Rhs,iterGMRESMax, precision, "GMRES","LU")
+ dUk=LS.solve();
+ residu = dUk.norm()
+ #stop
+ Uk+=dUk
+ if(not LS.getStatus()):
+ print("Linear system did not converge ", LS.getNumberOfIter(), " GMRES iterations")
+ raise ValueError("Pas de convergence du système linéaire");
+ k=k+1
+ Un = Uk.deepCopy()
+ dUn-=Un
+ for k in range(nbCells):
+ rho_field[k] = Un[k*(dim+1)+0]
+ q_field[k] = Un[k*(dim+1)+1]
+ v_field = [ q_field[i]/rho_field[i] for i in range(nx)]
+ p_field = [ rho_field[i]*(c0*c0) for i in range(nx)]
+
+ lineDensity.set_ydata(rho_field)
+ lineMomentum.set_ydata(q_field)
+ lineVelocity.set_ydata(v_field)
+ linePressure.set_ydata(p_field)
+ writer.grab_frame()
+
+ time=time+dt;
+ it=it+1;
+ #Sauvegardes
+ if(it==1 or it%output_freq==0 or it>=ntmax or isStationary or time >=tmax):
+ print("-- Iter: " + str(it) + ", Time: " + str(time) + ", dt: " + str(dt))
+ #print("Last linear system converged in ", LS.getNumberOfIter(), " GMRES iterations", ", residu final: ",residu)
+
+ plt.savefig("EulerSystem"+str(dim)+"D_ConservativeStaggered"+meshName+"_"+str(it)+".png")
+ print
+ print("-- Iter: " + str(it) + ", Time: " + str(time) + ", dt: " + str(dt))
+ if(it>=ntmax):
+ print("Nombre de pas de temps maximum ntmax= ", ntmax, " atteint")
+ return
+ elif(isStationary):
+ print("Régime stationnaire atteint au pas de temps ", it, ", t= ", time)
+ print("------------------------------------------------------------------------------------")
+ plt.savefig("EulerSystem"+str(dim)+"Staggered"+meshName+"_Stat.png")
+ return
+ else:
+ print("Temps maximum Tmax= ", tmax, " atteint")
+ return
+
+def solve( a,b,nx, meshName, meshType, cfl):
+ print("Resolution of the Euler system in dimension 1 on "+str(nx)+ " cells")
+ print("Initial data : ", "Riemann problem")
+ print("Boundary conditions : ", "Neumann")
+ print("Mesh name : ",meshName , ", ", nx, " cells")
+ # Problem data
+ tmax = 10.
+ ntmax = 50
+ output_freq = 10
+ EulerSystemStaggered(ntmax, tmax, cfl, a,b,nx, output_freq, meshName)
+ return
+
+if __name__ == """__main__""":
+ a=0.
+ b=1.
+ nx=100
+ cfl=0.99
+ solve( a,b,nx,"SquareRegularSquares","RegularSquares",cfl)
--- /dev/null
+# Isothermal Euler system
+# d rho/d t + d q/d x =0
+# d q/d t + d (q^2/rho+p)/d x = 0, where q = rho*u and p = c^2*rho
+# UU = (rho,q) : conservative variable
+# Scheme : Roe scheme (without any entropy correction)
+# Comment : the solution displays a non entropic (non physical) shock instead of a rarefaction wave
+# Date : November 2020
+
+from cdmath import *
+import numpy as np
+from math import sqrt
+
+import matplotlib
+matplotlib.use("Agg")
+import matplotlib.pyplot as plt
+import matplotlib.animation as manimation
+
+c0=330.#reference sound speed for water at 15 bars
+precision=1e-5
+
+def flux(rho,q):
+ fl = Vector(2);
+ fl[0] = q;
+ fl[1] = q*q/rho + c0*c0*rho;
+ return fl
+
+def compTimeStep(UU,dx,CFL):
+ M = UU.getMesh();
+ maxAbsEigVa = 0;
+ for i in range(M.getNumberOfCells()):
+ maxAbsEigVa = max(maxAbsEigVa,abs(UU[i,1]/UU[i,0]+c0),abs(UU[i,1]/UU[i,0]-c0));
+ pass
+ dt = CFL*dx/maxAbsEigVa;
+ return dt
+
+def AbsRoeMatrix(rho_l,q_l,rho_r,q_r):
+ AbsRoeMa = Matrix(2,2);
+ u_l = q_l/rho_l;
+ u_r = q_r/rho_r;
+ u = (u_l*sqrt(rho_l)+u_r*sqrt(rho_r))/(sqrt(rho_l)+sqrt(rho_r));
+ AbsRoeMa[0,0]=(abs(u-c0)*(u+c0)+abs(u+c0)*(c0-u))/(2*c0);
+ AbsRoeMa[0,1]= (abs(u+c0)-abs(u-c0))/(2*c0);
+ AbsRoeMa[1,0]=(abs(u-c0)*(u*u-c0*c0)+abs(u+c0)*(c0*c0-u*u))/(2*c0);
+ AbsRoeMa[1,1]=((u+c0)*abs(u+c0)-(u-c0)*abs(u-c0))/(2*c0);
+ return AbsRoeMa;
+
+def main(meshName):
+ # mesh definition
+ xmin=.0;
+ xmax= 1;
+ nx=100#1000;
+ ntmax=50;
+ dx = (xmax-xmin)/nx;
+ M=Mesh(xmin,xmax,nx);
+ dim=1
+
+ # initial time
+ time=0.;
+ # maximum time T
+ tmax=0.01;
+ it=0;
+ freqSortie=10;
+ # CFL condition
+ CFL=0.99;
+ # conservative variable (unknown)
+ UU=Field("Conservative variables",CELLS,M,2);
+
+ # initial condition for Riemann problem
+ print("Construction of the initial condition …")
+ rhoL=1
+ vL=-300
+ rhoR=1.45
+ vR=-300
+ for i in range(M.getNumberOfCells()):
+ x=M.getCell(i).x();
+ if x < (xmax-xmin)/2:
+ UU[i,0] = rhoL;
+ UU[i,1] = rhoL*vL;
+ else:
+ UU[i,0] = rhoR;
+ UU[i,1] = rhoR*vR;
+ pass
+
+ rho_field = [ UU[i,0] for i in range(nx)]
+ q_field = [ UU[i,1] for i in range(nx)]
+ v_field = [ UU[i,1]/UU[i,0] for i in range(nx)]
+ p_field = [ UU[i,0]*(c0*c0) for i in range(nx)]
+
+ max_initial_rho=max(rhoL,rhoR)
+ min_initial_rho=min(rhoL,rhoR)
+ max_initial_p=max_initial_rho*c0*c0
+ min_initial_p=min_initial_rho*c0*c0
+ max_initial_v=max(vL,vR)
+ min_initial_v=min(vL,vR)
+ max_initial_q=max(vL*rhoL,vR*rhoR)
+ min_initial_q=min(vL*rhoL,vR*rhoR)
+
+ print( "Numerical solution of the 1D Euler equation")
+
+ # prepare some memory
+ UU_limL = Vector(2);
+ UU_limR = Vector(2);
+ Del_UU_RL = Vector(2);
+ UU_new = UU;
+
+ # Picture settings
+ fig, ([axDensity, axMomentum],[axVelocity, axPressure]) = plt.subplots(2, 2,sharex=True, figsize=(10,10))
+ fig.suptitle('Explicit Upwind scheme')
+ lineDensity, = axDensity.plot([xmin+0.5*dx + i*dx for i in range(nx)], rho_field, label='Density') #new picture for video # Returns a tuple of line objects, thus the comma
+ axDensity.set(xlabel='x (m)', ylabel='Density')
+ axDensity.set_xlim(xmin,xmax)
+ axDensity.set_ylim(min_initial_rho - 0.1*(max_initial_rho-min_initial_rho), max_initial_rho + 0.1*(max_initial_rho-min_initial_rho) )
+ axDensity.legend()
+ lineMomentum, = axMomentum.plot([xmin+0.5*dx + i*dx for i in range(nx)], q_field, label='Momentum')
+ axMomentum.set(xlabel='x (m)', ylabel='Momentum')
+ axMomentum.set_xlim(xmin,xmax)
+ axMomentum.set_ylim(min_initial_q - 0.1*(max_initial_q-min_initial_q), max_initial_q + 0.1*(max_initial_q-min_initial_q) )
+ axMomentum.legend()
+ lineVelocity, = axVelocity.plot([xmin+0.5*dx + i*dx for i in range(nx)], v_field, label='Velocity')
+ axVelocity.set(xlabel='x (m)', ylabel='Velocity')
+ axVelocity.set_xlim(xmin,xmax)
+ axVelocity.set_ylim(min_initial_v - 0.25*abs(min_initial_v), max_initial_v + 0.05*abs(max_initial_v) )
+ axVelocity.legend()
+ linePressure, = axPressure.plot([xmin+0.5*dx + i*dx for i in range(nx)], p_field, label='Pressure')
+ axPressure.set(xlabel='x (m)', ylabel='Pressure')
+ axPressure.set_xlim(xmin,xmax)
+ axPressure.set_ylim(min_initial_p - 0.05*abs(min_initial_p), max_initial_p + 0.05*abs(max_initial_p) )
+ axPressure.ticklabel_format(axis='y', style='sci', scilimits=(0,0))
+ axPressure.legend()
+
+ # Video settings
+ FFMpegWriter = manimation.writers['ffmpeg']
+ metadata = dict(title="Upwind (Roe) scheme for the 1D isothermal Euler System", artist = "CEA Saclay", comment="Shock tube")
+ writer=FFMpegWriter(fps=10, metadata=metadata, codec='h264')
+ with writer.saving(fig, "1DEuler_System_Upwind"+".mp4", ntmax):
+ writer.grab_frame()
+ plt.savefig("EulerSystem"+str(dim)+"UpwindExplicit"+meshName+"_0"+".png")
+
+ # loop in time
+ while (it<ntmax and time <= tmax ):
+ # time step
+ dt=compTimeStep(UU,dx,CFL);
+ # Neumann boundary condition
+ UU_limL[0] = UU[0,0];
+ UU_limL[1] = UU[0,1];
+ UU_limR[0] = UU[M.getNumberOfCells()-1,0];
+ UU_limR[1] = UU[M.getNumberOfCells()-1,1];
+ flux_iminus = flux(UU_limL[0],UU_limL[1]);
+ #Main loop on cells
+ for i in range(0,M.getNumberOfCells()):
+ if (i<M.getNumberOfCells()-1):
+ Del_UU_RL[0]=UU[i+1,0]-UU[i,0];
+ Del_UU_RL[1]=UU[i+1,1]-UU[i,1];
+ AbsRoeMa = AbsRoeMatrix(UU[i,0],UU[i,1],UU[i+1,0],UU[i+1,1]);
+ flux_iplus = (flux(UU[i,0],UU[i,1])+flux(UU[i+1,0],UU[i+1,1]))-AbsRoeMa*Del_UU_RL;
+ flux_iplus[0]*=0.5;
+ flux_iplus[1]*=0.5;
+ else:
+ flux_iplus= flux(UU_limR[0],UU_limR[1]);
+ UU_new[i,0] = UU[i,0] - dt/dx*(flux_iplus[0]-flux_iminus[0]);
+ UU_new[i,1] = UU[i,1] - dt/dx*(flux_iplus[1]-flux_iminus[1]);
+ flux_iminus = flux_iplus;
+ pass
+ UU = UU_new;
+ time+=dt;
+ it+=1;
+
+ rho_field = [ UU[i,0] for i in range(nx)]
+ q_field = [ UU[i,1] for i in range(nx)]
+ v_field = [ UU[i,1]/UU[i,0] for i in range(nx)]
+ p_field = [ UU[i,0]*(c0*c0) for i in range(nx)]
+
+ lineDensity.set_ydata(rho_field)
+ lineMomentum.set_ydata(q_field)
+ lineVelocity.set_ydata(v_field)
+ linePressure.set_ydata(p_field)
+ writer.grab_frame()
+
+ if (it%freqSortie==0):
+ print( "-- Iter : ", it," Time : ",time," dt : ",dt)
+ plt.savefig("EulerSystem"+str(dim)+"UpwindExplicit"+meshName+"_"+str(it)+".png")
+ pass
+
+ print("-- Iter: " + str(it) + ", Time: " + str(time) + ", dt: " + str(dt))
+ if(it>=ntmax):
+ print("Nombre de pas de temps maximum ntmax= ", ntmax, " atteint")
+ return
+ elif(isStationary):
+ print("Régime stationnaire atteint au pas de temps ", it, ", t= ", time)
+ print("------------------------------------------------------------------------------------")
+ plt.savefig("EulerSystem"+str(dim)+"UpwindExplicit"+meshName+"_Stat.png")
+ return
+ else:
+ print("Temps maximum Tmax= ", tmax, " atteint")
+
+ return
+
+if __name__ == '__main__':
+ main("SquareRegularSquares")
--- /dev/null
+# Isothermal Euler system
+# d rho/d t + d q/d x =0
+# d q/d t + d (q^2/rho+p)/d x = 0, where q = rho*u and p = c^2*rho
+# UU = (rho,q) : conservative variable
+# Scheme : Roe scheme with entropy correction
+# Comment : If one removes the entropic correction the result displays a non entropic (non physical) shock instead of a rarefaction wave
+# Date : November 2020
+
+from cdmath import *
+import numpy as np
+from math import sqrt
+
+import matplotlib
+matplotlib.use("Agg")
+import matplotlib.pyplot as plt
+import matplotlib.animation as manimation
+
+c0=330.#reference sound speed for water at 15 bars
+precision=1e-5
+
+def flux(rho,q):
+ fl = Vector(2);
+ fl[0] = q;
+ fl[1] = q*q/rho + c0*c0*rho;
+ return fl
+
+def compTimeStep(UU,dx,CFL):
+ M = UU.getMesh();
+ maxAbsEigVa = 0;
+ for i in range(M.getNumberOfCells()):
+ maxAbsEigVa = max(maxAbsEigVa,abs(UU[i,1]/UU[i,0]+c0),abs(UU[i,1]/UU[i,0]-c0));
+ pass
+ dt = CFL*dx/maxAbsEigVa;
+ return dt
+
+def AbsRoeMatrix(rho_l,q_l,rho_r,q_r):
+ AbsRoeMa = Matrix(2,2);
+ u_l = q_l/rho_l;
+ u_r = q_r/rho_r;
+ u = (u_l*sqrt(rho_l)+u_r*sqrt(rho_r))/(sqrt(rho_l)+sqrt(rho_r));
+ AbsRoeMa[0,0]=(abs(u-c0)*(u+c0)+abs(u+c0)*(c0-u))/(2*c0);
+ AbsRoeMa[0,1]= (abs(u+c0)-abs(u-c0))/(2*c0);
+ AbsRoeMa[1,0]=(abs(u-c0)*(u*u-c0*c0)+abs(u+c0)*(c0*c0-u*u))/(2*c0);
+ AbsRoeMa[1,1]=((u+c0)*abs(u+c0)-(u-c0)*abs(u-c0))/(2*c0);
+ return AbsRoeMa;
+
+def main(meshName):
+ # mesh definition
+ xmin=.0;
+ xmax= 1;
+ nx=200#1000;
+ ntmax=100;
+ dx = (xmax-xmin)/nx;
+ M=Mesh(xmin,xmax,nx);
+ dim=1
+
+ # initial time
+ time=0.;
+ # maximum time T
+ tmax=0.01;
+ it=0;
+ freqSortie=10;
+ # CFL condition
+ CFL=0.9;
+ # conservative variable (unknown)
+ UU=Field("Conservative variables",CELLS,M,2);
+
+ # initial condition for Riemann problem
+ print("Construction of the initial condition …")
+ rhoL=1
+ vL=-300
+ rhoR=2
+ vR=-300
+ for i in range(M.getNumberOfCells()):
+ x=M.getCell(i).x();
+ if x < (xmax-xmin)/2:
+ UU[i,0] = rhoL;
+ UU[i,1] = rhoL*vL;
+ else:
+ UU[i,0] = rhoR;
+ UU[i,1] = rhoR*vR;
+ pass
+
+ rho_field = [ UU[i,0] for i in range(nx)]
+ q_field = [ UU[i,1] for i in range(nx)]
+ v_field = [ UU[i,1]/UU[i,0] for i in range(nx)]
+ p_field = [ UU[i,0]*(c0*c0) for i in range(nx)]
+
+ max_initial_rho=max(rhoL,rhoR)
+ min_initial_rho=min(rhoL,rhoR)
+ max_initial_p=max_initial_rho*c0*c0
+ min_initial_p=min_initial_rho*c0*c0
+ max_initial_v=max(vL,vR)
+ min_initial_v=min(vL,vR)
+ max_initial_q=max(vL*rhoL,vR*rhoR)
+ min_initial_q=min(vL*rhoL,vR*rhoR)
+
+ print( "Numerical solution of the 1D Euler equation")
+
+ # prepare some memory
+ UU_limL = Vector(2);
+ UU_limR = Vector(2);
+ Del_UU_RL = Vector(2);
+ UU_new = UU;
+
+ # Picture settings
+ fig, ([axDensity, axMomentum],[axVelocity, axPressure]) = plt.subplots(2, 2,sharex=True, figsize=(10,10))
+ fig.suptitle('Explicit Upwind scheme')
+ lineDensity, = axDensity.plot([xmin+0.5*dx + i*dx for i in range(nx)], rho_field, label='Density') #new picture for video # Returns a tuple of line objects, thus the comma
+ axDensity.set(xlabel='x (m)', ylabel='Density')
+ axDensity.set_xlim(xmin,xmax)
+ axDensity.set_ylim(min_initial_rho - 0.1*(max_initial_rho-min_initial_rho), max_initial_rho + 0.1*(max_initial_rho-min_initial_rho) )
+ axDensity.legend()
+ lineMomentum, = axMomentum.plot([xmin+0.5*dx + i*dx for i in range(nx)], q_field, label='Momentum')
+ axMomentum.set(xlabel='x (m)', ylabel='Momentum')
+ axMomentum.set_xlim(xmin,xmax)
+ axMomentum.set_ylim(min_initial_q - 0.1*(max_initial_q-min_initial_q), max_initial_q + 0.1*(max_initial_q-min_initial_q) )
+ axMomentum.legend()
+ lineVelocity, = axVelocity.plot([xmin+0.5*dx + i*dx for i in range(nx)], v_field, label='Velocity')
+ axVelocity.set(xlabel='x (m)', ylabel='Velocity')
+ axVelocity.set_xlim(xmin,xmax)
+ axVelocity.set_ylim(min_initial_v - 0.4*abs(min_initial_v), max_initial_v + 0.05*abs(max_initial_v) )
+ axVelocity.legend()
+ linePressure, = axPressure.plot([xmin+0.5*dx + i*dx for i in range(nx)], p_field, label='Pressure')
+ axPressure.set(xlabel='x (m)', ylabel='Pressure')
+ axPressure.set_xlim(xmin,xmax)
+ axPressure.set_ylim(min_initial_p - 0.05*abs(min_initial_p), max_initial_p + 0.05*abs(max_initial_p) )
+ axPressure.ticklabel_format(axis='y', style='sci', scilimits=(0,0))
+ axPressure.legend()
+
+ # Video settings
+ FFMpegWriter = manimation.writers['ffmpeg']
+ metadata = dict(title="Upwind (Roe) scheme with entropic correction for the 1D isothermal Euler System", artist = "CEA Saclay", comment="Shock tube")
+ writer=FFMpegWriter(fps=10, metadata=metadata, codec='h264')
+ with writer.saving(fig, "1DEuler_System_Upwind"+".mp4", ntmax):
+ writer.grab_frame()
+ plt.savefig("EulerSystem"+str(dim)+"UpwindExplicit"+meshName+"_0"+".png")
+
+ # loop in time
+ while (it<ntmax and time <= tmax ):
+ # time step
+ dt=compTimeStep(UU,dx,CFL);
+ # Neumann boundary condition
+ UU_limL[0] = UU[0,0];
+ UU_limL[1] = UU[0,1];
+ UU_limR[0] = UU[M.getNumberOfCells()-1,0];
+ UU_limR[1] = UU[M.getNumberOfCells()-1,1];
+ flux_iminus = flux(UU_limL[0],UU_limL[1]);
+ #Main loop on cells
+ for i in range(0,M.getNumberOfCells()):
+ if (i<M.getNumberOfCells()-1):
+ Del_UU_RL[0]=UU[i+1,0]-UU[i,0];
+ Del_UU_RL[1]=UU[i+1,1]-UU[i,1];
+ AbsRoeMa = AbsRoeMatrix(UU[i,0],UU[i,1],UU[i+1,0],UU[i+1,1]);
+ flux_iplus = (flux(UU[i,0],UU[i,1])+flux(UU[i+1,0],UU[i+1,1]))-AbsRoeMa*Del_UU_RL - Del_UU_RL*abs(UU[i+1,1]/UU[i+1,0]-UU[i,1]/UU[i,0])#The last term is the entropic correction
+ flux_iplus[0]*=0.5;
+ flux_iplus[1]*=0.5;
+ else:
+ flux_iplus= flux(UU_limR[0],UU_limR[1]);
+ UU_new[i,0] = UU[i,0] - dt/dx*(flux_iplus[0]-flux_iminus[0]);
+ UU_new[i,1] = UU[i,1] - dt/dx*(flux_iplus[1]-flux_iminus[1]);
+ flux_iminus = flux_iplus;
+ pass
+ UU = UU_new;
+ time+=dt;
+ it+=1;
+
+ rho_field = [ UU[i,0] for i in range(nx)]
+ q_field = [ UU[i,1] for i in range(nx)]
+ v_field = [ UU[i,1]/UU[i,0] for i in range(nx)]
+ p_field = [ UU[i,0]*(c0*c0) for i in range(nx)]
+
+ lineDensity.set_ydata(rho_field)
+ lineMomentum.set_ydata(q_field)
+ lineVelocity.set_ydata(v_field)
+ linePressure.set_ydata(p_field)
+ writer.grab_frame()
+
+ if (it%freqSortie==0):
+ print( "-- Iter : ", it," Time : ",time," dt : ",dt)
+ plt.savefig("EulerSystem"+str(dim)+"UpwindExplicit"+meshName+"_"+str(it)+".png")
+ pass
+
+ print("-- Iter: " + str(it) + ", Time: " + str(time) + ", dt: " + str(dt))
+ if(it>=ntmax):
+ print("Nombre de pas de temps maximum ntmax= ", ntmax, " atteint")
+ return
+ elif(isStationary):
+ print("Régime stationnaire atteint au pas de temps ", it, ", t= ", time)
+ print("------------------------------------------------------------------------------------")
+ plt.savefig("EulerSystem"+str(dim)+"UpwindExplicit"+meshName+"_Stat.png")
+ return
+ else:
+ print("Temps maximum Tmax= ", tmax, " atteint")
+
+ return
+
+if __name__ == '__main__':
+ main("SquareRegularSquares")
--- /dev/null
+#!/usr/bin/env python3\r
+# -*-coding:utf-8 -*\r
+\r
+#===============================================================================================================================\r
+# Name : Résolution VF des équations d'Euler isotherme 2D sans terme source\r
+# \partial_t rho + \div q = 0\r
+# \partial_t q + \div q \otimes q/rho \grad p = 0\r
+# Author : Michaël Ndjinga, Coraline Mounier\r
+# Copyright : CEA Saclay 2020\r
+# Description : Propagation d'une onde de choc droite\r
+# Utilisation du schéma upwind explicite ou implicite sur un maillage général\r
+# Initialisation par une discontinuité verticale\r
+# Conditions aux limites de Neumann\r
+# Création et sauvegarde du champ résultant et des figures\r
+#================================================================================================================================\r
+\r
+\r
+import cdmath\r
+import numpy as np\r
+import matplotlib\r
+matplotlib.use("Agg")\r
+import matplotlib.pyplot as plt\r
+import matplotlib.animation as manimation\r
+from math import sqrt\r
+import sys\r
+import PV_routines\r
+\r
+p0=155.e5 #reference pressure in a pressurised nuclear vessel\r
+c0=700. #reference sound speed for water at 155 bars\r
+rho0=p0/(c0*c0) #reference density\r
+precision=1e-5\r
+\r
+def initial_conditions_Riemann_problem(my_mesh):\r
+ print( "Initial data : Riemann problem" )\r
+ dim = my_mesh.getMeshDimension()\r
+ nbCells = my_mesh.getNumberOfCells()\r
+\r
+ xcentre = 0.5\r
+\r
+ density_field = cdmath.Field("Density", cdmath.CELLS, my_mesh, 1)\r
+ q_x_field = cdmath.Field("Momentum x", cdmath.CELLS, my_mesh, 1)\r
+ q_y_field = cdmath.Field("Momentum y", cdmath.CELLS, my_mesh, 1)\r
+ #Velocity field with 3 components should be created for streamlines to be drawn\r
+ velocity_field = cdmath.Field("Velocity", cdmath.CELLS, my_mesh, 3)\r
+\r
+ for i in range(nbCells):\r
+ x = my_mesh.getCell(i).x()\r
+\r
+ # Initial momentum is zero\r
+ q_x_field[i] = 0\r
+ q_y_field[i] = 0\r
+\r
+ # Initial velocity is zero\r
+ velocity_field[i,0] = 0\r
+ velocity_field[i,1] = 0\r
+ velocity_field[i,2] = 0\r
+\r
+ if x < xcentre:\r
+ density_field[i] = rho0\r
+ pass\r
+ else:\r
+ density_field[i] = rho0/2\r
+ pass\r
+ pass\r
+\r
+ return density_field, q_x_field, q_y_field, velocity_field\r
+\r
+\r
+def jacobianMatricesm_x(coeff,rho_l,q_lx,q_ly,rho_r,q_rx,q_ry):\r
+ RoeMatx = cdmath.Matrix(3,3);\r
+ Dmacx = cdmath.Matrix(3,3);\r
+\r
+ u_lx=q_lx/rho_l\r
+ u_rx=q_rx/rho_r\r
+ u_ly=q_ly/rho_l\r
+ u_ry=q_ry/rho_r\r
+ if rho_l<0 or rho_r<0:\r
+ print("rho_l=",rho_l, " rho_r= ",rho_r)\r
+ raise ValueError("Negative density")\r
+ ux = (u_lx*sqrt(rho_l)+u_rx*sqrt(rho_r))/(sqrt(rho_l)+sqrt(rho_r));\r
+ uy = (u_ly*sqrt(rho_l)+u_ry*sqrt(rho_r))/(sqrt(rho_l)+sqrt(rho_r));\r
+\r
+ RoeMatx[0,0] = 0\r
+ RoeMatx[0,1] = 1\r
+ RoeMatx[0,2] = 0\r
+ RoeMatx[1,0] = c0*c0 - ux*ux\r
+ RoeMatx[1,1] = 2*ux\r
+ RoeMatx[1,2] = 0\r
+ RoeMatx[2,0] = -ux*uy\r
+ RoeMatx[2,1] = uy\r
+ RoeMatx[2,2] = ux\r
+\r
+ Dmacx[0,0] = abs(ux)-ux\r
+ Dmacx[0,1] = 1\r
+ Dmacx[0,2] = 0\r
+ Dmacx[1,0] = -c0*c0-ux*ux\r
+ Dmacx[1,1] = abs(ux)+ux\r
+ Dmacx[1,2] = 0\r
+ Dmacx[2,0] = -ux*uy\r
+ Dmacx[2,1] = uy\r
+ Dmacx[2,2] = abs(ux)\r
+\r
+ return (RoeMatx-Dmacx)*coeff*0.5\r
+\r
+def jacobianMatricesp_x(coeff,rho_l,q_lx,q_ly,rho_r,q_rx,q_ry):\r
+ RoeMatx = cdmath.Matrix(3,3)\r
+ Dmacx = cdmath.Matrix(3,3)\r
+\r
+ u_lx=q_lx/rho_l\r
+ u_rx=q_rx/rho_r\r
+ u_ly=q_ly/rho_l\r
+ u_ry=q_ry/rho_r\r
+ if rho_l<0 or rho_r<0:\r
+ print("rho_l=",rho_l, " rho_r= ",rho_r)\r
+ raise ValueError("Negative density")\r
+ ux = (u_lx*sqrt(rho_l)+u_rx*sqrt(rho_r))/(sqrt(rho_l)+sqrt(rho_r));\r
+ uy = (u_ly*sqrt(rho_l)+u_ry*sqrt(rho_r))/(sqrt(rho_l)+sqrt(rho_r));\r
+\r
+ RoeMatx[0,0] = 0\r
+ RoeMatx[0,1] = 1\r
+ RoeMatx[0,2] = 0\r
+ RoeMatx[1,0] = c0*c0 - ux*ux\r
+ RoeMatx[1,1] = 2*ux\r
+ RoeMatx[1,2] = 0\r
+ RoeMatx[2,0] = -ux*uy\r
+ RoeMatx[2,1] = uy\r
+ RoeMatx[2,2] = ux\r
+\r
+ Dmacx[0,0] = abs(ux)-ux\r
+ Dmacx[0,1] = 1\r
+ Dmacx[0,2] = 0\r
+ Dmacx[1,0] = -c0*c0-ux*ux\r
+ Dmacx[1,1] = abs(ux)+ux\r
+ Dmacx[1,2] = 0\r
+ Dmacx[2,0] = -ux*uy\r
+ Dmacx[2,1] = uy\r
+ Dmacx[2,2] = abs(ux)\r
+\r
+ return (RoeMatx+Dmacx)*coeff*0.5\r
+\r
+def jacobianMatricesm_y(coeff,rho_l,q_lx,q_ly,rho_r,q_rx,q_ry):\r
+ RoeMaty = cdmath.Matrix(3,3);\r
+ Dmacy = cdmath.Matrix(3,3);\r
+\r
+ u_lx=q_lx/rho_l\r
+ u_rx=q_rx/rho_r\r
+ u_ly=q_ly/rho_l\r
+ u_ry=q_ry/rho_r\r
+ if rho_l<0 or rho_r<0:\r
+ print("rho_l=",rho_l, " rho_r= ",rho_r)\r
+ raise ValueError("Negative density")\r
+ ux = (u_lx*sqrt(rho_l)+u_rx*sqrt(rho_r))/(sqrt(rho_l)+sqrt(rho_r));\r
+ uy = (u_ly*sqrt(rho_l)+u_ry*sqrt(rho_r))/(sqrt(rho_l)+sqrt(rho_r));\r
+\r
+ RoeMaty[0,0] = 0\r
+ RoeMaty[0,1] = 0\r
+ RoeMaty[0,2] = 1\r
+ RoeMaty[1,0] = -ux*uy\r
+ RoeMaty[1,1] = uy\r
+ RoeMaty[1,2] = ux\r
+ RoeMaty[2,0] = c0*c0-uy*uy\r
+ RoeMaty[2,1] = 0\r
+ RoeMaty[2,2] = 2*uy\r
+\r
+ Dmacy[0,0] = abs(uy)-uy\r
+ Dmacy[0,1] = 0\r
+ Dmacy[0,2] = 1\r
+ Dmacy[1,0] = -ux*uy\r
+ Dmacy[1,1] = abs(uy)\r
+ Dmacy[1,2] = ux\r
+ Dmacy[2,0] = -c0*c0-uy*uy\r
+ Dmacy[2,1] = 0\r
+ Dmacy[2,2] = abs(uy)+uy\r
+\r
+ return (RoeMaty-Dmacy)*coeff*0.5\r
+\r
+def jacobianMatricesp_y(coeff,rho_l,q_lx,q_ly,rho_r,q_rx,q_ry):\r
+ RoeMaty = cdmath.Matrix(3,3);\r
+ Dmacy = cdmath.Matrix(3,3);\r
+\r
+ u_lx=q_lx/rho_l\r
+ u_rx=q_rx/rho_r\r
+ u_ly=q_ly/rho_l\r
+ u_ry=q_ry/rho_r\r
+ if rho_l<0 or rho_r<0:\r
+ print("rho_l=",rho_l, " rho_r= ",rho_r)\r
+ raise ValueError("Negative density")\r
+ ux = (u_lx*sqrt(rho_l)+u_rx*sqrt(rho_r))/(sqrt(rho_l)+sqrt(rho_r));\r
+ uy = (u_ly*sqrt(rho_l)+u_ry*sqrt(rho_r))/(sqrt(rho_l)+sqrt(rho_r));\r
+\r
+ RoeMaty[0,0] = 0\r
+ RoeMaty[0,1] = 0\r
+ RoeMaty[0,2] = 1\r
+ RoeMaty[1,0] = -ux*uy\r
+ RoeMaty[1,1] = uy\r
+ RoeMaty[1,2] = ux\r
+ RoeMaty[2,0] = c0*c0-uy*uy\r
+ RoeMaty[2,1] = 0\r
+ RoeMaty[2,2] = 2*uy\r
+\r
+ Dmacy[0,0] = abs(uy)-uy\r
+ Dmacy[0,1] = 0\r
+ Dmacy[0,2] = 1\r
+ Dmacy[1,0] = -ux*uy\r
+ Dmacy[1,1] = abs(uy)\r
+ Dmacy[1,2] = ux\r
+ Dmacy[2,0] = -c0*c0-uy*uy\r
+ Dmacy[2,1] = 0\r
+ Dmacy[2,2] = abs(uy)+uy\r
+\r
+ return (RoeMaty+Dmacy)*coeff*0.5\r
+\r
+def EulerSystemStaggered(ntmax, tmax, cfl,output_freq, my_mesh, meshName):\r
+\r
+ if not my_mesh.isStructured() :\r
+ raise ValueError("Euler_ConservativeStaggered2D_RiemannProblem.py expects a structured mesh")\r
+\r
+ dim=my_mesh.getMeshDimension()\r
+ if dim != 2 :\r
+ raise ValueError("Euler_ConservativeStaggered2D_RiemannProblem.py expects a 2D mesh")\r
+\r
+ nbComp=dim+1\r
+ # Mesh parameters\r
+ nbCells = my_mesh.getNumberOfCells()\r
+ nbCells_x = my_mesh.getNx()\r
+ nbCells_y = my_mesh.getNy()\r
+ dx,dy = my_mesh.getDXYZ()\r
+ nbVoisinsMax=my_mesh.getMaxNbNeighbours(cdmath.CELLS)\r
+ dx_min=my_mesh.minRatioVolSurf()\r
+ \r
+ # Time evolution parameters\r
+ time = 0.\r
+ it=0;\r
+ isStationary=False\r
+ iterGMRESMax = 50\r
+ newton_max = 50\r
+\r
+ #iteration vectors\r
+ Un =cdmath.Vector(nbCells*nbComp)\r
+ dUn=cdmath.Vector(nbCells*nbComp)\r
+ dUk=cdmath.Vector(nbCells*nbComp)\r
+ Rhs=cdmath.Vector(nbCells*nbComp)\r
+\r
+ dUi_x=cdmath.Vector(nbComp)\r
+ dUi_y=cdmath.Vector(nbComp)\r
+ dUi1_x=cdmath.Vector(nbComp)\r
+ dUi2_x=cdmath.Vector(nbComp)\r
+ dUi1_y=cdmath.Vector(nbComp)\r
+ dUi2_y=cdmath.Vector(nbComp)\r
+\r
+ # Initial conditions #\r
+ print("Construction of the initial condition …")\r
+ rho_field, q_x_field, q_y_field, velocity_field= initial_conditions_Riemann_problem(my_mesh)\r
+\r
+ for i in range(nbCells):\r
+ Un[nbComp*i] = rho_field[i]\r
+ Un[nbComp*i+1] = q_x_field[i]\r
+ Un[nbComp*i+2] = q_y_field[i]\r
+\r
+ #Sauvegarde de la donnée initiale\r
+ rho_field.setTime(time,it);\r
+ rho_field.writeVTK("EulerIsothermal_"+str(dim)+"DConservativeStaggered_"+meshName+"_density");\r
+ q_x_field.setTime(time,it);\r
+ q_x_field.writeVTK("EulerIsothermal_"+str(dim)+"DConservativeStaggered_"+meshName+"_momentumX");\r
+ q_y_field.setTime(time,it);\r
+ q_y_field.writeVTK("EulerIsothermal_"+str(dim)+"DConservativeStaggered_"+meshName+"_momentumY");\r
+ velocity_field.setTime(time,it);\r
+ velocity_field.writeVTK("EulerIsothermal_"+str(dim)+"DConservativeStaggered_"+meshName+"_velocity");\r
+\r
+ print("Starting computation of the isothermal Euler system with a conservative staggered scheme …")\r
+ divMat=cdmath.SparseMatrixPetsc(nbCells*nbComp,nbCells*nbComp,(nbVoisinsMax+1)*nbComp)\r
+\r
+ # Starting time loop\r
+ while (it<ntmax and time <= tmax and not isStationary):\r
+ dUn = Un.deepCopy()\r
+ Uk = Un.deepCopy()\r
+ residu = 1e10\r
+ k=0\r
+ while (k<newton_max and residu > 1/precision ):\r
+\r
+ 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
+ #DEBUT BOUCLE NEWTON\r
+ for j in range(nbCells_y):\r
+ for i in range(nbCells_x):\r
+ #Traitement des coins\r
+ if ( j==0 and i==0) :\r
+ #Calcul de Am_x\r
+ rho_l=Uk[3*nbCells_x*j+3*i]\r
+ qx_l =Uk[3*nbCells_x*j+3*i+1]\r
+ qy_l =Uk[3*nbCells_x*j+3*i+2]\r
+ rho_r=Uk[3*nbCells_x*j+3*(i+1)]\r
+ qx_r =Uk[3*nbCells_x*j+3*(i+1)+1]\r
+ qy_r =Uk[3*nbCells_x*j+3*(i+1)+2]\r
+ Am_x= jacobianMatricesm_x(dt/dx,rho_l,qx_l,qy_l,rho_r,qx_r,qy_r)\r
+\r
+ #calcul de Am_y\r
+ rho_l=Uk[3*nbCells_x*j+3*i]\r
+ qx_l =Uk[3*nbCells_x*j+3*i+1]\r
+ qy_l =Uk[3*nbCells_x*j+3*i+2]\r
+ rho_r=Uk[3*nbCells_x*(j+1)+3*i]\r
+ qx_r =Uk[3*nbCells_x*(j+1)+3*i+1]\r
+ qy_r =Uk[3*nbCells_x*(j+1)+3*i+2]\r
+ Am_y= jacobianMatricesm_y(dt/dx,rho_l,qx_l,qy_l,rho_r,qx_r,qy_r)\r
+\r
+ divMat.addValue(3*nbCells_x*j+3*i,3*nbCells_x*j+3*(i+1),Am_x)\r
+ divMat.addValue(3*nbCells_x*j+3*i,3*nbCells_x*(j+1)+3*i,Am_y)\r
+ divMat.addValue(3*nbCells_x*j+3*i,3*nbCells_x*j+3*i,Am_x*(-1.)-Am_y)\r
+\r
+ dUi_x[0] = Uk[3*nbCells_x*j+3*(i+1)]-Uk[3*nbCells_x*j+3*i]\r
+ dUi_x[1] = Uk[3*nbCells_x*j+3*(i+1)+1]-Uk[3*nbCells_x*j+3*i+1]\r
+ dUi_x[2] = Uk[3*nbCells_x*j+3*(i+1)+2]-Uk[3*nbCells_x*j+3*i+2]\r
+ dUi_y[0] = Uk[3*nbCells_x*(j+1)+3*i]-Uk[3*nbCells_x*j+3*i]\r
+ dUi_y[1] = Uk[3*nbCells_x*(j+1)+3*i+1]-Uk[3*nbCells_x*j+3*i+1]\r
+ dUi_y[2] = Uk[3*nbCells_x*(j+1)+3*i+2]-Uk[3*nbCells_x*j+3*i+2]\r
+\r
+ temp_x = Am_x*dUi_x\r
+ temp_y = Am_y*dUi_y\r
+ #print("Bloc 0 matrix : ", Am)\r
+ 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
+ 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
+ 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
+\r
+ elif(i==0 and j==nbCells_y-1):\r
+ #Calcul de Am_x\r
+ rho_l=Uk[3*nbCells_x*j+3*i]\r
+ qx_l =Uk[3*nbCells_x*j+3*i+1]\r
+ qy_l =Uk[3*nbCells_x*j+3*i+2]\r
+ rho_r=Uk[3*nbCells_x*j+3*(i+1)]\r
+ qx_r =Uk[3*nbCells_x*j+3*(i+1)+1]\r
+ qy_r =Uk[3*nbCells_x*j+3*(i+1)+2]\r
+ Am_x= jacobianMatricesm_x(dt/dx,rho_l,qx_l,qy_l,rho_r,qx_r,qy_r)\r
+\r
+ #calcul de Ap_y\r
+ rho_l=Uk[3*nbCells_x*(j-1)+3*i]\r
+ qx_l =Uk[3*nbCells_x*(j-1)+3*i+1]\r
+ qy_l =Uk[3*nbCells_x*(j-1)+3*i+2]\r
+ rho_r=Uk[3*nbCells_x*j+3*i]\r
+ qx_r =Uk[3*nbCells_x*j+3*i+1]\r
+ qy_r =Uk[3*nbCells_x*j+3*i+2]\r
+ Ap_y= jacobianMatricesp_y(dt/dx,rho_l,qx_l,qy_l,rho_r,qx_r,qy_r)\r
+\r
+ divMat.addValue(3*nbCells_x*j+3*i,3*nbCells_x*j+3*(i+1),Am_x)\r
+ divMat.addValue(3*nbCells_x*j+3*i,3*nbCells_x*(j-1)+3*i,Ap_y*(-1.))\r
+ divMat.addValue(3*nbCells_x*j+3*i,3*nbCells_x*j+3*i,Am_x*(-1.)+Ap_y)\r
+\r
+ dUi_x[0] = Uk[3*nbCells_x*j+3*(i+1)]-Uk[3*nbCells_x*j+3*i]\r
+ dUi_x[1] = Uk[3*nbCells_x*j+3*(i+1)+1]-Uk[3*nbCells_x*j+3*i+1]\r
+ dUi_x[2] = Uk[3*nbCells_x*j+3*(i+1)+2]-Uk[3*nbCells_x*j+3*i+2]\r
+ dUi_y[0] = Uk[3*nbCells_x*j+3*i]-Uk[3*nbCells_x*(j-1)+3*i]\r
+ dUi_y[1] = Uk[3*nbCells_x*j+3*i+1]-Uk[3*nbCells_x*(j-1)+3*i+1]\r
+ dUi_y[2] = Uk[3*nbCells_x*j+3*i+2]-Uk[3*nbCells_x*(j-1)+3*i+2]\r
+\r
+ temp_x = Am_x*dUi_x\r
+ temp_y = Ap_y*dUi_y\r
+ #print("Bloc 0 matrix : ", Am)\r
+ 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
+ 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
+ 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
+\r
+ elif(i==nbCells_x-1 and j==0):\r
+\r
+ #Calcul de Ap_x\r
+ rho_l=Uk[3*nbCells_x*j+3*(i-1)]\r
+ qx_l =Uk[3*nbCells_x*j+3*(i-1)+1]\r
+ qy_l =Uk[3*nbCells_x*j+3*(i-1)+2]\r
+ rho_r=Uk[3*nbCells_x*j+3*i]\r
+ qx_r =Uk[3*nbCells_x*j+3*i+1]\r
+ qy_r =Uk[3*nbCells_x*j+3*i+2]\r
+ Ap_x= jacobianMatricesp_x(dt/dx,rho_l,qx_l,qy_l,rho_r,qx_r,qy_r)\r
+\r
+ #calcul de Am_y\r
+ rho_l=Uk[3*nbCells_x*j+3*i]\r
+ qx_l =Uk[3*nbCells_x*j+3*i+1]\r
+ qy_l =Uk[3*nbCells_x*j+3*i+2]\r
+ rho_r=Uk[3*nbCells_x*(j+1)+3*i]\r
+ qx_r =Uk[3*nbCells_x*(j+1)+3*i+1]\r
+ qy_r =Uk[3*nbCells_x*(j+1)+3*i+2]\r
+ Am_y= jacobianMatricesm_y(dt/dx,rho_l,qx_l,qy_l,rho_r,qx_r,qy_r)\r
+\r
+ divMat.addValue(3*nbCells_x*j+3*i,3*nbCells_x*j+3*(i-1),Ap_x*(-1))\r
+ divMat.addValue(3*nbCells_x*j+3*i,3*nbCells_x*(j+1)+3*i,Am_y)\r
+ divMat.addValue(3*nbCells_x*j+3*i,3*nbCells_x*j+3*i,Ap_x-Am_y)\r
+\r
+ dUi_x[0] = Uk[3*nbCells_x*j+3*i]-Uk[3*nbCells_x*j+3*(i-1)]\r
+ dUi_x[1] = Uk[3*nbCells_x*j+3*i+1]-Uk[3*nbCells_x*j+3*(i-1)+1]\r
+ dUi_x[2] = Uk[3*nbCells_x*j+3*i+2]-Uk[3*nbCells_x*j+3*(i-1)+2]\r
+ dUi_y[0] = Uk[3*nbCells_x*(j+1)+3*i]-Uk[3*nbCells_x*j+3*i]\r
+ dUi_y[1] = Uk[3*nbCells_x*(j+1)+3*i+1]-Uk[3*nbCells_x*j+3*i+1]\r
+ dUi_y[2] = Uk[3*nbCells_x*(j+1)+3*i+2]-Uk[3*nbCells_x*j+3*i+2]\r
+\r
+ temp_x = Ap_x*dUi_x\r
+ temp_y= Am_y*dUi_y\r
+ #print("Bloc 0 matrix : ", Am)\r
+ 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
+ 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
+ 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
+\r
+ elif ( j==nbCells_y-1 and i==nbCells_x-1) :\r
+\r
+ #Calcul de Ap_x\r
+ rho_l=Uk[3*nbCells_x*j+3*(i-1)]\r
+ qx_l =Uk[3*nbCells_x*j+3*(i-1)+1]\r
+ qy_l =Uk[3*nbCells_x*j+3*(i-1)+2]\r
+ rho_r=Uk[3*nbCells_x*j+3*i]\r
+ qx_r =Uk[3*nbCells_x*j+3*i+1]\r
+ qy_r =Uk[3*nbCells_x*j+3*i+2]\r
+ Ap_x= jacobianMatricesp_x(dt/dx,rho_l,qx_l,qy_l,rho_r,qx_r,qy_r)\r
+\r
+ #calcul de Ap_y\r
+ rho_l=Uk[3*nbCells_x*(j-1)+3*i]\r
+ qx_l =Uk[3*nbCells_x*(j-1)+3*i+1]\r
+ qy_l =Uk[3*nbCells_x*(j-1)+3*i+2]\r
+ rho_r=Uk[3*nbCells_x*j+3*i]\r
+ qx_r =Uk[3*nbCells_x*j+3*i+1]\r
+ qy_r =Uk[3*nbCells_x*j+3*i+2]\r
+ Ap_y= jacobianMatricesp_y(dt/dx,rho_l,qx_l,qy_l,rho_r,qx_r,qy_r)\r
+\r
+ divMat.addValue(3*nbCells_x*j+3*i,3*nbCells_x*j+3*(i-1),Ap_x*(-1.))\r
+ divMat.addValue(3*nbCells_x*j+3*i,3*nbCells_x*(j-1)+3*i,Ap_y*(-1.))\r
+ divMat.addValue(3*nbCells_x*j+3*i,3*nbCells_x*j+3*i,Ap_x+Ap_y)\r
+\r
+ dUi_x[0] = Uk[3*nbCells_x*j+3*i]-Uk[3*nbCells_x*j+3*(i-1)]\r
+ dUi_x[1] = Uk[3*nbCells_x*j+3*i+1]-Uk[3*nbCells_x*j+3*(i-1)+1]\r
+ dUi_x[2] = Uk[3*nbCells_x*j+3*i+2]-Uk[3*nbCells_x*j+3*(i-1)+2]\r
+ dUi_y[0] = Uk[3*nbCells_x*j+3*i]-Uk[3*nbCells_x*(j-1)+3*i]\r
+ dUi_y[1] = Uk[3*nbCells_x*j+3*i+1]-Uk[3*nbCells_x*(j-1)+3*i+1]\r
+ dUi_y[2] = Uk[3*nbCells_x*j+3*i+2]-Uk[3*nbCells_x*(j-1)+3*i+2]\r
+\r
+ temp_x = Ap_x*dUi_x\r
+ temp_y = Ap_y*dUi_y\r
+ 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
+ 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
+ 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
+\r
+ #Traitement des bords restants\r
+ elif i==0 :\r
+\r
+ #Calcul de Am_x\r
+ rho_l=Uk[3*nbCells_x*j+3*i]\r
+ qx_l =Uk[3*nbCells_x*j+3*i+1]\r
+ qy_l =Uk[3*nbCells_x*j+3*i+2]\r
+ rho_r=Uk[3*nbCells_x*j+3*(i+1)]\r
+ qx_r =Uk[3*nbCells_x*j+3*(i+1)+1]\r
+ qy_r =Uk[3*nbCells_x*j+3*(i+1)+2]\r
+ Am_x= jacobianMatricesm_x(dt/dx,rho_l,qx_l,qy_l,rho_r,qx_r,qy_r)\r
+\r
+ #calcul de Am_y\r
+ rho_l=Uk[3*nbCells_x*j+3*i]\r
+ qx_l =Uk[3*nbCells_x*j+3*i+1]\r
+ qy_l =Uk[3*nbCells_x*j+3*i+2]\r
+ rho_r=Uk[3*nbCells_x*(j+1)+3*i]\r
+ qx_r =Uk[3*nbCells_x*(j+1)+3*i+1]\r
+ qy_r =Uk[3*nbCells_x*(j+1)+3*i+2]\r
+ Am_y= jacobianMatricesm_y(dt/dx,rho_l,qx_l,qy_l,rho_r,qx_r,qy_r)\r
+\r
+ #calcul de Ap_y\r
+ rho_l=Uk[3*nbCells_x*(j-1)+3*i]\r
+ qx_l =Uk[3*nbCells_x*(j-1)+3*i+1]\r
+ qy_l =Uk[3*nbCells_x*(j-1)+3*i+2]\r
+ rho_r=Uk[3*nbCells_x*j+3*i]\r
+ qx_r =Uk[3*nbCells_x*j+3*i+1]\r
+ qy_r =Uk[3*nbCells_x*j+3*i+2]\r
+ Ap_y= jacobianMatricesp_y(dt/dx,rho_l,qx_l,qy_l,rho_r,qx_r,qy_r)\r
+\r
+ #remplissage de la divMat\r
+ divMat.addValue(3*nbCells_x*j+3*i,3*nbCells_x*j+3*(i+1),Am_x)\r
+ divMat.addValue(3*nbCells_x*j+3*i,3*nbCells_x*(j+1)+3*i,Am_y)\r
+ divMat.addValue(3*nbCells_x*j+3*i,3*nbCells_x*(j-1)+3*i,Ap_y*(-1.))\r
+ divMat.addValue(3*nbCells_x*j+3*i,3*nbCells_x*j+3*i,Ap_y-Am_x-Am_y)\r
+\r
+ #remplissage du membre de droite\r
+ dUi_x[0] = Uk[3*nbCells_x*j+3*(i+1)]-Uk[3*nbCells_x*j+3*i]\r
+ dUi_x[1] = Uk[3*nbCells_x*j+3*(i+1)+1]-Uk[3*nbCells_x*j+3*i+1]\r
+ dUi_x[2] = Uk[3*nbCells_x*j+3*(i+1)+2]-Uk[3*nbCells_x*j+3*i+2]\r
+ dUi1_y[0] = Uk[3*nbCells_x*(j+1)+3*i]-Uk[3*nbCells_x*j+3*i]\r
+ dUi1_y[1] = Uk[3*nbCells_x*(j+1)+3*i+1]-Uk[3*nbCells_x*j+3*i+1]\r
+ dUi1_y[2] = Uk[3*nbCells_x*(j+1)+3*i+2]-Uk[3*nbCells_x*j+3*i+2]\r
+ dUi2_y[0] = Uk[3*nbCells_x*j+3*i]-Uk[3*nbCells_x*(j-1)+3*i]\r
+ dUi2_y[1] = Uk[3*nbCells_x*j+3*i+1]-Uk[3*nbCells_x*(j-1)+3*i+1]\r
+ dUi2_y[2] = Uk[3*nbCells_x*j+3*i+2]-Uk[3*nbCells_x*(j-1)+3*i+2]\r
+\r
+ temp_x = Am_x*dUi_x\r
+ temp1_y = Am_y*dUi1_y\r
+ temp2_y = Ap_y*dUi2_y\r
+ 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
+ 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
+ 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
+\r
+ elif i==nbCells_x-1:\r
+\r
+ #Calcul de Ap_x\r
+ rho_l=Uk[3*nbCells_x*j+3*(i-1)]\r
+ qx_l =Uk[3*nbCells_x*j+3*(i-1)+1]\r
+ qy_l =Uk[3*nbCells_x*j+3*(i-1)+2]\r
+ rho_r=Uk[3*nbCells_x*j+3*i]\r
+ qx_r =Uk[3*nbCells_x*j+3*i+1]\r
+ qy_r =Uk[3*nbCells_x*j+3*i+2]\r
+ Ap_x= jacobianMatricesp_x(dt/dx,rho_l,qx_l,qy_l,rho_r,qx_r,qy_r)\r
+\r
+ #calcul de Am_y\r
+ rho_l=Uk[3*nbCells_x*j+3*i]\r
+ qx_l =Uk[3*nbCells_x*j+3*i+1]\r
+ qy_l =Uk[3*nbCells_x*j+3*i+2]\r
+ rho_r=Uk[3*nbCells_x*(j+1)+3*i]\r
+ qx_r =Uk[3*nbCells_x*(j+1)+3*i+1]\r
+ qy_r =Uk[3*nbCells_x*(j+1)+3*i+2]\r
+ Am_y= jacobianMatricesm_y(dt/dx,rho_l,qx_l,qy_l,rho_r,qx_r,qy_r)\r
+\r
+ #calcul de Ap_y\r
+ rho_l=Uk[3*nbCells_x*(j-1)+3*i]\r
+ qx_l =Uk[3*nbCells_x*(j-1)+3*i+1]\r
+ qy_l =Uk[3*nbCells_x*(j-1)+3*i+2]\r
+ rho_r=Uk[3*nbCells_x*j+3*i]\r
+ qx_r =Uk[3*nbCells_x*j+3*i+1]\r
+ qy_r =Uk[3*nbCells_x*j+3*i+2]\r
+ Ap_y= jacobianMatricesp_y(dt/dx,rho_l,qx_l,qy_l,rho_r,qx_r,qy_r)\r
+\r
+ #remplissage de la divMat\r
+ divMat.addValue(3*nbCells_x*j+3*i,3*nbCells_x*j+3*(i-1),Ap_x*(-1.))\r
+ divMat.addValue(3*nbCells_x*j+3*i,3*nbCells_x*(j+1)+3*i,Am_y)\r
+ divMat.addValue(3*nbCells_x*j+3*i,3*nbCells_x*(j-1)+3*i,Ap_y*(-1.))\r
+ divMat.addValue(3*nbCells_x*j+3*i,3*nbCells_x*j+3*i,Ap_y+Ap_x-Am_y)\r
+\r
+ #remplissage du membre de droite\r
+ dUi_x[0] = Uk[3*nbCells_x*j+3*i]-Uk[3*nbCells_x*j+3*(i-1)]\r
+ dUi_x[1] = Uk[3*nbCells_x*j+3*i+1]-Uk[3*nbCells_x*j+3*(i-1)+1]\r
+ dUi_x[2] = Uk[3*nbCells_x*j+3*i+2]-Uk[3*nbCells_x*j+3*(i-1)+2]\r
+ dUi1_y[0] = Uk[3*nbCells_x*(j+1)+3*i]-Uk[3*nbCells_x*j+3*i]\r
+ dUi1_y[1] = Uk[3*nbCells_x*(j+1)+3*i+1]-Uk[3*nbCells_x*j+3*i+1]\r
+ dUi1_y[2] = Uk[3*nbCells_x*(j+1)+3*i+2]-Uk[3*nbCells_x*j+3*i+2]\r
+ dUi2_y[0] = Uk[3*nbCells_x*j+3*i]-Uk[3*nbCells_x*(j-1)+3*i]\r
+ dUi2_y[1] = Uk[3*nbCells_x*j+3*i+1]-Uk[3*nbCells_x*(j-1)+3*i+1]\r
+ dUi2_y[2] = Uk[3*nbCells_x*j+3*i+2]-Uk[3*nbCells_x*(j-1)+3*i+2]\r
+\r
+ temp_x = Ap_x*dUi_x\r
+ temp1_y = Am_y*dUi1_y\r
+ temp2_y = Ap_y*dUi2_y\r
+ 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
+ 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
+ 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
+\r
+ elif j==0:\r
+\r
+ #Calcul de Am_x\r
+ rho_l=Uk[3*nbCells_x*j+3*i]\r
+ qx_l =Uk[3*nbCells_x*j+3*i+1]\r
+ qy_l =Uk[3*nbCells_x*j+3*i+2]\r
+ rho_r=Uk[3*nbCells_x*j+3*(i+1)]\r
+ qx_r =Uk[3*nbCells_x*j+3*(i+1)+1]\r
+ qy_r =Uk[3*nbCells_x*j+3*(i+1)+2]\r
+ Am_x= jacobianMatricesm_x(dt/dx,rho_l,qx_l,qy_l,rho_r,qx_r,qy_r)\r
+\r
+ #Calcul de Ap_x\r
+ rho_l=Uk[3*nbCells_x*j+3*(i-1)]\r
+ qx_l =Uk[3*nbCells_x*j+3*(i-1)+1]\r
+ qy_l =Uk[3*nbCells_x*j+3*(i-1)+2]\r
+ rho_r=Uk[3*nbCells_x*j+3*i]\r
+ qx_r =Uk[3*nbCells_x*j+3*i+1]\r
+ qy_r =Uk[3*nbCells_x*j+3*i+2]\r
+ Ap_x= jacobianMatricesp_x(dt/dx,rho_l,qx_l,qy_l,rho_r,qx_r,qy_r)\r
+\r
+ #calcul de Am_y\r
+ rho_l=Uk[3*nbCells_x*j+3*i]\r
+ qx_l =Uk[3*nbCells_x*j+3*i+1]\r
+ qy_l =Uk[3*nbCells_x*j+3*i+2]\r
+ rho_r=Uk[3*nbCells_x*(j+1)+3*i]\r
+ qx_r =Uk[3*nbCells_x*(j+1)+3*i+1]\r
+ qy_r =Uk[3*nbCells_x*(j+1)+3*i+2]\r
+ Am_y= jacobianMatricesm_y(dt/dx,rho_l,qx_l,qy_l,rho_r,qx_r,qy_r)\r
+\r
+ #remplissage de la divMat\r
+ divMat.addValue(3*nbCells_x*j+3*i,3*nbCells_x*j+3*(i-1),Ap_x*(-1.))\r
+ divMat.addValue(3*nbCells_x*j+3*i,3*nbCells_x*j+3*(i+1),Am_x)\r
+ divMat.addValue(3*nbCells_x*j+3*i,3*nbCells_x*(j+1)+3*i,Am_y)\r
+ divMat.addValue(3*nbCells_x*j+3*i,3*nbCells_x*j+3*i,Ap_x-Am_x-Am_y)\r
+\r
+ #remplissage du membre de droite\r
+ dUi1_x[0] = Uk[3*nbCells_x*j+3*(i+1)]-Uk[3*nbCells_x*j+3*i]\r
+ dUi1_x[1] = Uk[3*nbCells_x*j+3*(i+1)+1]-Uk[3*nbCells_x*j+3*i+1]\r
+ dUi1_x[2] = Uk[3*nbCells_x*j+3*(i+1)+2]-Uk[3*nbCells_x*j+3*i+2]\r
+ dUi2_x[0] = Uk[3*nbCells_x*j+3*i]-Uk[3*nbCells_x*j+3*(i-1)]\r
+ dUi2_x[1] = Uk[3*nbCells_x*j+3*i+1]-Uk[3*nbCells_x*j+3*(i-1)+1]\r
+ dUi2_x[2] = Uk[3*nbCells_x*j+3*i+2]-Uk[3*nbCells_x*j+3*(i-1)+2]\r
+ dUi_y[0] = Uk[3*nbCells_x*(j+1)+3*i]-Uk[3*nbCells_x*j+3*i]\r
+ dUi_y[1] = Uk[3*nbCells_x*(j+1)+3*i+1]-Uk[3*nbCells_x*j+3*i+1]\r
+ dUi_y[2] = Uk[3*nbCells_x*(j+1)+3*i+2]-Uk[3*nbCells_x*j+3*i+2]\r
+\r
+ temp1_x = Am_x*dUi1_x\r
+ temp2_x = Ap_x*dUi2_x\r
+ temp_y = Am_y*dUi_y\r
+ 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
+ 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
+ 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
+\r
+ elif j==nbCells_y-1:\r
+\r
+ #Calcul de Am_x\r
+ rho_l=Uk[3*nbCells_x*j+3*i]\r
+ qx_l =Uk[3*nbCells_x*j+3*i+1]\r
+ qy_l =Uk[3*nbCells_x*j+3*i+2]\r
+ rho_r=Uk[3*nbCells_x*j+3*(i+1)]\r
+ qx_r =Uk[3*nbCells_x*j+3*(i+1)+1]\r
+ qy_r =Uk[3*nbCells_x*j+3*(i+1)+2]\r
+ Am_x= jacobianMatricesm_x(dt/dx,rho_l,qx_l,qy_l,rho_r,qx_r,qy_r)\r
+\r
+ #Calcul de Ap_x\r
+ rho_l=Uk[3*nbCells_x*j+3*(i-1)]\r
+ qx_l =Uk[3*nbCells_x*j+3*(i-1)+1]\r
+ qy_l =Uk[3*nbCells_x*j+3*(i-1)+2]\r
+ rho_r=Uk[3*nbCells_x*j+3*i]\r
+ qx_r =Uk[3*nbCells_x*j+3*i+1]\r
+ qy_r =Uk[3*nbCells_x*j+3*i+2]\r
+ Ap_x= jacobianMatricesp_x(dt/dx,rho_l,qx_l,qy_l,rho_r,qx_r,qy_r)\r
+\r
+ #calcul de Ap_y\r
+ rho_l=Uk[3*nbCells_x*(j-1)+3*i]\r
+ qx_l =Uk[3*nbCells_x*(j-1)+3*i+1]\r
+ qy_l =Uk[3*nbCells_x*(j-1)+3*i+2]\r
+ rho_r=Uk[3*nbCells_x*j+3*i]\r
+ qx_r =Uk[3*nbCells_x*j+3*i+1]\r
+ qy_r =Uk[3*nbCells_x*j+3*i+2]\r
+ Ap_y= jacobianMatricesp_y(dt/dx,rho_l,qx_l,qy_l,rho_r,qx_r,qy_r)\r
+\r
+ #remplissage de la divMat\r
+ divMat.addValue(3*nbCells_x*j+3*i,3*nbCells_x*j+3*(i-1),Ap_x*(-1.))\r
+ divMat.addValue(3*nbCells_x*j+3*i,3*nbCells_x*j+3*(i+1),Am_x)\r
+ divMat.addValue(3*nbCells_x*j+3*i,3*nbCells_x*(j-1)+3*i,Ap_y*(-1.))\r
+ divMat.addValue(3*nbCells_x*j+3*i,3*nbCells_x*j+3*i,Ap_x+Ap_y-Am_x)\r
+\r
+ #remplissage du membre de droite\r
+ dUi1_x[0] = Uk[3*nbCells_x*j+3*(i+1)]-Uk[3*nbCells_x*j+3*i]\r
+ dUi1_x[1] = Uk[3*nbCells_x*j+3*(i+1)+1]-Uk[3*nbCells_x*j+3*i+1]\r
+ dUi1_x[2] = Uk[3*nbCells_x*j+3*(i+1)+2]-Uk[3*nbCells_x*j+3*i+2]\r
+ dUi2_x[0] = Uk[3*nbCells_x*j+3*i]-Uk[3*nbCells_x*j+3*(i-1)]\r
+ dUi2_x[1] = Uk[3*nbCells_x*j+3*i+1]-Uk[3*nbCells_x*j+3*(i-1)+1]\r
+ dUi2_x[2] = Uk[3*nbCells_x*j+3*i+2]-Uk[3*nbCells_x*j+3*(i-1)+2]\r
+ dUi_y[0] = Uk[3*nbCells_x*j+3*i]-Uk[3*nbCells_x*(j-1)+3*i]\r
+ dUi_y[1] = Uk[3*nbCells_x*j+3*i+1]-Uk[3*nbCells_x*(j-1)+3*i+1]\r
+ dUi_y[2] = Uk[3*nbCells_x*j+3*i+2]-Uk[3*nbCells_x*(j-1)+3*i+2]\r
+\r
+ temp1_x = Am_x*dUi1_x\r
+ temp2_x = Ap_x*dUi2_x\r
+ temp_y = Ap_y*dUi_y\r
+ 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
+ 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
+ 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
+\r
+ #Traitement des autres cellules\r
+ else:\r
+\r
+ #Calcul de Am_x\r
+ rho_l=Uk[3*nbCells_x*j+3*i]\r
+ qx_l =Uk[3*nbCells_x*j+3*i+1]\r
+ qy_l =Uk[3*nbCells_x*j+3*i+2]\r
+ rho_r=Uk[3*nbCells_x*j+3*(i+1)]\r
+ qx_r =Uk[3*nbCells_x*j+3*(i+1)+1]\r
+ qy_r =Uk[3*nbCells_x*j+3*(i+1)+2]\r
+ Am_x= jacobianMatricesm_x(dt/dx,rho_l,qx_l,qy_l,rho_r,qx_r,qy_r)\r
+\r
+ #Calcul de Ap_x\r
+ rho_l=Uk[3*nbCells_x*j+3*(i-1)]\r
+ qx_l =Uk[3*nbCells_x*j+3*(i-1)+1]\r
+ qy_l =Uk[3*nbCells_x*j+3*(i-1)+2]\r
+ rho_r=Uk[3*nbCells_x*j+3*i]\r
+ qx_r =Uk[3*nbCells_x*j+3*i+1]\r
+ qy_r =Uk[3*nbCells_x*j+3*i+2]\r
+ Ap_x= jacobianMatricesp_x(dt/dx,rho_l,qx_l,qy_l,rho_r,qx_r,qy_r)\r
+\r
+ #calcul de Am_y\r
+ rho_l=Uk[3*nbCells_x*j+3*i]\r
+ qx_l =Uk[3*nbCells_x*j+3*i+1]\r
+ qy_l =Uk[3*nbCells_x*j+3*i+2]\r
+ rho_r=Uk[3*nbCells_x*(j+1)+3*i]\r
+ qx_r =Uk[3*nbCells_x*(j+1)+3*i+1]\r
+ qy_r =Uk[3*nbCells_x*(j+1)+3*i+2]\r
+ Am_y= jacobianMatricesm_y(dt/dx,rho_l,qx_l,qy_l,rho_r,qx_r,qy_r)\r
+\r
+ #calcul de Ap_y\r
+ rho_l=Uk[3*nbCells_x*(j-1)+3*i]\r
+ qx_l =Uk[3*nbCells_x*(j-1)+3*i+1]\r
+ qy_l =Uk[3*nbCells_x*(j-1)+3*i+2]\r
+ rho_r=Uk[3*nbCells_x*j+3*i]\r
+ qx_r =Uk[3*nbCells_x*j+3*i+1]\r
+ qy_r =Uk[3*nbCells_x*j+3*i+2]\r
+ Ap_y= jacobianMatricesp_y(dt/dx,rho_l,qx_l,qy_l,rho_r,qx_r,qy_r)\r
+\r
+ #remplissage de la divMat\r
+ divMat.addValue(3*nbCells_x*j+3*i,3*nbCells_x*j+3*(i-1),Ap_x*(-1.))\r
+ divMat.addValue(3*nbCells_x*j+3*i,3*nbCells_x*(j-1)+3*i,Ap_y*(-1.))\r
+ divMat.addValue(3*nbCells_x*j+3*i,3*nbCells_x*j+3*(i+1),Am_x)\r
+ divMat.addValue(3*nbCells_x*j+3*i,3*nbCells_x*(j+1)+3*i,Am_y)\r
+ divMat.addValue(3*nbCells_x*j+3*i,3*nbCells_x*j+3*i,Ap_x+Ap_y-Am_x-Am_y)\r
+\r
+ #remplissage du membre de droite\r
+\r
+ dUi1_x[0] = Uk[3*nbCells_x*j+3*(i+1)]-Uk[3*nbCells_x*j+3*i]\r
+ dUi1_x[1] = Uk[3*nbCells_x*j+3*(i+1)+1]-Uk[3*nbCells_x*j+3*i+1]\r
+ dUi1_x[2] = Uk[3*nbCells_x*j+3*(i+1)+2]-Uk[3*nbCells_x*j+3*i+2]\r
+\r
+ dUi2_x[0] = Uk[3*nbCells_x*j+3*i]-Uk[3*nbCells_x*j+3*(i-1)]\r
+ dUi2_x[1] = Uk[3*nbCells_x*j+3*i+1]-Uk[3*nbCells_x*j+3*(i-1)+1]\r
+ dUi2_x[2] = Uk[3*nbCells_x*j+3*i+2]-Uk[3*nbCells_x*j+3*(i-1)+2]\r
+\r
+ dUi1_y[0] = Uk[3*nbCells_x*(j+1)+3*i]-Uk[3*nbCells_x*j+3*i]\r
+ dUi1_y[1] = Uk[3*nbCells_x*(j+1)+3*i+1]-Uk[3*nbCells_x*j+3*i+1]\r
+ dUi1_y[2] = Uk[3*nbCells_x*(j+1)+3*i+2]-Uk[3*nbCells_x*j+3*i+2]\r
+\r
+ dUi2_y[0] = Uk[3*nbCells_x*j+3*i]-Uk[3*nbCells_x*(j-1)+3*i]\r
+ dUi2_y[1] = Uk[3*nbCells_x*j+3*i+1]-Uk[3*nbCells_x*(j-1)+3*i+1]\r
+ dUi2_y[2] = Uk[3*nbCells_x*j+3*i+2]-Uk[3*nbCells_x*(j-1)+3*i+2]\r
+\r
+ temp1_x = Am_x*dUi1_x\r
+ temp2_x = Ap_x*dUi2_x\r
+ temp1_y = Am_y*dUi1_y\r
+ temp2_y = Ap_y*dUi2_y\r
+ 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
+ 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
+ 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
+\r
+ divMat.diagonalShift(1) #only after filling all coefficients\r
+ LS=cdmath.LinearSolver(divMat,Rhs,iterGMRESMax, precision, "GMRES","LU")\r
+ dUk=LS.solve();\r
+ residu = dUk.norm()\r
+ Uk+=dUk\r
+ #print("Newton iteration number ",k, "residu = ",residu)\r
+ if(not LS.getStatus()):\r
+ print("Linear system did not converge ", LS.getNumberOfIter(), " GMRES iterations")\r
+ raise ValueError("Pas de convergence du système linéaire");\r
+ k=k+1\r
+ Un = Uk.deepCopy()\r
+ dUn-=Un\r
+ isStationary = dUn.norm()<precision\r
+ \r
+ for i in range(nbCells):\r
+ rho_field[i] = Un[nbComp*i]\r
+ q_x_field[i] = Un[nbComp*i+1]\r
+ q_y_field[i] = Un[nbComp*i+2]\r
+\r
+ time=time+dt;\r
+ it=it+1;\r
+ #Sauvegardes\r
+ if(it==1 or it%output_freq==0 or it>=ntmax or isStationary or time >=tmax):\r
+ print("-- Iter: " + str(it) + ", Time: " + str(time) + ", dt: " + str(dt) + ", Newton iterations: "+str(k)+", ||dUn|| = "+str(dUn.norm()))\r
+ print("Linear system converged in ", LS.getNumberOfIter(), " GMRES iterations, résidu = ", residu)\r
+ for k in range(nbCells):\r
+ rho =rho_field[k]\r
+ velocity_field[k,0]=q_x_field[i]/rho\r
+ if(dim>1):\r
+ velocity_field[k,1]=q_y_field[k]/rho\r
+ print\r
+ rho_field.setTime(time,it);\r
+ rho_field.writeVTK("EulerIsothermal_"+str(dim)+"DConservativeStaggered_"+meshName+"_density",False);\r
+ q_x_field.setTime(time,it);\r
+ q_x_field.writeVTK("EulerIsothermal_"+str(dim)+"DConservativeStaggered_"+meshName+"_momentumX",False);\r
+ q_y_field.setTime(time,it);\r
+ q_y_field.writeVTK("EulerIsothermal_"+str(dim)+"DConservativeStaggered_"+meshName+"_momentumY",False);\r
+ velocity_field.setTime(time,it);\r
+ velocity_field.writeVTK("EulerIsothermal_"+str(dim)+"DConservativeStaggered_"+meshName+"_velocity",False);\r
+ #Postprocessing : save 2D picture\r
+ 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
+ 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
+ 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
+ print("-- Iter: " + str(it) + ", Time: " + str(time) + ", dt: " + str(dt))\r
+ if(it>=ntmax):\r
+ print("Nombre de pas de temps maximum ntmax= ", ntmax, " atteint")\r
+ return\r
+ elif(isStationary):\r
+ print("Régime stationnaire atteint au pas de temps ", it, ", t= ", time)\r
+ print("------------------------------------------------------------------------------------")\r
+ return\r
+ else:\r
+ print("Temps maximum Tmax= ", tmax, " atteint")\r
+ return\r
+\r
+def solve( my_mesh,filename, resolution):\r
+ print("Resolution of the Isothermal Euler system in dimension 2 on "+str(my_mesh.getNumberOfCells())+ " cells")\r
+ print("Initial data : ", "Riemann problem")\r
+ print("Boundary conditions : ", "Neumann")\r
+ print("Mesh name : ",filename )\r
+ # Problem data\r
+ tmax = 10.\r
+ ntmax = 10\r
+ cfl=1\r
+ output_freq = 1\r
+ EulerSystemStaggered(ntmax, tmax, cfl, output_freq, my_mesh, filename)\r
+ return\r
+\r
+if __name__ == """__main__""":\r
+ if len(sys.argv) >1 :\r
+ filename=sys.argv[1]\r
+ my_mesh = cdmath.Mesh(filename)\r
+ else :\r
+ filename = "CartesianGrid"\r
+ ax=0;bx=1;nx=20;\r
+ ay=0;by=1;ny=10; \r
+ my_mesh = cdmath.Mesh(ax,bx,nx,ay,by,ny)\r
+\r
+ solve(my_mesh,filename,100)\r
