--- /dev/null
+import EulerEquations1D
+import exact_rs_stiffenedgas
+
+import matplotlib
+matplotlib.use("Agg")
+import matplotlib.pyplot as plt
+import numpy as np
+from math import log10, sqrt
+import sys
+import time, json
+
+
+def test_validation1DEulerEquationsToro(scheme,isImplicit):
+ start = time.time()
+ #### 1D regular grid
+ meshList=[50,100,200,400]
+ meshType="1D regular grid"
+ testColor="Orange : à debugger"
+ nbMeshes=len(meshList)
+ mesh_name='RegularGrid'
+ simulation_name='Riemann problem'
+ a=0. ; b=1.
+ error_rho_tab=[0]*nbMeshes
+ error_v_tab =[0]*nbMeshes
+ error_p_tab =[0]*nbMeshes
+ error_q_tab =[0]*nbMeshes
+ error_h_tab =[0]*nbMeshes
+ error_e_tab =[0]*nbMeshes
+ time_tab =[0]*nbMeshes
+
+ plt.close('all')
+
+ # Initial parameters for Riemann problem (p in Pa, v in m/s, rho in Kg/m**3)
+
+ p_L = 460.894
+ p_R = 46.095
+ v_L = 19.5975
+ v_R = -6.19633
+ rho_L = 5.99924
+ rho_R = 5.99242
+
+ WL=[rho_L, v_L, p_L]
+ WR=[rho_R, v_R, p_R]
+
+ pos_disc = 0.4
+
+ precision = 1.e-5
+
+ # Stifened gas equation of state
+ gamma = 1.4
+ p0 = 0.
+ q = 0
+
+ tmax=0.035
+ lam_max= max(abs(v_L)+sqrt(gamma*p_L/rho_L), abs(v_R)+sqrt(gamma*p_R/rho_R) )
+
+ ###### Compute exact solution #####
+ RS = exact_rs_stiffenedgas.exact_rs_stiffenedgas(gamma, gamma, p0, p0);# Set the problem EOS
+ RS.solve_RP(WL,WR);# Solve the problem structure
+
+ numsamples=max(meshList) #nombre d'échantillons = taille du plus grand maillage
+ dx=(b-a)/numsamples
+
+ rho_field_exact=[0]*numsamples
+ v_field_exact =[0]*numsamples
+ p_field_exact =[0]*numsamples
+ q_field_exact =[0]*numsamples
+ h_field_exact =[0]*numsamples
+ e_field_exact =[0]*numsamples
+
+ for i in range(numsamples):
+
+ soln = RS.sample_solution(WL, WR, ( a+(i+1/2)*dx - pos_disc)/tmax); # Determine the solution ar time t in cell number i
+
+ rho_field_exact[i]=soln[0]
+ v_field_exact[i] =soln[1]
+ p_field_exact[i] =soln[2]
+ q_field_exact[i] =rho_field_exact[i]*v_field_exact[i]
+ h_field_exact[i] =exact_rs_stiffenedgas.stiffenedgas_h(rho_field_exact[i], p_field_exact[i], gamma, p0)
+ e_field_exact[i] =exact_rs_stiffenedgas.stiffenedgas_e(rho_field_exact[i], p_field_exact[i], gamma, p0)
+
+ #Set convergence picture parameters
+ max_initial_rho = max(rho_L,rho_R)
+ min_initial_rho = min(rho_L,rho_R)
+ max_initial_p = max(p_L,p_R)
+ min_initial_p = min(p_L,p_R)
+ max_initial_v = max(v_L,v_R)
+ min_initial_v = min(v_L,v_R)
+ max_initial_q = max_initial_rho*max_initial_v
+ min_initial_q = min_initial_rho*min_initial_v
+
+ e_L=exact_rs_stiffenedgas.p_to_e_StiffenedGaz(p_L, rho_L, gamma, p0)
+ e_R=exact_rs_stiffenedgas.p_to_e_StiffenedGaz(p_R, rho_R, gamma, p0)
+ h_L=e_L+p_L/rho_L
+ h_R=e_R+p_R/rho_R
+ max_initial_e = max(e_L, e_R)
+ min_initial_e = min(e_L, e_R)
+ min_initial_h = min(h_L,h_R)
+ max_initial_h = max(h_L,h_R)
+
+ fig, ([axDensity, axMomentum, axEnthalpie],[axPressure, axVitesse, axEinterne]) = plt.subplots(2, 3,sharex=True, figsize=(14,10))
+ plt.gcf().subplots_adjust(wspace = 0.5)
+
+ axDensity.plot( [a+0.5*dx + i*dx for i in range(numsamples)], rho_field_exact, label='Exact solution')
+ axDensity.set(xlabel='x (m)', ylabel='Densité (Kg/m3)')
+ axDensity.set_xlim(a,b)
+ axDensity.set_ylim(0.9*min_initial_rho , 6*max_initial_rho )
+ axDensity.legend()
+
+ axMomentum.plot( [a+0.5*dx + i*dx for i in range(numsamples)], q_field_exact, label='Exact solution')
+ axMomentum.set(xlabel='x (m)', ylabel='Momentum (Kg/m2/s)')
+ axMomentum.set_xlim(a,b)
+ axMomentum.set_ylim(0.9*min_initial_q , 2.5*max_initial_q )
+ axMomentum.legend()
+
+ axEnthalpie.plot( [a+0.5*dx + i*dx for i in range(numsamples)], h_field_exact, label='Exact solution')
+ axEnthalpie.set(xlabel='x (m)', ylabel='Enthalpy (J/Kg)')
+ axEnthalpie.set_xlim(a,b)
+ axEnthalpie.set_ylim(0.9*min_initial_h , 1.75*max_initial_h )
+ axEnthalpie.legend()
+
+ axPressure.plot( [a+0.5*dx + i*dx for i in range(numsamples)], p_field_exact, label='Exact solution')
+ axPressure.set(xlabel='x (m)', ylabel='Pression (Pa)')
+ axPressure.set_xlim(a,b)
+ axPressure.set_ylim(0.9*min_initial_p , 4*max_initial_p)
+ axPressure.ticklabel_format(axis='y', style='sci', scilimits=(0,0))
+ axPressure.legend()
+
+ axVitesse.plot( [a+0.5*dx + i*dx for i in range(numsamples)], v_field_exact, label='Exact solution')
+ axVitesse.set(xlabel='x (m)', ylabel='Vitesse (m/s)')
+ axVitesse.set_xlim(a,b)
+ axVitesse.set_ylim(0.9*min_initial_v , 1.1*max_initial_v)
+ axVitesse.ticklabel_format(axis='y', style='sci', scilimits=(0,0))
+ axVitesse.legend()
+
+ axEinterne.plot( [a+0.5*dx + i*dx for i in range(numsamples)], e_field_exact, label='Exact solution')
+ axEinterne.set(xlabel='x (m)', ylabel='Energie interne (J/Kg)')
+ axEinterne.set_xlim(a,b)
+ axEinterne.set_ylim(0.9*min_initial_e , 1.75*max_initial_e)
+ axEinterne.ticklabel_format(axis='y', style='sci', scilimits=(0,0))
+ axEinterne.legend()
+
+ ### Main loop to Store of numerical errors, mesh sizes and solutions
+ j=0
+ for nx in meshList:
+ p_field, v_field, e_field, rho_field, q_field, h_field, time_tab[j] = EulerEquations1D.solve(a, b, nx, mesh_name, meshType, isImplicit, scheme, simulation_name, testColor,fig,p_L,v_L,rho_L,p_R,v_R,rho_R, pos_disc, tmax, lam_max, gamma, p0)
+
+ axDensity.plot( [a+0.5*dx + i*dx for i in range(nx)], rho_field, label=scheme+' scheme on '+str(nx)+' cells') #new picture for video # Returns a tuple of line objects, thus the comma
+ axMomentum.plot( [a+0.5*dx + i*dx for i in range(nx)], q_field, label=scheme+' scheme on '+str(nx)+' cells')
+ axEnthalpie.plot([a+0.5*dx + i*dx for i in range(nx)], h_field, label=scheme+' scheme on '+str(nx)+' cells')
+ axPressure.plot( [a+0.5*dx + i*dx for i in range(nx)], p_field, label=scheme+' scheme on '+str(nx)+' cells')
+ axVitesse.plot( [a+0.5*dx + i*dx for i in range(nx)], v_field, label=scheme+' scheme on '+str(nx)+' cells')
+ axEinterne.plot( [a+0.5*dx + i*dx for i in range(nx)], e_field, label=scheme+' scheme on '+str(nx)+' cells')
+
+ for i in range(nx):
+ soln = RS.sample_solution(WL, WR, ( a+(i+1/2)*dx - pos_disc)/tmax); # Determine the solution at time t in cell number i
+
+ error_rho_tab[j]=max(error_rho_tab[j],abs(rho_field[i]-soln[0]))
+ error_v_tab[j] =max(error_v_tab[j], abs( v_field[i]-soln[1]))
+ error_p_tab[j] =max(error_p_tab[j], abs( p_field[i]-soln[2]))
+
+ j=j+1
+
+ end = time.time()
+
+ for i in range(nbMeshes):
+ error_rho_tab[i]=log10(error_rho_tab[i])
+ error_v_tab[i] =log10(error_v_tab[i])
+ error_p_tab[i] =log10(error_p_tab[i])
+ time_tab[i] =log10(time_tab[i])
+
+ if(isImplicit):
+ ImplicitOrExplicit="Implicit"
+ else:
+ ImplicitOrExplicit="Explicit"
+
+ # Plot of the 6 results fields
+ plt.suptitle('Euler equations : Convergence of ' + ImplicitOrExplicit + "_" + scheme + ' scheme ')
+ plt.savefig(mesh_name+"_1DEulerEquations_Toro4_" + ImplicitOrExplicit + "_" + scheme + ".png")
+ plt.close()
+
+
+ # Least square linear regression
+ # Find the best a,b such that f(x)=ax+b best approximates the convergence curve
+ # The vector X=(a,b) solves a symmetric linear system AX=B with A=(a1,a2\\a2,a3), B=(b1,b2)
+ a1=np.dot(meshList,meshList)
+ a2=np.sum(meshList)
+ a3=nbMeshes
+
+ det=a1*a3-a2*a2
+ assert det!=0, 'test_validation1DEulerEquation'+ImplicitOrExplicit+scheme+' : Make sure you use distinct meshes and at least two meshes'
+
+ b1r=np.dot(error_rho_tab,meshList)
+ b2r=np.sum(error_rho_tab)
+ ar=( a3*b1r-a2*b2r)/det
+ br=(-a2*b1r+a1*b2r)/det
+
+ print(ImplicitOrExplicit + "_" + scheme + " scheme for Euler Equations on 1D regular grid : order for density is ", -ar)
+
+ b1u=np.dot(error_v_tab,meshList)
+ b2u=np.sum(error_v_tab)
+ au=( a3*b1u-a2*b2u)/det
+ bu=(-a2*b1u+a1*b2u)/det
+
+ print(ImplicitOrExplicit + "_" + scheme + " scheme for Euler Equations on 1D regular grid : order for velocity is ", -au)
+
+ b1p=np.dot(error_p_tab,meshList)
+ b2p=np.sum(error_p_tab)
+ ap=( a3*b1p-a2*b2p)/det
+ bp=(-a2*b1p+a1*b2p)/det
+
+ print(ImplicitOrExplicit + "_" + scheme + " scheme for Euler Equations on 1D regular grid : order for velocity is ", -ap)
+
+ # Plot of convergence curves
+ fig_convergence, ([axDensity, axVitesse, axPressure]) = plt.subplots(1, 3,sharex=True, figsize=(14,5))
+ plt.gcf().subplots_adjust(wspace = 0.5)
+
+ plt.close()
+ axDensity.plot(meshList, error_rho_tab, label='log(|error density|)')
+ axDensity.plot(meshList, a*np.array(meshList)+b,label='least square slope : '+'%.3f' % ar)
+ axDensity.legend()
+ axDensity.set(xlabel='log(Number of cells)', ylabel='log(|error density|)')
+ #axDensity.title('Mesh convergence of density')
+
+ axVitesse.plot(meshList, error_v_tab, label='log(|error velocity|)')
+ axVitesse.plot(meshList, a*np.array(meshList)+b,label='least square slope : '+'%.3f' % au)
+ axVitesse.legend()
+ axVitesse.set(xlabel='log(Number of cells)', ylabel='log(|error velocity|)')
+ #axVitesse.title('Mesh convergence of Velocity')
+
+ axPressure.plot(meshList, error_p_tab, label='log(|error Pressure|)')
+ axPressure.plot(meshList, a*np.array(meshList)+b,label='least square slope : '+'%.3f' % ap)
+ axPressure.legend()
+ axPressure.set(xlabel='log(Number of cells)', ylabel='log(|error Pressure|)')
+ #axPressure.title('Mesh convergence of Pressure')
+
+ plt.suptitle('Convergence of finite volumes for the Euler equations : Toro 4 \n Riemann problem with '+ImplicitOrExplicit+scheme+' scheme on a 1D regular grid')
+ plt.savefig(mesh_name+"1DEulerEquation_"+ImplicitOrExplicit+scheme+"_ConvergenceCurves.png")
+
+ # Plot of computational time
+ plt.close()
+ plt.plot(meshList, time_tab, label='log(cpu time)')
+ plt.legend()
+ plt.xlabel('log(Number of cells)')
+ plt.ylabel('log(cpu time)')
+ plt.title('Computational time of finite volumes for the Euler equations : Toro 4 \n Riemann problem with '+ImplicitOrExplicit+scheme+' scheme on a 1D regular grid')
+ plt.savefig(mesh_name+"1DEulerEquation_"+ImplicitOrExplicit+scheme+"_ComputationalTime.png")
+
+ plt.close('all')
+
+ convergence_synthesis={}
+
+ convergence_synthesis["PDE_model"]="Euler_Equations"
+ convergence_synthesis["PDE_is_stationary"]=False
+ convergence_synthesis["PDE_search_for_stationary_solution"]=True
+ convergence_synthesis["Numerical_method_name"]="Upwind scheme"
+ convergence_synthesis["Numerical_method_space_discretization"]="Finite volumes"
+ convergence_synthesis["Numerical_method_time_discretization"]=ImplicitOrExplicit
+ convergence_synthesis["Initial_data"]="Riemann Problem"
+ convergence_synthesis["Boundary_conditions"]="Neumann"
+ convergence_synthesis["Space_dimension"]=1
+ convergence_synthesis["Mesh_dimension"]=1
+ convergence_synthesis["Mesh_names"]=meshList
+ convergence_synthesis["Mesh_type"]=meshType
+ convergence_synthesis["Mesh_description"]=mesh_name
+ convergence_synthesis["Mesh_sizes"]=meshList
+ convergence_synthesis["Mesh_cell_type"]="1D regular grid"
+ convergence_synthesis["Scheme_order_density"]=-ar
+ convergence_synthesis["Scheme_order_velocity"]=-au
+ convergence_synthesis["Scheme_order_pressure"]=-ap
+ convergence_synthesis["Test_color"]=testColor
+ convergence_synthesis["Computational_time"]=end-start
+
+ with open('Convergence_1DEulerEquations_Toro4_'+ImplicitOrExplicit+scheme+mesh_name+'.json', 'w') as outfile:
+ json.dump(convergence_synthesis, outfile)
+
+if __name__ == """__main__""":
+ if len(sys.argv) >2 :
+ scheme = sys.argv[1]
+ isImplicit = bool(int(sys.argv[2]))
+ test_validation1DEulerEquationsToro(scheme,isImplicit)
+ else :
+ test_validation1DEulerEquationsToro('Roe',False)
+
--- /dev/null
+#!/usr/bin/env python3\r
+# -*-coding:utf-8 -*\r
+\r
+"""\r
+Created on Mon Aug 30 2021\r
+@author: Michael NDJINGA, Katia Ait Ameur, Coraline Mounier\r
+\r
+Euler system without heating source term (phi) in one dimension on regular domain [a,b]\r
+Ghost cell (Neumann) boundary condition\r
+Roe scheme or conservative MAC scheme\r
+Regular square mesh\r
+\r
+State law Stiffened gaz : p = (gamma - 1) * rho * (e - q) - gamma * pinf\r
+4 choices of parameters gamma and pinf are available : \r
+ - Lagrange interpolation (with q=0)\r
+ - Hermite interpolation with reference point at 575K (with q=0)\r
+ - Hermite interpolation with reference point at 590K (with q=0)\r
+ - Hermite interpolation with reference point at 617.94K (saturation at 155 bar) with q=0\r
+ \r
+Linearized enthalpy : h = h_sat + cp * (T - T_sat)\r
+Values for cp and T_sat parameters are taken at the reference point chosen for the state law\r
+\r
+To do correct the computation of the time step : lambda_max (maximum eigenvalue) should be computed first)\r
+"""\r
+\r
+import cdmath\r
+import numpy as np\r
+import matplotlib\r
+\r
+matplotlib.use("Agg")\r
+import matplotlib.pyplot as plt\r
+import sys, time, json\r
+from math import sqrt\r
+from numpy import sign\r
+\r
+\r
+precision = 1.e-5\r
+\r
+def p_to_e_StiffenedGaz(p, rho, gamma, pinf):\r
+ e_field = (p + gamma*pinf) / (gamma - 1.) / rho\r
+ return e_field\r
+\r
+def initial_conditions_Riemann_problem(a, b, nx,p_L,v_L,rho_L,p_R,v_R,rho_R,pos_disc, gamma, pinf):\r
+ print("Initial data Riemann problem")\r
+ dx = (b - a) / nx # space step\r
+ x = [a + 0.5 * dx + i * dx for i in range(nx)] # array of cell center (1D mesh) \r
+\r
+ p_initial = np.array([ (xi < pos_disc) * p_L + (xi >= pos_disc) * p_R for xi in x])\r
+ v_initial = np.array([ (xi < pos_disc) * v_L + (xi >= pos_disc) * v_R for xi in x])\r
+ rho_initial = np.array([ (xi < pos_disc) * rho_L + (xi >= pos_disc) * rho_R for xi in x])\r
+\r
+ e_initial = p_to_e_StiffenedGaz(p_initial, rho_initial, gamma, pinf)\r
+ q_initial = rho_initial * v_initial\r
+ rho_E_initial = rho_initial*e_initial + 0.5 * q_initial**2/rho_initial\r
+\r
+ return rho_initial, q_initial, rho_E_initial, p_initial, v_initial, e_initial\r
+\r
+def p_to_rho_StiffenedGaz(p_field, T_field, gamma, pinf):\r
+ rho_field = (p_field + pinf) * gamma / (gamma - 1) * 1. / (h_sat + cp * (T_field - T_sat))\r
+ return rho_field\r
+ \r
+def T_to_rhoE_StiffenedGaz(T_field, rho_field, q_field, gamma, pinf):\r
+ rho_E_field = 1. / 2. * (q_field) ** 2 / rho_field + pinf + rho_field / gamma * (h_sat + cp * (T_field- T_sat))\r
+ return rho_E_field\r
+\r
+def rho_to_p_StiffenedGaz(rho_field, q_field, rho_E_field, gamma, pinf):\r
+ p_field = (gamma - 1) * (rho_E_field - 1. / 2 * q_field ** 2 / rho_field) - gamma * pinf\r
+ return p_field\r
+\r
+def T_to_E_StiffenedGaz(p_field, T_field, v_field, gamma, pinf):\r
+ rho_field = p_to_rho_StiffenedGaz(p_field, T_field)\r
+ E_field = (p_field + gamma * pinf) / ((gamma-1) * rho_field) + 0.5 * v_field **2\r
+ return E_field\r
+ \r
+def dp_drho_e_const_StiffenedGaz( e , gamma, pinf):\r
+ return (gamma-1)*(e)\r
+\r
+def dp_de_rho_const_StiffenedGaz( rho, gamma, pinf ):\r
+ return (gamma-1)*rho\r
+\r
+def sound_speed_StiffenedGaz( h, gamma, pinf ):\r
+ return np.sqrt((gamma-1)* h)\r
+\r
+def rho_h_to_p_StiffenedGaz( rho, h, gamma, pinf ):\r
+ return (gamma - 1) * rho * h / gamma - pinf\r
+\r
+def rhoE_to_e_StiffenedGaz(rho_field, q_field, rho_E_field):\r
+ e_field = rho_E_field / rho_field - 1. / 2. * (q_field / rho_field)**2\r
+ return e_field\r
+\r
+def MatRoe_StiffenedGaz( rho_l, q_l, rho_E_l, rho_r, q_r, rho_E_r, gamma, pinf):\r
+ RoeMat = cdmath.Matrix(3, 3)\r
+\r
+ u_l = q_l / rho_l\r
+ u_r = q_r / rho_r\r
+ p_l = rho_to_p_StiffenedGaz(rho_l, q_l, rho_E_l, gamma, pinf)\r
+ p_r = rho_to_p_StiffenedGaz(rho_r, q_r, rho_E_r, gamma, pinf)\r
+ H_l = rho_E_l / rho_l + p_l / rho_l\r
+ H_r = rho_E_r / rho_r + p_r / rho_r\r
+\r
+ # Roe averages\r
+ rho = np.sqrt(rho_l * rho_r)\r
+ u = (u_l * sqrt(rho_l) + u_r * sqrt(rho_r)) / (sqrt(rho_l) + sqrt(rho_r))\r
+ H = (H_l * sqrt(rho_l) + H_r * sqrt(rho_r)) / (sqrt(rho_l) + sqrt(rho_r))\r
+\r
+ p = rho_h_to_p_StiffenedGaz( rho, H - u**2/2., gamma, pinf )\r
+ e = H - p / rho - 1./2 * u**2\r
+ dp_drho = dp_drho_e_const_StiffenedGaz( e, gamma, pinf )\r
+ dp_de = dp_de_rho_const_StiffenedGaz( rho, gamma, pinf )\r
+\r
+ RoeMat[0, 0] = 0\r
+ RoeMat[0, 1] = 1\r
+ RoeMat[0, 2] = 0\r
+ RoeMat[1, 0] = dp_drho - u ** 2 + dp_de / rho * (u**2/2 - e)\r
+ RoeMat[1, 1] = 2 * u - u * dp_de / rho\r
+ RoeMat[1, 2] = dp_de / rho\r
+ RoeMat[2, 0] = -u * ( -dp_drho + H - dp_de / rho * (u**2/2 - e) )\r
+ RoeMat[2, 1] = H - dp_de / rho * u ** 2\r
+ RoeMat[2, 2] = (dp_de / rho +1) * u\r
+ \r
+ return(RoeMat)\r
+\r
+ \r
+def Dmac_StiffenedGaz( rho_l, q_l, rho_E_l, rho_r, q_r, rho_E_r, gamma, pinf):\r
+ Dmac = cdmath.Matrix(3, 3)\r
+\r
+ u_l = q_l / rho_l\r
+ u_r = q_r / rho_r\r
+ p_l = rho_to_p_StiffenedGaz(rho_l, q_l, rho_E_l, gamma, pinf)\r
+ p_r = rho_to_p_StiffenedGaz(rho_r, q_r, rho_E_r, gamma, pinf)\r
+ H_l = rho_E_l / rho_l + p_l / rho_l\r
+ H_r = rho_E_r / rho_r + p_r / rho_r\r
+\r
+ # Roe averages\r
+ rho = np.sqrt(rho_l * rho_r)\r
+ u = (u_l * sqrt(rho_l) + u_r * sqrt(rho_r)) / (sqrt(rho_l) + sqrt(rho_r))\r
+ H = (H_l * sqrt(rho_l) + H_r * sqrt(rho_r)) / (sqrt(rho_l) + sqrt(rho_r))\r
+\r
+ p = rho_h_to_p_StiffenedGaz( rho, H - u**2/2., gamma, pinf )\r
+ e = H - p / rho - 1./2 * u**2\r
+ dp_drho = dp_drho_e_const_StiffenedGaz( e, gamma, pinf )\r
+ dp_de = dp_de_rho_const_StiffenedGaz( rho, gamma, pinf )\r
+ \r
+ #Third choice for Dstag\r
+ Dmac[0, 0] = 0\r
+ Dmac[0, 1] = 1\r
+ Dmac[0, 2] = 0\r
+ Dmac[1, 0] = -dp_drho - u ** 2 - dp_de / rho * (u**2/2 - e)\r
+ Dmac[1, 1] = 2 * u + u * dp_de / rho\r
+ Dmac[1, 2] = -dp_de / rho\r
+ Dmac[2, 0] = -u * ( dp_drho + H + dp_de / rho * (u**2/2 - e) )\r
+ Dmac[2, 1] = H + dp_de / rho * u ** 2\r
+ Dmac[2, 2] = (-dp_de / rho +1) * u\r
+ \r
+ #return Dmac * sign(u)\r
+ \r
+ #Fifth choice for Dstag\r
+ Dmac[0, 0] = abs(u)-u\r
+ Dmac[0, 1] = 1\r
+ Dmac[0, 2] = 0\r
+ Dmac[1, 0] = -dp_drho - u ** 2 - dp_de / rho * (u**2/2 - e)\r
+ Dmac[1, 1] = abs(u) + u + u * dp_de / rho\r
+ Dmac[1, 2] = -dp_de / rho\r
+ Dmac[2, 0] = -u * ( dp_drho + H + u*(u-abs(u)) + dp_de / rho * (u**2/2 - e) )\r
+ Dmac[2, 1] = H + dp_de / rho * u ** 2\r
+ Dmac[2, 2] = -u*dp_de / rho + abs(u)\r
+ \r
+ return Dmac \r
+\r
+def Droe_StiffenedGaz( rho_l, q_l, rho_E_l, rho_r, q_r, rho_E_r, gamma, pinf):\r
+ Droe = cdmath.Matrix(3, 3)\r
+\r
+ u_l = q_l / rho_l\r
+ u_r = q_r / rho_r\r
+ p_l = rho_to_p_StiffenedGaz(rho_l, q_l, rho_E_l, gamma, pinf)\r
+ p_r = rho_to_p_StiffenedGaz(rho_r, q_r, rho_E_r, gamma, pinf)\r
+ H_l = rho_E_l / rho_l + p_l / rho_l\r
+ H_r = rho_E_r / rho_r + p_r / rho_r\r
+\r
+ # Roe averages\r
+ rho = np.sqrt(rho_l * rho_r)\r
+ u = (u_l * sqrt(rho_l) + u_r * sqrt(rho_r)) / (sqrt(rho_l) + sqrt(rho_r))\r
+ H = (H_l * sqrt(rho_l) + H_r * sqrt(rho_r)) / (sqrt(rho_l) + sqrt(rho_r))\r
+\r
+ c = sound_speed_StiffenedGaz( H - u**2/2., gamma, pinf )\r
+ \r
+ lamb = cdmath.Vector(3)\r
+ lamb[0] = u-c\r
+ lamb[1] = u\r
+ lamb[2] = u+c \r
+\r
+ r = cdmath.Matrix(3, 3)\r
+ r[0,0] = 1.\r
+ r[1,0] = u-c\r
+ r[2,0] = H-u*c \r
+ r[0,1] = 1.\r
+ r[1,1] = u \r
+ r[2,1] = H-c**2/(gamma-1) \r
+ r[0,2] = 1.\r
+ r[1,2] = u+c\r
+ r[2,2] = H+u*c \r
+\r
+ l = cdmath.Matrix(3, 3)\r
+ l[0,0] = (1./(2*c**2))*(0.5*(gamma-1)*u**2+u*c)\r
+ l[1,0] = (1./(2*c**2))*(-u*(gamma-1)-c)\r
+ l[2,0] = (1./(2*c**2))*(gamma-1)\r
+ l[0,1] = ((gamma-1)/c**2)*(H-u**2)\r
+ l[1,1] = ((gamma-1)/c**2)*u \r
+ l[2,1] = -((gamma-1)/c**2) \r
+ l[0,2] = (1./(2*c**2))*(0.5*(gamma-1)*u**2-u*c)\r
+ l[1,2] = (1./(2*c**2))*(c-u*(gamma-1))\r
+ l[2,2] = (1./(2*c**2))*(gamma-1)\r
+\r
+ M1 = cdmath.Matrix(3, 3) #abs(lamb[0])*r[:,0].tensProduct(l[:,0])\r
+ M2 = cdmath.Matrix(3, 3) #abs(lamb[1])*r[:,1].tensProduct(l[:,1]) \r
+ M3 = cdmath.Matrix(3, 3) #abs(lamb[2])*r[:,2].tensProduct(l[:,2])\r
+ for i in range(3):\r
+ for j in range(3):\r
+ M1[i,j] = abs(lamb[0])*r[i,0]*l[j,0]\r
+ M2[i,j] = abs(lamb[1])*r[i,1]*l[j,1] \r
+ M3[i,j] = abs(lamb[2])*r[i,2]*l[j,2]\r
+ \r
+ Droe = M1+M2+M3 \r
+ \r
+ return(Droe) \r
+\r
+\r
+def jacobianMatricesm(coeff, rho_l, q_l, rho_E_l, rho_r, q_r, rho_E_r, scheme, gamma, pinf):\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
+ if rho_E_l < 0 or rho_E_r < 0:\r
+ print("rho_E_l=", rho_E_l, " rho_E_r= ", rho_E_r)\r
+ raise ValueError("Negative total energy")\r
+\r
+ RoeMat = MatRoe_StiffenedGaz( rho_l, q_l, rho_E_l, rho_r, q_r, rho_E_r, gamma, pinf)\r
+ \r
+ if scheme == 'Stag':\r
+ Dmac = Dmac_StiffenedGaz( rho_l, q_l, rho_E_l, rho_r, q_r, rho_E_r, gamma, pinf)\r
+ return (RoeMat - Dmac) * coeff * 0.5\r
+ \r
+ elif scheme == 'Roe':\r
+ Droe = Droe_StiffenedGaz( rho_l, q_l, rho_E_l, rho_r, q_r, rho_E_r, gamma, pinf) \r
+ return (RoeMat - Droe) * coeff * 0.5\r
+ else:\r
+ raise ValueError("Wrong scheme given : ", scheme, ". Accepted scheme are Roe and Stag")\r
+\r
+\r
+def jacobianMatricesp(coeff, rho_l, q_l, rho_E_l, rho_r, q_r, rho_E_r, scheme, gamma, pinf):\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
+ if rho_E_l < 0 or rho_E_r < 0:\r
+ print("rho_E_l=", rho_E_l, " rho_E_r= ", rho_E_r)\r
+ raise ValueError("Negative total energy")\r
+\r
+ RoeMat = MatRoe_StiffenedGaz( rho_l, q_l, rho_E_l, rho_r, q_r, rho_E_r, gamma, pinf)\r
+\r
+ if scheme == 'Stag':\r
+ Dmac = Dmac_StiffenedGaz( rho_l, q_l, rho_E_l, rho_r, q_r, rho_E_r, gamma, pinf)\r
+ return (RoeMat - Dmac) * coeff * 0.5\r
+ \r
+ elif scheme == 'Roe':\r
+ Droe = Droe_StiffenedGaz( rho_l, q_l, rho_E_l, rho_r, q_r, rho_E_r, gamma, pinf) \r
+ return (RoeMat - Droe) * coeff * 0.5\r
+\r
+ else:\r
+ raise ValueError("Wrong scheme given : ", scheme, ". Accepted scheme are Roe and Stag")\r
+\r
+def FillEdges(nx, j, Uk, nbComp, divMat, Rhs, Un, dt, dx, isImplicit, scheme, gamma, pinf):\r
+ dUi1 = cdmath.Vector(3)\r
+ dUi2 = cdmath.Vector(3)\r
+ temp1 = cdmath.Vector(3)\r
+ temp2 = cdmath.Vector(3)\r
+\r
+ if (j == 0):\r
+ rho_l = Uk[j * nbComp + 0]\r
+ q_l = Uk[j * nbComp + 1]\r
+ rho_E_l = Uk[j * nbComp + 2]\r
+ rho_r = Uk[(j + 1) * nbComp + 0]\r
+ q_r = Uk[(j + 1) * nbComp + 1]\r
+ rho_E_r = Uk[(j + 1) * nbComp + 2] \r
+\r
+ Am = jacobianMatricesm(dt / dx, rho_l, q_l, rho_E_l, rho_r, q_r, rho_E_r, scheme, gamma, pinf)\r
+ divMat.addValue(j * nbComp, (j + 1) * nbComp, Am)\r
+ divMat.addValue(j * nbComp, j * nbComp, Am * (-1.))\r
+\r
+ Ap = jacobianMatricesp(dt / dx, rho_l, q_l, rho_E_l, rho_l, q_l, rho_E_l, scheme, gamma, pinf)\r
+ divMat.addValue(j * nbComp, j * nbComp, Ap)\r
+ \r
+ if(isImplicit):\r
+ dUi1[0] = Uk[(j + 1) * nbComp + 0] - Uk[j * nbComp + 0]\r
+ dUi1[1] = Uk[(j + 1) * nbComp + 1] - Uk[j * nbComp + 1]\r
+ dUi1[2] = Uk[(j + 1) * nbComp + 2] - Uk[j * nbComp + 2]\r
+ temp1 = Am * dUi1\r
+ else:\r
+ dUi2[0] = rho_l \r
+ dUi2[1] = q_l \r
+ dUi2[2] = rho_E_l \r
+ temp2 = Ap * dUi2\r
+\r
+ elif (j == nx - 1):\r
+ rho_l = Uk[(j - 1) * nbComp + 0]\r
+ q_l = Uk[(j - 1) * nbComp + 1]\r
+ rho_E_l = Uk[(j - 1) * nbComp + 2]\r
+ rho_r = Uk[j * nbComp + 0]\r
+ q_r = Uk[j * nbComp + 1]\r
+ rho_E_r = Uk[j * nbComp + 2]\r
+\r
+ Ap = jacobianMatricesp(dt / dx, rho_l, q_l, rho_E_l, rho_r, q_r, rho_E_r, scheme, gamma, pinf)\r
+ divMat.addValue(j * nbComp, j * nbComp, Ap)\r
+ divMat.addValue(j * nbComp, (j - 1) * nbComp, Ap * (-1.))\r
+\r
+ Am = jacobianMatricesm(dt / dx, rho_r, q_r, rho_E_r, rho_r, q_r, rho_E_r, scheme, gamma, pinf) \r
+ divMat.addValue(j * nbComp, j * nbComp, Am * (-1.))\r
+\r
+ if(isImplicit):\r
+ dUi2[0] = Uk[(j - 1) * nbComp + 0] - Uk[j * nbComp + 0]\r
+ dUi2[1] = Uk[(j - 1) * nbComp + 1] - Uk[j * nbComp + 1]\r
+ dUi2[2] = Uk[(j - 1) * nbComp + 2] - Uk[j * nbComp + 2]\r
+ temp2 = Ap * dUi2\r
+ else:\r
+ dUi1[0] = rho_r \r
+ dUi1[1] = q_r \r
+ dUi1[2] = rho_E_r \r
+ temp1 = Am * dUi1\r
+ \r
+ if(isImplicit):#implicit scheme, contribution from the Newton scheme residual\r
+ Rhs[j * nbComp + 0] = -temp1[0] + temp2[0] - (Uk[j * nbComp + 0] - Un[j * nbComp + 0])\r
+ Rhs[j * nbComp + 1] = -temp1[1] + temp2[1] - (Uk[j * nbComp + 1] - Un[j * nbComp + 1])\r
+ Rhs[j * nbComp + 2] = -temp1[2] + temp2[2] - (Uk[j * nbComp + 2] - Un[j * nbComp + 2])\r
+ else:#explicit scheme, contribution from the boundary data the right hand side\r
+ Rhs[j * nbComp + 0] = -temp1[0] + temp2[0] \r
+ Rhs[j * nbComp + 1] = -temp1[1] + temp2[1] \r
+ Rhs[j * nbComp + 2] = -temp1[2] + temp2[2] \r
+\r
+def FillInnerCell(j, Uk, nbComp, divMat, Rhs, Un, dt, dx, isImplicit, scheme, gamma, pinf):\r
+\r
+ rho_l = Uk[(j - 1) * nbComp + 0]\r
+ q_l = Uk[(j - 1) * nbComp + 1]\r
+ rho_E_l = Uk[(j - 1) * nbComp + 2]\r
+ rho_r = Uk[j * nbComp + 0]\r
+ q_r = Uk[j * nbComp + 1]\r
+ rho_E_r = Uk[j * nbComp + 2]\r
+ Ap = jacobianMatricesp(dt / dx, rho_l, q_l, rho_E_l, rho_r, q_r, rho_E_r, scheme, gamma, pinf)\r
+ \r
+ rho_l = Uk[j * nbComp + 0]\r
+ q_l = Uk[j * nbComp + 1]\r
+ rho_E_l = Uk[j * nbComp + 2]\r
+ rho_r = Uk[(j + 1) * nbComp + 0]\r
+ q_r = Uk[(j + 1) * nbComp + 1]\r
+ rho_E_r = Uk[(j + 1) * nbComp + 2]\r
+ Am = jacobianMatricesm(dt / dx, rho_l, q_l, rho_E_l, rho_r, q_r, rho_E_r, scheme, gamma, pinf)\r
+\r
+ divMat.addValue(j * nbComp, (j + 1) * nbComp, Am)\r
+ divMat.addValue(j * nbComp, j * nbComp, Am * (-1.))\r
+ divMat.addValue(j * nbComp, j * nbComp, Ap)\r
+ divMat.addValue(j * nbComp, (j - 1) * nbComp, Ap * (-1.))\r
+\r
+ if(isImplicit):\r
+ dUi1 = cdmath.Vector(3)\r
+ dUi2 = cdmath.Vector(3)\r
+ \r
+ dUi1[0] = Uk[(j + 1) * nbComp + 0] - Uk[j * nbComp + 0]\r
+ dUi1[1] = Uk[(j + 1) * nbComp + 1] - Uk[j * nbComp + 1]\r
+ dUi1[2] = Uk[(j + 1) * nbComp + 2] - Uk[j * nbComp + 2]\r
+ \r
+ dUi2[0] = Uk[(j - 1) * nbComp + 0] - Uk[j * nbComp + 0]\r
+ dUi2[1] = Uk[(j - 1) * nbComp + 1] - Uk[j * nbComp + 1]\r
+ dUi2[2] = Uk[(j - 1) * nbComp + 2] - Uk[j * nbComp + 2]\r
+ \r
+ temp1 = Am * dUi1\r
+ temp2 = Ap * dUi2\r
+\r
+ Rhs[j * nbComp + 0] = -temp1[0] + temp2[0] - (Uk[j * nbComp + 0] - Un[j * nbComp + 0])\r
+ Rhs[j * nbComp + 1] = -temp1[1] + temp2[1] - (Uk[j * nbComp + 1] - Un[j * nbComp + 1])\r
+ Rhs[j * nbComp + 2] = -temp1[2] + temp2[2] - (Uk[j * nbComp + 2] - Un[j * nbComp + 2])\r
+ else:\r
+ Rhs[j * nbComp + 0] = 0\r
+ Rhs[j * nbComp + 1] = 0\r
+ Rhs[j * nbComp + 2] = 0\r
+\r
+\r
+def EulerSystem(ntmax, tmax, cfl, a, b, nbCells, output_freq, meshName, isImplicit, scheme,fig,p_L,v_L,rho_L,p_R,v_R,rho_R, lam_max, gamma, pinf,pos_disc):\r
+ dim = 1\r
+ nbComp = 3\r
+ dt = 0.\r
+ time = 0.\r
+ it = 0\r
+ isStationary = False\r
+ dx = (b - a) / nbCells\r
+ dt = cfl * dx / lam_max\r
+ nbVoisinsMax = 2\r
+ \r
+ e_L=p_to_e_StiffenedGaz(p_L, rho_L, gamma, pinf)\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
+ # Initial conditions\r
+ print("Construction of the initial condition …")\r
+\r
+ rho_field, q_field, rho_E_field, p_field, v_field, e_field = initial_conditions_Riemann_problem(a, b, nbCells, p_L,v_L,rho_L,p_R,v_R,rho_R,pos_disc, gamma, pinf)\r
+ h_field = (rho_E_field + p_field) / rho_field - 0.5 * (q_field / rho_field) **2\r
+ p_field = p_field * 10 ** (-5)\r
+ \r
+\r
+ for k in range(nbCells):\r
+ Un[k * nbComp + 0] = rho_field[k]\r
+ Un[k * nbComp + 1] = q_field[k]\r
+ Un[k * nbComp + 2] = rho_E_field[k]\r
+\r
+ divMat = cdmath.SparseMatrixPetsc(nbCells * nbComp, nbCells * nbComp, (nbVoisinsMax + 1) * nbComp)\r
+ \r
+ iterGMRESMax = 50\r
+ newton_max = 100\r
+\r
+ print("Starting computation of the non linear Euler non isentropic system with scheme : ", scheme)\r
+ # STARTING TIME LOOP\r
+ while (it < ntmax and time <= tmax and not isStationary):\r
+ dUn = Un.deepCopy()\r
+ Uk = Un.deepCopy()\r
+ residu = 1.\r
+ \r
+ k = 0\r
+ while (k < newton_max and residu > precision):\r
+ # STARTING NEWTON LOOP\r
+ divMat.zeroEntries() #sets the matrix coefficients to zero\r
+ for j in range(nbCells):\r
+ \r
+ # traitements des bords\r
+ if (j == 0):\r
+ FillEdges(nbCells,j, Uk, nbComp, divMat, Rhs, Un, dt, dx, isImplicit, scheme, gamma, pinf)\r
+ elif (j == nbCells - 1):\r
+ FillEdges(nbCells,j, Uk, nbComp, divMat, Rhs, Un, dt, dx, isImplicit, scheme, gamma, pinf)\r
+\r
+ # traitement des cellules internes\r
+ else:\r
+ FillInnerCell( j, Uk, nbComp, divMat, Rhs, Un, dt, dx, isImplicit, scheme, gamma, pinf)\r
+ \r
+ if(isImplicit):\r
+ #solving the linear system divMat * dUk = Rhs\r
+ divMat.diagonalShift(1.)\r
+ LS = cdmath.LinearSolver(divMat, Rhs, iterGMRESMax, precision, "GMRES", "LU")\r
+ dUk = LS.solve()\r
+ vector_residu = dUk.maxVector(nbComp)\r
+ residu = max(abs(vector_residu[0])/rho_L, abs(vector_residu[1])/(rho_L*v_L), abs(vector_residu[2])/rho_L*e_L )\r
+ else:\r
+ dUk=Rhs - divMat*Un\r
+ residu = 0.#Convergence schéma Newton\r
+ \r
+ if (isImplicit and (it == 1 or it % output_freq == 0 or it >= ntmax or isStationary or time >= tmax)):\r
+ print("Residu Newton at iteration ",k, " : ", residu)\r
+ print("Linear system converged in ", LS.getNumberOfIter(), " GMRES iterations")\r
+\r
+ #updates for Newton loop\r
+ Uk += dUk\r
+ k = k + 1\r
+ if (isImplicit and not LS.getStatus()):\r
+ print("Linear system did not converge ", LS.getNumberOfIter(), " GMRES iterations")\r
+ raise ValueError("No convergence of the linear system")\r
+ \r
+ if k == newton_max:\r
+ raise ValueError("No convergence of Newton algorithm for " + scheme + " scheme")\r
+\r
+ #updating fields\r
+ Un = Uk.deepCopy()\r
+ dUn -= Un\r
+\r
+ #Testing stationarity\r
+ residu_stat = dUn.maxVector(nbComp)#On prend le max de chaque composante\r
+ if (it % output_freq == 0 ):\r
+ print("Test de stationarité : Un+1-Un= ", max(abs(residu_stat[0])/rho_L, abs(residu_stat[1])/(rho_L*v_L), abs(residu_stat[2])/rho_L*e_L ))\r
+\r
+ if ( it>1 and abs(residu_stat[0])/rho_L<precision and abs(residu_stat[1])/(rho_L*v_L)<precision and abs(residu_stat[2])/rho_L*e_L<precision):\r
+ isStationary = True\r
+ \r
+ for k in range(nbCells):\r
+ rho_field[k] = Un[k * nbComp + 0]\r
+ q_field[k] = Un[k * nbComp + 1]\r
+ rho_E_field[k] = Un[k * nbComp + 2]\r
+\r
+ v_field = q_field / rho_field\r
+ p_field = rho_to_p_StiffenedGaz(rho_field, q_field, rho_E_field, gamma, pinf)\r
+ e_field = rhoE_to_e_StiffenedGaz(rho_field, q_field, rho_E_field )\r
+ h_field = (rho_E_field + p_field) / rho_field - 0.5 * (q_field / rho_field) **2\r
+ p_field = p_field * 10 ** (-5)\r
+ \r
+ if( min(p_field)<0) :\r
+ raise ValueError("Negative pressure, stopping calculation")\r
+\r
+ time = time + dt\r
+ it = it + 1\r
+\r
+ # Savings\r
+ if (it == 1 or it % output_freq == 0 or it >= ntmax or isStationary or time >= tmax):\r
+ print("-- Time step : " + str(it) + ", Time: " + str(time) + ", dt: " + str(dt))\r
+\r
+ print("-- Iter: " + str(it) + ", Time: " + str(time) + ", dt: " + str(dt)+"\n")\r
+ \r
+ if(isImplicit):\r
+ test_desc["Linear_solver_algorithm"]=LS.getNameOfMethod()\r
+ test_desc["Linear_solver_preconditioner"]=LS.getNameOfPc()\r
+ test_desc["Linear_solver_precision"]=LS.getTolerance()\r
+ test_desc["Linear_solver_maximum_iterations"]=LS.getNumberMaxOfIter()\r
+ test_desc["Linear_system_max_actual_iterations_number"]=LS.getNumberOfIter()\r
+ test_desc["Linear_system_max_actual_error"]=LS.getResidu()\r
+ test_desc["Linear_system_max_actual_condition number"]=LS.getConditionNumber()\r
+\r
+ if (it >= ntmax):\r
+ print("Maximum number of time steps ntmax= ", ntmax, " reached")\r
+\r
+ elif (isStationary):\r
+ print("Stationary regime reached at time step ", it, ", t= ", time)\r
+ print("------------------------------------------------------------------------------------")\r
+ np.savetxt("EulerComplet_" + str(dim) + "_" + scheme + meshName + "_rho_Stat.txt", rho_field, delimiter="\n")\r
+ np.savetxt("EulerComplet_" + str(dim) + "_" + scheme + meshName + "_q_Stat.txt", q_field, delimiter="\n")\r
+ np.savetxt("EulerComplet_" + str(dim) + "_" + scheme + meshName + "_rhoE_Stat.txt", rho_E_field, delimiter="\n")\r
+ np.savetxt("EulerComplet_" + str(dim) + "_" + scheme + meshName + "_p_Stat.txt", p_field, delimiter="\n")\r
+ else:\r
+ print("Maximum time Tmax= ", tmax, " reached")\r
+\r
+ return p_field, v_field, e_field, rho_field, q_field, h_field\r
+\r
+\r
+def solve(a, b, nx, meshName, meshType, isImplicit, scheme, simulation_name, testColor,fig,p_L,v_L,rho_L,p_R,v_R,rho_R, pos_disc, tmax, lam_max, gamma, pinf):\r
+\r
+ print("Simulation of Euler equations in dimension 1 on " + str(nx) + " cells")\r
+ print("Problem name ", simulation_name)\r
+ print("State Law Stiffened Gaz, gamma= " , gamma, " pinf= ", pinf)\r
+ print("Initial data : ", "Riemann data")\r
+ print("Boundary conditions : ", "Neumann")\r
+ print("Mesh name : ", meshName, ", ", nx, " cells")\r
+\r
+ # Problem data\r
+ ntmax = 100000\r
+ output_freq = 1000\r
+ cfl=0.95\r
+ \r
+ test_desc={}\r
+ test_desc["Initial_data"]=simulation_name\r
+ test_desc["Boundary_conditions"]="Neumann"\r
+ test_desc["Global_name"]="FV simulation of "+simulation_name\r
+ test_desc["Global_comment"]="Regular 1D mesh"\r
+ test_desc["PDE_model"]="Euler equations"\r
+ test_desc["PDE_is_stationary"]=False\r
+ test_desc["PDE_search_for_stationary_solution"]=False\r
+ test_desc["Numerical_method_name"]=scheme\r
+ test_desc["Numerical_method_space_discretization"]="Finite volumes"\r
+ test_desc["Numerical_method_time_discretization"]="Explicit"\r
+ test_desc["Mesh_is_unstructured"]=False\r
+ test_desc["Mesh_type"]=meshType\r
+ test_desc["Test_color"]=testColor\r
+ test_desc["Geometry"]="Segment"\r
+ test_desc["Part_of_mesh_convergence_analysis"]=True\r
+ test_desc["Space_dimension"]=1\r
+ test_desc["Mesh_dimension"]=1\r
+ test_desc["Mesh_number_of_elements"]=nx\r
+ test_desc["Mesh_cell_type"]='SEG'\r
+ test_desc["Mesh_max_number_of_neighbours"]=2\r
+\r
+ start = time.time()\r
+ p_field, v_field, e_field, rho_field, q_field, h_field = EulerSystem(ntmax, tmax, cfl, a, b, nx, output_freq, meshName, isImplicit,scheme,fig,p_L,v_L,rho_L,p_R,v_R,rho_R, pos_disc,lam_max, gamma, pinf)\r
+ end = time.time()\r
+\r
+ test_desc["Computational_time_taken_by_run"]=end-start\r
+ with open('test_EulerEquations_1D_VF_'+simulation_name+"_"+str(nx)+ "Cells.json", 'w') as outfile: \r
+ json.dump(test_desc, outfile)\r
+\r
+ \r
+ return p_field, v_field, e_field, rho_field, q_field, h_field, end - start\r
+ \r
+\r
+if __name__ == """__main__""":\r
+ nbComp=3 # number of equations \r
+ a = 0.# domain is interval [a,b]\r
+ b = 4.2# domain is interval [a,b]\r
+ nx = 10# number of cells\r
+ dx = (b - a) / nx # space step\r
+ x = [a + 0.5 * dx + i * dx for i in range(nx)] # array of cell center (1D mesh)\r
+ rho0=p_to_rho_StiffenedGaz(p_0, T_0)\r
+ rhoE0=T_to_rhoE_StiffenedGaz(T_0, rho0, rho0*v_0)\r
+\r
+\r
+#### initial condition (T in K, v in m/s, p in Pa)\r
+ p_initial = np.array([ p_outlet for xi in x])\r
+ v_initial = np.array([ v_inlet for xi in x])\r
+ T_initial = np.array([ T_inlet for xi in x])\r
+ \r
+ rho_field = p_to_rho_StiffenedGaz(p_initial, T_initial)\r
+ q_field = rho_field * v_initial\r
+ rho_E_field = rho_field * T_to_E_StiffenedGaz(p_initial, T_initial, v_initial)\r
+\r
+ U = cdmath.Vector(nx * (nbComp))#Inutile à terme mais nécessaire pour le moment\r
+\r
+ for k in range(nx):\r
+ U[k * nbComp + 0] = rho_field[k]\r
+ U[k * nbComp + 1] = q_field[k]\r
+ U[k * nbComp + 2] = rho_E_field[k]\r
+\r
+ print("\n ############# Testing Roe scheme #############\n")\r
+ scheme='Roe'\r
+\r
+ print("\n Testing function solve (Implicit Roe scheme) \n")\r
+ isImplicit=True\r
+ cfl = 1000.\r
+ solve(a, b, nx, "RegularSquares", "", cfl, isImplicit, scheme, gamma, pinf) #Implicit Roe simulation\r
+ \r
+ print("\n Testing function solve (Explicit Roe scheme) \n")\r
+ isImplicit=False\r
+ cfl = 0.5\r
+ solve(a, b, nx, "RegularSquares", "", cfl, isImplicit, scheme, gamma, pinf) #Explicit Roe simulation\r
+\r
+ print("\n ############# Testing MAC scheme #############\n")\r
+ scheme='Mac'\r
+ isImplicit=True #Scheme must be implicit (explicit version is unstable)\r
+\r
+ print("\n Testing function solve (Implicit MAC scheme) \n")\r
+ isImplicit=True\r
+ cfl = 1000.\r
+ solve(a, b, nx, "RegularSquares", "", cfl, isImplicit, scheme, gamma, pinf) #Implicit MAC simulation\r
+ \r
