From c77ac8b1d7b223e9cee5efc23d1862e22da30c9f Mon Sep 17 00:00:00 2001 From: michael Date: Sun, 10 Oct 2021 00:00:57 +0200 Subject: [PATCH] Added validation tests for the Full Euler system --- .../1DToroBenchmark/CMakeLists.txt | 36 + .../test_validation1DEulerEquationsToro.py | 284 ++++++++ .../1DEulerEquations/EulerEquations1D.py | 626 ++++++++++++++++++ 3 files changed, 946 insertions(+) create mode 100644 CDMATH/tests/validation/EulerEquations/1DToroBenchmark/CMakeLists.txt create mode 100755 CDMATH/tests/validation/EulerEquations/1DToroBenchmark/test_validation1DEulerEquationsToro.py create mode 100644 CDMATH/tests/validation/scripts/1DEulerEquations/EulerEquations1D.py diff --git a/CDMATH/tests/validation/EulerEquations/1DToroBenchmark/CMakeLists.txt b/CDMATH/tests/validation/EulerEquations/1DToroBenchmark/CMakeLists.txt new file mode 100644 index 0000000..ff097a2 --- /dev/null +++ b/CDMATH/tests/validation/EulerEquations/1DToroBenchmark/CMakeLists.txt @@ -0,0 +1,36 @@ + +SET(SCRIPT + ../../scripts/1DEulerEquations/EulerEquations1D.py + ../../scripts/ReferenceSolutions/exact_rs_stiffenedgas.py + ./test_validation1DEulerEquationsToro.py + ) + +file(COPY ${SCRIPT} DESTINATION ${CMAKE_CURRENT_BINARY_DIR}) +install(FILES ${SCRIPT} DESTINATION share/validation/test_validation1DEulerEquationsToro) + +if (CDMATH_WITH_PYTHON ) + + SET(SCHEME Roe )#Roe numerical method + + SET(ISIMPLICIT 0 )#Explicit scheme + + ADD_TEST(validation1DEulerEquationsToro_Roe_Explicit ${PYTHON_EXECUTABLE} ${CMAKE_CURRENT_SOURCE_DIR}/test_validation1DEulerEquationsToro.py ${SCHEME} ${ISIMPLICIT} ) + + SET(ISIMPLICIT 1 )#Implicit scheme + + ADD_TEST(validation1DEulerEquationsToro_Roe_Implicit ${PYTHON_EXECUTABLE} ${CMAKE_CURRENT_SOURCE_DIR}/test_validation1DEulerEquationsToro.py ${SCHEME} ${ISIMPLICIT} ) + + SET(SCHEME Stag )#Conservative staggered numerical method + + SET(ISIMPLICIT 1 )#Implicit scheme + + ADD_TEST(validation1DEulerEquationsToro_Stag_Implicit ${PYTHON_EXECUTABLE} ${CMAKE_CURRENT_SOURCE_DIR}/test_validation1DEulerEquationsToro.py ${SCHEME} ${ISIMPLICIT} ) + + +endif (CDMATH_WITH_PYTHON ) + +install( DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}/ DESTINATION share/validation/test_validation1DEulerEquationsToro + FILES_MATCHING PATTERN "*.py" + PATTERN "*.png" +) + diff --git a/CDMATH/tests/validation/EulerEquations/1DToroBenchmark/test_validation1DEulerEquationsToro.py b/CDMATH/tests/validation/EulerEquations/1DToroBenchmark/test_validation1DEulerEquationsToro.py new file mode 100755 index 0000000..b49c902 --- /dev/null +++ b/CDMATH/tests/validation/EulerEquations/1DToroBenchmark/test_validation1DEulerEquationsToro.py @@ -0,0 +1,284 @@ +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) + diff --git a/CDMATH/tests/validation/scripts/1DEulerEquations/EulerEquations1D.py b/CDMATH/tests/validation/scripts/1DEulerEquations/EulerEquations1D.py new file mode 100644 index 0000000..2effdf8 --- /dev/null +++ b/CDMATH/tests/validation/scripts/1DEulerEquations/EulerEquations1D.py @@ -0,0 +1,626 @@ +#!/usr/bin/env python3 +# -*-coding:utf-8 -* + +""" +Created on Mon Aug 30 2021 +@author: Michael NDJINGA, Katia Ait Ameur, Coraline Mounier + +Euler system without heating source term (phi) in one dimension on regular domain [a,b] +Ghost cell (Neumann) boundary condition +Roe scheme or conservative MAC scheme +Regular square mesh + +State law Stiffened gaz : p = (gamma - 1) * rho * (e - q) - gamma * pinf +4 choices of parameters gamma and pinf are available : + - Lagrange interpolation (with q=0) + - Hermite interpolation with reference point at 575K (with q=0) + - Hermite interpolation with reference point at 590K (with q=0) + - Hermite interpolation with reference point at 617.94K (saturation at 155 bar) with q=0 + +Linearized enthalpy : h = h_sat + cp * (T - T_sat) +Values for cp and T_sat parameters are taken at the reference point chosen for the state law + +To do correct the computation of the time step : lambda_max (maximum eigenvalue) should be computed first) +""" + +import cdmath +import numpy as np +import matplotlib + +matplotlib.use("Agg") +import matplotlib.pyplot as plt +import sys, time, json +from math import sqrt +from numpy import sign + + +precision = 1.e-5 + +def p_to_e_StiffenedGaz(p, rho, gamma, pinf): + e_field = (p + gamma*pinf) / (gamma - 1.) / rho + return e_field + +def initial_conditions_Riemann_problem(a, b, nx,p_L,v_L,rho_L,p_R,v_R,rho_R,pos_disc, gamma, pinf): + 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) + + p_initial = np.array([ (xi < pos_disc) * p_L + (xi >= pos_disc) * p_R for xi in x]) + v_initial = np.array([ (xi < pos_disc) * v_L + (xi >= pos_disc) * v_R for xi in x]) + rho_initial = np.array([ (xi < pos_disc) * rho_L + (xi >= pos_disc) * rho_R for xi in x]) + + e_initial = p_to_e_StiffenedGaz(p_initial, rho_initial, gamma, pinf) + q_initial = rho_initial * v_initial + rho_E_initial = rho_initial*e_initial + 0.5 * q_initial**2/rho_initial + + return rho_initial, q_initial, rho_E_initial, p_initial, v_initial, e_initial + +def p_to_rho_StiffenedGaz(p_field, T_field, gamma, pinf): + rho_field = (p_field + pinf) * gamma / (gamma - 1) * 1. / (h_sat + cp * (T_field - T_sat)) + return rho_field + +def T_to_rhoE_StiffenedGaz(T_field, rho_field, q_field, gamma, pinf): + rho_E_field = 1. / 2. * (q_field) ** 2 / rho_field + pinf + rho_field / gamma * (h_sat + cp * (T_field- T_sat)) + return rho_E_field + +def rho_to_p_StiffenedGaz(rho_field, q_field, rho_E_field, gamma, pinf): + p_field = (gamma - 1) * (rho_E_field - 1. / 2 * q_field ** 2 / rho_field) - gamma * pinf + return p_field + +def T_to_E_StiffenedGaz(p_field, T_field, v_field, gamma, pinf): + rho_field = p_to_rho_StiffenedGaz(p_field, T_field) + E_field = (p_field + gamma * pinf) / ((gamma-1) * rho_field) + 0.5 * v_field **2 + return E_field + +def dp_drho_e_const_StiffenedGaz( e , gamma, pinf): + return (gamma-1)*(e) + +def dp_de_rho_const_StiffenedGaz( rho, gamma, pinf ): + return (gamma-1)*rho + +def sound_speed_StiffenedGaz( h, gamma, pinf ): + return np.sqrt((gamma-1)* h) + +def rho_h_to_p_StiffenedGaz( rho, h, gamma, pinf ): + return (gamma - 1) * rho * h / gamma - pinf + +def rhoE_to_e_StiffenedGaz(rho_field, q_field, rho_E_field): + e_field = rho_E_field / rho_field - 1. / 2. * (q_field / rho_field)**2 + return e_field + +def MatRoe_StiffenedGaz( rho_l, q_l, rho_E_l, rho_r, q_r, rho_E_r, gamma, pinf): + RoeMat = cdmath.Matrix(3, 3) + + u_l = q_l / rho_l + u_r = q_r / rho_r + p_l = rho_to_p_StiffenedGaz(rho_l, q_l, rho_E_l, gamma, pinf) + p_r = rho_to_p_StiffenedGaz(rho_r, q_r, rho_E_r, gamma, pinf) + H_l = rho_E_l / rho_l + p_l / rho_l + H_r = rho_E_r / rho_r + p_r / rho_r + + # Roe averages + rho = np.sqrt(rho_l * rho_r) + u = (u_l * sqrt(rho_l) + u_r * sqrt(rho_r)) / (sqrt(rho_l) + sqrt(rho_r)) + H = (H_l * sqrt(rho_l) + H_r * sqrt(rho_r)) / (sqrt(rho_l) + sqrt(rho_r)) + + p = rho_h_to_p_StiffenedGaz( rho, H - u**2/2., gamma, pinf ) + e = H - p / rho - 1./2 * u**2 + dp_drho = dp_drho_e_const_StiffenedGaz( e, gamma, pinf ) + dp_de = dp_de_rho_const_StiffenedGaz( rho, gamma, pinf ) + + RoeMat[0, 0] = 0 + RoeMat[0, 1] = 1 + RoeMat[0, 2] = 0 + RoeMat[1, 0] = dp_drho - u ** 2 + dp_de / rho * (u**2/2 - e) + RoeMat[1, 1] = 2 * u - u * dp_de / rho + RoeMat[1, 2] = dp_de / rho + RoeMat[2, 0] = -u * ( -dp_drho + H - dp_de / rho * (u**2/2 - e) ) + RoeMat[2, 1] = H - dp_de / rho * u ** 2 + RoeMat[2, 2] = (dp_de / rho +1) * u + + return(RoeMat) + + +def Dmac_StiffenedGaz( rho_l, q_l, rho_E_l, rho_r, q_r, rho_E_r, gamma, pinf): + Dmac = cdmath.Matrix(3, 3) + + u_l = q_l / rho_l + u_r = q_r / rho_r + p_l = rho_to_p_StiffenedGaz(rho_l, q_l, rho_E_l, gamma, pinf) + p_r = rho_to_p_StiffenedGaz(rho_r, q_r, rho_E_r, gamma, pinf) + H_l = rho_E_l / rho_l + p_l / rho_l + H_r = rho_E_r / rho_r + p_r / rho_r + + # Roe averages + rho = np.sqrt(rho_l * rho_r) + u = (u_l * sqrt(rho_l) + u_r * sqrt(rho_r)) / (sqrt(rho_l) + sqrt(rho_r)) + H = (H_l * sqrt(rho_l) + H_r * sqrt(rho_r)) / (sqrt(rho_l) + sqrt(rho_r)) + + p = rho_h_to_p_StiffenedGaz( rho, H - u**2/2., gamma, pinf ) + e = H - p / rho - 1./2 * u**2 + dp_drho = dp_drho_e_const_StiffenedGaz( e, gamma, pinf ) + dp_de = dp_de_rho_const_StiffenedGaz( rho, gamma, pinf ) + + #Third choice for Dstag + Dmac[0, 0] = 0 + Dmac[0, 1] = 1 + Dmac[0, 2] = 0 + Dmac[1, 0] = -dp_drho - u ** 2 - dp_de / rho * (u**2/2 - e) + Dmac[1, 1] = 2 * u + u * dp_de / rho + Dmac[1, 2] = -dp_de / rho + Dmac[2, 0] = -u * ( dp_drho + H + dp_de / rho * (u**2/2 - e) ) + Dmac[2, 1] = H + dp_de / rho * u ** 2 + Dmac[2, 2] = (-dp_de / rho +1) * u + + #return Dmac * sign(u) + + #Fifth choice for Dstag + Dmac[0, 0] = abs(u)-u + Dmac[0, 1] = 1 + Dmac[0, 2] = 0 + Dmac[1, 0] = -dp_drho - u ** 2 - dp_de / rho * (u**2/2 - e) + Dmac[1, 1] = abs(u) + u + u * dp_de / rho + Dmac[1, 2] = -dp_de / rho + Dmac[2, 0] = -u * ( dp_drho + H + u*(u-abs(u)) + dp_de / rho * (u**2/2 - e) ) + Dmac[2, 1] = H + dp_de / rho * u ** 2 + Dmac[2, 2] = -u*dp_de / rho + abs(u) + + return Dmac + +def Droe_StiffenedGaz( rho_l, q_l, rho_E_l, rho_r, q_r, rho_E_r, gamma, pinf): + Droe = cdmath.Matrix(3, 3) + + u_l = q_l / rho_l + u_r = q_r / rho_r + p_l = rho_to_p_StiffenedGaz(rho_l, q_l, rho_E_l, gamma, pinf) + p_r = rho_to_p_StiffenedGaz(rho_r, q_r, rho_E_r, gamma, pinf) + H_l = rho_E_l / rho_l + p_l / rho_l + H_r = rho_E_r / rho_r + p_r / rho_r + + # Roe averages + rho = np.sqrt(rho_l * rho_r) + u = (u_l * sqrt(rho_l) + u_r * sqrt(rho_r)) / (sqrt(rho_l) + sqrt(rho_r)) + H = (H_l * sqrt(rho_l) + H_r * sqrt(rho_r)) / (sqrt(rho_l) + sqrt(rho_r)) + + c = sound_speed_StiffenedGaz( H - u**2/2., gamma, pinf ) + + lamb = cdmath.Vector(3) + lamb[0] = u-c + lamb[1] = u + lamb[2] = u+c + + r = cdmath.Matrix(3, 3) + r[0,0] = 1. + r[1,0] = u-c + r[2,0] = H-u*c + r[0,1] = 1. + r[1,1] = u + r[2,1] = H-c**2/(gamma-1) + r[0,2] = 1. + r[1,2] = u+c + r[2,2] = H+u*c + + l = cdmath.Matrix(3, 3) + l[0,0] = (1./(2*c**2))*(0.5*(gamma-1)*u**2+u*c) + l[1,0] = (1./(2*c**2))*(-u*(gamma-1)-c) + l[2,0] = (1./(2*c**2))*(gamma-1) + l[0,1] = ((gamma-1)/c**2)*(H-u**2) + l[1,1] = ((gamma-1)/c**2)*u + l[2,1] = -((gamma-1)/c**2) + l[0,2] = (1./(2*c**2))*(0.5*(gamma-1)*u**2-u*c) + l[1,2] = (1./(2*c**2))*(c-u*(gamma-1)) + l[2,2] = (1./(2*c**2))*(gamma-1) + + M1 = cdmath.Matrix(3, 3) #abs(lamb[0])*r[:,0].tensProduct(l[:,0]) + M2 = cdmath.Matrix(3, 3) #abs(lamb[1])*r[:,1].tensProduct(l[:,1]) + M3 = cdmath.Matrix(3, 3) #abs(lamb[2])*r[:,2].tensProduct(l[:,2]) + for i in range(3): + for j in range(3): + M1[i,j] = abs(lamb[0])*r[i,0]*l[j,0] + M2[i,j] = abs(lamb[1])*r[i,1]*l[j,1] + M3[i,j] = abs(lamb[2])*r[i,2]*l[j,2] + + Droe = M1+M2+M3 + + return(Droe) + + +def jacobianMatricesm(coeff, rho_l, q_l, rho_E_l, rho_r, q_r, rho_E_r, scheme, gamma, pinf): + + if rho_l < 0 or rho_r < 0: + print("rho_l=", rho_l, " rho_r= ", rho_r) + raise ValueError("Negative density") + if rho_E_l < 0 or rho_E_r < 0: + print("rho_E_l=", rho_E_l, " rho_E_r= ", rho_E_r) + raise ValueError("Negative total energy") + + RoeMat = MatRoe_StiffenedGaz( rho_l, q_l, rho_E_l, rho_r, q_r, rho_E_r, gamma, pinf) + + if scheme == 'Stag': + Dmac = Dmac_StiffenedGaz( rho_l, q_l, rho_E_l, rho_r, q_r, rho_E_r, gamma, pinf) + return (RoeMat - Dmac) * coeff * 0.5 + + elif scheme == 'Roe': + Droe = Droe_StiffenedGaz( rho_l, q_l, rho_E_l, rho_r, q_r, rho_E_r, gamma, pinf) + return (RoeMat - Droe) * coeff * 0.5 + else: + raise ValueError("Wrong scheme given : ", scheme, ". Accepted scheme are Roe and Stag") + + +def jacobianMatricesp(coeff, rho_l, q_l, rho_E_l, rho_r, q_r, rho_E_r, scheme, gamma, pinf): + if rho_l < 0 or rho_r < 0: + print("rho_l=", rho_l, " rho_r= ", rho_r) + raise ValueError("Negative density") + if rho_E_l < 0 or rho_E_r < 0: + print("rho_E_l=", rho_E_l, " rho_E_r= ", rho_E_r) + raise ValueError("Negative total energy") + + RoeMat = MatRoe_StiffenedGaz( rho_l, q_l, rho_E_l, rho_r, q_r, rho_E_r, gamma, pinf) + + if scheme == 'Stag': + Dmac = Dmac_StiffenedGaz( rho_l, q_l, rho_E_l, rho_r, q_r, rho_E_r, gamma, pinf) + return (RoeMat - Dmac) * coeff * 0.5 + + elif scheme == 'Roe': + Droe = Droe_StiffenedGaz( rho_l, q_l, rho_E_l, rho_r, q_r, rho_E_r, gamma, pinf) + return (RoeMat - Droe) * coeff * 0.5 + + else: + raise ValueError("Wrong scheme given : ", scheme, ". Accepted scheme are Roe and Stag") + +def FillEdges(nx, j, Uk, nbComp, divMat, Rhs, Un, dt, dx, isImplicit, scheme, gamma, pinf): + dUi1 = cdmath.Vector(3) + dUi2 = cdmath.Vector(3) + temp1 = cdmath.Vector(3) + temp2 = cdmath.Vector(3) + + if (j == 0): + rho_l = Uk[j * nbComp + 0] + q_l = Uk[j * nbComp + 1] + rho_E_l = Uk[j * nbComp + 2] + rho_r = Uk[(j + 1) * nbComp + 0] + q_r = Uk[(j + 1) * nbComp + 1] + rho_E_r = Uk[(j + 1) * nbComp + 2] + + Am = jacobianMatricesm(dt / dx, rho_l, q_l, rho_E_l, rho_r, q_r, rho_E_r, scheme, gamma, pinf) + divMat.addValue(j * nbComp, (j + 1) * nbComp, Am) + divMat.addValue(j * nbComp, j * nbComp, Am * (-1.)) + + Ap = jacobianMatricesp(dt / dx, rho_l, q_l, rho_E_l, rho_l, q_l, rho_E_l, scheme, gamma, pinf) + divMat.addValue(j * nbComp, j * nbComp, Ap) + + if(isImplicit): + dUi1[0] = Uk[(j + 1) * nbComp + 0] - Uk[j * nbComp + 0] + dUi1[1] = Uk[(j + 1) * nbComp + 1] - Uk[j * nbComp + 1] + dUi1[2] = Uk[(j + 1) * nbComp + 2] - Uk[j * nbComp + 2] + temp1 = Am * dUi1 + else: + dUi2[0] = rho_l + dUi2[1] = q_l + dUi2[2] = rho_E_l + temp2 = Ap * dUi2 + + elif (j == nx - 1): + rho_l = Uk[(j - 1) * nbComp + 0] + q_l = Uk[(j - 1) * nbComp + 1] + rho_E_l = Uk[(j - 1) * nbComp + 2] + rho_r = Uk[j * nbComp + 0] + q_r = Uk[j * nbComp + 1] + rho_E_r = Uk[j * nbComp + 2] + + Ap = jacobianMatricesp(dt / dx, rho_l, q_l, rho_E_l, rho_r, q_r, rho_E_r, scheme, gamma, pinf) + divMat.addValue(j * nbComp, j * nbComp, Ap) + divMat.addValue(j * nbComp, (j - 1) * nbComp, Ap * (-1.)) + + Am = jacobianMatricesm(dt / dx, rho_r, q_r, rho_E_r, rho_r, q_r, rho_E_r, scheme, gamma, pinf) + divMat.addValue(j * nbComp, j * nbComp, Am * (-1.)) + + if(isImplicit): + dUi2[0] = Uk[(j - 1) * nbComp + 0] - Uk[j * nbComp + 0] + dUi2[1] = Uk[(j - 1) * nbComp + 1] - Uk[j * nbComp + 1] + dUi2[2] = Uk[(j - 1) * nbComp + 2] - Uk[j * nbComp + 2] + temp2 = Ap * dUi2 + else: + dUi1[0] = rho_r + dUi1[1] = q_r + dUi1[2] = rho_E_r + temp1 = Am * dUi1 + + if(isImplicit):#implicit scheme, contribution from the Newton scheme residual + 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]) + Rhs[j * nbComp + 2] = -temp1[2] + temp2[2] - (Uk[j * nbComp + 2] - Un[j * nbComp + 2]) + else:#explicit scheme, contribution from the boundary data the right hand side + Rhs[j * nbComp + 0] = -temp1[0] + temp2[0] + Rhs[j * nbComp + 1] = -temp1[1] + temp2[1] + Rhs[j * nbComp + 2] = -temp1[2] + temp2[2] + +def FillInnerCell(j, Uk, nbComp, divMat, Rhs, Un, dt, dx, isImplicit, scheme, gamma, pinf): + + rho_l = Uk[(j - 1) * nbComp + 0] + q_l = Uk[(j - 1) * nbComp + 1] + rho_E_l = Uk[(j - 1) * nbComp + 2] + rho_r = Uk[j * nbComp + 0] + q_r = Uk[j * nbComp + 1] + rho_E_r = Uk[j * nbComp + 2] + Ap = jacobianMatricesp(dt / dx, rho_l, q_l, rho_E_l, rho_r, q_r, rho_E_r, scheme, gamma, pinf) + + rho_l = Uk[j * nbComp + 0] + q_l = Uk[j * nbComp + 1] + rho_E_l = Uk[j * nbComp + 2] + rho_r = Uk[(j + 1) * nbComp + 0] + q_r = Uk[(j + 1) * nbComp + 1] + rho_E_r = Uk[(j + 1) * nbComp + 2] + Am = jacobianMatricesm(dt / dx, rho_l, q_l, rho_E_l, rho_r, q_r, rho_E_r, scheme, gamma, pinf) + + 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.)) + + if(isImplicit): + dUi1 = cdmath.Vector(3) + dUi2 = cdmath.Vector(3) + + dUi1[0] = Uk[(j + 1) * nbComp + 0] - Uk[j * nbComp + 0] + dUi1[1] = Uk[(j + 1) * nbComp + 1] - Uk[j * nbComp + 1] + dUi1[2] = Uk[(j + 1) * nbComp + 2] - Uk[j * nbComp + 2] + + dUi2[0] = Uk[(j - 1) * nbComp + 0] - Uk[j * nbComp + 0] + dUi2[1] = Uk[(j - 1) * nbComp + 1] - Uk[j * nbComp + 1] + dUi2[2] = Uk[(j - 1) * nbComp + 2] - Uk[j * nbComp + 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]) + Rhs[j * nbComp + 2] = -temp1[2] + temp2[2] - (Uk[j * nbComp + 2] - Un[j * nbComp + 2]) + else: + Rhs[j * nbComp + 0] = 0 + Rhs[j * nbComp + 1] = 0 + Rhs[j * nbComp + 2] = 0 + + +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): + dim = 1 + nbComp = 3 + dt = 0. + time = 0. + it = 0 + isStationary = False + dx = (b - a) / nbCells + dt = cfl * dx / lam_max + nbVoisinsMax = 2 + + e_L=p_to_e_StiffenedGaz(p_L, rho_L, gamma, pinf) + + # iteration vectors + Un = cdmath.Vector(nbCells * (nbComp)) + dUn = cdmath.Vector(nbCells * (nbComp)) + dUk = cdmath.Vector(nbCells * (nbComp)) + Rhs = cdmath.Vector(nbCells * (nbComp)) + + # Initial conditions + print("Construction of the initial condition …") + + 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) + h_field = (rho_E_field + p_field) / rho_field - 0.5 * (q_field / rho_field) **2 + p_field = p_field * 10 ** (-5) + + + for k in range(nbCells): + Un[k * nbComp + 0] = rho_field[k] + Un[k * nbComp + 1] = q_field[k] + Un[k * nbComp + 2] = rho_E_field[k] + + divMat = cdmath.SparseMatrixPetsc(nbCells * nbComp, nbCells * nbComp, (nbVoisinsMax + 1) * nbComp) + + iterGMRESMax = 50 + newton_max = 100 + + print("Starting computation of the non linear Euler non isentropic system with scheme : ", scheme) + # 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): + # STARTING NEWTON LOOP + divMat.zeroEntries() #sets the matrix coefficients to zero + for j in range(nbCells): + + # traitements des bords + if (j == 0): + FillEdges(nbCells,j, Uk, nbComp, divMat, Rhs, Un, dt, dx, isImplicit, scheme, gamma, pinf) + elif (j == nbCells - 1): + FillEdges(nbCells,j, Uk, nbComp, divMat, Rhs, Un, dt, dx, isImplicit, scheme, gamma, pinf) + + # traitement des cellules internes + else: + FillInnerCell( j, Uk, nbComp, divMat, Rhs, Un, dt, dx, isImplicit, scheme, gamma, pinf) + + if(isImplicit): + #solving the linear system divMat * dUk = Rhs + divMat.diagonalShift(1.) + LS = cdmath.LinearSolver(divMat, Rhs, iterGMRESMax, precision, "GMRES", "LU") + dUk = LS.solve() + vector_residu = dUk.maxVector(nbComp) + residu = max(abs(vector_residu[0])/rho_L, abs(vector_residu[1])/(rho_L*v_L), abs(vector_residu[2])/rho_L*e_L ) + else: + dUk=Rhs - divMat*Un + residu = 0.#Convergence schéma Newton + + if (isImplicit and (it == 1 or it % output_freq == 0 or it >= ntmax or isStationary or time >= tmax)): + print("Residu Newton at iteration ",k, " : ", residu) + print("Linear system converged in ", LS.getNumberOfIter(), " GMRES iterations") + + #updates for Newton loop + Uk += dUk + k = k + 1 + if (isImplicit and not LS.getStatus()): + print("Linear system did not converge ", LS.getNumberOfIter(), " GMRES iterations") + raise ValueError("No convergence of the linear system") + + if k == newton_max: + raise ValueError("No convergence of Newton algorithm for " + scheme + " scheme") + + #updating fields + Un = Uk.deepCopy() + dUn -= Un + + #Testing stationarity + residu_stat = dUn.maxVector(nbComp)#On prend le max de chaque composante + if (it % output_freq == 0 ): + 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 )) + + if ( it>1 and abs(residu_stat[0])/rho_L= ntmax or isStationary or time >= tmax): + print("-- Time step : " + str(it) + ", Time: " + str(time) + ", dt: " + str(dt)) + + print("-- Iter: " + str(it) + ", Time: " + str(time) + ", dt: " + str(dt)+"\n") + + if(isImplicit): + test_desc["Linear_solver_algorithm"]=LS.getNameOfMethod() + test_desc["Linear_solver_preconditioner"]=LS.getNameOfPc() + test_desc["Linear_solver_precision"]=LS.getTolerance() + test_desc["Linear_solver_maximum_iterations"]=LS.getNumberMaxOfIter() + test_desc["Linear_system_max_actual_iterations_number"]=LS.getNumberOfIter() + test_desc["Linear_system_max_actual_error"]=LS.getResidu() + test_desc["Linear_system_max_actual_condition number"]=LS.getConditionNumber() + + if (it >= ntmax): + print("Maximum number of time steps ntmax= ", ntmax, " reached") + + elif (isStationary): + print("Stationary regime reached at time step ", it, ", t= ", time) + print("------------------------------------------------------------------------------------") + np.savetxt("EulerComplet_" + str(dim) + "_" + scheme + meshName + "_rho_Stat.txt", rho_field, delimiter="\n") + np.savetxt("EulerComplet_" + str(dim) + "_" + scheme + meshName + "_q_Stat.txt", q_field, delimiter="\n") + np.savetxt("EulerComplet_" + str(dim) + "_" + scheme + meshName + "_rhoE_Stat.txt", rho_E_field, delimiter="\n") + np.savetxt("EulerComplet_" + str(dim) + "_" + scheme + meshName + "_p_Stat.txt", p_field, delimiter="\n") + else: + print("Maximum time Tmax= ", tmax, " reached") + + return p_field, v_field, e_field, rho_field, q_field, h_field + + +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): + + print("Simulation of Euler equations in dimension 1 on " + str(nx) + " cells") + print("Problem name ", simulation_name) + print("State Law Stiffened Gaz, gamma= " , gamma, " pinf= ", pinf) + print("Initial data : ", "Riemann data") + print("Boundary conditions : ", "Neumann") + print("Mesh name : ", meshName, ", ", nx, " cells") + + # Problem data + ntmax = 100000 + output_freq = 1000 + cfl=0.95 + + test_desc={} + test_desc["Initial_data"]=simulation_name + test_desc["Boundary_conditions"]="Neumann" + test_desc["Global_name"]="FV simulation of "+simulation_name + test_desc["Global_comment"]="Regular 1D mesh" + test_desc["PDE_model"]="Euler equations" + test_desc["PDE_is_stationary"]=False + test_desc["PDE_search_for_stationary_solution"]=False + test_desc["Numerical_method_name"]=scheme + test_desc["Numerical_method_space_discretization"]="Finite volumes" + test_desc["Numerical_method_time_discretization"]="Explicit" + test_desc["Mesh_is_unstructured"]=False + test_desc["Mesh_type"]=meshType + test_desc["Test_color"]=testColor + test_desc["Geometry"]="Segment" + test_desc["Part_of_mesh_convergence_analysis"]=True + test_desc["Space_dimension"]=1 + test_desc["Mesh_dimension"]=1 + test_desc["Mesh_number_of_elements"]=nx + test_desc["Mesh_cell_type"]='SEG' + test_desc["Mesh_max_number_of_neighbours"]=2 + + start = time.time() + 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) + end = time.time() + + test_desc["Computational_time_taken_by_run"]=end-start + with open('test_EulerEquations_1D_VF_'+simulation_name+"_"+str(nx)+ "Cells.json", 'w') as outfile: + json.dump(test_desc, outfile) + + + return p_field, v_field, e_field, rho_field, q_field, h_field, end - start + + +if __name__ == """__main__""": + nbComp=3 # number of equations + a = 0.# domain is interval [a,b] + b = 4.2# domain is interval [a,b] + nx = 10# number of cells + dx = (b - a) / nx # space step + x = [a + 0.5 * dx + i * dx for i in range(nx)] # array of cell center (1D mesh) + rho0=p_to_rho_StiffenedGaz(p_0, T_0) + rhoE0=T_to_rhoE_StiffenedGaz(T_0, rho0, rho0*v_0) + + +#### initial condition (T in K, v in m/s, p in Pa) + p_initial = np.array([ p_outlet for xi in x]) + v_initial = np.array([ v_inlet for xi in x]) + T_initial = np.array([ T_inlet for xi in x]) + + rho_field = p_to_rho_StiffenedGaz(p_initial, T_initial) + q_field = rho_field * v_initial + rho_E_field = rho_field * T_to_E_StiffenedGaz(p_initial, T_initial, v_initial) + + U = cdmath.Vector(nx * (nbComp))#Inutile à terme mais nécessaire pour le moment + + for k in range(nx): + U[k * nbComp + 0] = rho_field[k] + U[k * nbComp + 1] = q_field[k] + U[k * nbComp + 2] = rho_E_field[k] + + print("\n ############# Testing Roe scheme #############\n") + scheme='Roe' + + print("\n Testing function solve (Implicit Roe scheme) \n") + isImplicit=True + cfl = 1000. + solve(a, b, nx, "RegularSquares", "", cfl, isImplicit, scheme, gamma, pinf) #Implicit Roe simulation + + print("\n Testing function solve (Explicit Roe scheme) \n") + isImplicit=False + cfl = 0.5 + solve(a, b, nx, "RegularSquares", "", cfl, isImplicit, scheme, gamma, pinf) #Explicit Roe simulation + + print("\n ############# Testing MAC scheme #############\n") + scheme='Mac' + isImplicit=True #Scheme must be implicit (explicit version is unstable) + + print("\n Testing function solve (Implicit MAC scheme) \n") + isImplicit=True + cfl = 1000. + solve(a, b, nx, "RegularSquares", "", cfl, isImplicit, scheme, gamma, pinf) #Implicit MAC simulation + -- 2.39.2