From: michael Date: Sat, 12 Dec 2020 12:32:07 +0000 (+0100) Subject: Added new examples for 1D single phase flows X-Git-Tag: V9_7_0~104 X-Git-Url: http://git.salome-platform.org/gitweb/?a=commitdiff_plain;h=3101e62a85cfd865175ca3fbe3fa2624f79bfb9d;p=tools%2Fsolverlab.git Added new examples for 1D single phase flows --- diff --git a/CDMATH/tests/examples/CMakeLists.txt b/CDMATH/tests/examples/CMakeLists.txt index 7bb78e6..8543fb9 100755 --- a/CDMATH/tests/examples/CMakeLists.txt +++ b/CDMATH/tests/examples/CMakeLists.txt @@ -9,6 +9,7 @@ file(GLOB NICE_EXAMPLES_TO_INSTALL WaveSystem1DUpwind WaveSystem1DUpwind_RiemannProblem WaveSystem1Staggered_RiemannProblem WaveSystem2DUpwind_RiemannProblem WaveSystem_Stationary/WaveSystemUpwind WaveSystem_Stationary/WaveSystemCentered WaveSystem_Stationary/WaveSystemPStag WaveSystem_Stationary/WaveSystemStaggered WaveSystem_Shock/WaveSystemUpwind WaveSystem_Shock/WaveSystemCentered WaveSystem_Shock/WaveSystemPStag WaveSystem_Shock/WaveSystemStaggered +EulerSystem1D/EulerSystem1DUpwind EulerSystem1D/EulerSystem1DUpwindEntrCorr EulerSystem1D/EulerSystem1DConservativeStaggered EulerSystem_Shock/EulerSystemStaggered EulerSystem_Shock/EulerSystemUpwind ) install(DIRECTORY ${NICE_EXAMPLES_TO_INSTALL} DESTINATION share/examples) @@ -54,7 +55,11 @@ IF (CDMATH_WITH_PYTHON AND CDMATH_WITH_PETSC AND CDMATH_WITH_POSTPRO) ADD_SUBDIRECTORY(WaveSystem_Shock/WaveSystemPStag) ADD_SUBDIRECTORY(WaveSystem_Shock/WaveSystemStaggered) ADD_SUBDIRECTORY(WaveSystem_Shock/WaveSystemCentered) + ADD_SUBDIRECTORY(EulerSystem1D/EulerSystem1DUpwind) + ADD_SUBDIRECTORY(EulerSystem1D/EulerSystem1DUpwindEntrCorr) + ADD_SUBDIRECTORY(EulerSystem1D/EulerSystem1DConservativeStaggered) ADD_SUBDIRECTORY(EulerSystem_Shock/EulerSystemUpwind) + ADD_SUBDIRECTORY(EulerSystem_Shock/EulerSystemStaggered) ENDIF (CDMATH_WITH_PYTHON AND CDMATH_WITH_PETSC AND CDMATH_WITH_POSTPRO) diff --git a/CDMATH/tests/examples/EulerSystem1D/EulerSystem1DConservativeStaggered/CMakeLists.txt b/CDMATH/tests/examples/EulerSystem1D/EulerSystem1DConservativeStaggered/CMakeLists.txt new file mode 100755 index 0000000..a4899a3 --- /dev/null +++ b/CDMATH/tests/examples/EulerSystem1D/EulerSystem1DConservativeStaggered/CMakeLists.txt @@ -0,0 +1,8 @@ + +if (CDMATH_WITH_PYTHON AND CDMATH_WITH_PETSC AND CDMATH_WITH_POSTPRO) + + ADD_TEST(ExampleEulerSystem_1DShock_ConservativeStaggered ${PYTHON_EXECUTABLE} ${CMAKE_CURRENT_SOURCE_DIR}/EulerSystem1DConservativeStaggered.py ) + +endif (CDMATH_WITH_PYTHON AND CDMATH_WITH_PETSC AND CDMATH_WITH_POSTPRO) + + diff --git a/CDMATH/tests/examples/EulerSystem1D/EulerSystem1DConservativeStaggered/EulerSystem1DConservativeStaggered.py b/CDMATH/tests/examples/EulerSystem1D/EulerSystem1DConservativeStaggered/EulerSystem1DConservativeStaggered.py new file mode 100644 index 0000000..c1158ac --- /dev/null +++ b/CDMATH/tests/examples/EulerSystem1D/EulerSystem1DConservativeStaggered/EulerSystem1DConservativeStaggered.py @@ -0,0 +1,291 @@ +#!/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 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) diff --git a/CDMATH/tests/examples/EulerSystem1D/EulerSystem1DUpwind/CMakeLists.txt b/CDMATH/tests/examples/EulerSystem1D/EulerSystem1DUpwind/CMakeLists.txt new file mode 100755 index 0000000..8b41b92 --- /dev/null +++ b/CDMATH/tests/examples/EulerSystem1D/EulerSystem1DUpwind/CMakeLists.txt @@ -0,0 +1,8 @@ + +if (CDMATH_WITH_PYTHON AND CDMATH_WITH_PETSC AND CDMATH_WITH_POSTPRO) + + ADD_TEST(ExampleEulerSystem_1DShock_UpwindExplicit ${PYTHON_EXECUTABLE} ${CMAKE_CURRENT_SOURCE_DIR}/EulerSystem1DUpwind.py ) + +endif (CDMATH_WITH_PYTHON AND CDMATH_WITH_PETSC AND CDMATH_WITH_POSTPRO) + + diff --git a/CDMATH/tests/examples/EulerSystem1D/EulerSystem1DUpwind/EulerSystem1DUpwind.py b/CDMATH/tests/examples/EulerSystem1D/EulerSystem1DUpwind/EulerSystem1DUpwind.py new file mode 100644 index 0000000..ee66cd1 --- /dev/null +++ b/CDMATH/tests/examples/EulerSystem1D/EulerSystem1DUpwind/EulerSystem1DUpwind.py @@ -0,0 +1,199 @@ +# 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): + 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") diff --git a/CDMATH/tests/examples/EulerSystem1D/EulerSystem1DUpwindEntrCorr/CMakeLists.txt b/CDMATH/tests/examples/EulerSystem1D/EulerSystem1DUpwindEntrCorr/CMakeLists.txt new file mode 100755 index 0000000..655c40f --- /dev/null +++ b/CDMATH/tests/examples/EulerSystem1D/EulerSystem1DUpwindEntrCorr/CMakeLists.txt @@ -0,0 +1,8 @@ + +if (CDMATH_WITH_PYTHON AND CDMATH_WITH_PETSC AND CDMATH_WITH_POSTPRO) + + ADD_TEST(ExampleEulerSystem_1DShock_UpwindExplicit_EntrCorr ${PYTHON_EXECUTABLE} ${CMAKE_CURRENT_SOURCE_DIR}/EulerSystem1DUpwindEntrCorr.py ) + +endif (CDMATH_WITH_PYTHON AND CDMATH_WITH_PETSC AND CDMATH_WITH_POSTPRO) + + diff --git a/CDMATH/tests/examples/EulerSystem1D/EulerSystem1DUpwindEntrCorr/EulerSystem1DUpwindEntrCorr.py b/CDMATH/tests/examples/EulerSystem1D/EulerSystem1DUpwindEntrCorr/EulerSystem1DUpwindEntrCorr.py new file mode 100644 index 0000000..e1ae7f0 --- /dev/null +++ b/CDMATH/tests/examples/EulerSystem1D/EulerSystem1DUpwindEntrCorr/EulerSystem1DUpwindEntrCorr.py @@ -0,0 +1,199 @@ +# 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): + 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") diff --git a/CDMATH/tests/examples/EulerSystem_Shock/EulerSystemStaggered/CMakeLists.txt b/CDMATH/tests/examples/EulerSystem_Shock/EulerSystemStaggered/CMakeLists.txt new file mode 100755 index 0000000..211c08d --- /dev/null +++ b/CDMATH/tests/examples/EulerSystem_Shock/EulerSystemStaggered/CMakeLists.txt @@ -0,0 +1,8 @@ + +if (CDMATH_WITH_PYTHON AND CDMATH_WITH_PETSC AND CDMATH_WITH_POSTPRO) + + ADD_TEST(ExampleEulerSystem_2DShock_ConservativeStaggered ${PYTHON_EXECUTABLE} ${CMAKE_CURRENT_SOURCE_DIR}/EulerIsothermal_2DConservativeStaggered.py ) + +endif (CDMATH_WITH_PYTHON AND CDMATH_WITH_PETSC AND CDMATH_WITH_POSTPRO) + + diff --git a/CDMATH/tests/examples/EulerSystem_Shock/EulerSystemStaggered/EulerIsothermal_2DConservativeStaggered.py b/CDMATH/tests/examples/EulerSystem_Shock/EulerSystemStaggered/EulerIsothermal_2DConservativeStaggered.py new file mode 100644 index 0000000..5480b7e --- /dev/null +++ b/CDMATH/tests/examples/EulerSystem_Shock/EulerSystemStaggered/EulerIsothermal_2DConservativeStaggered.py @@ -0,0 +1,798 @@ +#!/usr/bin/env python3 +# -*-coding:utf-8 -* + +#=============================================================================================================================== +# Name : Résolution VF des équations d'Euler isotherme 2D sans terme source +# \partial_t rho + \div q = 0 +# \partial_t q + \div q \otimes q/rho \grad p = 0 +# Author : Michaël Ndjinga, Coraline Mounier +# Copyright : CEA Saclay 2020 +# Description : Propagation d'une onde de choc droite +# Utilisation du schéma upwind explicite ou implicite sur un maillage général +# Initialisation par une discontinuité verticale +# Conditions aux limites de Neumann +# Création et sauvegarde du champ résultant et des figures +#================================================================================================================================ + + +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 +import sys +import PV_routines + +p0=155.e5 #reference pressure in a pressurised nuclear vessel +c0=700. #reference sound speed for water at 155 bars +rho0=p0/(c0*c0) #reference density +precision=1e-5 + +def initial_conditions_Riemann_problem(my_mesh): + print( "Initial data : Riemann problem" ) + dim = my_mesh.getMeshDimension() + nbCells = my_mesh.getNumberOfCells() + + xcentre = 0.5 + + density_field = cdmath.Field("Density", cdmath.CELLS, my_mesh, 1) + q_x_field = cdmath.Field("Momentum x", cdmath.CELLS, my_mesh, 1) + q_y_field = cdmath.Field("Momentum y", cdmath.CELLS, my_mesh, 1) + #Velocity field with 3 components should be created for streamlines to be drawn + velocity_field = cdmath.Field("Velocity", cdmath.CELLS, my_mesh, 3) + + for i in range(nbCells): + x = my_mesh.getCell(i).x() + + # Initial momentum is zero + q_x_field[i] = 0 + q_y_field[i] = 0 + + # Initial velocity is zero + velocity_field[i,0] = 0 + velocity_field[i,1] = 0 + velocity_field[i,2] = 0 + + if x < xcentre: + density_field[i] = rho0 + pass + else: + density_field[i] = rho0/2 + pass + pass + + return density_field, q_x_field, q_y_field, velocity_field + + +def jacobianMatricesm_x(coeff,rho_l,q_lx,q_ly,rho_r,q_rx,q_ry): + RoeMatx = cdmath.Matrix(3,3); + Dmacx = cdmath.Matrix(3,3); + + u_lx=q_lx/rho_l + u_rx=q_rx/rho_r + u_ly=q_ly/rho_l + u_ry=q_ry/rho_r + if rho_l<0 or rho_r<0: + print("rho_l=",rho_l, " rho_r= ",rho_r) + raise ValueError("Negative density") + ux = (u_lx*sqrt(rho_l)+u_rx*sqrt(rho_r))/(sqrt(rho_l)+sqrt(rho_r)); + uy = (u_ly*sqrt(rho_l)+u_ry*sqrt(rho_r))/(sqrt(rho_l)+sqrt(rho_r)); + + RoeMatx[0,0] = 0 + RoeMatx[0,1] = 1 + RoeMatx[0,2] = 0 + RoeMatx[1,0] = c0*c0 - ux*ux + RoeMatx[1,1] = 2*ux + RoeMatx[1,2] = 0 + RoeMatx[2,0] = -ux*uy + RoeMatx[2,1] = uy + RoeMatx[2,2] = ux + + Dmacx[0,0] = abs(ux)-ux + Dmacx[0,1] = 1 + Dmacx[0,2] = 0 + Dmacx[1,0] = -c0*c0-ux*ux + Dmacx[1,1] = abs(ux)+ux + Dmacx[1,2] = 0 + Dmacx[2,0] = -ux*uy + Dmacx[2,1] = uy + Dmacx[2,2] = abs(ux) + + return (RoeMatx-Dmacx)*coeff*0.5 + +def jacobianMatricesp_x(coeff,rho_l,q_lx,q_ly,rho_r,q_rx,q_ry): + RoeMatx = cdmath.Matrix(3,3) + Dmacx = cdmath.Matrix(3,3) + + u_lx=q_lx/rho_l + u_rx=q_rx/rho_r + u_ly=q_ly/rho_l + u_ry=q_ry/rho_r + if rho_l<0 or rho_r<0: + print("rho_l=",rho_l, " rho_r= ",rho_r) + raise ValueError("Negative density") + ux = (u_lx*sqrt(rho_l)+u_rx*sqrt(rho_r))/(sqrt(rho_l)+sqrt(rho_r)); + uy = (u_ly*sqrt(rho_l)+u_ry*sqrt(rho_r))/(sqrt(rho_l)+sqrt(rho_r)); + + RoeMatx[0,0] = 0 + RoeMatx[0,1] = 1 + RoeMatx[0,2] = 0 + RoeMatx[1,0] = c0*c0 - ux*ux + RoeMatx[1,1] = 2*ux + RoeMatx[1,2] = 0 + RoeMatx[2,0] = -ux*uy + RoeMatx[2,1] = uy + RoeMatx[2,2] = ux + + Dmacx[0,0] = abs(ux)-ux + Dmacx[0,1] = 1 + Dmacx[0,2] = 0 + Dmacx[1,0] = -c0*c0-ux*ux + Dmacx[1,1] = abs(ux)+ux + Dmacx[1,2] = 0 + Dmacx[2,0] = -ux*uy + Dmacx[2,1] = uy + Dmacx[2,2] = abs(ux) + + return (RoeMatx+Dmacx)*coeff*0.5 + +def jacobianMatricesm_y(coeff,rho_l,q_lx,q_ly,rho_r,q_rx,q_ry): + RoeMaty = cdmath.Matrix(3,3); + Dmacy = cdmath.Matrix(3,3); + + u_lx=q_lx/rho_l + u_rx=q_rx/rho_r + u_ly=q_ly/rho_l + u_ry=q_ry/rho_r + if rho_l<0 or rho_r<0: + print("rho_l=",rho_l, " rho_r= ",rho_r) + raise ValueError("Negative density") + ux = (u_lx*sqrt(rho_l)+u_rx*sqrt(rho_r))/(sqrt(rho_l)+sqrt(rho_r)); + uy = (u_ly*sqrt(rho_l)+u_ry*sqrt(rho_r))/(sqrt(rho_l)+sqrt(rho_r)); + + RoeMaty[0,0] = 0 + RoeMaty[0,1] = 0 + RoeMaty[0,2] = 1 + RoeMaty[1,0] = -ux*uy + RoeMaty[1,1] = uy + RoeMaty[1,2] = ux + RoeMaty[2,0] = c0*c0-uy*uy + RoeMaty[2,1] = 0 + RoeMaty[2,2] = 2*uy + + Dmacy[0,0] = abs(uy)-uy + Dmacy[0,1] = 0 + Dmacy[0,2] = 1 + Dmacy[1,0] = -ux*uy + Dmacy[1,1] = abs(uy) + Dmacy[1,2] = ux + Dmacy[2,0] = -c0*c0-uy*uy + Dmacy[2,1] = 0 + Dmacy[2,2] = abs(uy)+uy + + return (RoeMaty-Dmacy)*coeff*0.5 + +def jacobianMatricesp_y(coeff,rho_l,q_lx,q_ly,rho_r,q_rx,q_ry): + RoeMaty = cdmath.Matrix(3,3); + Dmacy = cdmath.Matrix(3,3); + + u_lx=q_lx/rho_l + u_rx=q_rx/rho_r + u_ly=q_ly/rho_l + u_ry=q_ry/rho_r + if rho_l<0 or rho_r<0: + print("rho_l=",rho_l, " rho_r= ",rho_r) + raise ValueError("Negative density") + ux = (u_lx*sqrt(rho_l)+u_rx*sqrt(rho_r))/(sqrt(rho_l)+sqrt(rho_r)); + uy = (u_ly*sqrt(rho_l)+u_ry*sqrt(rho_r))/(sqrt(rho_l)+sqrt(rho_r)); + + RoeMaty[0,0] = 0 + RoeMaty[0,1] = 0 + RoeMaty[0,2] = 1 + RoeMaty[1,0] = -ux*uy + RoeMaty[1,1] = uy + RoeMaty[1,2] = ux + RoeMaty[2,0] = c0*c0-uy*uy + RoeMaty[2,1] = 0 + RoeMaty[2,2] = 2*uy + + Dmacy[0,0] = abs(uy)-uy + Dmacy[0,1] = 0 + Dmacy[0,2] = 1 + Dmacy[1,0] = -ux*uy + Dmacy[1,1] = abs(uy) + Dmacy[1,2] = ux + Dmacy[2,0] = -c0*c0-uy*uy + Dmacy[2,1] = 0 + Dmacy[2,2] = abs(uy)+uy + + return (RoeMaty+Dmacy)*coeff*0.5 + +def EulerSystemStaggered(ntmax, tmax, cfl,output_freq, my_mesh, meshName): + + if not my_mesh.isStructured() : + raise ValueError("Euler_ConservativeStaggered2D_RiemannProblem.py expects a structured mesh") + + dim=my_mesh.getMeshDimension() + if dim != 2 : + raise ValueError("Euler_ConservativeStaggered2D_RiemannProblem.py expects a 2D mesh") + + nbComp=dim+1 + # Mesh parameters + nbCells = my_mesh.getNumberOfCells() + nbCells_x = my_mesh.getNx() + nbCells_y = my_mesh.getNy() + dx,dy = my_mesh.getDXYZ() + nbVoisinsMax=my_mesh.getMaxNbNeighbours(cdmath.CELLS) + dx_min=my_mesh.minRatioVolSurf() + + # Time evolution parameters + time = 0. + it=0; + isStationary=False + iterGMRESMax = 50 + newton_max = 50 + + #iteration vectors + Un =cdmath.Vector(nbCells*nbComp) + dUn=cdmath.Vector(nbCells*nbComp) + dUk=cdmath.Vector(nbCells*nbComp) + Rhs=cdmath.Vector(nbCells*nbComp) + + dUi_x=cdmath.Vector(nbComp) + dUi_y=cdmath.Vector(nbComp) + dUi1_x=cdmath.Vector(nbComp) + dUi2_x=cdmath.Vector(nbComp) + dUi1_y=cdmath.Vector(nbComp) + dUi2_y=cdmath.Vector(nbComp) + + # Initial conditions # + print("Construction of the initial condition …") + rho_field, q_x_field, q_y_field, velocity_field= initial_conditions_Riemann_problem(my_mesh) + + for i in range(nbCells): + Un[nbComp*i] = rho_field[i] + Un[nbComp*i+1] = q_x_field[i] + Un[nbComp*i+2] = q_y_field[i] + + #Sauvegarde de la donnée initiale + rho_field.setTime(time,it); + rho_field.writeVTK("EulerIsothermal_"+str(dim)+"DConservativeStaggered_"+meshName+"_density"); + q_x_field.setTime(time,it); + q_x_field.writeVTK("EulerIsothermal_"+str(dim)+"DConservativeStaggered_"+meshName+"_momentumX"); + q_y_field.setTime(time,it); + q_y_field.writeVTK("EulerIsothermal_"+str(dim)+"DConservativeStaggered_"+meshName+"_momentumY"); + velocity_field.setTime(time,it); + velocity_field.writeVTK("EulerIsothermal_"+str(dim)+"DConservativeStaggered_"+meshName+"_velocity"); + + print("Starting computation of the isothermal Euler system with a conservative staggered scheme …") + divMat=cdmath.SparseMatrixPetsc(nbCells*nbComp,nbCells*nbComp,(nbVoisinsMax+1)*nbComp) + + # Starting time loop + while (it 1/precision ): + + 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 + #DEBUT BOUCLE NEWTON + for j in range(nbCells_y): + for i in range(nbCells_x): + #Traitement des coins + if ( j==0 and i==0) : + #Calcul de Am_x + rho_l=Uk[3*nbCells_x*j+3*i] + qx_l =Uk[3*nbCells_x*j+3*i+1] + qy_l =Uk[3*nbCells_x*j+3*i+2] + rho_r=Uk[3*nbCells_x*j+3*(i+1)] + qx_r =Uk[3*nbCells_x*j+3*(i+1)+1] + qy_r =Uk[3*nbCells_x*j+3*(i+1)+2] + Am_x= jacobianMatricesm_x(dt/dx,rho_l,qx_l,qy_l,rho_r,qx_r,qy_r) + + #calcul de Am_y + rho_l=Uk[3*nbCells_x*j+3*i] + qx_l =Uk[3*nbCells_x*j+3*i+1] + qy_l =Uk[3*nbCells_x*j+3*i+2] + rho_r=Uk[3*nbCells_x*(j+1)+3*i] + qx_r =Uk[3*nbCells_x*(j+1)+3*i+1] + qy_r =Uk[3*nbCells_x*(j+1)+3*i+2] + Am_y= jacobianMatricesm_y(dt/dx,rho_l,qx_l,qy_l,rho_r,qx_r,qy_r) + + divMat.addValue(3*nbCells_x*j+3*i,3*nbCells_x*j+3*(i+1),Am_x) + divMat.addValue(3*nbCells_x*j+3*i,3*nbCells_x*(j+1)+3*i,Am_y) + divMat.addValue(3*nbCells_x*j+3*i,3*nbCells_x*j+3*i,Am_x*(-1.)-Am_y) + + dUi_x[0] = Uk[3*nbCells_x*j+3*(i+1)]-Uk[3*nbCells_x*j+3*i] + dUi_x[1] = Uk[3*nbCells_x*j+3*(i+1)+1]-Uk[3*nbCells_x*j+3*i+1] + dUi_x[2] = Uk[3*nbCells_x*j+3*(i+1)+2]-Uk[3*nbCells_x*j+3*i+2] + dUi_y[0] = Uk[3*nbCells_x*(j+1)+3*i]-Uk[3*nbCells_x*j+3*i] + dUi_y[1] = Uk[3*nbCells_x*(j+1)+3*i+1]-Uk[3*nbCells_x*j+3*i+1] + dUi_y[2] = Uk[3*nbCells_x*(j+1)+3*i+2]-Uk[3*nbCells_x*j+3*i+2] + + temp_x = Am_x*dUi_x + temp_y = Am_y*dUi_y + #print("Bloc 0 matrix : ", Am) + 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]) + 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]) + 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]) + + elif(i==0 and j==nbCells_y-1): + #Calcul de Am_x + rho_l=Uk[3*nbCells_x*j+3*i] + qx_l =Uk[3*nbCells_x*j+3*i+1] + qy_l =Uk[3*nbCells_x*j+3*i+2] + rho_r=Uk[3*nbCells_x*j+3*(i+1)] + qx_r =Uk[3*nbCells_x*j+3*(i+1)+1] + qy_r =Uk[3*nbCells_x*j+3*(i+1)+2] + Am_x= jacobianMatricesm_x(dt/dx,rho_l,qx_l,qy_l,rho_r,qx_r,qy_r) + + #calcul de Ap_y + rho_l=Uk[3*nbCells_x*(j-1)+3*i] + qx_l =Uk[3*nbCells_x*(j-1)+3*i+1] + qy_l =Uk[3*nbCells_x*(j-1)+3*i+2] + rho_r=Uk[3*nbCells_x*j+3*i] + qx_r =Uk[3*nbCells_x*j+3*i+1] + qy_r =Uk[3*nbCells_x*j+3*i+2] + Ap_y= jacobianMatricesp_y(dt/dx,rho_l,qx_l,qy_l,rho_r,qx_r,qy_r) + + divMat.addValue(3*nbCells_x*j+3*i,3*nbCells_x*j+3*(i+1),Am_x) + divMat.addValue(3*nbCells_x*j+3*i,3*nbCells_x*(j-1)+3*i,Ap_y*(-1.)) + divMat.addValue(3*nbCells_x*j+3*i,3*nbCells_x*j+3*i,Am_x*(-1.)+Ap_y) + + dUi_x[0] = Uk[3*nbCells_x*j+3*(i+1)]-Uk[3*nbCells_x*j+3*i] + dUi_x[1] = Uk[3*nbCells_x*j+3*(i+1)+1]-Uk[3*nbCells_x*j+3*i+1] + dUi_x[2] = Uk[3*nbCells_x*j+3*(i+1)+2]-Uk[3*nbCells_x*j+3*i+2] + dUi_y[0] = Uk[3*nbCells_x*j+3*i]-Uk[3*nbCells_x*(j-1)+3*i] + dUi_y[1] = Uk[3*nbCells_x*j+3*i+1]-Uk[3*nbCells_x*(j-1)+3*i+1] + dUi_y[2] = Uk[3*nbCells_x*j+3*i+2]-Uk[3*nbCells_x*(j-1)+3*i+2] + + temp_x = Am_x*dUi_x + temp_y = Ap_y*dUi_y + #print("Bloc 0 matrix : ", Am) + 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]) + 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]) + 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]) + + elif(i==nbCells_x-1 and j==0): + + #Calcul de Ap_x + rho_l=Uk[3*nbCells_x*j+3*(i-1)] + qx_l =Uk[3*nbCells_x*j+3*(i-1)+1] + qy_l =Uk[3*nbCells_x*j+3*(i-1)+2] + rho_r=Uk[3*nbCells_x*j+3*i] + qx_r =Uk[3*nbCells_x*j+3*i+1] + qy_r =Uk[3*nbCells_x*j+3*i+2] + Ap_x= jacobianMatricesp_x(dt/dx,rho_l,qx_l,qy_l,rho_r,qx_r,qy_r) + + #calcul de Am_y + rho_l=Uk[3*nbCells_x*j+3*i] + qx_l =Uk[3*nbCells_x*j+3*i+1] + qy_l =Uk[3*nbCells_x*j+3*i+2] + rho_r=Uk[3*nbCells_x*(j+1)+3*i] + qx_r =Uk[3*nbCells_x*(j+1)+3*i+1] + qy_r =Uk[3*nbCells_x*(j+1)+3*i+2] + Am_y= jacobianMatricesm_y(dt/dx,rho_l,qx_l,qy_l,rho_r,qx_r,qy_r) + + divMat.addValue(3*nbCells_x*j+3*i,3*nbCells_x*j+3*(i-1),Ap_x*(-1)) + divMat.addValue(3*nbCells_x*j+3*i,3*nbCells_x*(j+1)+3*i,Am_y) + divMat.addValue(3*nbCells_x*j+3*i,3*nbCells_x*j+3*i,Ap_x-Am_y) + + dUi_x[0] = Uk[3*nbCells_x*j+3*i]-Uk[3*nbCells_x*j+3*(i-1)] + dUi_x[1] = Uk[3*nbCells_x*j+3*i+1]-Uk[3*nbCells_x*j+3*(i-1)+1] + dUi_x[2] = Uk[3*nbCells_x*j+3*i+2]-Uk[3*nbCells_x*j+3*(i-1)+2] + dUi_y[0] = Uk[3*nbCells_x*(j+1)+3*i]-Uk[3*nbCells_x*j+3*i] + dUi_y[1] = Uk[3*nbCells_x*(j+1)+3*i+1]-Uk[3*nbCells_x*j+3*i+1] + dUi_y[2] = Uk[3*nbCells_x*(j+1)+3*i+2]-Uk[3*nbCells_x*j+3*i+2] + + temp_x = Ap_x*dUi_x + temp_y= Am_y*dUi_y + #print("Bloc 0 matrix : ", Am) + 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]) + 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]) + 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]) + + elif ( j==nbCells_y-1 and i==nbCells_x-1) : + + #Calcul de Ap_x + rho_l=Uk[3*nbCells_x*j+3*(i-1)] + qx_l =Uk[3*nbCells_x*j+3*(i-1)+1] + qy_l =Uk[3*nbCells_x*j+3*(i-1)+2] + rho_r=Uk[3*nbCells_x*j+3*i] + qx_r =Uk[3*nbCells_x*j+3*i+1] + qy_r =Uk[3*nbCells_x*j+3*i+2] + Ap_x= jacobianMatricesp_x(dt/dx,rho_l,qx_l,qy_l,rho_r,qx_r,qy_r) + + #calcul de Ap_y + rho_l=Uk[3*nbCells_x*(j-1)+3*i] + qx_l =Uk[3*nbCells_x*(j-1)+3*i+1] + qy_l =Uk[3*nbCells_x*(j-1)+3*i+2] + rho_r=Uk[3*nbCells_x*j+3*i] + qx_r =Uk[3*nbCells_x*j+3*i+1] + qy_r =Uk[3*nbCells_x*j+3*i+2] + Ap_y= jacobianMatricesp_y(dt/dx,rho_l,qx_l,qy_l,rho_r,qx_r,qy_r) + + divMat.addValue(3*nbCells_x*j+3*i,3*nbCells_x*j+3*(i-1),Ap_x*(-1.)) + divMat.addValue(3*nbCells_x*j+3*i,3*nbCells_x*(j-1)+3*i,Ap_y*(-1.)) + divMat.addValue(3*nbCells_x*j+3*i,3*nbCells_x*j+3*i,Ap_x+Ap_y) + + dUi_x[0] = Uk[3*nbCells_x*j+3*i]-Uk[3*nbCells_x*j+3*(i-1)] + dUi_x[1] = Uk[3*nbCells_x*j+3*i+1]-Uk[3*nbCells_x*j+3*(i-1)+1] + dUi_x[2] = Uk[3*nbCells_x*j+3*i+2]-Uk[3*nbCells_x*j+3*(i-1)+2] + dUi_y[0] = Uk[3*nbCells_x*j+3*i]-Uk[3*nbCells_x*(j-1)+3*i] + dUi_y[1] = Uk[3*nbCells_x*j+3*i+1]-Uk[3*nbCells_x*(j-1)+3*i+1] + dUi_y[2] = Uk[3*nbCells_x*j+3*i+2]-Uk[3*nbCells_x*(j-1)+3*i+2] + + temp_x = Ap_x*dUi_x + temp_y = Ap_y*dUi_y + 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]) + 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]) + 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]) + + #Traitement des bords restants + elif i==0 : + + #Calcul de Am_x + rho_l=Uk[3*nbCells_x*j+3*i] + qx_l =Uk[3*nbCells_x*j+3*i+1] + qy_l =Uk[3*nbCells_x*j+3*i+2] + rho_r=Uk[3*nbCells_x*j+3*(i+1)] + qx_r =Uk[3*nbCells_x*j+3*(i+1)+1] + qy_r =Uk[3*nbCells_x*j+3*(i+1)+2] + Am_x= jacobianMatricesm_x(dt/dx,rho_l,qx_l,qy_l,rho_r,qx_r,qy_r) + + #calcul de Am_y + rho_l=Uk[3*nbCells_x*j+3*i] + qx_l =Uk[3*nbCells_x*j+3*i+1] + qy_l =Uk[3*nbCells_x*j+3*i+2] + rho_r=Uk[3*nbCells_x*(j+1)+3*i] + qx_r =Uk[3*nbCells_x*(j+1)+3*i+1] + qy_r =Uk[3*nbCells_x*(j+1)+3*i+2] + Am_y= jacobianMatricesm_y(dt/dx,rho_l,qx_l,qy_l,rho_r,qx_r,qy_r) + + #calcul de Ap_y + rho_l=Uk[3*nbCells_x*(j-1)+3*i] + qx_l =Uk[3*nbCells_x*(j-1)+3*i+1] + qy_l =Uk[3*nbCells_x*(j-1)+3*i+2] + rho_r=Uk[3*nbCells_x*j+3*i] + qx_r =Uk[3*nbCells_x*j+3*i+1] + qy_r =Uk[3*nbCells_x*j+3*i+2] + Ap_y= jacobianMatricesp_y(dt/dx,rho_l,qx_l,qy_l,rho_r,qx_r,qy_r) + + #remplissage de la divMat + divMat.addValue(3*nbCells_x*j+3*i,3*nbCells_x*j+3*(i+1),Am_x) + divMat.addValue(3*nbCells_x*j+3*i,3*nbCells_x*(j+1)+3*i,Am_y) + divMat.addValue(3*nbCells_x*j+3*i,3*nbCells_x*(j-1)+3*i,Ap_y*(-1.)) + divMat.addValue(3*nbCells_x*j+3*i,3*nbCells_x*j+3*i,Ap_y-Am_x-Am_y) + + #remplissage du membre de droite + dUi_x[0] = Uk[3*nbCells_x*j+3*(i+1)]-Uk[3*nbCells_x*j+3*i] + dUi_x[1] = Uk[3*nbCells_x*j+3*(i+1)+1]-Uk[3*nbCells_x*j+3*i+1] + dUi_x[2] = Uk[3*nbCells_x*j+3*(i+1)+2]-Uk[3*nbCells_x*j+3*i+2] + dUi1_y[0] = Uk[3*nbCells_x*(j+1)+3*i]-Uk[3*nbCells_x*j+3*i] + dUi1_y[1] = Uk[3*nbCells_x*(j+1)+3*i+1]-Uk[3*nbCells_x*j+3*i+1] + dUi1_y[2] = Uk[3*nbCells_x*(j+1)+3*i+2]-Uk[3*nbCells_x*j+3*i+2] + dUi2_y[0] = Uk[3*nbCells_x*j+3*i]-Uk[3*nbCells_x*(j-1)+3*i] + dUi2_y[1] = Uk[3*nbCells_x*j+3*i+1]-Uk[3*nbCells_x*(j-1)+3*i+1] + dUi2_y[2] = Uk[3*nbCells_x*j+3*i+2]-Uk[3*nbCells_x*(j-1)+3*i+2] + + temp_x = Am_x*dUi_x + temp1_y = Am_y*dUi1_y + temp2_y = Ap_y*dUi2_y + 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]) + 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]) + 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]) + + elif i==nbCells_x-1: + + #Calcul de Ap_x + rho_l=Uk[3*nbCells_x*j+3*(i-1)] + qx_l =Uk[3*nbCells_x*j+3*(i-1)+1] + qy_l =Uk[3*nbCells_x*j+3*(i-1)+2] + rho_r=Uk[3*nbCells_x*j+3*i] + qx_r =Uk[3*nbCells_x*j+3*i+1] + qy_r =Uk[3*nbCells_x*j+3*i+2] + Ap_x= jacobianMatricesp_x(dt/dx,rho_l,qx_l,qy_l,rho_r,qx_r,qy_r) + + #calcul de Am_y + rho_l=Uk[3*nbCells_x*j+3*i] + qx_l =Uk[3*nbCells_x*j+3*i+1] + qy_l =Uk[3*nbCells_x*j+3*i+2] + rho_r=Uk[3*nbCells_x*(j+1)+3*i] + qx_r =Uk[3*nbCells_x*(j+1)+3*i+1] + qy_r =Uk[3*nbCells_x*(j+1)+3*i+2] + Am_y= jacobianMatricesm_y(dt/dx,rho_l,qx_l,qy_l,rho_r,qx_r,qy_r) + + #calcul de Ap_y + rho_l=Uk[3*nbCells_x*(j-1)+3*i] + qx_l =Uk[3*nbCells_x*(j-1)+3*i+1] + qy_l =Uk[3*nbCells_x*(j-1)+3*i+2] + rho_r=Uk[3*nbCells_x*j+3*i] + qx_r =Uk[3*nbCells_x*j+3*i+1] + qy_r =Uk[3*nbCells_x*j+3*i+2] + Ap_y= jacobianMatricesp_y(dt/dx,rho_l,qx_l,qy_l,rho_r,qx_r,qy_r) + + #remplissage de la divMat + divMat.addValue(3*nbCells_x*j+3*i,3*nbCells_x*j+3*(i-1),Ap_x*(-1.)) + divMat.addValue(3*nbCells_x*j+3*i,3*nbCells_x*(j+1)+3*i,Am_y) + divMat.addValue(3*nbCells_x*j+3*i,3*nbCells_x*(j-1)+3*i,Ap_y*(-1.)) + divMat.addValue(3*nbCells_x*j+3*i,3*nbCells_x*j+3*i,Ap_y+Ap_x-Am_y) + + #remplissage du membre de droite + dUi_x[0] = Uk[3*nbCells_x*j+3*i]-Uk[3*nbCells_x*j+3*(i-1)] + dUi_x[1] = Uk[3*nbCells_x*j+3*i+1]-Uk[3*nbCells_x*j+3*(i-1)+1] + dUi_x[2] = Uk[3*nbCells_x*j+3*i+2]-Uk[3*nbCells_x*j+3*(i-1)+2] + dUi1_y[0] = Uk[3*nbCells_x*(j+1)+3*i]-Uk[3*nbCells_x*j+3*i] + dUi1_y[1] = Uk[3*nbCells_x*(j+1)+3*i+1]-Uk[3*nbCells_x*j+3*i+1] + dUi1_y[2] = Uk[3*nbCells_x*(j+1)+3*i+2]-Uk[3*nbCells_x*j+3*i+2] + dUi2_y[0] = Uk[3*nbCells_x*j+3*i]-Uk[3*nbCells_x*(j-1)+3*i] + dUi2_y[1] = Uk[3*nbCells_x*j+3*i+1]-Uk[3*nbCells_x*(j-1)+3*i+1] + dUi2_y[2] = Uk[3*nbCells_x*j+3*i+2]-Uk[3*nbCells_x*(j-1)+3*i+2] + + temp_x = Ap_x*dUi_x + temp1_y = Am_y*dUi1_y + temp2_y = Ap_y*dUi2_y + 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]) + 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]) + 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]) + + elif j==0: + + #Calcul de Am_x + rho_l=Uk[3*nbCells_x*j+3*i] + qx_l =Uk[3*nbCells_x*j+3*i+1] + qy_l =Uk[3*nbCells_x*j+3*i+2] + rho_r=Uk[3*nbCells_x*j+3*(i+1)] + qx_r =Uk[3*nbCells_x*j+3*(i+1)+1] + qy_r =Uk[3*nbCells_x*j+3*(i+1)+2] + Am_x= jacobianMatricesm_x(dt/dx,rho_l,qx_l,qy_l,rho_r,qx_r,qy_r) + + #Calcul de Ap_x + rho_l=Uk[3*nbCells_x*j+3*(i-1)] + qx_l =Uk[3*nbCells_x*j+3*(i-1)+1] + qy_l =Uk[3*nbCells_x*j+3*(i-1)+2] + rho_r=Uk[3*nbCells_x*j+3*i] + qx_r =Uk[3*nbCells_x*j+3*i+1] + qy_r =Uk[3*nbCells_x*j+3*i+2] + Ap_x= jacobianMatricesp_x(dt/dx,rho_l,qx_l,qy_l,rho_r,qx_r,qy_r) + + #calcul de Am_y + rho_l=Uk[3*nbCells_x*j+3*i] + qx_l =Uk[3*nbCells_x*j+3*i+1] + qy_l =Uk[3*nbCells_x*j+3*i+2] + rho_r=Uk[3*nbCells_x*(j+1)+3*i] + qx_r =Uk[3*nbCells_x*(j+1)+3*i+1] + qy_r =Uk[3*nbCells_x*(j+1)+3*i+2] + Am_y= jacobianMatricesm_y(dt/dx,rho_l,qx_l,qy_l,rho_r,qx_r,qy_r) + + #remplissage de la divMat + divMat.addValue(3*nbCells_x*j+3*i,3*nbCells_x*j+3*(i-1),Ap_x*(-1.)) + divMat.addValue(3*nbCells_x*j+3*i,3*nbCells_x*j+3*(i+1),Am_x) + divMat.addValue(3*nbCells_x*j+3*i,3*nbCells_x*(j+1)+3*i,Am_y) + divMat.addValue(3*nbCells_x*j+3*i,3*nbCells_x*j+3*i,Ap_x-Am_x-Am_y) + + #remplissage du membre de droite + dUi1_x[0] = Uk[3*nbCells_x*j+3*(i+1)]-Uk[3*nbCells_x*j+3*i] + dUi1_x[1] = Uk[3*nbCells_x*j+3*(i+1)+1]-Uk[3*nbCells_x*j+3*i+1] + dUi1_x[2] = Uk[3*nbCells_x*j+3*(i+1)+2]-Uk[3*nbCells_x*j+3*i+2] + dUi2_x[0] = Uk[3*nbCells_x*j+3*i]-Uk[3*nbCells_x*j+3*(i-1)] + dUi2_x[1] = Uk[3*nbCells_x*j+3*i+1]-Uk[3*nbCells_x*j+3*(i-1)+1] + dUi2_x[2] = Uk[3*nbCells_x*j+3*i+2]-Uk[3*nbCells_x*j+3*(i-1)+2] + dUi_y[0] = Uk[3*nbCells_x*(j+1)+3*i]-Uk[3*nbCells_x*j+3*i] + dUi_y[1] = Uk[3*nbCells_x*(j+1)+3*i+1]-Uk[3*nbCells_x*j+3*i+1] + dUi_y[2] = Uk[3*nbCells_x*(j+1)+3*i+2]-Uk[3*nbCells_x*j+3*i+2] + + temp1_x = Am_x*dUi1_x + temp2_x = Ap_x*dUi2_x + temp_y = Am_y*dUi_y + 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]) + 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]) + 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]) + + elif j==nbCells_y-1: + + #Calcul de Am_x + rho_l=Uk[3*nbCells_x*j+3*i] + qx_l =Uk[3*nbCells_x*j+3*i+1] + qy_l =Uk[3*nbCells_x*j+3*i+2] + rho_r=Uk[3*nbCells_x*j+3*(i+1)] + qx_r =Uk[3*nbCells_x*j+3*(i+1)+1] + qy_r =Uk[3*nbCells_x*j+3*(i+1)+2] + Am_x= jacobianMatricesm_x(dt/dx,rho_l,qx_l,qy_l,rho_r,qx_r,qy_r) + + #Calcul de Ap_x + rho_l=Uk[3*nbCells_x*j+3*(i-1)] + qx_l =Uk[3*nbCells_x*j+3*(i-1)+1] + qy_l =Uk[3*nbCells_x*j+3*(i-1)+2] + rho_r=Uk[3*nbCells_x*j+3*i] + qx_r =Uk[3*nbCells_x*j+3*i+1] + qy_r =Uk[3*nbCells_x*j+3*i+2] + Ap_x= jacobianMatricesp_x(dt/dx,rho_l,qx_l,qy_l,rho_r,qx_r,qy_r) + + #calcul de Ap_y + rho_l=Uk[3*nbCells_x*(j-1)+3*i] + qx_l =Uk[3*nbCells_x*(j-1)+3*i+1] + qy_l =Uk[3*nbCells_x*(j-1)+3*i+2] + rho_r=Uk[3*nbCells_x*j+3*i] + qx_r =Uk[3*nbCells_x*j+3*i+1] + qy_r =Uk[3*nbCells_x*j+3*i+2] + Ap_y= jacobianMatricesp_y(dt/dx,rho_l,qx_l,qy_l,rho_r,qx_r,qy_r) + + #remplissage de la divMat + divMat.addValue(3*nbCells_x*j+3*i,3*nbCells_x*j+3*(i-1),Ap_x*(-1.)) + divMat.addValue(3*nbCells_x*j+3*i,3*nbCells_x*j+3*(i+1),Am_x) + divMat.addValue(3*nbCells_x*j+3*i,3*nbCells_x*(j-1)+3*i,Ap_y*(-1.)) + divMat.addValue(3*nbCells_x*j+3*i,3*nbCells_x*j+3*i,Ap_x+Ap_y-Am_x) + + #remplissage du membre de droite + dUi1_x[0] = Uk[3*nbCells_x*j+3*(i+1)]-Uk[3*nbCells_x*j+3*i] + dUi1_x[1] = Uk[3*nbCells_x*j+3*(i+1)+1]-Uk[3*nbCells_x*j+3*i+1] + dUi1_x[2] = Uk[3*nbCells_x*j+3*(i+1)+2]-Uk[3*nbCells_x*j+3*i+2] + dUi2_x[0] = Uk[3*nbCells_x*j+3*i]-Uk[3*nbCells_x*j+3*(i-1)] + dUi2_x[1] = Uk[3*nbCells_x*j+3*i+1]-Uk[3*nbCells_x*j+3*(i-1)+1] + dUi2_x[2] = Uk[3*nbCells_x*j+3*i+2]-Uk[3*nbCells_x*j+3*(i-1)+2] + dUi_y[0] = Uk[3*nbCells_x*j+3*i]-Uk[3*nbCells_x*(j-1)+3*i] + dUi_y[1] = Uk[3*nbCells_x*j+3*i+1]-Uk[3*nbCells_x*(j-1)+3*i+1] + dUi_y[2] = Uk[3*nbCells_x*j+3*i+2]-Uk[3*nbCells_x*(j-1)+3*i+2] + + temp1_x = Am_x*dUi1_x + temp2_x = Ap_x*dUi2_x + temp_y = Ap_y*dUi_y + 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]) + 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]) + 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]) + + #Traitement des autres cellules + else: + + #Calcul de Am_x + rho_l=Uk[3*nbCells_x*j+3*i] + qx_l =Uk[3*nbCells_x*j+3*i+1] + qy_l =Uk[3*nbCells_x*j+3*i+2] + rho_r=Uk[3*nbCells_x*j+3*(i+1)] + qx_r =Uk[3*nbCells_x*j+3*(i+1)+1] + qy_r =Uk[3*nbCells_x*j+3*(i+1)+2] + Am_x= jacobianMatricesm_x(dt/dx,rho_l,qx_l,qy_l,rho_r,qx_r,qy_r) + + #Calcul de Ap_x + rho_l=Uk[3*nbCells_x*j+3*(i-1)] + qx_l =Uk[3*nbCells_x*j+3*(i-1)+1] + qy_l =Uk[3*nbCells_x*j+3*(i-1)+2] + rho_r=Uk[3*nbCells_x*j+3*i] + qx_r =Uk[3*nbCells_x*j+3*i+1] + qy_r =Uk[3*nbCells_x*j+3*i+2] + Ap_x= jacobianMatricesp_x(dt/dx,rho_l,qx_l,qy_l,rho_r,qx_r,qy_r) + + #calcul de Am_y + rho_l=Uk[3*nbCells_x*j+3*i] + qx_l =Uk[3*nbCells_x*j+3*i+1] + qy_l =Uk[3*nbCells_x*j+3*i+2] + rho_r=Uk[3*nbCells_x*(j+1)+3*i] + qx_r =Uk[3*nbCells_x*(j+1)+3*i+1] + qy_r =Uk[3*nbCells_x*(j+1)+3*i+2] + Am_y= jacobianMatricesm_y(dt/dx,rho_l,qx_l,qy_l,rho_r,qx_r,qy_r) + + #calcul de Ap_y + rho_l=Uk[3*nbCells_x*(j-1)+3*i] + qx_l =Uk[3*nbCells_x*(j-1)+3*i+1] + qy_l =Uk[3*nbCells_x*(j-1)+3*i+2] + rho_r=Uk[3*nbCells_x*j+3*i] + qx_r =Uk[3*nbCells_x*j+3*i+1] + qy_r =Uk[3*nbCells_x*j+3*i+2] + Ap_y= jacobianMatricesp_y(dt/dx,rho_l,qx_l,qy_l,rho_r,qx_r,qy_r) + + #remplissage de la divMat + divMat.addValue(3*nbCells_x*j+3*i,3*nbCells_x*j+3*(i-1),Ap_x*(-1.)) + divMat.addValue(3*nbCells_x*j+3*i,3*nbCells_x*(j-1)+3*i,Ap_y*(-1.)) + divMat.addValue(3*nbCells_x*j+3*i,3*nbCells_x*j+3*(i+1),Am_x) + divMat.addValue(3*nbCells_x*j+3*i,3*nbCells_x*(j+1)+3*i,Am_y) + divMat.addValue(3*nbCells_x*j+3*i,3*nbCells_x*j+3*i,Ap_x+Ap_y-Am_x-Am_y) + + #remplissage du membre de droite + + dUi1_x[0] = Uk[3*nbCells_x*j+3*(i+1)]-Uk[3*nbCells_x*j+3*i] + dUi1_x[1] = Uk[3*nbCells_x*j+3*(i+1)+1]-Uk[3*nbCells_x*j+3*i+1] + dUi1_x[2] = Uk[3*nbCells_x*j+3*(i+1)+2]-Uk[3*nbCells_x*j+3*i+2] + + dUi2_x[0] = Uk[3*nbCells_x*j+3*i]-Uk[3*nbCells_x*j+3*(i-1)] + dUi2_x[1] = Uk[3*nbCells_x*j+3*i+1]-Uk[3*nbCells_x*j+3*(i-1)+1] + dUi2_x[2] = Uk[3*nbCells_x*j+3*i+2]-Uk[3*nbCells_x*j+3*(i-1)+2] + + dUi1_y[0] = Uk[3*nbCells_x*(j+1)+3*i]-Uk[3*nbCells_x*j+3*i] + dUi1_y[1] = Uk[3*nbCells_x*(j+1)+3*i+1]-Uk[3*nbCells_x*j+3*i+1] + dUi1_y[2] = Uk[3*nbCells_x*(j+1)+3*i+2]-Uk[3*nbCells_x*j+3*i+2] + + dUi2_y[0] = Uk[3*nbCells_x*j+3*i]-Uk[3*nbCells_x*(j-1)+3*i] + dUi2_y[1] = Uk[3*nbCells_x*j+3*i+1]-Uk[3*nbCells_x*(j-1)+3*i+1] + dUi2_y[2] = Uk[3*nbCells_x*j+3*i+2]-Uk[3*nbCells_x*(j-1)+3*i+2] + + temp1_x = Am_x*dUi1_x + temp2_x = Ap_x*dUi2_x + temp1_y = Am_y*dUi1_y + temp2_y = Ap_y*dUi2_y + 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]) + 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]) + 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]) + + divMat.diagonalShift(1) #only after filling all coefficients + LS=cdmath.LinearSolver(divMat,Rhs,iterGMRESMax, precision, "GMRES","LU") + dUk=LS.solve(); + residu = dUk.norm() + Uk+=dUk + #print("Newton iteration number ",k, "residu = ",residu) + 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 + isStationary = dUn.norm()=ntmax or isStationary or time >=tmax): + print("-- Iter: " + str(it) + ", Time: " + str(time) + ", dt: " + str(dt) + ", Newton iterations: "+str(k)+", ||dUn|| = "+str(dUn.norm())) + print("Linear system converged in ", LS.getNumberOfIter(), " GMRES iterations, résidu = ", residu) + for k in range(nbCells): + rho =rho_field[k] + velocity_field[k,0]=q_x_field[i]/rho + if(dim>1): + velocity_field[k,1]=q_y_field[k]/rho + print + rho_field.setTime(time,it); + rho_field.writeVTK("EulerIsothermal_"+str(dim)+"DConservativeStaggered_"+meshName+"_density",False); + q_x_field.setTime(time,it); + q_x_field.writeVTK("EulerIsothermal_"+str(dim)+"DConservativeStaggered_"+meshName+"_momentumX",False); + q_y_field.setTime(time,it); + q_y_field.writeVTK("EulerIsothermal_"+str(dim)+"DConservativeStaggered_"+meshName+"_momentumY",False); + velocity_field.setTime(time,it); + velocity_field.writeVTK("EulerIsothermal_"+str(dim)+"DConservativeStaggered_"+meshName+"_velocity",False); + #Postprocessing : save 2D picture + 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)) + 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)) + 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)) + 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("------------------------------------------------------------------------------------") + return + else: + print("Temps maximum Tmax= ", tmax, " atteint") + return + +def solve( my_mesh,filename, resolution): + print("Resolution of the Isothermal Euler system in dimension 2 on "+str(my_mesh.getNumberOfCells())+ " cells") + print("Initial data : ", "Riemann problem") + print("Boundary conditions : ", "Neumann") + print("Mesh name : ",filename ) + # Problem data + tmax = 10. + ntmax = 10 + cfl=1 + output_freq = 1 + EulerSystemStaggered(ntmax, tmax, cfl, output_freq, my_mesh, filename) + return + +if __name__ == """__main__""": + if len(sys.argv) >1 : + filename=sys.argv[1] + my_mesh = cdmath.Mesh(filename) + else : + filename = "CartesianGrid" + ax=0;bx=1;nx=20; + ay=0;by=1;ny=10; + my_mesh = cdmath.Mesh(ax,bx,nx,ay,by,ny) + + solve(my_mesh,filename,100)