From: michael Date: Sat, 9 Oct 2021 21:59:25 +0000 (+0200) Subject: Added examples for the Full Euler system X-Git-Tag: V9_8_0~88 X-Git-Url: http://git.salome-platform.org/gitweb/?a=commitdiff_plain;h=386965127ce8969481a837335e9fe38e21351ab7;p=tools%2Fsolverlab.git Added examples for the Full Euler system --- diff --git a/CDMATH/tests/examples/EulerEquations/FullEulerSystem/CMakeLists.txt b/CDMATH/tests/examples/EulerEquations/FullEulerSystem/CMakeLists.txt index cb95639..7db8212 100755 --- a/CDMATH/tests/examples/EulerEquations/FullEulerSystem/CMakeLists.txt +++ b/CDMATH/tests/examples/EulerEquations/FullEulerSystem/CMakeLists.txt @@ -1,12 +1,14 @@ file(GLOB FullEulerEquations_EXAMPLES_TO_INSTALL -EulerSystem1D/EulerSystem1DUpwind EulerSystem1D/EulerSystem1DUpwindEntrCorr EulerSystem1D/EulerSystem1DConservativeStaggered EulerSystem_Shock/EulerSystemStaggered EulerSystem_Shock/EulerSystemUpwind +EulerSystem1D_HeatedChannel_Roe EulerSystem1D_HeatedChannel_MAC EulerSystem1D_RiemannProblem_Roe EulerSystem1D_RiemannProblem_MAC ) install(DIRECTORY ${FullEulerEquations_EXAMPLES_TO_INSTALL} DESTINATION share/examples/EulerEquations/FullEulerEquations) IF (CDMATH_WITH_PYTHON AND CDMATH_WITH_PETSC AND CDMATH_WITH_POSTPRO) - ADD_SUBDIRECTORY(EulerSystem1D_RiemannProblem) - ADD_SUBDIRECTORY(EulerSystem1D_HeatedChannel) + ADD_SUBDIRECTORY(EulerSystem1D_RiemannProblem_Roe) + ADD_SUBDIRECTORY(EulerSystem1D_RiemannProblem_MAC) + ADD_SUBDIRECTORY(EulerSystem1D_HeatedChannel_Roe) + ADD_SUBDIRECTORY(EulerSystem1D_HeatedChannel_MAC) ENDIF (CDMATH_WITH_PYTHON AND CDMATH_WITH_PETSC AND CDMATH_WITH_POSTPRO) diff --git a/CDMATH/tests/examples/EulerEquations/FullEulerSystem/EulerSystem1D_HeatedChannel/CMakeLists.txt b/CDMATH/tests/examples/EulerEquations/FullEulerSystem/EulerSystem1D_HeatedChannel/CMakeLists.txt deleted file mode 100755 index 8415053..0000000 --- a/CDMATH/tests/examples/EulerEquations/FullEulerSystem/EulerSystem1D_HeatedChannel/CMakeLists.txt +++ /dev/null @@ -1,8 +0,0 @@ - -if (CDMATH_WITH_PYTHON AND CDMATH_WITH_PETSC AND CDMATH_WITH_POSTPRO) - - ADD_TEST(ExampleFullEulerSystem_1DHeatedChannel_Roe ${PYTHON_EXECUTABLE} ${CMAKE_CURRENT_SOURCE_DIR}/Euler_complet_HeatedChanel.py ) - -endif (CDMATH_WITH_PYTHON AND CDMATH_WITH_PETSC AND CDMATH_WITH_POSTPRO) - - diff --git a/CDMATH/tests/examples/EulerEquations/FullEulerSystem/EulerSystem1D_HeatedChannel/Euler_complet_HeatedChanel.py b/CDMATH/tests/examples/EulerEquations/FullEulerSystem/EulerSystem1D_HeatedChannel/Euler_complet_HeatedChanel.py deleted file mode 100644 index e63c38c..0000000 --- a/CDMATH/tests/examples/EulerEquations/FullEulerSystem/EulerSystem1D_HeatedChannel/Euler_complet_HeatedChanel.py +++ /dev/null @@ -1,981 +0,0 @@ -#!/usr/bin/env python3 -# -*-coding:utf-8 -* - -""" -Created on Mon Aug 30 2021 -@author: Michael NDJINGA, Katia Ait Ameur, Coraline Mounier - -Euler system with heating source term (phi) in one dimension on regular domain [a,b] -Riemann problemn with ghost cell boundary condition -Left : Inlet boundary condition (velocity and temperature imposed) -Right : Outlet boundary condition (pressure imposed) -Roe scheme -Regular square mesh - -State law Stiffened gaz : p = (gamma - 1) * rho * (e - q) - gamma * p0 -4 choices of parameters gamma and p0 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 matplotlib.animation as manimation -import sys -from math import sqrt, atan, pi -from numpy import sign - - -#### Initial and boundary condition (T in K, v in m/s, p in Pa) -T_inlet = 565. -v_inlet = 5. -p_outlet = 155.0 * 10**5 - -#initial parameters are determined from boundary conditions -p_0 = p_outlet #initial pressure -v_0 = v_inlet #initial velocity -T_0 = T_inlet #initial temperature -### Heating source term -phi=1.e8 - -## Numerical parameter -precision = 1e-6 - -#state law parameter : can be 'Lagrange', 'Hermite590K', 'Hermite617K', or 'FLICA' -state_law = "Hermite575K" - -def state_law_parameters(state_law): - #state law Stiffened Gaz : p = (gamma - 1) * rho * e - gamma * p0 - global gamma - global p0 - global q - global c0 - global cp - global h_sat - global T_sat - - if state_law == "Lagrange": - # reference values for Lagrange interpolation - p_ref = 155. * 10**5 #Reference pressure in a REP 900 nuclear power plant - p1 = 153. * 10**5 # value of pressure at inlet of a 900 MWe PWR vessel - rho_ref = 594.38 #density of water at saturation temperature of 617.94K and 155 bars - rho1 = 742.36 # value of density at inlet of a 900 MWe PWR vessel (T1 = 565K) - e_ref = 1603.8 * 10**3 #internal energy of water at saturation temperature of 617.94K and 155 bars - e1 = 1273.6 * 10**3 # value of internal energy at inlet of a 900 MWe PWR vessel - - gamma = (p1 - p_ref) / (rho1 * e1 - rho_ref *e_ref) + 1. - p0 = - 1. / gamma * ( - (gamma - 1) * rho_ref * e_ref + p_ref) - q=0. - c0 = sqrt((gamma - 1) * (e_ref + p_ref / rho_ref)) - - cp = 8950. - h_sat = 1.63 * 10 ** 6 # saturation enthalpy of water at 155 bars - T_sat = 617.94 - - elif state_law == "Hermite617K": - # reference values for Hermite interpolation - p_ref = 155. * 10**5 #Reference pressure in a REP 900 nuclear power plant - T_ref = 617.94 #Reference temperature for interpolation at 617.94K - rho_ref = 594.38 #density of water at saturation temperature of 617.94K and 155 bars - e_ref = 1603.8 * 10**3 #internal energy of water at saturation temperature of 617.94K and 155 bars - h_ref = e_ref + p_ref / rho_ref - c_ref = 621.43 #sound speed for water at 155 bars and 617.94K - - gamma = 1. + c_ref * c_ref / (e_ref + p_ref / rho_ref) # From the sound speed formula - p0 = 1. / gamma * ( (gamma - 1) * rho_ref * e_ref - p_ref) - q=0. - c0 = sqrt((gamma - 1) * (e_ref + p_ref / rho_ref)) - - cp = 8950. # value at 155 bar and 617.94K - h_sat = 1.63 * 10 ** 6 # saturation enthalpy of water at 155 bars - T_sat = 617.94 - - elif state_law == 'Hermite590K': - # reference values for Hermite interpolation - p_ref = 155. * 10**5 #Reference pressure in a REP 900 nuclear power plant - T_ref = 590. #Reference temperature for interpolation at 590K - rho_ref = 688.3 #density of water at 590K and 155 bars - e_ref = 1411.4 * 10**3 #internal energy of water at 590K and 155 bars - h_ref = e_ref + p_ref / rho_ref - c_ref = 866.29 #sound speed for water at 155 bars and 590K - - gamma = 1. + c_ref * c_ref / (e_ref + p_ref / rho_ref) # From the sound speed formula - p0 = 1. / gamma * ( (gamma - 1) * rho_ref * e_ref - p_ref) - q=0. - c0 = sqrt((gamma - 1) * (e_ref + p_ref / rho_ref)) - - cp = 5996.8 # value at 155 bar and 590K - h_sat = 1433.9 * 10 ** 3 # saturation enthalpy of water at 155 bars - T_sat = 590. - - elif state_law == 'Hermite575K': - # reference values for Hermite interpolation - p_ref = 155 * 10**5 #Reference pressure in a REP 900 nuclear power plant - T_ref = 575 #Reference temperature at inlet in a REP 900 nuclear power plant - #Remaining values determined using iapws python package - rho_ref = 722.66 #density of water at 575K and 155 bars - e_ref = 1326552.66 #internal energy of water at 575K and 155 bars - h_ref = e_ref + p_ref / rho_ref - c_ref = 959.28 #sound speed for water at 155 bars and 575K - - gamma = 1 + c_ref * c_ref / (e_ref + p_ref / rho_ref) # From the sound speed formula - p0 = 1 / gamma * ( (gamma - 1) * rho_ref * e_ref - p_ref) - q=0. - c0 = sqrt((gamma - 1) * (e_ref + p_ref / rho_ref))#This is actually c_ref - - cp = 5504.05 # value at 155 bar and 590K - h_sat = h_ref # saturation enthalpy of water at 155 bars - T_sat = T_ref - else: - raise ValueError("Incorrect value for parameter state_law") - -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) - - p_initial = np.array([ p_0 for xi in x]) - v_initial = np.array([ v_0 for xi in x]) - T_initial = np.array([ T_0 for xi in x]) - - rho_initial = p_to_rho_StiffenedGaz(p_initial, T_initial) - q_initial = rho_initial * v_initial - rho_E_initial = T_to_rhoE_StiffenedGaz(T_initial, rho_initial, q_initial) - - return rho_initial, q_initial, rho_E_initial, p_initial, v_initial, T_initial - -def p_to_rho_StiffenedGaz(p_field, T_field): - rho_field = (p_field + p0) * 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): - rho_E_field = 1. / 2. * (q_field) ** 2 / rho_field + p0 + rho_field / gamma * (h_sat + cp * (T_field- T_sat)) - return rho_E_field - -def rhoE_to_T_StiffenedGaz(rho_field, q_field, rho_E_field): - T_field = T_sat + 1 / cp * (gamma * (rho_E_field / rho_field - 1 / 2 * (q_field / rho_field) ** 2) - gamma * p0 / rho_field - h_sat) - return T_field - -def rho_to_p_StiffenedGaz(rho_field, q_field, rho_E_field): - p_field = (gamma - 1) * (rho_E_field - 1. / 2 * q_field ** 2 / rho_field) - gamma * p0 - return p_field - -def T_to_E_StiffenedGaz(p_field, T_field, v_field): - rho_field = p_to_rho_StiffenedGaz(p_field, T_field) - E_field = (p_field + gamma * p0) / ((gamma-1) * rho_field) + 0.5 * v_field **2 - return E_field - -def dp_drho_e_const_StiffenedGaz( e ): - return (gamma-1)*(e-q) - -def dp_de_rho_const_StiffenedGaz( rho ): - return (gamma-1)*rho - -def sound_speed_StiffenedGaz( h ): - return np.sqrt((gamma-1)*(h-q)) - -def rho_h_to_p_StiffenedGaz( rho, h ): - return (gamma - 1) * rho * ( h - q ) / gamma - p0 - -def MatRoe_StiffenedGaz( rho_l, q_l, rho_E_l, rho_r, q_r, rho_E_r): - 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) - p_r = rho_to_p_StiffenedGaz(rho_r, q_r, rho_E_r) - 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. ) - e = H - p / rho - 1./2 * u**2 - dp_drho = dp_drho_e_const_StiffenedGaz( e ) - dp_de = dp_de_rho_const_StiffenedGaz( rho ) - - 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 Droe_StiffenedGaz( rho_l, q_l, rho_E_l, rho_r, q_r, rho_E_r): - 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) - p_r = rho_to_p_StiffenedGaz(rho_r, q_r, rho_E_r) - 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. ) - - 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): - - 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) - - Droe = Droe_StiffenedGaz( rho_l, q_l, rho_E_l, rho_r, q_r, rho_E_r) - return (RoeMat - Droe) * coeff * 0.5 - - -def jacobianMatricesp(coeff, rho_l, q_l, rho_E_l, rho_r, q_r, rho_E_r): - 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) - Droe = Droe_StiffenedGaz( rho_l, q_l, rho_E_l, rho_r, q_r, rho_E_r) - - return (RoeMat + Droe) * coeff * 0.5 - - -def FillEdges(j, Uk, nbComp, divMat, Rhs, Un, dt, dx, isImplicit): - 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) - divMat.addValue(j * nbComp, (j + 1) * nbComp, Am) - divMat.addValue(j * nbComp, j * nbComp, Am * (-1.)) - - p_inlet = rho_to_p_StiffenedGaz(Uk[j * nbComp + 0], Uk[j * nbComp + 1], Uk[j * nbComp + 2])# We take p from inside the domain - rho_l=p_to_rho_StiffenedGaz(p_inlet, T_inlet) # rho is computed from the temperature BC and the inner pressure - #rho_l = Uk[j * nbComp + 0] # We take rho from inside the domain - q_l = rho_l * v_inlet # q is imposed by the boundary condition v_inlet - rho_E_l = T_to_rhoE_StiffenedGaz(T_inlet, rho_l, q_l) #rhoE is obtained using the two boundary conditions v_inlet and e_inlet - 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) - 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 - - dUi2[0] = rho_l - Uk[(j ) * nbComp + 0] - dUi2[1] = q_l - Uk[(j ) * nbComp + 1] - dUi2[2] = rho_E_l - Uk[(j ) * nbComp + 2] - temp2 = Ap * dUi2 - 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) - divMat.addValue(j * nbComp, j * nbComp, Ap) - divMat.addValue(j * nbComp, (j - 1) * nbComp, Ap * (-1.)) - - rho_l = Uk[j * nbComp + 0] - q_l = Uk[j * nbComp + 1] - rho_E_l = Uk[j * nbComp + 2] - rho_r = Uk[(j ) * nbComp + 0] # We take rho inside the domain - q_r = Uk[(j ) * nbComp + 1] # We take q from inside the domain - rho_E_r = (p_outlet+gamma*p0)/(gamma-1) + 0.5*q_r**2/rho_r # rhoE is obtained using the boundary condition p_outlet - - Am = jacobianMatricesm(dt / dx, rho_l, q_l, rho_E_l, rho_r, q_r, rho_E_r) - divMat.addValue(j * nbComp, j * nbComp, Am * (-1.)) - - if(isImplicit): - dUi1[0] = rho_r - Uk[j * nbComp + 0] - dUi1[1] = q_r - Uk[j * nbComp + 1] - dUi1[2] = rho_E_r - Uk[j * nbComp + 2] - temp1 = Am * dUi1 - - 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): - - 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) - - 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) - - 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 SetPicture(rho_field_Roe, q_field_Roe, h_field_Roe, p_field_Roe, v_field_Roe, T_field_Roe, dx): - max_initial_rho = max(rho_field_Roe) - min_initial_rho = min(rho_field_Roe) - max_initial_q = max(q_field_Roe) - min_initial_q = min(q_field_Roe) - min_initial_h = min(h_field_Roe) - max_initial_h = max(h_field_Roe) - max_initial_p = max(p_field_Roe) - min_initial_p = min(p_field_Roe) - max_initial_v = max(v_field_Roe) - min_initial_v = min(v_field_Roe) - max_initial_T = max(T_field_Roe) - min_initial_T = min(T_field_Roe) - - fig, ([axDensity, axMomentum, axh],[axPressure, axVitesse, axTemperature]) = plt.subplots(2, 3,sharex=True, figsize=(14,10)) - plt.gcf().subplots_adjust(wspace = 0.5) - - lineDensity_Roe, = axDensity.plot([a+0.5*dx + i*dx for i in range(nx)], rho_field_Roe, label='Roe') - axDensity.set(xlabel='x (m)', ylabel='Densité (kg/m3)') - axDensity.set_xlim(a,b) - axDensity.set_ylim(680, 800) - axDensity.legend() - - lineMomentum_Roe, = axMomentum.plot([a+0.5*dx + i*dx for i in range(nx)], q_field_Roe, label='Roe') - axMomentum.set(xlabel='x (m)', ylabel='Momentum (kg/(m2.s))') - axMomentum.set_xlim(a,b) - axMomentum.set_ylim(3500, 4000) - axMomentum.legend() - - lineh_Roe, = axh.plot([a+0.5*dx + i*dx for i in range(nx)], h_field_Roe, label='Roe') - axh.set(xlabel='x (m)', ylabel='h (J/m3)') - axh.set_xlim(a,b) - axh.set_ylim(1.2 * 10**6, 1.5*10**6) - axh.legend() - - linePressure_Roe, = axPressure.plot([a+0.5*dx + i*dx for i in range(nx)], p_field_Roe, label='Roe') - axPressure.set(xlabel='x (m)', ylabel='Pression (bar)') - axPressure.set_xlim(a,b) - axPressure.set_ylim(155, 155.5) - axPressure.ticklabel_format(axis='y', style='sci', scilimits=(0,0)) - axPressure.legend() - - lineVitesse_Roe, = axVitesse.plot([a+0.5*dx + i*dx for i in range(nx)], v_field_Roe, label='Roe') - axVitesse.set(xlabel='x (m)', ylabel='Vitesse (m/s)') - axVitesse.set_xlim(a,b) - axVitesse.set_ylim(v_0-1, v_0+1) - axVitesse.ticklabel_format(axis='y', style='sci', scilimits=(0,0)) - axVitesse.legend() - - lineTemperature_Roe, = axTemperature.plot([a+0.5*dx + i*dx for i in range(nx)], T_field_Roe, label='Roe') - axTemperature.set(xlabel='x (m)', ylabel='Température (K)') - axTemperature.set_xlim(a,b) - axTemperature.set_ylim(T_0-10, T_0+30) - axTemperature.ticklabel_format(axis='y', style='sci', scilimits=(0,0)) - axTemperature.legend() - - return(fig, lineDensity_Roe, lineMomentum_Roe, lineh_Roe, linePressure_Roe, lineVitesse_Roe, lineTemperature_Roe) - - -def EulerSystemRoe(ntmax, tmax, cfl, a, b, nbCells, output_freq, meshName, state_law, isImplicit): - state_law_parameters(state_law) - dim = 1 - nbComp = 3 - dt = 0. - time = 0. - it = 0 - isStationary = False - dx = (b - a) / nx - dt = cfl * dx / c0 - #dt = 5*10**(-6) - nbVoisinsMax = 2 - - # iteration vectors - Un_Roe = cdmath.Vector(nbCells * (nbComp)) - dUn_Roe = cdmath.Vector(nbCells * (nbComp)) - dUk_Roe = cdmath.Vector(nbCells * (nbComp)) - Rhs_Roe = cdmath.Vector(nbCells * (nbComp)) - - # Initial conditions - print("Construction of the initial condition …") - - rho_field_Roe, q_field_Roe, rho_E_field_Roe, p_field_Roe, v_field_Roe, T_field_Roe = initial_conditions_Riemann_problem(a, b, nx) - h_field_Roe = (rho_E_field_Roe + p_field_Roe) / rho_field_Roe - 0.5 * (q_field_Roe / rho_field_Roe) **2 - p_field_Roe = p_field_Roe * 10 ** (-5) - - - for k in range(nbCells): - Un_Roe[k * nbComp + 0] = rho_field_Roe[k] - Un_Roe[k * nbComp + 1] = q_field_Roe[k] - Un_Roe[k * nbComp + 2] = rho_E_field_Roe[k] - - divMat_Roe = cdmath.SparseMatrixPetsc(nbCells * nbComp, nbCells * nbComp, (nbVoisinsMax + 1) * nbComp) - - # Picture settings - fig, lineDensity_Roe, lineMomentum_Roe, lineRhoE_Roe, linePressure_Roe, lineVitesse_Roe, lineTemperature_Roe = SetPicture( rho_field_Roe, q_field_Roe, h_field_Roe, p_field_Roe, v_field_Roe, T_field_Roe, dx) - - plt.savefig("EulerComplet_HeatedChannel_" + str(dim) + "D_Roe" + meshName + "0" + ".png") - iterGMRESMax = 50 - newton_max = 100 - - print("Starting computation of the non linear Euler non isentropic system with Roe scheme …") - # STARTING TIME LOOP - while (it < ntmax and time <= tmax and not isStationary): - dUn_Roe = Un_Roe.deepCopy() - Uk_Roe = Un_Roe.deepCopy() - residu_Roe = 1. - - k_Roe = 0 - while (k_Roe < newton_max and residu_Roe > precision): - # STARTING NEWTON LOOP - divMat_Roe.zeroEntries() #sets the matrix coefficients to zero - for j in range(nbCells): - - # traitements des bords - if (j == 0): - FillEdges(j, Uk_Roe, nbComp, divMat_Roe, Rhs_Roe, Un_Roe, dt, dx, isImplicit) - elif (j == nbCells - 1): - FillEdges(j, Uk_Roe, nbComp, divMat_Roe, Rhs_Roe, Un_Roe, dt, dx, isImplicit) - - # traitement des cellules internes - else: - FillInnerCell(j, Uk_Roe, nbComp, divMat_Roe, Rhs_Roe, Un_Roe, dt, dx, isImplicit) - - Rhs_Roe[j * nbComp + 2] += phi*dt - - if(isImplicit): - #solving the linear system divMat * dUk = Rhs - divMat_Roe.diagonalShift(1.) - LS_Roe = cdmath.LinearSolver(divMat_Roe, Rhs_Roe, iterGMRESMax, precision, "GMRES", "LU") - dUk_Roe = LS_Roe.solve() - vector_residu_Roe = dUk_Roe.maxVector(nbComp) - residu_Roe = max(abs(vector_residu_Roe[0])/rho0, abs(vector_residu_Roe[1])/(rho0*v_0), abs(vector_residu_Roe[2])/rhoE0 ) - else: - dUk_Roe=Rhs_Roe - divMat_Roe*Un_Roe - residu_Roe = 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_Roe, " : ", residu_Roe) - print("Linear system converged in ", LS_Roe.getNumberOfIter(), " GMRES iterations") - - #updates for Newton loop - Uk_Roe += dUk_Roe - k_Roe = k_Roe + 1 - if (isImplicit and not LS_Roe.getStatus()): - print("Linear system did not converge ", LS_Roe.getNumberOfIter(), " GMRES iterations") - raise ValueError("No convergence of the linear system") - - if k_Roe == newton_max: - raise ValueError("No convergence of Newton Roe Scheme") - - #updating fields - Un_Roe = Uk_Roe.deepCopy() - dUn_Roe -= Un_Roe - - #Testing stationarity - residu_stat = dUn_Roe.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])/rho0, abs(residu_stat[1])/(rho0*v_0), abs(residu_stat[2])/rhoE0 )) - - if ( it>1 and abs(residu_stat[0])/rho0= ntmax or isStationary or time >= tmax): - - print("-- Time step : " + str(it) + ", Time: " + str(time) + ", dt: " + str(dt)) - - print("Temperature gain between inlet and outlet is ", T_field_Roe[nbCells-1]-T_field_Roe[0],"\n") - - plt.savefig("EulerComplet_HeatedChannel_" + str(dim) + "D_Roe" + meshName + "Implicit"+str(isImplicit)+ str(it) + '_time' + str(time) + ".png") - - print("-- Iter: " + str(it) + ", Time: " + str(time) + ", dt: " + str(dt)+"\n") - - if (it >= ntmax): - print("Maximum number of time steps ntmax= ", ntmax, " reached") - return - - elif (isStationary): - print("Stationary regime reached at time step ", it, ", t= ", time) - print("------------------------------------------------------------------------------------") - np.savetxt("EulerComplet_HeatedChannel_" + str(dim) + "D_Roe" + meshName + "_rho_Stat.txt", rho_field_Roe, delimiter="\n") - np.savetxt("EulerComplet_HeatedChannel_" + str(dim) + "D_Roe" + meshName + "_q_Stat.txt", q_field_Roe, delimiter="\n") - np.savetxt("EulerComplet_HeatedChannel_" + str(dim) + "D_Roe" + meshName + "_rhoE_Stat.txt", rho_E_field_Roe, delimiter="\n") - np.savetxt("EulerComplet_HeatedChannel_" + str(dim) + "D_Roe" + meshName + "_p_Stat.txt", p_field_Roe, delimiter="\n") - plt.savefig("EulerComplet_HeatedChannel_" + str(dim) + "D_Roe" + meshName + "Implicit"+str(isImplicit)+"_Stat.png") - return - else: - print("Maximum time Tmax= ", tmax, " reached") - return - - -def solve(a, b, nx, meshName, meshType, cfl, state_law, isImplicit): - print("Simulation of a heated channel in dimension 1 on " + str(nx) + " cells") - print("State Law Stiffened Gaz, " + state_law) - print("Initial data : ", "constant fields") - print("Boundary conditions : ", "Inlet (Left), Outlet (Right)") - print("Mesh name : ", meshName, ", ", nx, " cells") - # Problem data - tmax = 10. - ntmax = 100000 - output_freq = 1000 - EulerSystemRoe(ntmax, tmax, cfl, a, b, nx, output_freq, meshName, state_law, isImplicit) - return - -def FillMatrixFromEdges(j, Uk, nbComp, divMat, dt, dx): - - 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) - divMat.addValue(j * nbComp, (j + 1) * nbComp, Am) - divMat.addValue(j * nbComp, j * nbComp, Am * (-1.)) - - p_inlet = rho_to_p_StiffenedGaz(Uk[j * nbComp + 0], Uk[j * nbComp + 1], Uk[j * nbComp + 2])# We take p from inside the domain - rho_l=p_to_rho_StiffenedGaz(p_inlet, T_inlet) # rho is computed from the temperature BC and the inner pressure - #rho_l = Uk[j * nbComp + 0] # We take rho from inside the domain - q_l = rho_l * v_inlet # q is imposed by the boundary condition v_inlet - rho_E_l = T_to_rhoE_StiffenedGaz(T_inlet, rho_l, q_l) #rhoE is obtained using the two boundary conditions v_inlet and e_inlet - 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) - divMat.addValue(j * nbComp, j * nbComp, Ap) - - elif (j == nx - 1): - rho_l = Uk[j * nbComp + 0] - q_l = Uk[j * nbComp + 1] - rho_E_l = Uk[j * nbComp + 2] - rho_r = Uk[(j ) * nbComp + 0] # We take rho inside the domain - q_r = Uk[(j ) * nbComp + 1] # We take q from inside the domain - rho_E_r = (p_outlet+gamma*p0)/(gamma-1) + 0.5*q_r**2/rho_r # rhoE is obtained using the boundary condition p_outlet - - Am = jacobianMatricesm(dt / dx, rho_l, q_l, rho_E_l, rho_r, q_r, rho_E_r) - divMat.addValue(j * nbComp, j * nbComp, Am * (-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) - divMat.addValue(j * nbComp, j * nbComp, Ap) - divMat.addValue(j * nbComp, (j - 1) * nbComp, Ap * (-1.)) - - -def FillMatrixFromInnerCell(j, Uk, nbComp, divMat, dt, dx): - - 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) - - 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) - - 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.)) - -def FillRHSFromEdges(j, Uk, nbComp, Rhs, Un, dt, dx, isImplicit): - 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) - - p_inlet = rho_to_p_StiffenedGaz(Uk[j * nbComp + 0], Uk[j * nbComp + 1], Uk[j * nbComp + 2])# We take p from inside the domain - rho_l=p_to_rho_StiffenedGaz(p_inlet, T_inlet) # rho is computed from the temperature BC and the inner pressure - #rho_l = Uk[j * nbComp + 0] # We take rho from inside the domain - q_l = rho_l * v_inlet # q is imposed by the boundary condition v_inlet - rho_E_l = T_to_rhoE_StiffenedGaz(T_inlet, rho_l, q_l) #rhoE is obtained using the two boundary conditions v_inlet and e_inlet - 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) - - 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 - - dUi2[0] = rho_l - Uk[(j ) * nbComp + 0] - dUi2[1] = q_l - Uk[(j ) * nbComp + 1] - dUi2[2] = rho_E_l - Uk[(j ) * nbComp + 2] - temp2 = Ap * dUi2 - 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) - - rho_l = Uk[j * nbComp + 0] - q_l = Uk[j * nbComp + 1] - rho_E_l = Uk[j * nbComp + 2] - rho_r = Uk[(j ) * nbComp + 0] # We take rho inside the domain - q_r = Uk[(j ) * nbComp + 1] # We take q from inside the domain - rho_E_r = (p_outlet+gamma*p0)/(gamma-1) + 0.5*q_r**2/rho_r # rhoE is obtained using the boundary condition p_outlet - - Am = jacobianMatricesm(dt / dx, rho_l, q_l, rho_E_l, rho_r, q_r, rho_E_r) - - if(isImplicit): - dUi1[0] = rho_r - Uk[j * nbComp + 0] - dUi1[1] = q_r - Uk[j * nbComp + 1] - dUi1[2] = rho_E_r - Uk[j * nbComp + 2] - temp1 = Am * dUi1 - - 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): - 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 FillRHSFromInnerCell(j, Uk, nbComp, Rhs, Un, dt, dx, isImplicit): - - 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) - - 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) - - if(isImplicit):#Contribution to the right hand side if te scheme is implicit - 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 computeSystemMatrix(a,b,nx, cfl, Uk, isImplicit): - dim = 1 - nbComp = 3 - dx = (b - a) / nx - dt = cfl * dx / c0 - nbVoisinsMax = 2 - - nbCells = nx - divMat = cdmath.SparseMatrixPetsc(nbCells * nbComp, nbCells * nbComp, (nbVoisinsMax + 1) * nbComp) - - divMat.zeroEntries() #sets the matrix coefficients to zero - for j in range(nbCells): - - # traitements des bords - if (j == 0): - FillMatrixFromEdges(j, Uk, nbComp, divMat, dt, dx) - elif (j == nbCells - 1): - FillMatrixFromEdges(j, Uk, nbComp, divMat, dt, dx) - # traitement des cellules internes - else: - FillMatrixFromInnerCell(j, Uk, nbComp, divMat, dt, dx) - - if(isImplicit): - divMat.diagonalShift(1.) # add one on the diagonal - else: - divMat*=-1. - divMat.diagonalShift(1.) # add one on the diagonal - - return divMat - -def computeRHSVector(a,b,nx, cfl, Uk, Un, isImplicit): - dim = 1 - nbComp = 3 - dx = (b - a) / nx - dt = cfl * dx / c0 - nbVoisinsMax = 2 - - nbCells = nx - Rhs = cdmath.Vector(nbCells * (nbComp)) - - for j in range(nbCells): - - # traitements des bords - if (j == 0): - FillRHSFromEdges(j, Uk, nbComp, Rhs, Un, dt, dx, isImplicit) - elif (j == nbCells - 1): - FillRHSFromEdges(j, Uk, nbComp, Rhs, Un, dt, dx, isImplicit) - # traitement des cellules internes - else: - FillRHSFromInnerCell(j, Uk, nbComp, Rhs, Un, dt, dx, isImplicit) - - return Rhs - - -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) - state_law = "Hermite575K" - state_law_parameters(state_law) - 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 function computeSystemMatrix \n") - cfl = 0.5 - computeSystemMatrix(a, b, nx, cfl, U,True) #Implicit matrix - computeSystemMatrix(a, b, nx, cfl, U,False) #Explicit matrix - - print("\n Testing function computeRHSVector \n") - cfl = 0.5 - computeRHSVector(a, b, nx, cfl, U, U,True) #Implicit RHS - computeRHSVector(a, b, nx, cfl, U, U,False) #Explicit RHS - - print("\n Testing function solve (Implicit scheme) \n") - isImplicit=True - cfl = 1000. - solve(a, b, nx, "RegularSquares", "", cfl, state_law, isImplicit) - - print("\n Testing function solve (Explicit scheme) \n") - isImplicit=False - cfl = 0.5 - solve(a, b, nx, "RegularSquares", "", cfl, state_law, isImplicit) - diff --git a/CDMATH/tests/examples/EulerEquations/FullEulerSystem/EulerSystem1D_HeatedChannel_MAC/CMakeLists.txt b/CDMATH/tests/examples/EulerEquations/FullEulerSystem/EulerSystem1D_HeatedChannel_MAC/CMakeLists.txt new file mode 100755 index 0000000..6cc4002 --- /dev/null +++ b/CDMATH/tests/examples/EulerEquations/FullEulerSystem/EulerSystem1D_HeatedChannel_MAC/CMakeLists.txt @@ -0,0 +1,8 @@ + +if (CDMATH_WITH_PYTHON AND CDMATH_WITH_PETSC AND CDMATH_WITH_POSTPRO) + + ADD_TEST(ExampleFullEulerSystem_1DHeatedChannel_MAC ${PYTHON_EXECUTABLE} ${CMAKE_CURRENT_SOURCE_DIR}/Euler_complet_HeatedChanel_MAC.py ) + +endif (CDMATH_WITH_PYTHON AND CDMATH_WITH_PETSC AND CDMATH_WITH_POSTPRO) + + diff --git a/CDMATH/tests/examples/EulerEquations/FullEulerSystem/EulerSystem1D_HeatedChannel_MAC/Euler_complet_HeatedChanel_MAC.py b/CDMATH/tests/examples/EulerEquations/FullEulerSystem/EulerSystem1D_HeatedChannel_MAC/Euler_complet_HeatedChanel_MAC.py new file mode 100644 index 0000000..a2a0a23 --- /dev/null +++ b/CDMATH/tests/examples/EulerEquations/FullEulerSystem/EulerSystem1D_HeatedChannel_MAC/Euler_complet_HeatedChanel_MAC.py @@ -0,0 +1,893 @@ +#!/usr/bin/env python3 +# -*-coding:utf-8 -* + +""" +Created on Mon Sep 27 2021 +@author: Michael NDJINGA, Katia Ait Ameur, Coraline Mounier + +Euler system with heating source term (phi) in one dimension on regular domain [a,b] +Riemann problemn with ghost cell boundary condition +Left : Inlet boundary condition (velocity and temperature imposed) +Right : Outlet boundary condition (pressure imposed) +Staggered scheme +Regular square mesh + +State law Stiffened gaz : p = (gamma - 1) * rho * (e - q) - gamma * p0 +4 choices of parameters gamma and p0 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 +from math import sqrt +from numpy import sign + + +#### Initial and boundary condition (T in K, v in m/s, p in Pa) +T_inlet = 565. +v_inlet = 5. +p_outlet = 155.0 * 10**5 + +#initial parameters are determined from boundary conditions +p_0 = p_outlet #initial pressure +v_0 = v_inlet #initial velocity +T_0 = T_inlet #initial temperature +### Heating source term +phi=1.e8 + +## Numerical parameter +precision = 1e-6 + +#state law parameter : can be 'Lagrange', 'Hermite590K', 'Hermite617K', or 'FLICA' +state_law = "Hermite575K" + +def state_law_parameters(state_law): + #state law Stiffened Gaz : p = (gamma - 1) * rho * e - gamma * p0 + global gamma + global p0 + global q + global c0 + global cp + global h_sat + global T_sat + + if state_law == "Lagrange": + # reference values for Lagrange interpolation + p_ref = 155. * 10**5 #Reference pressure in a REP 900 nuclear power plant + p1 = 153. * 10**5 # value of pressure at inlet of a 900 MWe PWR vessel + rho_ref = 594.38 #density of water at saturation temperature of 617.94K and 155 bars + rho1 = 742.36 # value of density at inlet of a 900 MWe PWR vessel (T1 = 565K) + e_ref = 1603.8 * 10**3 #internal energy of water at saturation temperature of 617.94K and 155 bars + e1 = 1273.6 * 10**3 # value of internal energy at inlet of a 900 MWe PWR vessel + + gamma = (p1 - p_ref) / (rho1 * e1 - rho_ref *e_ref) + 1. + p0 = - 1. / gamma * ( - (gamma - 1) * rho_ref * e_ref + p_ref) + q=0. + c0 = sqrt((gamma - 1) * (e_ref + p_ref / rho_ref)) + + cp = 8950. + h_sat = 1.63 * 10 ** 6 # saturation enthalpy of water at 155 bars + T_sat = 617.94 + + elif state_law == "Hermite617K": + # reference values for Hermite interpolation + p_ref = 155. * 10**5 #Reference pressure in a REP 900 nuclear power plant + T_ref = 617.94 #Reference temperature for interpolation at 617.94K + rho_ref = 594.38 #density of water at saturation temperature of 617.94K and 155 bars + e_ref = 1603.8 * 10**3 #internal energy of water at saturation temperature of 617.94K and 155 bars + h_ref = e_ref + p_ref / rho_ref + c_ref = 621.43 #sound speed for water at 155 bars and 617.94K + + gamma = 1. + c_ref * c_ref / (e_ref + p_ref / rho_ref) # From the sound speed formula + p0 = 1. / gamma * ( (gamma - 1) * rho_ref * e_ref - p_ref) + q=0. + c0 = sqrt((gamma - 1) * (e_ref + p_ref / rho_ref)) + + cp = 8950. # value at 155 bar and 617.94K + h_sat = 1.63 * 10 ** 6 # saturation enthalpy of water at 155 bars + T_sat = 617.94 + + elif state_law == 'Hermite590K': + # reference values for Hermite interpolation + p_ref = 155. * 10**5 #Reference pressure in a REP 900 nuclear power plant + T_ref = 590. #Reference temperature for interpolation at 590K + rho_ref = 688.3 #density of water at 590K and 155 bars + e_ref = 1411.4 * 10**3 #internal energy of water at 590K and 155 bars + h_ref = e_ref + p_ref / rho_ref + c_ref = 866.29 #sound speed for water at 155 bars and 590K + + gamma = 1. + c_ref * c_ref / (e_ref + p_ref / rho_ref) # From the sound speed formula + p0 = 1. / gamma * ( (gamma - 1) * rho_ref * e_ref - p_ref) + q=0. + c0 = sqrt((gamma - 1) * (e_ref + p_ref / rho_ref)) + + cp = 5996.8 # value at 155 bar and 590K + h_sat = 1433.9 * 10 ** 3 # saturation enthalpy of water at 155 bars + T_sat = 590. + + elif state_law == 'Hermite575K': + # reference values for Hermite interpolation + p_ref = 155 * 10**5 #Reference pressure in a REP 900 nuclear power plant + T_ref = 575 #Reference temperature at inlet in a REP 900 nuclear power plant + #Remaining values determined using iapws python package + rho_ref = 722.66 #density of water at 575K and 155 bars + e_ref = 1326552.66 #internal energy of water at 575K and 155 bars + h_ref = e_ref + p_ref / rho_ref + c_ref = 959.28 #sound speed for water at 155 bars and 575K + + gamma = 1 + c_ref * c_ref / (e_ref + p_ref / rho_ref) # From the sound speed formula + p0 = 1 / gamma * ( (gamma - 1) * rho_ref * e_ref - p_ref) + q=0. + c0 = sqrt((gamma - 1) * (e_ref + p_ref / rho_ref))#This is actually c_ref + + cp = 5504.05 # value at 155 bar and 590K + h_sat = h_ref # saturation enthalpy of water at 155 bars + T_sat = T_ref + else: + raise ValueError("Incorrect value for parameter state_law") + +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) + + p_initial = np.array([ p_0 for xi in x]) + v_initial = np.array([ v_0 for xi in x]) + T_initial = np.array([ T_0 for xi in x]) + + rho_initial = p_to_rho_StiffenedGaz(p_initial, T_initial) + q_initial = rho_initial * v_initial + rho_E_initial = T_to_rhoE_StiffenedGaz(T_initial, rho_initial, q_initial) + + return rho_initial, q_initial, rho_E_initial, p_initial, v_initial, T_initial + +def p_to_rho_StiffenedGaz(p_field, T_field): + rho_field = (p_field + p0) * 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): + rho_E_field = 1. / 2. * (q_field) ** 2 / rho_field + p0 + rho_field / gamma * (h_sat + cp * (T_field- T_sat)) + return rho_E_field + +def rhoE_to_T_StiffenedGaz(rho_field, q_field, rho_E_field): + T_field = T_sat + 1 / cp * (gamma * (rho_E_field / rho_field - 1 / 2 * (q_field / rho_field) ** 2) - gamma * p0 / rho_field - h_sat) + return T_field + +def rho_to_p_StiffenedGaz(rho_field, q_field, rho_E_field): + p_field = (gamma - 1) * (rho_E_field - 1. / 2 * q_field ** 2 / rho_field) - gamma * p0 + return p_field + +def T_to_E_StiffenedGaz(p_field, T_field, v_field): + rho_field = p_to_rho_StiffenedGaz(p_field, T_field) + E_field = (p_field + gamma * p0) / ((gamma-1) * rho_field) + 0.5 * v_field **2 + return E_field + +def dp_drho_e_const_StiffenedGaz( e ): + return (gamma-1)*(e-q) + +def dp_de_rho_const_StiffenedGaz( rho ): + return (gamma-1)*rho + +def sound_speed_StiffenedGaz( h ): + return np.sqrt((gamma-1)*(h-q)) + +def rho_h_to_p_StiffenedGaz( rho, h ): + return (gamma - 1) * rho * ( h - q ) / gamma - p0 + +def MatRoe_StiffenedGaz( rho_l, q_l, rho_E_l, rho_r, q_r, rho_E_r): + 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) + p_r = rho_to_p_StiffenedGaz(rho_r, q_r, rho_E_r) + 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. ) + e = H - p / rho - 1./2 * u**2 + dp_drho = dp_drho_e_const_StiffenedGaz( e ) + dp_de = dp_de_rho_const_StiffenedGaz( rho ) + + 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): + 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) + p_r = rho_to_p_StiffenedGaz(rho_r, q_r, rho_E_r) + 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. ) + e = H - p / rho - 1./2 * u**2 + dp_drho = dp_drho_e_const_StiffenedGaz( e ) + dp_de = dp_de_rho_const_StiffenedGaz( rho ) + + #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) + +def jacobianMatricesm(coeff, rho_l, q_l, rho_E_l, rho_r, q_r, rho_E_r): + + 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) + + Dmac = Dmac_StiffenedGaz( rho_l, q_l, rho_E_l, rho_r, q_r, rho_E_r) + return (RoeMat - Dmac) * coeff * 0.5 + + +def jacobianMatricesp(coeff, rho_l, q_l, rho_E_l, rho_r, q_r, rho_E_r): + 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) + Dmac = Dmac_StiffenedGaz( rho_l, q_l, rho_E_l, rho_r, q_r, rho_E_r) + + return (RoeMat + Dmac) * coeff * 0.5 + + +def FillEdges(j, Uk, nbComp, divMat, Rhs, Un, dt, dx): + 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) + divMat.addValue(j * nbComp, (j + 1) * nbComp, Am) + divMat.addValue(j * nbComp, j * nbComp, Am * (-1.)) + + p_inlet = rho_to_p_StiffenedGaz(Uk[j * nbComp + 0], Uk[j * nbComp + 1], Uk[j * nbComp + 2])# We take p from inside the domain + rho_l=p_to_rho_StiffenedGaz(p_inlet, T_inlet) # rho is computed from the temperature BC and the inner pressure + q_l = rho_l * v_inlet # q is imposed by the boundary condition v_inlet + rho_E_l = T_to_rhoE_StiffenedGaz(T_inlet, rho_l, q_l) #rhoE is obtained using the two boundary conditions v_inlet and e_inlet + 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) + divMat.addValue(j * nbComp, j * nbComp, Ap) + + 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 + + dUi2[0] = rho_l - Uk[(j ) * nbComp + 0] + dUi2[1] = q_l - Uk[(j ) * nbComp + 1] + dUi2[2] = rho_E_l - Uk[(j ) * nbComp + 2] + 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) + divMat.addValue(j * nbComp, j * nbComp, Ap) + divMat.addValue(j * nbComp, (j - 1) * nbComp, Ap * (-1.)) + + rho_l = Uk[j * nbComp + 0] + q_l = Uk[j * nbComp + 1] + rho_E_l = Uk[j * nbComp + 2] + rho_r = Uk[(j ) * nbComp + 0] # We take rho inside the domain + q_r = Uk[(j ) * nbComp + 1] # We take q from inside the domain + rho_E_r = (p_outlet+gamma*p0)/(gamma-1) + 0.5*q_r**2/rho_r # rhoE is obtained using the boundary condition p_outlet + + Am = jacobianMatricesm(dt / dx, rho_l, q_l, rho_E_l, rho_r, q_r, rho_E_r) + divMat.addValue(j * nbComp, j * nbComp, Am * (-1.)) + + dUi1[0] = rho_r - Uk[j * nbComp + 0] + dUi1[1] = q_r - Uk[j * nbComp + 1] + dUi1[2] = rho_E_r - Uk[j * nbComp + 2] + temp1 = Am * dUi1 + + 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 + + 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]) + +def FillInnerCell(j, Uk, nbComp, divMat, Rhs, Un, dt, dx): + + 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) + + 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) + + 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(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]) + +def SetPicture(rho_field, q_field, h_field, p_field, v_field, T_field, dx): + max_initial_rho = max(rho_field) + min_initial_rho = min(rho_field) + max_initial_q = max(q_field) + min_initial_q = min(q_field) + min_initial_h = min(h_field) + max_initial_h = max(h_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) + max_initial_T = max(T_field) + min_initial_T = min(T_field) + + fig, ([axDensity, axMomentum, axh],[axPressure, axVitesse, axTemperature]) = plt.subplots(2, 3,sharex=True, figsize=(14,10)) + plt.gcf().subplots_adjust(wspace = 0.5) + + lineDensity, = axDensity.plot([a+0.5*dx + i*dx for i in range(nx)], rho_field, label='MAC scheme') + axDensity.set(xlabel='x (m)', ylabel='Densité (kg/m3)') + axDensity.set_xlim(a,b) + axDensity.set_ylim(680, 800) + axDensity.legend() + + lineMomentum, = axMomentum.plot([a+0.5*dx + i*dx for i in range(nx)], q_field, label='MAC scheme') + axMomentum.set(xlabel='x (m)', ylabel='Momentum (kg/(m2.s))') + axMomentum.set_xlim(a,b) + axMomentum.set_ylim(3760, 3780) + axMomentum.legend() + + lineh, = axh.plot([a+0.5*dx + i*dx for i in range(nx)], h_field, label='MAC scheme') + axh.set(xlabel='x (m)', ylabel='h (J/m3)') + axh.set_xlim(a,b) + axh.set_ylim(1.25 * 10**6, 1.45*10**6) + axh.legend() + + linePressure, = axPressure.plot([a+0.5*dx + i*dx for i in range(nx)], p_field, label='MAC scheme') + axPressure.set(xlabel='x (m)', ylabel='Pression (bar)') + axPressure.set_xlim(a,b) + axPressure.set_ylim(154.99, 155.02) + axPressure.ticklabel_format(axis='y', style='sci', scilimits=(0,0)) + axPressure.legend() + + lineVitesse, = axVitesse.plot([a+0.5*dx + i*dx for i in range(nx)], v_field, label='MAC scheme') + axVitesse.set(xlabel='x (m)', ylabel='Vitesse (m/s)') + axVitesse.set_xlim(a,b) + axVitesse.set_ylim(v_0-1, v_0+1) + axVitesse.ticklabel_format(axis='y', style='sci', scilimits=(0,0)) + axVitesse.legend() + + lineTemperature, = axTemperature.plot([a+0.5*dx + i*dx for i in range(nx)], T_field, label='MAC scheme') + axTemperature.set(xlabel='x (m)', ylabel='Température (K)') + axTemperature.set_xlim(a,b) + axTemperature.set_ylim(T_0-10, T_0+30) + axTemperature.ticklabel_format(axis='y', style='sci', scilimits=(0,0)) + axTemperature.legend() + + return(fig, lineDensity, lineMomentum, lineh, linePressure, lineVitesse, lineTemperature) + + +def EulerSystemMAC(ntmax, tmax, cfl, a, b, nbCells, output_freq, meshName, state_law): + state_law_parameters(state_law) + dim = 1 + nbComp = 3 + dt = 0. + time = 0. + it = 0 + isStationary = False + dx = (b - a) / nx + dt = cfl * dx / c0 + #dt = 5*10**(-6) + nbVoisinsMax = 2 + + # 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, T_field = initial_conditions_Riemann_problem(a, b, nx) + 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) + + # Picture settings + fig, lineDensity, lineMomentum, lineRhoE, linePressure, lineVitesse, lineTemperature = SetPicture( rho_field, q_field, h_field, p_field, v_field, T_field, dx) + + plt.savefig("EulerComplet_HeatedChannel_" + str(dim) + "D_MAC" + meshName + "0" + ".png") + iterGMRESMax = 50 + newton_max = 100 + + print("Starting computation of the non linear Euler non isentropic system with MAC 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(j, Uk, nbComp, divMat, Rhs, Un, dt, dx) + elif (j == nbCells - 1): + FillEdges(j, Uk, nbComp, divMat, Rhs, Un, dt, dx) + + # traitement des cellules internes + else: + FillInnerCell(j, Uk, nbComp, divMat, Rhs, Un, dt, dx) + + Rhs[j * nbComp + 2] += phi*dt + + #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])/rho0, abs(vector_residu[1])/(rho0*v_0), abs(vector_residu[2])/rhoE0 ) + + if (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 (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 MAC 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])/rho0, abs(residu_stat[1])/(rho0*v_0), abs(residu_stat[2])/rhoE0 )) + + if ( it>1 and abs(residu_stat[0])/rho0= ntmax or isStationary or time >= tmax): + + print("-- Time step : " + str(it) + ", Time: " + str(time) + ", dt: " + str(dt)) + + print("Temperature gain between inlet and outlet is ", T_field[nbCells-1]-T_field[0],"\n") + + plt.savefig("EulerComplet_HeatedChannel_" + str(dim) + "D_MAC" + meshName + str(it) + '_time' + str(time) + ".png") + + print("-- Iter: " + str(it) + ", Time: " + str(time) + ", dt: " + str(dt)+"\n") + + if (it >= ntmax): + print("Maximum number of time steps ntmax= ", ntmax, " reached") + return + + elif (isStationary): + print("Stationary regime reached at time step ", it, ", t= ", time) + print("------------------------------------------------------------------------------------") + np.savetxt("EulerComplet_HeatedChannel_" + str(dim) + "D_MAC" + meshName + "_rho_Stat.txt", rho_field, delimiter="\n") + np.savetxt("EulerComplet_HeatedChannel_" + str(dim) + "D_MAC" + meshName + "_q_Stat.txt", q_field, delimiter="\n") + np.savetxt("EulerComplet_HeatedChannel_" + str(dim) + "D_MAC" + meshName + "_rhoE_Stat.txt", rho_E_field, delimiter="\n") + np.savetxt("EulerComplet_HeatedChannel_" + str(dim) + "D_MAC" + meshName + "_p_Stat.txt", p_field, delimiter="\n") + plt.savefig("EulerComplet_HeatedChannel_" + str(dim) + "D_MAC" + meshName +"_Stat.png") + return + else: + print("Maximum time Tmax= ", tmax, " reached") + return + + +def solve(a, b, nx, meshName, meshType, cfl, state_law): + print("Simulation of a heated channel in dimension 1 on " + str(nx) + " cells") + print("State Law Stiffened Gaz, " + state_law) + print("Initial data : ", "constant fields") + print("Boundary conditions : ", "Inlet (Left), Outlet (Right)") + print("Mesh name : ", meshName, ", ", nx, " cells") + # Problem data + tmax = 10. + ntmax = 100000 + output_freq = 1000 + EulerSystemMAC(ntmax, tmax, cfl, a, b, nx, output_freq, meshName, state_law) + return + +def FillMatrixFromEdges(j, Uk, nbComp, divMat, dt, dx): + + 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) + divMat.addValue(j * nbComp, (j + 1) * nbComp, Am) + divMat.addValue(j * nbComp, j * nbComp, Am * (-1.)) + + p_inlet = rho_to_p_StiffenedGaz(Uk[j * nbComp + 0], Uk[j * nbComp + 1], Uk[j * nbComp + 2])# We take p from inside the domain + rho_l=p_to_rho_StiffenedGaz(p_inlet, T_inlet) # rho is computed from the temperature BC and the inner pressure + #rho_l = Uk[j * nbComp + 0] # We take rho from inside the domain + q_l = rho_l * v_inlet # q is imposed by the boundary condition v_inlet + rho_E_l = T_to_rhoE_StiffenedGaz(T_inlet, rho_l, q_l) #rhoE is obtained using the two boundary conditions v_inlet and e_inlet + 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) + divMat.addValue(j * nbComp, j * nbComp, Ap) + + elif (j == nx - 1): + rho_l = Uk[j * nbComp + 0] + q_l = Uk[j * nbComp + 1] + rho_E_l = Uk[j * nbComp + 2] + rho_r = Uk[(j ) * nbComp + 0] # We take rho inside the domain + q_r = Uk[(j ) * nbComp + 1] # We take q from inside the domain + rho_E_r = (p_outlet+gamma*p0)/(gamma-1) + 0.5*q_r**2/rho_r # rhoE is obtained using the boundary condition p_outlet + + Am = jacobianMatricesm(dt / dx, rho_l, q_l, rho_E_l, rho_r, q_r, rho_E_r) + divMat.addValue(j * nbComp, j * nbComp, Am * (-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) + divMat.addValue(j * nbComp, j * nbComp, Ap) + divMat.addValue(j * nbComp, (j - 1) * nbComp, Ap * (-1.)) + + +def FillMatrixFromInnerCell(j, Uk, nbComp, divMat, dt, dx): + + 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) + + 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) + + 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.)) + +def FillRHSFromEdges(j, Uk, nbComp, Rhs, Un, dt, dx): + 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) + + p_inlet = rho_to_p_StiffenedGaz(Uk[j * nbComp + 0], Uk[j * nbComp + 1], Uk[j * nbComp + 2])# We take p from inside the domain + rho_l=p_to_rho_StiffenedGaz(p_inlet, T_inlet) # rho is computed from the temperature BC and the inner pressure + #rho_l = Uk[j * nbComp + 0] # We take rho from inside the domain + q_l = rho_l * v_inlet # q is imposed by the boundary condition v_inlet + rho_E_l = T_to_rhoE_StiffenedGaz(T_inlet, rho_l, q_l) #rhoE is obtained using the two boundary conditions v_inlet and e_inlet + 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) + + 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 + + dUi2[0] = rho_l - Uk[(j ) * nbComp + 0] + dUi2[1] = q_l - Uk[(j ) * nbComp + 1] + dUi2[2] = rho_E_l - Uk[(j ) * nbComp + 2] + 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) + + rho_l = Uk[j * nbComp + 0] + q_l = Uk[j * nbComp + 1] + rho_E_l = Uk[j * nbComp + 2] + rho_r = Uk[(j ) * nbComp + 0] # We take rho inside the domain + q_r = Uk[(j ) * nbComp + 1] # We take q from inside the domain + rho_E_r = (p_outlet+gamma*p0)/(gamma-1) + 0.5*q_r**2/rho_r # rhoE is obtained using the boundary condition p_outlet + + Am = jacobianMatricesm(dt / dx, rho_l, q_l, rho_E_l, rho_r, q_r, rho_E_r) + + dUi1[0] = rho_r - Uk[j * nbComp + 0] + dUi1[1] = q_r - Uk[j * nbComp + 1] + dUi1[2] = rho_E_r - Uk[j * nbComp + 2] + temp1 = Am * dUi1 + + 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 + + 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]) + +def FillRHSFromInnerCell(j, Uk, nbComp, Rhs, Un, dt, dx): + + 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) + + 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) + + 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]) + + +def computeSystemMatrix(a,b,nx, cfl, Uk): + dim = 1 + nbComp = 3 + dx = (b - a) / nx + dt = cfl * dx / c0 + nbVoisinsMax = 2 + + nbCells = nx + divMat = cdmath.SparseMatrixPetsc(nbCells * nbComp, nbCells * nbComp, (nbVoisinsMax + 1) * nbComp) + + divMat.zeroEntries() #sets the matrix coefficients to zero + for j in range(nbCells): + + # traitements des bords + if (j == 0): + FillMatrixFromEdges(j, Uk, nbComp, divMat, dt, dx) + elif (j == nbCells - 1): + FillMatrixFromEdges(j, Uk, nbComp, divMat, dt, dx) + # traitement des cellules internes + else: + FillMatrixFromInnerCell(j, Uk, nbComp, divMat, dt, dx) + + divMat.diagonalShift(1.) # add one on the diagonal + + return divMat + +def computeRHSVector(a,b,nx, cfl, Uk, Un): + dim = 1 + nbComp = 3 + dx = (b - a) / nx + dt = cfl * dx / c0 + nbVoisinsMax = 2 + + nbCells = nx + Rhs = cdmath.Vector(nbCells * (nbComp)) + + for j in range(nbCells): + + # traitements des bords + if (j == 0): + FillRHSFromEdges(j, Uk, nbComp, Rhs, Un, dt, dx) + elif (j == nbCells - 1): + FillRHSFromEdges(j, Uk, nbComp, Rhs, Un, dt, dx) + # traitement des cellules internes + else: + FillRHSFromInnerCell(j, Uk, nbComp, Rhs, Un, dt, dx) + + return Rhs + + +if __name__ == """__main__""": + nbComp=3 # number of equations + a = 0.# domain is interval [a,b] + b = 4.2# domain is interval [a,b] + nx = 50# 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) + state_law = "Hermite575K" + state_law_parameters(state_law) + 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 function computeSystemMatrix \n") + cfl = 0.5 + computeSystemMatrix(a, b, nx, cfl, U) + + print("\n Testing function computeRHSVector \n") + cfl = 0.5 + computeRHSVector(a, b, nx, cfl, U, U) + + print("\n Testing function solve \n") + cfl = 1000. + solve(a, b, nx, "RegularSquares", "", cfl, state_law) + diff --git a/CDMATH/tests/examples/EulerEquations/FullEulerSystem/EulerSystem1D_HeatedChannel_Roe/CMakeLists.txt b/CDMATH/tests/examples/EulerEquations/FullEulerSystem/EulerSystem1D_HeatedChannel_Roe/CMakeLists.txt new file mode 100755 index 0000000..bf750c8 --- /dev/null +++ b/CDMATH/tests/examples/EulerEquations/FullEulerSystem/EulerSystem1D_HeatedChannel_Roe/CMakeLists.txt @@ -0,0 +1,8 @@ + +if (CDMATH_WITH_PYTHON AND CDMATH_WITH_PETSC AND CDMATH_WITH_POSTPRO) + + ADD_TEST(ExampleFullEulerSystem_1DHeatedChannel_Roe ${PYTHON_EXECUTABLE} ${CMAKE_CURRENT_SOURCE_DIR}/Euler_complet_HeatedChanel_Roe.py ) + +endif (CDMATH_WITH_PYTHON AND CDMATH_WITH_PETSC AND CDMATH_WITH_POSTPRO) + + diff --git a/CDMATH/tests/examples/EulerEquations/FullEulerSystem/EulerSystem1D_HeatedChannel_Roe/Euler_complet_HeatedChanel_Roe.py b/CDMATH/tests/examples/EulerEquations/FullEulerSystem/EulerSystem1D_HeatedChannel_Roe/Euler_complet_HeatedChanel_Roe.py new file mode 100644 index 0000000..859cde8 --- /dev/null +++ b/CDMATH/tests/examples/EulerEquations/FullEulerSystem/EulerSystem1D_HeatedChannel_Roe/Euler_complet_HeatedChanel_Roe.py @@ -0,0 +1,979 @@ +#!/usr/bin/env python3 +# -*-coding:utf-8 -* + +""" +Created on Mon Aug 30 2021 +@author: Michael NDJINGA, Katia Ait Ameur, Coraline Mounier + +Euler system with heating source term (phi) in one dimension on regular domain [a,b] +Riemann problemn with ghost cell boundary condition +Left : Inlet boundary condition (velocity and temperature imposed) +Right : Outlet boundary condition (pressure imposed) +Roe scheme +Regular square mesh + +State law Stiffened gaz : p = (gamma - 1) * rho * (e - q) - gamma * p0 +4 choices of parameters gamma and p0 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 +from math import sqrt +from numpy import sign + + +#### Initial and boundary condition (T in K, v in m/s, p in Pa) +T_inlet = 565. +v_inlet = 5. +p_outlet = 155.0 * 10**5 + +#initial parameters are determined from boundary conditions +p_0 = p_outlet #initial pressure +v_0 = v_inlet #initial velocity +T_0 = T_inlet #initial temperature +### Heating source term +phi=1.e8 + +## Numerical parameter +precision = 1e-6 + +#state law parameter : can be 'Lagrange', 'Hermite590K', 'Hermite617K', or 'FLICA' +state_law = "Hermite575K" + +def state_law_parameters(state_law): + #state law Stiffened Gaz : p = (gamma - 1) * rho * e - gamma * p0 + global gamma + global p0 + global q + global c0 + global cp + global h_sat + global T_sat + + if state_law == "Lagrange": + # reference values for Lagrange interpolation + p_ref = 155. * 10**5 #Reference pressure in a REP 900 nuclear power plant + p1 = 153. * 10**5 # value of pressure at inlet of a 900 MWe PWR vessel + rho_ref = 594.38 #density of water at saturation temperature of 617.94K and 155 bars + rho1 = 742.36 # value of density at inlet of a 900 MWe PWR vessel (T1 = 565K) + e_ref = 1603.8 * 10**3 #internal energy of water at saturation temperature of 617.94K and 155 bars + e1 = 1273.6 * 10**3 # value of internal energy at inlet of a 900 MWe PWR vessel + + gamma = (p1 - p_ref) / (rho1 * e1 - rho_ref *e_ref) + 1. + p0 = - 1. / gamma * ( - (gamma - 1) * rho_ref * e_ref + p_ref) + q=0. + c0 = sqrt((gamma - 1) * (e_ref + p_ref / rho_ref)) + + cp = 8950. + h_sat = 1.63 * 10 ** 6 # saturation enthalpy of water at 155 bars + T_sat = 617.94 + + elif state_law == "Hermite617K": + # reference values for Hermite interpolation + p_ref = 155. * 10**5 #Reference pressure in a REP 900 nuclear power plant + T_ref = 617.94 #Reference temperature for interpolation at 617.94K + rho_ref = 594.38 #density of water at saturation temperature of 617.94K and 155 bars + e_ref = 1603.8 * 10**3 #internal energy of water at saturation temperature of 617.94K and 155 bars + h_ref = e_ref + p_ref / rho_ref + c_ref = 621.43 #sound speed for water at 155 bars and 617.94K + + gamma = 1. + c_ref * c_ref / (e_ref + p_ref / rho_ref) # From the sound speed formula + p0 = 1. / gamma * ( (gamma - 1) * rho_ref * e_ref - p_ref) + q=0. + c0 = sqrt((gamma - 1) * (e_ref + p_ref / rho_ref)) + + cp = 8950. # value at 155 bar and 617.94K + h_sat = 1.63 * 10 ** 6 # saturation enthalpy of water at 155 bars + T_sat = 617.94 + + elif state_law == 'Hermite590K': + # reference values for Hermite interpolation + p_ref = 155. * 10**5 #Reference pressure in a REP 900 nuclear power plant + T_ref = 590. #Reference temperature for interpolation at 590K + rho_ref = 688.3 #density of water at 590K and 155 bars + e_ref = 1411.4 * 10**3 #internal energy of water at 590K and 155 bars + h_ref = e_ref + p_ref / rho_ref + c_ref = 866.29 #sound speed for water at 155 bars and 590K + + gamma = 1. + c_ref * c_ref / (e_ref + p_ref / rho_ref) # From the sound speed formula + p0 = 1. / gamma * ( (gamma - 1) * rho_ref * e_ref - p_ref) + q=0. + c0 = sqrt((gamma - 1) * (e_ref + p_ref / rho_ref)) + + cp = 5996.8 # value at 155 bar and 590K + h_sat = 1433.9 * 10 ** 3 # saturation enthalpy of water at 155 bars + T_sat = 590. + + elif state_law == 'Hermite575K': + # reference values for Hermite interpolation + p_ref = 155 * 10**5 #Reference pressure in a REP 900 nuclear power plant + T_ref = 575 #Reference temperature at inlet in a REP 900 nuclear power plant + #Remaining values determined using iapws python package + rho_ref = 722.66 #density of water at 575K and 155 bars + e_ref = 1326552.66 #internal energy of water at 575K and 155 bars + h_ref = e_ref + p_ref / rho_ref + c_ref = 959.28 #sound speed for water at 155 bars and 575K + + gamma = 1 + c_ref * c_ref / (e_ref + p_ref / rho_ref) # From the sound speed formula + p0 = 1 / gamma * ( (gamma - 1) * rho_ref * e_ref - p_ref) + q=0. + c0 = sqrt((gamma - 1) * (e_ref + p_ref / rho_ref))#This is actually c_ref + + cp = 5504.05 # value at 155 bar and 590K + h_sat = h_ref # saturation enthalpy of water at 155 bars + T_sat = T_ref + else: + raise ValueError("Incorrect value for parameter state_law") + +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) + + p_initial = np.array([ p_0 for xi in x]) + v_initial = np.array([ v_0 for xi in x]) + T_initial = np.array([ T_0 for xi in x]) + + rho_initial = p_to_rho_StiffenedGaz(p_initial, T_initial) + q_initial = rho_initial * v_initial + rho_E_initial = T_to_rhoE_StiffenedGaz(T_initial, rho_initial, q_initial) + + return rho_initial, q_initial, rho_E_initial, p_initial, v_initial, T_initial + +def p_to_rho_StiffenedGaz(p_field, T_field): + rho_field = (p_field + p0) * 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): + rho_E_field = 1. / 2. * (q_field) ** 2 / rho_field + p0 + rho_field / gamma * (h_sat + cp * (T_field- T_sat)) + return rho_E_field + +def rhoE_to_T_StiffenedGaz(rho_field, q_field, rho_E_field): + T_field = T_sat + 1 / cp * (gamma * (rho_E_field / rho_field - 1 / 2 * (q_field / rho_field) ** 2) - gamma * p0 / rho_field - h_sat) + return T_field + +def rho_to_p_StiffenedGaz(rho_field, q_field, rho_E_field): + p_field = (gamma - 1) * (rho_E_field - 1. / 2 * q_field ** 2 / rho_field) - gamma * p0 + return p_field + +def T_to_E_StiffenedGaz(p_field, T_field, v_field): + rho_field = p_to_rho_StiffenedGaz(p_field, T_field) + E_field = (p_field + gamma * p0) / ((gamma-1) * rho_field) + 0.5 * v_field **2 + return E_field + +def dp_drho_e_const_StiffenedGaz( e ): + return (gamma-1)*(e-q) + +def dp_de_rho_const_StiffenedGaz( rho ): + return (gamma-1)*rho + +def sound_speed_StiffenedGaz( h ): + return np.sqrt((gamma-1)*(h-q)) + +def rho_h_to_p_StiffenedGaz( rho, h ): + return (gamma - 1) * rho * ( h - q ) / gamma - p0 + +def MatRoe_StiffenedGaz( rho_l, q_l, rho_E_l, rho_r, q_r, rho_E_r): + 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) + p_r = rho_to_p_StiffenedGaz(rho_r, q_r, rho_E_r) + 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. ) + e = H - p / rho - 1./2 * u**2 + dp_drho = dp_drho_e_const_StiffenedGaz( e ) + dp_de = dp_de_rho_const_StiffenedGaz( rho ) + + 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 Droe_StiffenedGaz( rho_l, q_l, rho_E_l, rho_r, q_r, rho_E_r): + 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) + p_r = rho_to_p_StiffenedGaz(rho_r, q_r, rho_E_r) + 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. ) + + 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): + + 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) + + Droe = Droe_StiffenedGaz( rho_l, q_l, rho_E_l, rho_r, q_r, rho_E_r) + return (RoeMat - Droe) * coeff * 0.5 + + +def jacobianMatricesp(coeff, rho_l, q_l, rho_E_l, rho_r, q_r, rho_E_r): + 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) + Droe = Droe_StiffenedGaz( rho_l, q_l, rho_E_l, rho_r, q_r, rho_E_r) + + return (RoeMat + Droe) * coeff * 0.5 + + +def FillEdges(j, Uk, nbComp, divMat, Rhs, Un, dt, dx, isImplicit): + 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) + divMat.addValue(j * nbComp, (j + 1) * nbComp, Am) + divMat.addValue(j * nbComp, j * nbComp, Am * (-1.)) + + p_inlet = rho_to_p_StiffenedGaz(Uk[j * nbComp + 0], Uk[j * nbComp + 1], Uk[j * nbComp + 2])# We take p from inside the domain + rho_l=p_to_rho_StiffenedGaz(p_inlet, T_inlet) # rho is computed from the temperature BC and the inner pressure + #rho_l = Uk[j * nbComp + 0] # We take rho from inside the domain + q_l = rho_l * v_inlet # q is imposed by the boundary condition v_inlet + rho_E_l = T_to_rhoE_StiffenedGaz(T_inlet, rho_l, q_l) #rhoE is obtained using the two boundary conditions v_inlet and e_inlet + 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) + 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 + + dUi2[0] = rho_l - Uk[(j ) * nbComp + 0] + dUi2[1] = q_l - Uk[(j ) * nbComp + 1] + dUi2[2] = rho_E_l - Uk[(j ) * nbComp + 2] + temp2 = Ap * dUi2 + 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) + divMat.addValue(j * nbComp, j * nbComp, Ap) + divMat.addValue(j * nbComp, (j - 1) * nbComp, Ap * (-1.)) + + rho_l = Uk[j * nbComp + 0] + q_l = Uk[j * nbComp + 1] + rho_E_l = Uk[j * nbComp + 2] + rho_r = Uk[(j ) * nbComp + 0] # We take rho inside the domain + q_r = Uk[(j ) * nbComp + 1] # We take q from inside the domain + rho_E_r = (p_outlet+gamma*p0)/(gamma-1) + 0.5*q_r**2/rho_r # rhoE is obtained using the boundary condition p_outlet + + Am = jacobianMatricesm(dt / dx, rho_l, q_l, rho_E_l, rho_r, q_r, rho_E_r) + divMat.addValue(j * nbComp, j * nbComp, Am * (-1.)) + + if(isImplicit): + dUi1[0] = rho_r - Uk[j * nbComp + 0] + dUi1[1] = q_r - Uk[j * nbComp + 1] + dUi1[2] = rho_E_r - Uk[j * nbComp + 2] + temp1 = Am * dUi1 + + 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): + + 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) + + 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) + + 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 SetPicture(rho_field_Roe, q_field_Roe, h_field_Roe, p_field_Roe, v_field_Roe, T_field_Roe, dx): + max_initial_rho = max(rho_field_Roe) + min_initial_rho = min(rho_field_Roe) + max_initial_q = max(q_field_Roe) + min_initial_q = min(q_field_Roe) + min_initial_h = min(h_field_Roe) + max_initial_h = max(h_field_Roe) + max_initial_p = max(p_field_Roe) + min_initial_p = min(p_field_Roe) + max_initial_v = max(v_field_Roe) + min_initial_v = min(v_field_Roe) + max_initial_T = max(T_field_Roe) + min_initial_T = min(T_field_Roe) + + fig, ([axDensity, axMomentum, axh],[axPressure, axVitesse, axTemperature]) = plt.subplots(2, 3,sharex=True, figsize=(14,10)) + plt.gcf().subplots_adjust(wspace = 0.5) + + lineDensity_Roe, = axDensity.plot([a+0.5*dx + i*dx for i in range(nx)], rho_field_Roe, label='Roe') + axDensity.set(xlabel='x (m)', ylabel='Densité (kg/m3)') + axDensity.set_xlim(a,b) + axDensity.set_ylim(680, 800) + axDensity.legend() + + lineMomentum_Roe, = axMomentum.plot([a+0.5*dx + i*dx for i in range(nx)], q_field_Roe, label='Roe') + axMomentum.set(xlabel='x (m)', ylabel='Momentum (kg/(m2.s))') + axMomentum.set_xlim(a,b) + axMomentum.set_ylim(3500, 4000) + axMomentum.legend() + + lineh_Roe, = axh.plot([a+0.5*dx + i*dx for i in range(nx)], h_field_Roe, label='Roe') + axh.set(xlabel='x (m)', ylabel='h (J/m3)') + axh.set_xlim(a,b) + axh.set_ylim(1.2 * 10**6, 1.5*10**6) + axh.legend() + + linePressure_Roe, = axPressure.plot([a+0.5*dx + i*dx for i in range(nx)], p_field_Roe, label='Roe') + axPressure.set(xlabel='x (m)', ylabel='Pression (bar)') + axPressure.set_xlim(a,b) + axPressure.set_ylim(155, 155.5) + axPressure.ticklabel_format(axis='y', style='sci', scilimits=(0,0)) + axPressure.legend() + + lineVitesse_Roe, = axVitesse.plot([a+0.5*dx + i*dx for i in range(nx)], v_field_Roe, label='Roe') + axVitesse.set(xlabel='x (m)', ylabel='Vitesse (m/s)') + axVitesse.set_xlim(a,b) + axVitesse.set_ylim(v_0-1, v_0+1) + axVitesse.ticklabel_format(axis='y', style='sci', scilimits=(0,0)) + axVitesse.legend() + + lineTemperature_Roe, = axTemperature.plot([a+0.5*dx + i*dx for i in range(nx)], T_field_Roe, label='Roe') + axTemperature.set(xlabel='x (m)', ylabel='Température (K)') + axTemperature.set_xlim(a,b) + axTemperature.set_ylim(T_0-10, T_0+30) + axTemperature.ticklabel_format(axis='y', style='sci', scilimits=(0,0)) + axTemperature.legend() + + return(fig, lineDensity_Roe, lineMomentum_Roe, lineh_Roe, linePressure_Roe, lineVitesse_Roe, lineTemperature_Roe) + + +def EulerSystemRoe(ntmax, tmax, cfl, a, b, nbCells, output_freq, meshName, state_law, isImplicit): + state_law_parameters(state_law) + dim = 1 + nbComp = 3 + dt = 0. + time = 0. + it = 0 + isStationary = False + dx = (b - a) / nx + dt = cfl * dx / c0 + #dt = 5*10**(-6) + nbVoisinsMax = 2 + + # iteration vectors + Un_Roe = cdmath.Vector(nbCells * (nbComp)) + dUn_Roe = cdmath.Vector(nbCells * (nbComp)) + dUk_Roe = cdmath.Vector(nbCells * (nbComp)) + Rhs_Roe = cdmath.Vector(nbCells * (nbComp)) + + # Initial conditions + print("Construction of the initial condition …") + + rho_field_Roe, q_field_Roe, rho_E_field_Roe, p_field_Roe, v_field_Roe, T_field_Roe = initial_conditions_Riemann_problem(a, b, nx) + h_field_Roe = (rho_E_field_Roe + p_field_Roe) / rho_field_Roe - 0.5 * (q_field_Roe / rho_field_Roe) **2 + p_field_Roe = p_field_Roe * 10 ** (-5) + + + for k in range(nbCells): + Un_Roe[k * nbComp + 0] = rho_field_Roe[k] + Un_Roe[k * nbComp + 1] = q_field_Roe[k] + Un_Roe[k * nbComp + 2] = rho_E_field_Roe[k] + + divMat_Roe = cdmath.SparseMatrixPetsc(nbCells * nbComp, nbCells * nbComp, (nbVoisinsMax + 1) * nbComp) + + # Picture settings + fig, lineDensity_Roe, lineMomentum_Roe, lineRhoE_Roe, linePressure_Roe, lineVitesse_Roe, lineTemperature_Roe = SetPicture( rho_field_Roe, q_field_Roe, h_field_Roe, p_field_Roe, v_field_Roe, T_field_Roe, dx) + + plt.savefig("EulerComplet_HeatedChannel_" + str(dim) + "D_Roe" + meshName + "0" + ".png") + iterGMRESMax = 50 + newton_max = 100 + + print("Starting computation of the non linear Euler non isentropic system with Roe scheme …") + # STARTING TIME LOOP + while (it < ntmax and time <= tmax and not isStationary): + dUn_Roe = Un_Roe.deepCopy() + Uk_Roe = Un_Roe.deepCopy() + residu_Roe = 1. + + k_Roe = 0 + while (k_Roe < newton_max and residu_Roe > precision): + # STARTING NEWTON LOOP + divMat_Roe.zeroEntries() #sets the matrix coefficients to zero + for j in range(nbCells): + + # traitements des bords + if (j == 0): + FillEdges(j, Uk_Roe, nbComp, divMat_Roe, Rhs_Roe, Un_Roe, dt, dx, isImplicit) + elif (j == nbCells - 1): + FillEdges(j, Uk_Roe, nbComp, divMat_Roe, Rhs_Roe, Un_Roe, dt, dx, isImplicit) + + # traitement des cellules internes + else: + FillInnerCell(j, Uk_Roe, nbComp, divMat_Roe, Rhs_Roe, Un_Roe, dt, dx, isImplicit) + + Rhs_Roe[j * nbComp + 2] += phi*dt + + if(isImplicit): + #solving the linear system divMat * dUk = Rhs + divMat_Roe.diagonalShift(1.) + LS_Roe = cdmath.LinearSolver(divMat_Roe, Rhs_Roe, iterGMRESMax, precision, "GMRES", "LU") + dUk_Roe = LS_Roe.solve() + vector_residu_Roe = dUk_Roe.maxVector(nbComp) + residu_Roe = max(abs(vector_residu_Roe[0])/rho0, abs(vector_residu_Roe[1])/(rho0*v_0), abs(vector_residu_Roe[2])/rhoE0 ) + else: + dUk_Roe=Rhs_Roe - divMat_Roe*Un_Roe + residu_Roe = 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_Roe, " : ", residu_Roe) + print("Linear system converged in ", LS_Roe.getNumberOfIter(), " GMRES iterations") + + #updates for Newton loop + Uk_Roe += dUk_Roe + k_Roe = k_Roe + 1 + if (isImplicit and not LS_Roe.getStatus()): + print("Linear system did not converge ", LS_Roe.getNumberOfIter(), " GMRES iterations") + raise ValueError("No convergence of the linear system") + + if k_Roe == newton_max: + raise ValueError("No convergence of Newton Roe Scheme") + + #updating fields + Un_Roe = Uk_Roe.deepCopy() + dUn_Roe -= Un_Roe + + #Testing stationarity + residu_stat = dUn_Roe.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])/rho0, abs(residu_stat[1])/(rho0*v_0), abs(residu_stat[2])/rhoE0 )) + + if ( it>1 and abs(residu_stat[0])/rho0= ntmax or isStationary or time >= tmax): + + print("-- Time step : " + str(it) + ", Time: " + str(time) + ", dt: " + str(dt)) + + print("Temperature gain between inlet and outlet is ", T_field_Roe[nbCells-1]-T_field_Roe[0],"\n") + + plt.savefig("EulerComplet_HeatedChannel_" + str(dim) + "D_Roe" + meshName + "Implicit"+str(isImplicit)+ str(it) + '_time' + str(time) + ".png") + + print("-- Iter: " + str(it) + ", Time: " + str(time) + ", dt: " + str(dt)+"\n") + + if (it >= ntmax): + print("Maximum number of time steps ntmax= ", ntmax, " reached") + return + + elif (isStationary): + print("Stationary regime reached at time step ", it, ", t= ", time) + print("------------------------------------------------------------------------------------") + np.savetxt("EulerComplet_HeatedChannel_" + str(dim) + "D_Roe" + meshName + "_rho_Stat.txt", rho_field_Roe, delimiter="\n") + np.savetxt("EulerComplet_HeatedChannel_" + str(dim) + "D_Roe" + meshName + "_q_Stat.txt", q_field_Roe, delimiter="\n") + np.savetxt("EulerComplet_HeatedChannel_" + str(dim) + "D_Roe" + meshName + "_rhoE_Stat.txt", rho_E_field_Roe, delimiter="\n") + np.savetxt("EulerComplet_HeatedChannel_" + str(dim) + "D_Roe" + meshName + "_p_Stat.txt", p_field_Roe, delimiter="\n") + plt.savefig("EulerComplet_HeatedChannel_" + str(dim) + "D_Roe" + meshName + "Implicit"+str(isImplicit)+"_Stat.png") + return + else: + print("Maximum time Tmax= ", tmax, " reached") + return + + +def solve(a, b, nx, meshName, meshType, cfl, state_law, isImplicit): + print("Simulation of a heated channel in dimension 1 on " + str(nx) + " cells") + print("State Law Stiffened Gaz, " + state_law) + print("Initial data : ", "constant fields") + print("Boundary conditions : ", "Inlet (Left), Outlet (Right)") + print("Mesh name : ", meshName, ", ", nx, " cells") + # Problem data + tmax = 10. + ntmax = 100000 + output_freq = 1000 + EulerSystemRoe(ntmax, tmax, cfl, a, b, nx, output_freq, meshName, state_law, isImplicit) + return + +def FillMatrixFromEdges(j, Uk, nbComp, divMat, dt, dx): + + 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) + divMat.addValue(j * nbComp, (j + 1) * nbComp, Am) + divMat.addValue(j * nbComp, j * nbComp, Am * (-1.)) + + p_inlet = rho_to_p_StiffenedGaz(Uk[j * nbComp + 0], Uk[j * nbComp + 1], Uk[j * nbComp + 2])# We take p from inside the domain + rho_l=p_to_rho_StiffenedGaz(p_inlet, T_inlet) # rho is computed from the temperature BC and the inner pressure + #rho_l = Uk[j * nbComp + 0] # We take rho from inside the domain + q_l = rho_l * v_inlet # q is imposed by the boundary condition v_inlet + rho_E_l = T_to_rhoE_StiffenedGaz(T_inlet, rho_l, q_l) #rhoE is obtained using the two boundary conditions v_inlet and e_inlet + 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) + divMat.addValue(j * nbComp, j * nbComp, Ap) + + elif (j == nx - 1): + rho_l = Uk[j * nbComp + 0] + q_l = Uk[j * nbComp + 1] + rho_E_l = Uk[j * nbComp + 2] + rho_r = Uk[(j ) * nbComp + 0] # We take rho inside the domain + q_r = Uk[(j ) * nbComp + 1] # We take q from inside the domain + rho_E_r = (p_outlet+gamma*p0)/(gamma-1) + 0.5*q_r**2/rho_r # rhoE is obtained using the boundary condition p_outlet + + Am = jacobianMatricesm(dt / dx, rho_l, q_l, rho_E_l, rho_r, q_r, rho_E_r) + divMat.addValue(j * nbComp, j * nbComp, Am * (-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) + divMat.addValue(j * nbComp, j * nbComp, Ap) + divMat.addValue(j * nbComp, (j - 1) * nbComp, Ap * (-1.)) + + +def FillMatrixFromInnerCell(j, Uk, nbComp, divMat, dt, dx): + + 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) + + 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) + + 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.)) + +def FillRHSFromEdges(j, Uk, nbComp, Rhs, Un, dt, dx, isImplicit): + 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) + + p_inlet = rho_to_p_StiffenedGaz(Uk[j * nbComp + 0], Uk[j * nbComp + 1], Uk[j * nbComp + 2])# We take p from inside the domain + rho_l=p_to_rho_StiffenedGaz(p_inlet, T_inlet) # rho is computed from the temperature BC and the inner pressure + #rho_l = Uk[j * nbComp + 0] # We take rho from inside the domain + q_l = rho_l * v_inlet # q is imposed by the boundary condition v_inlet + rho_E_l = T_to_rhoE_StiffenedGaz(T_inlet, rho_l, q_l) #rhoE is obtained using the two boundary conditions v_inlet and e_inlet + 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) + + 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 + + dUi2[0] = rho_l - Uk[(j ) * nbComp + 0] + dUi2[1] = q_l - Uk[(j ) * nbComp + 1] + dUi2[2] = rho_E_l - Uk[(j ) * nbComp + 2] + temp2 = Ap * dUi2 + 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) + + rho_l = Uk[j * nbComp + 0] + q_l = Uk[j * nbComp + 1] + rho_E_l = Uk[j * nbComp + 2] + rho_r = Uk[(j ) * nbComp + 0] # We take rho inside the domain + q_r = Uk[(j ) * nbComp + 1] # We take q from inside the domain + rho_E_r = (p_outlet+gamma*p0)/(gamma-1) + 0.5*q_r**2/rho_r # rhoE is obtained using the boundary condition p_outlet + + Am = jacobianMatricesm(dt / dx, rho_l, q_l, rho_E_l, rho_r, q_r, rho_E_r) + + if(isImplicit): + dUi1[0] = rho_r - Uk[j * nbComp + 0] + dUi1[1] = q_r - Uk[j * nbComp + 1] + dUi1[2] = rho_E_r - Uk[j * nbComp + 2] + temp1 = Am * dUi1 + + 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): + 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 FillRHSFromInnerCell(j, Uk, nbComp, Rhs, Un, dt, dx, isImplicit): + + 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) + + 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) + + if(isImplicit):#Contribution to the right hand side if te scheme is implicit + 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 computeSystemMatrix(a,b,nx, cfl, Uk, isImplicit): + dim = 1 + nbComp = 3 + dx = (b - a) / nx + dt = cfl * dx / c0 + nbVoisinsMax = 2 + + nbCells = nx + divMat = cdmath.SparseMatrixPetsc(nbCells * nbComp, nbCells * nbComp, (nbVoisinsMax + 1) * nbComp) + + divMat.zeroEntries() #sets the matrix coefficients to zero + for j in range(nbCells): + + # traitements des bords + if (j == 0): + FillMatrixFromEdges(j, Uk, nbComp, divMat, dt, dx) + elif (j == nbCells - 1): + FillMatrixFromEdges(j, Uk, nbComp, divMat, dt, dx) + # traitement des cellules internes + else: + FillMatrixFromInnerCell(j, Uk, nbComp, divMat, dt, dx) + + if(isImplicit): + divMat.diagonalShift(1.) # add one on the diagonal + else: + divMat*=-1. + divMat.diagonalShift(1.) # add one on the diagonal + + return divMat + +def computeRHSVector(a,b,nx, cfl, Uk, Un, isImplicit): + dim = 1 + nbComp = 3 + dx = (b - a) / nx + dt = cfl * dx / c0 + nbVoisinsMax = 2 + + nbCells = nx + Rhs = cdmath.Vector(nbCells * (nbComp)) + + for j in range(nbCells): + + # traitements des bords + if (j == 0): + FillRHSFromEdges(j, Uk, nbComp, Rhs, Un, dt, dx, isImplicit) + elif (j == nbCells - 1): + FillRHSFromEdges(j, Uk, nbComp, Rhs, Un, dt, dx, isImplicit) + # traitement des cellules internes + else: + FillRHSFromInnerCell(j, Uk, nbComp, Rhs, Un, dt, dx, isImplicit) + + return Rhs + + +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) + state_law = "Hermite575K" + state_law_parameters(state_law) + 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 function computeSystemMatrix \n") + cfl = 0.5 + computeSystemMatrix(a, b, nx, cfl, U,True) #Implicit matrix + computeSystemMatrix(a, b, nx, cfl, U,False) #Explicit matrix + + print("\n Testing function computeRHSVector \n") + cfl = 0.5 + computeRHSVector(a, b, nx, cfl, U, U,True) #Implicit RHS + computeRHSVector(a, b, nx, cfl, U, U,False) #Explicit RHS + + print("\n Testing function solve (Implicit scheme) \n") + isImplicit=True + cfl = 1000. + solve(a, b, nx, "RegularSquares", "", cfl, state_law, isImplicit) + + print("\n Testing function solve (Explicit scheme) \n") + isImplicit=False + cfl = 0.5 + solve(a, b, nx, "RegularSquares", "", cfl, state_law, isImplicit) + diff --git a/CDMATH/tests/examples/EulerEquations/FullEulerSystem/EulerSystem1D_RiemannProblem/CMakeLists.txt b/CDMATH/tests/examples/EulerEquations/FullEulerSystem/EulerSystem1D_RiemannProblem/CMakeLists.txt deleted file mode 100755 index 895c803..0000000 --- a/CDMATH/tests/examples/EulerEquations/FullEulerSystem/EulerSystem1D_RiemannProblem/CMakeLists.txt +++ /dev/null @@ -1,8 +0,0 @@ - -if (CDMATH_WITH_PYTHON AND CDMATH_WITH_PETSC AND CDMATH_WITH_POSTPRO) - - ADD_TEST(ExampleFullEulerSystem_1DRiemammProblem_Roe ${PYTHON_EXECUTABLE} ${CMAKE_CURRENT_SOURCE_DIR}/Euler_complet_RP.py ) - -endif (CDMATH_WITH_PYTHON AND CDMATH_WITH_PETSC AND CDMATH_WITH_POSTPRO) - - diff --git a/CDMATH/tests/examples/EulerEquations/FullEulerSystem/EulerSystem1D_RiemannProblem/Euler_complet_RP.py b/CDMATH/tests/examples/EulerEquations/FullEulerSystem/EulerSystem1D_RiemannProblem/Euler_complet_RP.py deleted file mode 100644 index 8785225..0000000 --- a/CDMATH/tests/examples/EulerEquations/FullEulerSystem/EulerSystem1D_RiemannProblem/Euler_complet_RP.py +++ /dev/null @@ -1,642 +0,0 @@ -#!/usr/bin/env python3 -# -*-coding:utf-8 -* - -""" -Created on Mon Aug 30 2021 -@author: Michael NDJINGA, Katia Ait Ameur, Coraline Mounier - -Euler system without source term in one dimension on regular domain [a,b] -Riemann problemn with ghost cell boundary condition -Left and right: Neumann boundary condition -Roe scheme -Regular square mesh - -State law Stiffened gaz : p = (gamma - 1) * rho * (e - q) - gamma * p0 -4 choices of parameters gamma and p0 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 matplotlib.animation as manimation -import sys -from math import sqrt, atan, pi -from numpy import sign - - -## Numerical parameter -precision = 1e-5 - -#state law parameter : can be 'Lagrange', 'Hermite590K', 'Hermite617K', or 'FLICA' -state_law = "Hermite590K" - -#indicates with test case is simulated to compare with FLICA5 results -#num_test = 0 means there are no FLICA5 results given here -num_test = 0 - -#def state_law_parameters(state_law): -#state law Stiffened Gaz : p = (gamma - 1) * rho * e - gamma * p0 -global gamma -global p0 -global q -global c0 -global cp -global h_sat -global T_sat - -if state_law == "Lagrange": - # reference values for Lagrange interpolation - p_ref = 155. * 10**5 #Reference pressure in a REP 900 nuclear power plant - p1 = 153. * 10**5 # value of pressure at inlet of a 900 MWe PWR vessel - rho_ref = 594.38 #density of water at saturation temperature of 617.94K and 155 bars - rho1 = 742.36 # value of density at inlet of a 900 MWe PWR vessel (T1 = 565K) - e_ref = 1603.8 * 10**3 #internal energy of water at saturation temperature of 617.94K and 155 bars - e1 = 1273.6 * 10**3 # value of internal energy at inlet of a 900 MWe PWR vessel - - gamma = (p1 - p_ref) / (rho1 * e1 - rho_ref *e_ref) + 1. - p0 = - 1. / gamma * ( - (gamma - 1) * rho_ref * e_ref + p_ref) - q=0. - c0 = sqrt((gamma - 1) * (e_ref + p_ref / rho_ref)) - - cp = 8950. - h_sat = 1.63 * 10 ** 6 # saturation enthalpy of water at 155 bars - T_sat = 617.94 - -elif state_law == "Hermite617K": - # reference values for Hermite interpolation - p_ref = 155. * 10**5 #Reference pressure in a REP 900 nuclear power plant - T_ref = 617.94 #Reference temperature for interpolation at 617.94K - rho_ref = 594.38 #density of water at saturation temperature of 617.94K and 155 bars - e_ref = 1603.8 * 10**3 #internal energy of water at saturation temperature of 617.94K and 155 bars - h_ref = e_ref + p_ref / rho_ref - c_ref = 621.43 #sound speed for water at 155 bars and 617.94K - - gamma = 1. + c_ref * c_ref / (e_ref + p_ref / rho_ref) # From the sound speed formula - p0 = 1. / gamma * ( (gamma - 1) * rho_ref * e_ref - p_ref) - q=0. - c0 = sqrt((gamma - 1) * (e_ref + p_ref / rho_ref)) - - cp = 8950. # value at 155 bar and 617.94K - h_sat = 1.63 * 10 ** 6 # saturation enthalpy of water at 155 bars - T_sat = 617.94 - -elif state_law == 'Hermite590K': - # reference values for Hermite interpolation - p_ref = 155. * 10**5 #Reference pressure in a REP 900 nuclear power plant - T_ref = 590. #Reference temperature for interpolation at 590K - rho_ref = 688.3 #density of water at 590K and 155 bars - e_ref = 1411.4 * 10**3 #internal energy of water at 590K and 155 bars - h_ref = e_ref + p_ref / rho_ref - c_ref = 866.29 #sound speed for water at 155 bars and 590K - - gamma = 1. + c_ref * c_ref / (e_ref + p_ref / rho_ref) # From the sound speed formula - p0 = 1. / gamma * ( (gamma - 1) * rho_ref * e_ref - p_ref) - q=0. - c0 = sqrt((gamma - 1) * (e_ref + p_ref / rho_ref)) - - cp = 5996.8 # value at 155 bar and 590K - h_sat = 1433.9 * 10 ** 3 # saturation enthalpy of water at 155 bars - T_sat = 590. - -elif state_law == 'Hermite575K': - # reference values for Hermite interpolation - p_ref = 155 * 10**5 #Reference pressure in a REP 900 nuclear power plant - T_ref = 575 #Reference temperature at inlet in a REP 900 nuclear power plant - #Remaining values determined using iapws python package - rho_ref = 722.66 #density of water at 575K and 155 bars - e_ref = 1326552.66 #internal energy of water at 575K and 155 bars - h_ref = e_ref + p_ref / rho_ref - c_ref = 959.28 #sound speed for water at 155 bars and 575K - - gamma = 1 + c_ref * c_ref / (e_ref + p_ref / rho_ref) # From the sound speed formula - p0 = 1 / gamma * ( (gamma - 1) * rho_ref * e_ref - p_ref) - q=0. - c0 = sqrt((gamma - 1) * (e_ref + p_ref / rho_ref))#This is actually c_ref - - cp = 5504.05 # value at 155 bar and 590K - h_sat = h_ref # saturation enthalpy of water at 155 bars - T_sat = T_ref -else: - raise ValueError("Incorrect value for parameter state_law") - - -#initial parameters for Riemann problem (p in Pa, v in m/s, T in K) -p_L = 155. * 10**5 -p_R = 150. * 10**5 -v_L = 0. -v_R = 0. -h_L = 1.4963*10**6 -h_R = 1.4963*10**6 - -T_L = (h_L - h_sat ) / cp + T_sat -T_R = (h_R - h_sat ) / cp + T_sat - -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) - - p_initial = np.array([ (xi < (a + b) / 2) * p_L + (xi >= (a + b) / 2) * p_R for xi in x]) - v_initial = np.array([ (xi < (a + b) / 2) * v_L + (xi >= (a + b) / 2) * v_R for xi in x]) - T_initial = np.array([ (xi < (a + b) / 2) * T_L + (xi >= (a + b) / 2) * T_R for xi in x]) - - rho_initial = p_to_rho_StiffenedGaz(p_initial, T_initial) - q_initial = rho_initial * v_initial - rho_E_initial = T_to_rhoE_StiffenedGaz(T_initial, rho_initial, q_initial) - - return rho_initial, q_initial, rho_E_initial, p_initial, v_initial, T_initial - -def rho_to_p_StiffenedGaz(rho_field, q_field, rho_E_field): - p_field = (gamma - 1) * ( rho_E_field - 1. / 2 * q_field ** 2 / rho_field - rho_field * q) - gamma * p0 - return p_field - - -def p_to_rho_StiffenedGaz(p_field, T_field): - rho_field = (p_field + p0) * gamma / (gamma - 1) * 1 / (h_sat + cp * (T_field - T_sat) - q) - return rho_field - - -def T_to_rhoE_StiffenedGaz(T_field, rho_field, q_field): - rho_E_field = 1 / 2 * (q_field) ** 2 / rho_field + p0 + rho_field / gamma * (h_sat + cp * (T_field- T_sat) + (gamma - 1) * q) - return rho_E_field - - -def rhoE_to_T_StiffenedGaz(rho_field, q_field, rho_E_field): - T_field = T_sat + 1 / cp * (gamma * (rho_E_field / rho_field - 1 / 2 * (q_field / rho_field) ** 2) - gamma * p0 / rho_field - (gamma - 1) * q - h_sat) - return T_field - -def dp_drho_e_const_StiffenedGaz( e ): - return (gamma-1)*(e-q) - -def dp_de_rho_const_StiffenedGaz( rho ): - return (gamma-1)*rho - -def sound_speed_StiffenedGaz( h ): - return np.sqrt((gamma-1)*(h-q)) - -def rho_h_to_p_StiffenedGaz( rho, h ): - return (gamma - 1) * rho * ( h - q ) / gamma - p0 - -def MatRoe_StiffenedGaz( rho_l, q_l, rho_E_l, rho_r, q_r, rho_E_r): - 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) - p_r = rho_to_p_StiffenedGaz(rho_r, q_r, rho_E_r) - 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. ) - e = H - p / rho - 1./2 * u**2 - dp_drho = dp_drho_e_const_StiffenedGaz( e ) - dp_de = dp_de_rho_const_StiffenedGaz( rho ) - - 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 Droe_StiffenedGaz( rho_l, q_l, rho_E_l, rho_r, q_r, rho_E_r): - 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) - p_r = rho_to_p_StiffenedGaz(rho_r, q_r, rho_E_r) - 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. ) - - 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): - - 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) - - Droe = Droe_StiffenedGaz( rho_l, q_l, rho_E_l, rho_r, q_r, rho_E_r) - return (RoeMat - Droe) * coeff * 0.5 - - -def jacobianMatricesp(coeff, rho_l, q_l, rho_E_l, rho_r, q_r, rho_E_r): - 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) - Droe = Droe_StiffenedGaz( rho_l, q_l, rho_E_l, rho_r, q_r, rho_E_r) - - return (RoeMat + Droe) * coeff * 0.5 - - -def FillEdges(j, Uk, nbComp, divMat, Rhs, Un, dt, dx, isImplicit): - dUi = 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) - divMat.addValue(j * nbComp, (j + 1) * nbComp, Am) - divMat.addValue(j * nbComp, j * nbComp, Am * (-1.)) - - dUi[0] = Uk[(j + 1) * nbComp + 0] - Uk[j * nbComp + 0] - dUi[1] = Uk[(j + 1) * nbComp + 1] - Uk[j * nbComp + 1] - dUi[2] = Uk[(j + 1) * nbComp + 2] - Uk[j * nbComp + 2] - temp = Am * dUi - - if(isImplicit): - 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]) - Rhs[j * nbComp + 2] = -temp[2] - (Uk[j * nbComp + 2] - Un[j * nbComp + 2]) - - 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) - divMat.addValue(j * nbComp, j * nbComp, Ap) - divMat.addValue(j * nbComp, (j - 1) * nbComp, Ap * (-1.)) - - dUi[0] = Uk[(j - 1) * nbComp + 0] - Uk[j * nbComp + 0] - dUi[1] = Uk[(j - 1) * nbComp + 1] - Uk[j * nbComp + 1] - dUi[2] = Uk[(j - 1) * nbComp + 2] - Uk[j * nbComp + 2] - - temp = Ap * dUi - - if(isImplicit): - 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]) - Rhs[j * nbComp + 2] = temp[2] - (Uk[j * nbComp + 2] - Un[j * nbComp + 2]) - -def FillInnerCell(j, Uk, nbComp, divMat, Rhs, Un, dt, dx, isImplicit): - - 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) - - 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) - - 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(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 - - if(isImplicit): - 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]) - -def SetPicture(rho_field_Roe, q_field_Roe, h_field_Roe, p_field_Roe, v_field_Roe, T_field_Roe, dx): - max_initial_rho = max(rho_field_Roe) - min_initial_rho = min(rho_field_Roe) - max_initial_q = max(q_field_Roe) - min_initial_q = min(q_field_Roe) - min_initial_h = min(h_field_Roe) - max_initial_h = max(h_field_Roe) - max_initial_p = max(p_field_Roe) - min_initial_p = min(p_field_Roe) - max_initial_v = max(v_field_Roe) - min_initial_v = min(v_field_Roe) - max_initial_T = max(T_field_Roe) - min_initial_T = min(T_field_Roe) - - fig, ([axDensity, axMomentum, axRhoE],[axPressure, axVitesse, axTemperature]) = plt.subplots(2, 3,sharex=True, figsize=(14,10)) - plt.gcf().subplots_adjust(wspace = 0.5) - - lineDensity_Roe, = axDensity.plot([a+0.5*dx + i*dx for i in range(nx)], rho_field_Roe, label='Roe') - axDensity.set(xlabel='x (m)', ylabel='Densité (kg/m3)') - axDensity.set_xlim(a,b) - #axDensity.set_ylim(min_initial_rho - 0.1 * (max_initial_rho - min_initial_rho), 700.) - axDensity.set_ylim(657, 660.5) - axDensity.legend() - - lineMomentum_Roe, = axMomentum.plot([a+0.5*dx + i*dx for i in range(nx)], q_field_Roe, label='Roe') - axMomentum.set(xlabel='x (m)', ylabel='Momentum (kg/(m2.s))') - axMomentum.set_xlim(a,b) - #axMomentum.set_ylim(min_initial_q - 0.1*(max_initial_q-min_initial_q), max_initial_q * 2.5) - axMomentum.set_ylim(-50, 500) - axMomentum.legend() - - lineRhoE_Roe, = axRhoE.plot([a+0.5*dx + i*dx for i in range(nx)], h_field_Roe, label='Roe') - axRhoE.set(xlabel='x (m)', ylabel='h (J/m3)') - axRhoE.set_xlim(a,b) - #axRhoE.set_ylim(min_initial_h - 0.05*abs(min_initial_h), max_initial_h + 0.5*(max_initial_h-min_initial_h)) - axRhoE.set_ylim(1.495 * 10**6, 1.5*10**6) - axRhoE.legend() - - linePressure_Roe, = axPressure.plot([a+0.5*dx + i*dx for i in range(nx)], p_field_Roe, label='Roe') - axPressure.set(xlabel='x (m)', ylabel='Pression (bar)') - axPressure.set_xlim(a,b) - #axPressure.set_ylim(min_initial_p - 0.05*abs(min_initial_p), max_initial_p + 0.5*(max_initial_p-min_initial_p)) - axPressure.set_ylim(149, 156) - axPressure.ticklabel_format(axis='y', style='sci', scilimits=(0,0)) - axPressure.legend() - - lineVitesse_Roe, = axVitesse.plot([a+0.5*dx + i*dx for i in range(nx)], v_field_Roe, label='Roe') - axVitesse.set(xlabel='x (m)', ylabel='Vitesse (m/s)') - axVitesse.set_xlim(a,b) - #axVitesse.set_ylim(min_initial_v - 0.05*abs(min_initial_v), 15) - axVitesse.set_ylim(-0.5, 1) - axVitesse.ticklabel_format(axis='y', style='sci', scilimits=(0,0)) - axVitesse.legend() - - lineTemperature_Roe, = axTemperature.plot([a+0.5*dx + i*dx for i in range(nx)], T_field_Roe, label='Roe') - axTemperature.set(xlabel='x (m)', ylabel='Température (K)') - axTemperature.set_xlim(a,b) - #axTemperature.set_ylim(min_initial_T - 0.005*abs(min_initial_T), 604) - axTemperature.set_ylim(600, 601) - axTemperature.ticklabel_format(axis='y', style='sci', scilimits=(0,0)) - axTemperature.legend() - - return(fig, lineDensity_Roe, lineMomentum_Roe, lineRhoE_Roe, linePressure_Roe, lineVitesse_Roe, lineTemperature_Roe, lineDensity_Roe, lineMomentum_Roe, lineRhoE_Roe, linePressure_Roe, lineVitesse_Roe, lineTemperature_Roe) - - -def EulerSystemRoe(ntmax, tmax, cfl, a, b, nbCells, output_freq, meshName, state_law): - #state_law_parameters(state_law) - dim = 1 - nbComp = 3 - dt = 0. - time = 0. - it = 0 - isStationary = False - isImplicit = False - dx = (b - a) / nx - dt = cfl * dx / c0 - #dt = 5*10**(-6) - nbVoisinsMax = 2 - - # iteration vectors - Un_Roe = cdmath.Vector(nbCells * (nbComp)) - dUn_Roe = cdmath.Vector(nbCells * (nbComp)) - dUk_Roe = cdmath.Vector(nbCells * (nbComp)) - Rhs_Roe = cdmath.Vector(nbCells * (nbComp)) - - # Initial conditions - print("Construction of the initial condition …") - - rho_field_Roe, q_field_Roe, rho_E_field_Roe, p_field_Roe, v_field_Roe, T_field_Roe = initial_conditions_Riemann_problem(a, b, nx) - h_field_Roe = (rho_E_field_Roe + p_field_Roe) / rho_field_Roe - 0.5 * (q_field_Roe / rho_field_Roe) **2 - p_field_Roe = p_field_Roe * 10 ** (-5) - - - for k in range(nbCells): - Un_Roe[k * nbComp + 0] = rho_field_Roe[k] - Un_Roe[k * nbComp + 1] = q_field_Roe[k] - Un_Roe[k * nbComp + 2] = rho_E_field_Roe[k] - - divMat_Roe = cdmath.SparseMatrixPetsc(nbCells * nbComp, nbCells * nbComp, (nbVoisinsMax + 1) * nbComp) - - # Picture settings - fig, lineDensity, lineMomentum, lineRhoE, linePressure, lineVitesse, lineTemperature, lineDensity_Roe, lineMomentum_Roe, lineRhoE_Roe, linePressure_Roe, lineVitesse_Roe, lineTemperature_Roe = SetPicture( rho_field_Roe, q_field_Roe, h_field_Roe, p_field_Roe, v_field_Roe, T_field_Roe, dx) - - # Video settings - FFMpegWriter = manimation.writers['ffmpeg'] - metadata = dict(title="Roe scheme for the 1D Euler system", artist="CEA Saclay", comment="Shock tube") - writer = FFMpegWriter(fps=10, metadata=metadata, codec='h264') - with writer.saving(fig, "1DEuler_complet_RP_Roe" + ".mp4", ntmax): - writer.grab_frame() - plt.savefig("EulerComplet_RP_" + str(dim) + "D_Roe" + meshName + "0" + ".png") - np.savetxt( "EulerComplet_RP_" + str(dim) + "D_Roe" + meshName + "_rho" + "0" + ".txt", rho_field_Roe, delimiter="\n") - np.savetxt( "EulerComplet_RP_" + str(dim) + "D_Roe" + meshName + "_q" + "0" + ".txt", q_field_Roe, delimiter="\n") - iterGMRESMax = 50 - newton_max = 100 - - print("Starting computation of the non linear Euler non isentropic system with Roe scheme …") - # STARTING TIME LOOP - while (it < ntmax and time <= tmax and not isStationary): - dUn_Roe = Un_Roe.deepCopy() - Uk_Roe = Un_Roe.deepCopy() - residu_Roe = 1. - - k_Roe = 0 - while (k_Roe < newton_max and residu_Roe > precision): - # STARTING NEWTON LOOP - divMat_Roe.zeroEntries() #sets the matrix coefficients to zero - for j in range(nbCells): - - # traitements des bords - if (j == 0): - FillEdges(j, Uk_Roe, nbComp, divMat_Roe, Rhs_Roe, Un_Roe, dt, dx, isImplicit) - elif (j == nbCells - 1): - FillEdges(j, Uk_Roe, nbComp, divMat_Roe, Rhs_Roe, Un_Roe, dt, dx, isImplicit) - - # traitement des cellules internes - else: - FillInnerCell(j, Uk_Roe, nbComp, divMat_Roe, Rhs_Roe, Un_Roe, dt, dx, isImplicit) - - #solving the linear system divMat * dUk = Rhs - - if(isImplicit): - divMat_Roe.diagonalShift(1.) - LS_Roe = cdmath.LinearSolver(divMat_Roe, Rhs_Roe, iterGMRESMax, precision, "GMRES", "LU") - dUk_Roe = LS_Roe.solve() - residu_Roe = dUk_Roe.norm() - else: - dUk_Roe=Rhs_Roe - divMat_Roe*Un_Roe - residu_Roe = 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 : ", residu_Roe) - print("Linear system converged in ", LS_Roe.getNumberOfIter(), " GMRES iterations") - - #updates for Newton loop - Uk_Roe += dUk_Roe - k_Roe = k_Roe + 1 - if (isImplicit and not LS_Roe.getStatus()): - print("Linear system did not converge ", LS_Roe.getNumberOfIter(), " GMRES iterations") - raise ValueError("No convergence of the linear system") - - if k_Roe == newton_max: - raise ValueError("No convergence of Newton Roe Scheme") - - #updating fields - Un_Roe = Uk_Roe.deepCopy() - dUn_Roe -= Un_Roe - if (dUn_Roe.norm()= ntmax or isStationary or time >= tmax): - - print("-- Time step : " + str(it) + ", Time: " + str(time) + ", dt: " + str(dt)) - - plt.savefig("EulerComplet_RP_" + str(dim) + "D_Roe" + meshName + str(it) + '_time' + str(time) + ".png") - - print("-- Iter: " + str(it) + ", Time: " + str(time) + ", dt: " + str(dt)) - if (it >= ntmax): - print("Maximum number of time steps ntmax= ", ntmax, " reached") - return - - elif (isStationary): - print("Stationary regime reached at time step ", it, ", t= ", time) - print("------------------------------------------------------------------------------------") - np.savetxt("EulerComplet_RP_" + str(dim) + "D_Roe" + meshName + "_rho_Stat.txt", rho_field_Roe, delimiter="\n") - np.savetxt("EulerComplet_RP_" + str(dim) + "D_Roe" + meshName + "_q_Stat.txt", q_field_Roe, delimiter="\n") - np.savetxt("EulerComplet_RP_" + str(dim) + "D_Roe" + meshName + "_rhoE_Stat.txt", rho_E_field_Roe, delimiter="\n") - np.savetxt("EulerComplet_RP_" + str(dim) + "D_Roe" + meshName + "_p_Stat.txt", p_field_Roe, delimiter="\n") - plt.savefig("EulerComplet_RP_" + str(dim) + "D_Roe" + meshName + "_Stat.png") - return - else: - print("Maximum time Tmax= ", tmax, " reached") - return - - -def solve(a, b, nx, meshName, meshType, cfl, state_law): - print("Resolution of the Euler System in dimension 1 on " + str(nx) + " cells") - print("State Law Stiffened Gaz, " + state_law) - print("Initial data : ", "Riemann problem") - print("Boundary conditions : ", "Neumann") - print("Mesh name : ", meshName, ", ", nx, " cells") - # Problem data - tmax = 10. - ntmax = 25 - output_freq = 1 - EulerSystemRoe(ntmax, tmax, cfl, a, b, nx, output_freq, meshName, state_law) - return - - -if __name__ == """__main__""": - a = 0. - b = 1. - nx = 50 - cfl = 0.5 - #state_law = "Hermite590K" - #state_law_parameters(state_law) - solve(a, b, nx, "RegularSquares", "", cfl, state_law) diff --git a/CDMATH/tests/examples/EulerEquations/FullEulerSystem/EulerSystem1D_RiemannProblem_MAC/CMakeLists.txt b/CDMATH/tests/examples/EulerEquations/FullEulerSystem/EulerSystem1D_RiemannProblem_MAC/CMakeLists.txt new file mode 100755 index 0000000..b832057 --- /dev/null +++ b/CDMATH/tests/examples/EulerEquations/FullEulerSystem/EulerSystem1D_RiemannProblem_MAC/CMakeLists.txt @@ -0,0 +1,9 @@ + +if (CDMATH_WITH_PYTHON AND CDMATH_WITH_PETSC AND CDMATH_WITH_POSTPRO) + file(COPY TTC4.dat DESTINATION ${CMAKE_CURRENT_BINARY_DIR}) + + ADD_TEST(ExampleFullEulerSystem_1DRiemammProblem_MAC ${PYTHON_EXECUTABLE} ${CMAKE_CURRENT_SOURCE_DIR}/Euler_complet_Toro4_convergence.py ) + +endif (CDMATH_WITH_PYTHON AND CDMATH_WITH_PETSC AND CDMATH_WITH_POSTPRO) + + diff --git a/CDMATH/tests/examples/EulerEquations/FullEulerSystem/EulerSystem1D_RiemannProblem_MAC/Euler_complet_Toro4_convergence.py b/CDMATH/tests/examples/EulerEquations/FullEulerSystem/EulerSystem1D_RiemannProblem_MAC/Euler_complet_Toro4_convergence.py new file mode 100644 index 0000000..cc02953 --- /dev/null +++ b/CDMATH/tests/examples/EulerEquations/FullEulerSystem/EulerSystem1D_RiemannProblem_MAC/Euler_complet_Toro4_convergence.py @@ -0,0 +1,548 @@ +#!/usr/bin/env python3 +# -*-coding:utf-8 -* + + + +""" +Created on Mon Apr 26 2021 +@author: Michael NDJINGA + +Toro 4 shock tube benchmarck from book 'Riemann Solvers and Numerical Methods for Fluid Dynamics, Third edition, Springer, 2009' (chapter 10) + +Euler system in one dimension on domain [a,b] +Riemann problemn with ghost cell boundary condition +Stag scheme compared with exact solution +Regular mesh + + +""" + +import cdmath +import numpy as np +import matplotlib + +matplotlib.use("Agg") +import matplotlib.pyplot as plt +import os +from math import sqrt +from numpy import sign + +#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 + +pos_disc = 0.4 + +precision = 1.e-5 + +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) ) + +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) + + + 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) + 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 rho_to_p_StiffenedGaz(rho_field, q_field, rho_E_field): + p_field = (gamma - 1.) * ( rho_E_field - 1. / 2 * q_field ** 2 / rho_field - rho_field * q) - gamma * p0 + return p_field + +def p_to_e_StiffenedGaz(p_field, rho_field): + e_field = (p_field + gamma*p0) / (gamma - 1.) / rho_field + q + return e_field + +#Does not involve stiffened gas +def e_to_rhoE_StiffenedGaz(e_field, rho_field, q_field): + rho_E_field = 1 / 2 * (q_field) ** 2 / rho_field + rho_field * e_field + return rho_E_field + +#Does not involve stiffened gas +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 dp_drho_e_const_StiffenedGaz( e ): + return (gamma-1)*(e-q) + +def dp_de_rho_const_StiffenedGaz( rho ): + return (gamma-1)*rho + +def sound_speed_StiffenedGaz( h ): + return np.sqrt((gamma-1)*(h-q)) + +def rho_h_to_p_StiffenedGaz( rho, h ): + return (gamma - 1) * rho * ( h - q ) / gamma - p0 + +def MatRoe_StiffenedGaz( rho_l, q_l, rho_E_l, rho_r, q_r, rho_E_r): + 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) + p_r = rho_to_p_StiffenedGaz(rho_r, q_r, rho_E_r) + 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. ) + e = H - p / rho - 1./2 * u**2 + dp_drho = dp_drho_e_const_StiffenedGaz( e ) + dp_de = dp_de_rho_const_StiffenedGaz( rho ) + + 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): + 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) + p_r = rho_to_p_StiffenedGaz(rho_r, q_r, rho_E_r) + 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. ) + e = H - p / rho - 1./2 * u**2 + dp_drho = dp_drho_e_const_StiffenedGaz( e ) + dp_de = dp_de_rho_const_StiffenedGaz( rho ) + + #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 jacobianMatricesm(coeff, rho_l, q_l, rho_E_l, rho_r, q_r, rho_E_r, scheme): + + 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) + + if scheme == 'Stag': + Dmac = Dmac_StiffenedGaz( rho_l, q_l, rho_E_l, rho_r, q_r, rho_E_r) + 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) + return (RoeMat - Droe) * coeff * 0.5 + + +def jacobianMatricesp(coeff, rho_l, q_l, rho_E_l, rho_r, q_r, rho_E_r, scheme): + 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) + if scheme == 'Stag': + Dmac = Dmac_StiffenedGaz( rho_l, q_l, rho_E_l, rho_r, q_r, rho_E_r) + 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) + return (RoeMat + Droe) * coeff * 0.5 + + +def FillEdges(j, Uk, nbComp, divMat, Rhs, Un, dt, dx, scheme): + 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) + divMat.addValue(j * nbComp, (j + 1) * nbComp, Am) + divMat.addValue(j * nbComp, j * nbComp, Am * (-1.)) + + dUi = cdmath.Vector(3) + dUi[0] = Uk[(j + 1) * nbComp + 0] - Uk[j * nbComp + 0] + dUi[1] = Uk[(j + 1) * nbComp + 1] - Uk[j * nbComp + 1] + dUi[2] = Uk[(j + 1) * nbComp + 2] - Uk[j * nbComp + 2] + temp = Am * 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]) + Rhs[j * nbComp + 2] = -temp[2] - (Uk[j * nbComp + 2] - Un[j * nbComp + 2]) + + elif (j == Un.size()/nbComp - 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) + divMat.addValue(j * nbComp, j * nbComp, Ap) + divMat.addValue(j * nbComp, (j - 1) * nbComp, Ap * (-1.)) + + dUi = cdmath.Vector(3) + dUi[0] = Uk[(j - 1) * nbComp + 0] - Uk[j * nbComp + 0] + dUi[1] = Uk[(j - 1) * nbComp + 1] - Uk[j * nbComp + 1] + dUi[2] = Uk[(j - 1) * nbComp + 2] - Uk[j * nbComp + 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]) + Rhs[j * nbComp + 2] = temp[2] - (Uk[j * nbComp + 2] - Un[j * nbComp + 2]) + +def FillInnerCell(j, Uk, nbComp, divMat, Rhs, Un, dt, dx, scheme): + + 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) + + 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) + + 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(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]) + +def SetPicture(): + 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=p_to_e_StiffenedGaz(p_L, rho_L) + e_R=p_to_e_StiffenedGaz(p_R, rho_R) + 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)) + #fig.suptitle('Staggered scheme') + plt.gcf().subplots_adjust(wspace = 0.5) + + 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.set_ylim(0.125, 1.01) + axDensity.legend() + + 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.set_ylim(0., 1.) + axMomentum.legend() + + axEnthalpie.set(xlabel='x (m)', ylabel='h (J/Kg)') + axEnthalpie.set_xlim(a,b) + axEnthalpie.set_ylim(0.9*min_initial_h , 1.75*max_initial_h ) + #axEnthalpie.set_ylim(2.5, 5.) + axEnthalpie.legend() + + 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.set_ylim(0.1, 1) + axPressure.ticklabel_format(axis='y', style='sci', scilimits=(0,0)) + axPressure.legend() + + 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.set_ylim(0., 1.5) + axVitesse.ticklabel_format(axis='y', style='sci', scilimits=(0,0)) + axVitesse.legend() + + 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.set_ylim(2., 3.5) + axEinterne.ticklabel_format(axis='y', style='sci', scilimits=(0,0)) + axEinterne.legend() + + return(fig) + +def addCurves(rho_field_Stag, q_field_Stag, h_field_Stag, p_field_Stag, v_field_Stag, e_field_Stag, dx, nbCells, fig): + + [axDensity, axMomentum, axEnthalpie,axPressure, axVitesse, axEinterne] = fig.get_axes() + + lineDensity_Stag, = axDensity.plot([a+0.5*dx + i*dx for i in range(nbCells)], rho_field_Stag, label='Stag on '+str(nbCells)+' cells') #new picture for video # Returns a tuple of line objects, thus the comma + + lineMomentum_Stag, = axMomentum.plot([a+0.5*dx + i*dx for i in range(nbCells)], q_field_Stag, label='Stag on '+str(nbCells)+' cells') + + lineEnthalpie_Stag, = axEnthalpie.plot([a+0.5*dx + i*dx for i in range(nbCells)], h_field_Stag, label='Stag on '+str(nbCells)+' cells') + + linePressure_Stag, = axPressure.plot([a+0.5*dx + i*dx for i in range(nbCells)], p_field_Stag, label='Stag on '+str(nbCells)+' cells') + + lineVitesse_Stag, = axVitesse.plot([a+0.5*dx + i*dx for i in range(nbCells)], v_field_Stag, label='Stag on '+str(nbCells)+' cells') + + lineEinterne_Stag, = axEinterne.plot([a+0.5*dx + i*dx for i in range(nbCells)], e_field_Stag, label='Stag on '+str(nbCells)+' cells') + + return + +def addSolutionToro(fig) : + + x_exact, rho_exact, v_exact, p_exact, e_exact, h_exact, Mach_exact, left_or_right = np.loadtxt(os.path.dirname(__file__)+'/TTC'+str(4) + '.dat', unpack = True ) + + [axDensity, axMomentum, axEnthalpie,axPressure, axVitesse, axEinterne] = fig.get_axes() + lineRhoToro, = axDensity.plot(x_exact, rho_exact , label='Solution exacte') + axDensity.legend() + lineVToro, = axVitesse.plot(x_exact, v_exact, label='Solution exacte') + axVitesse.legend() + lineQToro, = axMomentum.plot(x_exact, v_exact*rho_exact, label='Solution exacte') + axMomentum.legend() + lineEToro, = axEinterne.plot(x_exact, e_exact, label='Solution exacte') + axEinterne.legend() + linePToro, = axPressure.plot(x_exact, p_exact, label='Solution exacte') + axPressure.legend() + lineEnthalpieToro, = axEnthalpie.plot(x_exact, h_exact, label='Solution exacte') + axEnthalpie.legend() + + +def EulerSystemStaggered(ntmax, tmax, cfl, a, b, nbCells, output_freq, meshName,fig): + nbComp = 3 + time = 0. + it = 0 + isStationary = False + dx = (b - a) / nbCells + dt = cfl * dx / lam_max + nbVoisinsMax = 2 + + # iteration vectors + Un_Stag = cdmath.Vector(nbCells * (nbComp)) + dUn_Stag = cdmath.Vector(nbCells * (nbComp)) + dUk_Stag = cdmath.Vector(nbCells * (nbComp)) + Rhs_Stag = cdmath.Vector(nbCells * (nbComp)) + + # Initial conditions + print("Construction of the initial condition …") + rho_field_Stag, q_field_Stag, rho_E_field_Stag, p_field_Stag, v_field_Stag, e_field_Stag = initial_conditions_Riemann_problem(a, b, nbCells) + h_field_Stag = (rho_E_field_Stag + p_field_Stag) / rho_field_Stag - 0.5 * (q_field_Stag / rho_field_Stag) **2 + + + for k in range(nbCells): + Un_Stag[k * nbComp + 0] = rho_field_Stag[k] + Un_Stag[k * nbComp + 1] = q_field_Stag[k] + Un_Stag[k * nbComp + 2] = rho_E_field_Stag[k] + + divMat_Stag = 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 staggered scheme …") + # STARTING TIME LOOP + while (it < ntmax and time <= tmax and not isStationary): + dUn_Stag = Un_Stag.deepCopy() + Uk_Stag = Un_Stag.deepCopy() + residu_Stag = 1. + + k_Stag = 0 + while (k_Stag < newton_max and residu_Stag > precision): + # STARTING NEWTON LOOP + #print("Iteration k=", k_Stag, " residu = ",residu_Stag) + divMat_Stag.zeroEntries() #sets the matrix coefficients to zero + for j in range(nbCells): + # traitements des bords + if (j == 0): + FillEdges(j, Uk_Stag, nbComp, divMat_Stag, Rhs_Stag, Un_Stag, dt, dx, 'Stag') + + elif (j == nbCells - 1): + FillEdges(j, Uk_Stag, nbComp, divMat_Stag, Rhs_Stag, Un_Stag, dt, dx, 'Stag') + + # traitement des cellules internes + else: + FillInnerCell(j, Uk_Stag, nbComp, divMat_Stag, Rhs_Stag, Un_Stag, dt, dx, 'Stag') + + #print(divMat_Stag) + divMat_Stag.diagonalShift(1) # only after filling all coefficients + + #solving the linear system divMat * dUk = Rhs + LS_Stag = cdmath.LinearSolver(divMat_Stag, Rhs_Stag, iterGMRESMax, precision, "GMRES", "LU") + dUk_Stag = LS_Stag.solve() + residu_Stag = dUk_Stag.norm() + + #updates for Newton loop + Uk_Stag += dUk_Stag + k_Stag = k_Stag + 1 + if (not LS_Stag.getStatus()): + print("Linear system did not converge ", LS_Stag.getNumberOfIter(), " GMRES iterations") + raise ValueError("No convergence of the linear system") + + if k_Stag == newton_max: + raise ValueError("No convergence of Newton Staggered Scheme") + + + #updating fields + Un_Stag = Uk_Stag.deepCopy() + dUn_Stag -= Un_Stag + + if (dUn_Stag.norm()= ntmax or isStationary or time >= tmax): + + print("-- Time step : " + str(it) + ", Time: " + str(time) + ", dt: " + str(dt)) + print("Linear system converged in ", LS_Stag.getNumberOfIter(), " GMRES iterations") + + plt.savefig("EulerCompletStaggeredToro4" + meshName + str(it) + '_time' + str(time) + ".png") + + print("##################### Last time step: " + str(it) + ", Time: " + str(time) + ", dt: " + str(dt)) + + # Postprocessing + for k in range(nbCells): + rho_field_Stag[k] = Un_Stag[k * nbComp + 0] + q_field_Stag[k] = Un_Stag[k * nbComp + 1] + rho_E_field_Stag[k] = Un_Stag[k * nbComp + 2] + + v_field_Stag = q_field_Stag / rho_field_Stag + p_field_Stag = rho_to_p_StiffenedGaz(rho_field_Stag, q_field_Stag, rho_E_field_Stag) + e_field_Stag = rhoE_to_e_StiffenedGaz(rho_field_Stag, q_field_Stag, rho_E_field_Stag) + h_field_Stag = (rho_E_field_Stag + p_field_Stag) / rho_field_Stag - 0.5 * (q_field_Stag / rho_field_Stag) **2 + + # Picture settings + addCurves(rho_field_Stag, q_field_Stag, h_field_Stag, p_field_Stag, v_field_Stag, e_field_Stag, dx, nbCells,fig) + + np.savetxt("EulerCompletStaggeredToro4" + meshName + "_rho_" + str(nbCells) + "_cells.txt", rho_field_Stag, delimiter="\n") + np.savetxt("EulerCompletStaggeredToro4" + meshName + "_q_" + str(nbCells) + "_cells.txt", q_field_Stag, delimiter="\n") + np.savetxt("EulerCompletStaggeredToro4" + meshName + "_h_" + str(nbCells) + "._cellstxt", h_field_Stag, delimiter="\n") + np.savetxt("EulerCompletStaggeredToro4" + meshName + "_p_" + str(nbCells) + "_cells.txt", p_field_Stag, delimiter="\n") + np.savetxt("EulerCompletStaggeredToro4" + meshName + "_v_" + str(nbCells) + "_cells.txt", v_field_Stag, delimiter="\n") + np.savetxt("EulerCompletStaggeredToro4" + meshName + "_e_" + str(nbCells) + "_cells.txt", e_field_Stag, delimiter="\n") + + if (it >= ntmax): + print("Maximum number of time steps ntmax= ", ntmax, " reached on "+ str(nbCells) + " cells") + + elif (time >= tmax): + print("Maximum time tmax= ", tmax, " reached on "+ str(nbCells) + " cells") + elif (isStationary): + print("Stationary regime reached at time step ", it, ", t= ", time, ", on "+ str(nbCells) + " cells") + else: + raise ValueError("Unexpected end of simulation") + + return + +def solve(a, b, meshName, meshType, cfl): + print("Convergence study for the Euler System in dimension 1") + print("State Law Stiffened Gaz, gamma = ", gamma, " p0= ", p0 ) + print("Initial data : ", "Riemann problem") + print("Boundary conditions : ", "Neumann") + print("Mesh name : ", meshName) + # Problem data + ntmax = 1000000 + output_freq = 100 + fig=SetPicture() + for nx in [20, 50, 100]: + print("Resolution of the Euler System in dimension 1 on " + str(nx) + " cells") + EulerSystemStaggered(ntmax, tmax, cfl, a, b, nx, output_freq, meshName,fig) + + addSolutionToro(fig) + plt.savefig("EulerCompletStaggeredToro4" + meshName + "_convergence.png") + + return + + +if __name__ == """__main__""": + a = 0. + b = 1. + cfl = 0.99 + solve(a, b, "RegularGrid", "", cfl) diff --git a/CDMATH/tests/examples/EulerEquations/FullEulerSystem/EulerSystem1D_RiemannProblem_MAC/TTC4.dat b/CDMATH/tests/examples/EulerEquations/FullEulerSystem/EulerSystem1D_RiemannProblem_MAC/TTC4.dat new file mode 100644 index 0000000..163ea97 --- /dev/null +++ b/CDMATH/tests/examples/EulerEquations/FullEulerSystem/EulerSystem1D_RiemannProblem_MAC/TTC4.dat @@ -0,0 +1,100 @@ +0.0 5.99924 19.5975 460.894 192.063494709 268.888892593 1.88966250666 1.0 +0.01 5.99924 19.5975 460.894 192.063494709 268.888892593 1.88966250666 1.0 +0.02 5.99924 19.5975 460.894 192.063494709 268.888892593 1.88966250666 1.0 +0.03 5.99924 19.5975 460.894 192.063494709 268.888892593 1.88966250666 1.0 +0.04 5.99924 19.5975 460.894 192.063494709 268.888892593 1.88966250666 1.0 +0.05 5.99924 19.5975 460.894 192.063494709 268.888892593 1.88966250666 1.0 +0.06 5.99924 19.5975 460.894 192.063494709 268.888892593 1.88966250666 1.0 +0.07 5.99924 19.5975 460.894 192.063494709 268.888892593 1.88966250666 1.0 +0.08 5.99924 19.5975 460.894 192.063494709 268.888892593 1.88966250666 1.0 +0.09 5.99924 19.5975 460.894 192.063494709 268.888892593 1.88966250666 1.0 +0.1 5.99924 19.5975 460.894 192.063494709 268.888892593 1.88966250666 1.0 +0.11 5.99924 19.5975 460.894 192.063494709 268.888892593 1.88966250666 1.0 +0.12 5.99924 19.5975 460.894 192.063494709 268.888892593 1.88966250666 1.0 +0.13 5.99924 19.5975 460.894 192.063494709 268.888892593 1.88966250666 1.0 +0.14 5.99924 19.5975 460.894 192.063494709 268.888892593 1.88966250666 1.0 +0.15 5.99924 19.5975 460.894 192.063494709 268.888892593 1.88966250666 1.0 +0.16 5.99924 19.5975 460.894 192.063494709 268.888892593 1.88966250666 1.0 +0.17 5.99924 19.5975 460.894 192.063494709 268.888892593 1.88966250666 1.0 +0.18 5.99924 19.5975 460.894 192.063494709 268.888892593 1.88966250666 1.0 +0.19 5.99924 19.5975 460.894 192.063494709 268.888892593 1.88966250666 1.0 +0.2 5.99924 19.5975 460.894 192.063494709 268.888892593 1.88966250666 1.0 +0.21 5.99924 19.5975 460.894 192.063494709 268.888892593 1.88966250666 1.0 +0.22 5.99924 19.5975 460.894 192.063494709 268.888892593 1.88966250666 1.0 +0.23 5.99924 19.5975 460.894 192.063494709 268.888892593 1.88966250666 1.0 +0.24 5.99924 19.5975 460.894 192.063494709 268.888892593 1.88966250666 1.0 +0.25 5.99924 19.5975 460.894 192.063494709 268.888892593 1.88966250666 1.0 +0.26 5.99924 19.5975 460.894 192.063494709 268.888892593 1.88966250666 1.0 +0.27 5.99924 19.5975 460.894 192.063494709 268.888892593 1.88966250666 1.0 +0.28 5.99924 19.5975 460.894 192.063494709 268.888892593 1.88966250666 1.0 +0.29 5.99924 19.5975 460.894 192.063494709 268.888892593 1.88966250666 1.0 +0.3 5.99924 19.5975 460.894 192.063494709 268.888892593 1.88966250666 1.0 +0.31 5.99924 19.5975 460.894 192.063494709 268.888892593 1.88966250666 1.0 +0.32 5.99924 19.5975 460.894 192.063494709 268.888892593 1.88966250666 1.0 +0.33 5.99924 19.5975 460.894 192.063494709 268.888892593 1.88966250666 1.0 +0.34 5.99924 19.5975 460.894 192.063494709 268.888892593 1.88966250666 1.0 +0.35 5.99924 19.5975 460.894 192.063494709 268.888892593 1.88966250666 1.0 +0.36 5.99924 19.5975 460.894 192.063494709 268.888892593 1.88966250666 1.0 +0.37 5.99924 19.5975 460.894 192.063494709 268.888892593 1.88966250666 1.0 +0.38 5.99924 19.5975 460.894 192.063494709 268.888892593 1.88966250666 1.0 +0.39 5.99924 19.5975 460.894 192.063494709 268.888892593 1.88966250666 1.0 +0.4 5.99924 19.5975 460.894 192.063494709 268.888892593 1.88966250666 1.0 +0.41 5.99924 19.5975 460.894 192.063494709 268.888892593 1.88966250666 1.0 +0.42 5.99924 19.5975 460.894 192.063494709 268.888892593 1.88966250666 1.0 +0.43 14.282349952 8.68977441163 1691.6469554 296.107951613 414.551132258 0.674822343487 1.0 +0.44 14.282349952 8.68977441163 1691.6469554 296.107951613 414.551132258 0.674822343487 1.0 +0.45 14.282349952 8.68977441163 1691.6469554 296.107951613 414.551132258 0.674822343487 1.0 +0.46 14.282349952 8.68977441163 1691.6469554 296.107951613 414.551132258 0.674822343487 1.0 +0.47 14.282349952 8.68977441163 1691.6469554 296.107951613 414.551132258 0.674822343487 1.0 +0.48 14.282349952 8.68977441163 1691.6469554 296.107951613 414.551132258 0.674822343487 1.0 +0.49 14.282349952 8.68977441163 1691.6469554 296.107951613 414.551132258 0.674822343487 1.0 +0.5 14.282349952 8.68977441163 1691.6469554 296.107951613 414.551132258 0.674822343487 1.0 +0.51 14.282349952 8.68977441163 1691.6469554 296.107951613 414.551132258 0.674822343487 1.0 +0.52 14.282349952 8.68977441163 1691.6469554 296.107951613 414.551132258 0.674822343487 1.0 +0.53 14.282349952 8.68977441163 1691.6469554 296.107951613 414.551132258 0.674822343487 1.0 +0.54 14.282349952 8.68977441163 1691.6469554 296.107951613 414.551132258 0.674822343487 1.0 +0.55 14.282349952 8.68977441163 1691.6469554 296.107951613 414.551132258 0.674822343487 1.0 +0.56 14.282349952 8.68977441163 1691.6469554 296.107951613 414.551132258 0.674822343487 1.0 +0.57 14.282349952 8.68977441163 1691.6469554 296.107951613 414.551132258 0.674822343487 1.0 +0.58 14.282349952 8.68977441163 1691.6469554 296.107951613 414.551132258 0.674822343487 1.0 +0.59 14.282349952 8.68977441163 1691.6469554 296.107951613 414.551132258 0.674822343487 1.0 +0.6 14.282349952 8.68977441163 1691.6469554 296.107951613 414.551132258 0.674822343487 1.0 +0.61 14.282349952 8.68977441163 1691.6469554 296.107951613 414.551132258 0.674822343487 1.0 +0.62 14.282349952 8.68977441163 1691.6469554 296.107951613 414.551132258 0.674822343487 1.0 +0.63 14.282349952 8.68977441163 1691.6469554 296.107951613 414.551132258 0.674822343487 1.0 +0.64 14.282349952 8.68977441163 1691.6469554 296.107951613 414.551132258 0.674822343487 1.0 +0.65 14.282349952 8.68977441163 1691.6469554 296.107951613 414.551132258 0.674822343487 1.0 +0.66 14.282349952 8.68977441163 1691.6469554 296.107951613 414.551132258 0.674822343487 1.0 +0.67 14.282349952 8.68977441163 1691.6469554 296.107951613 414.551132258 0.674822343487 1.0 +0.68 14.282349952 8.68977441163 1691.6469554 296.107951613 414.551132258 0.674822343487 1.0 +0.69 14.282349952 8.68977441163 1691.6469554 296.107951613 414.551132258 0.674822343487 1.0 +0.7 14.282349952 8.68977441163 1691.6469554 296.107951613 414.551132258 0.674822343487 1.0 +0.71 31.0426016416 8.68977441163 1691.6469554 136.235919828 190.730287759 0.994875359291 0.0 +0.72 31.0426016416 8.68977441163 1691.6469554 136.235919828 190.730287759 0.994875359291 0.0 +0.73 31.0426016416 8.68977441163 1691.6469554 136.235919828 190.730287759 0.994875359291 0.0 +0.74 31.0426016416 8.68977441163 1691.6469554 136.235919828 190.730287759 0.994875359291 0.0 +0.75 31.0426016416 8.68977441163 1691.6469554 136.235919828 190.730287759 0.994875359291 0.0 +0.76 31.0426016416 8.68977441163 1691.6469554 136.235919828 190.730287759 0.994875359291 0.0 +0.77 31.0426016416 8.68977441163 1691.6469554 136.235919828 190.730287759 0.994875359291 0.0 +0.78 31.0426016416 8.68977441163 1691.6469554 136.235919828 190.730287759 0.994875359291 0.0 +0.79 31.0426016416 8.68977441163 1691.6469554 136.235919828 190.730287759 0.994875359291 0.0 +0.8 31.0426016416 8.68977441163 1691.6469554 136.235919828 190.730287759 0.994875359291 0.0 +0.81 31.0426016416 8.68977441163 1691.6469554 136.235919828 190.730287759 0.994875359291 0.0 +0.82 31.0426016416 8.68977441163 1691.6469554 136.235919828 190.730287759 0.994875359291 0.0 +0.83 5.99242 -6.19633 46.095 19.230544588 26.9227624232 1.88818582941 0.0 +0.84 5.99242 -6.19633 46.095 19.230544588 26.9227624232 1.88818582941 0.0 +0.85 5.99242 -6.19633 46.095 19.230544588 26.9227624232 1.88818582941 0.0 +0.86 5.99242 -6.19633 46.095 19.230544588 26.9227624232 1.88818582941 0.0 +0.87 5.99242 -6.19633 46.095 19.230544588 26.9227624232 1.88818582941 0.0 +0.88 5.99242 -6.19633 46.095 19.230544588 26.9227624232 1.88818582941 0.0 +0.89 5.99242 -6.19633 46.095 19.230544588 26.9227624232 1.88818582941 0.0 +0.9 5.99242 -6.19633 46.095 19.230544588 26.9227624232 1.88818582941 0.0 +0.91 5.99242 -6.19633 46.095 19.230544588 26.9227624232 1.88818582941 0.0 +0.92 5.99242 -6.19633 46.095 19.230544588 26.9227624232 1.88818582941 0.0 +0.93 5.99242 -6.19633 46.095 19.230544588 26.9227624232 1.88818582941 0.0 +0.94 5.99242 -6.19633 46.095 19.230544588 26.9227624232 1.88818582941 0.0 +0.95 5.99242 -6.19633 46.095 19.230544588 26.9227624232 1.88818582941 0.0 +0.96 5.99242 -6.19633 46.095 19.230544588 26.9227624232 1.88818582941 0.0 +0.97 5.99242 -6.19633 46.095 19.230544588 26.9227624232 1.88818582941 0.0 +0.98 5.99242 -6.19633 46.095 19.230544588 26.9227624232 1.88818582941 0.0 +0.99 5.99242 -6.19633 46.095 19.230544588 26.9227624232 1.88818582941 0.0 diff --git a/CDMATH/tests/examples/EulerEquations/FullEulerSystem/EulerSystem1D_RiemannProblem_Roe/CMakeLists.txt b/CDMATH/tests/examples/EulerEquations/FullEulerSystem/EulerSystem1D_RiemannProblem_Roe/CMakeLists.txt new file mode 100755 index 0000000..c21587b --- /dev/null +++ b/CDMATH/tests/examples/EulerEquations/FullEulerSystem/EulerSystem1D_RiemannProblem_Roe/CMakeLists.txt @@ -0,0 +1,12 @@ + +if (CDMATH_WITH_PYTHON AND CDMATH_WITH_PETSC AND CDMATH_WITH_POSTPRO) + + SET(ISIMPLICIT 0 )#Explicit scheme + ADD_TEST(ExampleEulerSystem_1DRiemammProblem_ExplicitRoe ${PYTHON_EXECUTABLE} ${CMAKE_CURRENT_SOURCE_DIR}/EulerEquations_RiemannProblem_Roe.py ${ISIMPLICIT} ) + + SET(ISIMPLICIT 1 )#Implicit scheme + ADD_TEST(ExampleEulerSystem_1DRiemammProblem_ImplicitRoe ${PYTHON_EXECUTABLE} ${CMAKE_CURRENT_SOURCE_DIR}/EulerEquations_RiemannProblem_Roe.py ${ISIMPLICIT} ) + +endif (CDMATH_WITH_PYTHON AND CDMATH_WITH_PETSC AND CDMATH_WITH_POSTPRO) + + diff --git a/CDMATH/tests/examples/EulerEquations/FullEulerSystem/EulerSystem1D_RiemannProblem_Roe/EulerEquations_RiemannProblem_Roe.py b/CDMATH/tests/examples/EulerEquations/FullEulerSystem/EulerSystem1D_RiemannProblem_Roe/EulerEquations_RiemannProblem_Roe.py new file mode 100644 index 0000000..2fd3d69 --- /dev/null +++ b/CDMATH/tests/examples/EulerEquations/FullEulerSystem/EulerSystem1D_RiemannProblem_Roe/EulerEquations_RiemannProblem_Roe.py @@ -0,0 +1,644 @@ +#!/usr/bin/env python3 +# -*-coding:utf-8 -* + +""" +Created on Mon Aug 30 2021 +@author: Michael NDJINGA, Katia Ait Ameur, Coraline Mounier + +Euler system without source term in one dimension on regular domain [a,b] +Riemann problemn with ghost cell boundary condition +Left and right: Neumann boundary condition +Roe scheme +Regular square mesh + +State law Stiffened gaz : p = (gamma - 1) * rho * (e - q) - gamma * p0 +4 choices of parameters gamma and p0 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 +import matplotlib.animation as manimation +matplotlib.use("Agg") +import matplotlib.pyplot as plt +from math import sqrt +from numpy import sign +import sys + +## Numerical parameter +precision = 1e-5 + +#state law parameter : can be 'Lagrange', 'Hermite590K', 'Hermite617K', or 'FLICA' +state_law = "Hermite590K" + +#indicates with test case is simulated to compare with FLICA5 results +#num_test = 0 means there are no FLICA5 results given here +num_test = 0 + +#def state_law_parameters(state_law): +#state law Stiffened Gaz : p = (gamma - 1) * rho * e - gamma * p0 +global gamma +global p0 +global q +global c0 +global cp +global h_sat +global T_sat + +if state_law == "Lagrange": + # reference values for Lagrange interpolation + p_ref = 155. * 10**5 #Reference pressure in a REP 900 nuclear power plant + p1 = 153. * 10**5 # value of pressure at inlet of a 900 MWe PWR vessel + rho_ref = 594.38 #density of water at saturation temperature of 617.94K and 155 bars + rho1 = 742.36 # value of density at inlet of a 900 MWe PWR vessel (T1 = 565K) + e_ref = 1603.8 * 10**3 #internal energy of water at saturation temperature of 617.94K and 155 bars + e1 = 1273.6 * 10**3 # value of internal energy at inlet of a 900 MWe PWR vessel + + gamma = (p1 - p_ref) / (rho1 * e1 - rho_ref *e_ref) + 1. + p0 = - 1. / gamma * ( - (gamma - 1) * rho_ref * e_ref + p_ref) + q=0. + c0 = sqrt((gamma - 1) * (e_ref + p_ref / rho_ref)) + + cp = 8950. + h_sat = 1.63 * 10 ** 6 # saturation enthalpy of water at 155 bars + T_sat = 617.94 + +elif state_law == "Hermite617K": + # reference values for Hermite interpolation + p_ref = 155. * 10**5 #Reference pressure in a REP 900 nuclear power plant + T_ref = 617.94 #Reference temperature for interpolation at 617.94K + rho_ref = 594.38 #density of water at saturation temperature of 617.94K and 155 bars + e_ref = 1603.8 * 10**3 #internal energy of water at saturation temperature of 617.94K and 155 bars + h_ref = e_ref + p_ref / rho_ref + c_ref = 621.43 #sound speed for water at 155 bars and 617.94K + + gamma = 1. + c_ref * c_ref / (e_ref + p_ref / rho_ref) # From the sound speed formula + p0 = 1. / gamma * ( (gamma - 1) * rho_ref * e_ref - p_ref) + q=0. + c0 = sqrt((gamma - 1) * (e_ref + p_ref / rho_ref)) + + cp = 8950. # value at 155 bar and 617.94K + h_sat = 1.63 * 10 ** 6 # saturation enthalpy of water at 155 bars + T_sat = 617.94 + +elif state_law == 'Hermite590K': + # reference values for Hermite interpolation + p_ref = 155. * 10**5 #Reference pressure in a REP 900 nuclear power plant + T_ref = 590. #Reference temperature for interpolation at 590K + rho_ref = 688.3 #density of water at 590K and 155 bars + e_ref = 1411.4 * 10**3 #internal energy of water at 590K and 155 bars + h_ref = e_ref + p_ref / rho_ref + c_ref = 866.29 #sound speed for water at 155 bars and 590K + + gamma = 1. + c_ref * c_ref / (e_ref + p_ref / rho_ref) # From the sound speed formula + p0 = 1. / gamma * ( (gamma - 1) * rho_ref * e_ref - p_ref) + q=0. + c0 = sqrt((gamma - 1) * (e_ref + p_ref / rho_ref)) + + cp = 5996.8 # value at 155 bar and 590K + h_sat = 1433.9 * 10 ** 3 # saturation enthalpy of water at 155 bars + T_sat = 590. + +elif state_law == 'Hermite575K': + # reference values for Hermite interpolation + p_ref = 155 * 10**5 #Reference pressure in a REP 900 nuclear power plant + T_ref = 575 #Reference temperature at inlet in a REP 900 nuclear power plant + #Remaining values determined using iapws python package + rho_ref = 722.66 #density of water at 575K and 155 bars + e_ref = 1326552.66 #internal energy of water at 575K and 155 bars + h_ref = e_ref + p_ref / rho_ref + c_ref = 959.28 #sound speed for water at 155 bars and 575K + + gamma = 1 + c_ref * c_ref / (e_ref + p_ref / rho_ref) # From the sound speed formula + p0 = 1 / gamma * ( (gamma - 1) * rho_ref * e_ref - p_ref) + q=0. + c0 = sqrt((gamma - 1) * (e_ref + p_ref / rho_ref))#This is actually c_ref + + cp = 5504.05 # value at 155 bar and 590K + h_sat = h_ref # saturation enthalpy of water at 155 bars + T_sat = T_ref +else: + raise ValueError("Incorrect value for parameter state_law") + + +#initial parameters for Riemann problem (p in Pa, v in m/s, T in K) +p_L = 155. * 10**5 +p_R = 150. * 10**5 +v_L = 0. +v_R = 0. +h_L = 1.4963*10**6 +h_R = 1.4963*10**6 + +T_L = (h_L - h_sat ) / cp + T_sat +T_R = (h_R - h_sat ) / cp + T_sat + +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) + + p_initial = np.array([ (xi < (a + b) / 2) * p_L + (xi >= (a + b) / 2) * p_R for xi in x]) + v_initial = np.array([ (xi < (a + b) / 2) * v_L + (xi >= (a + b) / 2) * v_R for xi in x]) + T_initial = np.array([ (xi < (a + b) / 2) * T_L + (xi >= (a + b) / 2) * T_R for xi in x]) + + rho_initial = p_to_rho_StiffenedGaz(p_initial, T_initial) + q_initial = rho_initial * v_initial + rho_E_initial = T_to_rhoE_StiffenedGaz(T_initial, rho_initial, q_initial) + + return rho_initial, q_initial, rho_E_initial, p_initial, v_initial, T_initial + +def rho_to_p_StiffenedGaz(rho_field, q_field, rho_E_field): + p_field = (gamma - 1) * ( rho_E_field - 1. / 2 * q_field ** 2 / rho_field - rho_field * q) - gamma * p0 + return p_field + + +def p_to_rho_StiffenedGaz(p_field, T_field): + rho_field = (p_field + p0) * gamma / (gamma - 1) * 1 / (h_sat + cp * (T_field - T_sat) - q) + return rho_field + + +def T_to_rhoE_StiffenedGaz(T_field, rho_field, q_field): + rho_E_field = 1 / 2 * (q_field) ** 2 / rho_field + p0 + rho_field / gamma * (h_sat + cp * (T_field- T_sat) + (gamma - 1) * q) + return rho_E_field + + +def rhoE_to_T_StiffenedGaz(rho_field, q_field, rho_E_field): + T_field = T_sat + 1 / cp * (gamma * (rho_E_field / rho_field - 1 / 2 * (q_field / rho_field) ** 2) - gamma * p0 / rho_field - (gamma - 1) * q - h_sat) + return T_field + +def dp_drho_e_const_StiffenedGaz( e ): + return (gamma-1)*(e-q) + +def dp_de_rho_const_StiffenedGaz( rho ): + return (gamma-1)*rho + +def sound_speed_StiffenedGaz( h ): + return np.sqrt((gamma-1)*(h-q)) + +def rho_h_to_p_StiffenedGaz( rho, h ): + return (gamma - 1) * rho * ( h - q ) / gamma - p0 + +def MatRoe_StiffenedGaz( rho_l, q_l, rho_E_l, rho_r, q_r, rho_E_r): + 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) + p_r = rho_to_p_StiffenedGaz(rho_r, q_r, rho_E_r) + 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. ) + e = H - p / rho - 1./2 * u**2 + dp_drho = dp_drho_e_const_StiffenedGaz( e ) + dp_de = dp_de_rho_const_StiffenedGaz( rho ) + + 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 Droe_StiffenedGaz( rho_l, q_l, rho_E_l, rho_r, q_r, rho_E_r): + 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) + p_r = rho_to_p_StiffenedGaz(rho_r, q_r, rho_E_r) + 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. ) + + 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): + + 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) + + Droe = Droe_StiffenedGaz( rho_l, q_l, rho_E_l, rho_r, q_r, rho_E_r) + return (RoeMat - Droe) * coeff * 0.5 + + +def jacobianMatricesp(coeff, rho_l, q_l, rho_E_l, rho_r, q_r, rho_E_r): + 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) + Droe = Droe_StiffenedGaz( rho_l, q_l, rho_E_l, rho_r, q_r, rho_E_r) + + return (RoeMat + Droe) * coeff * 0.5 + + +def FillEdges(j, Uk, nbComp, divMat, Rhs, Un, dt, dx, isImplicit): + dUi = 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) + divMat.addValue(j * nbComp, (j + 1) * nbComp, Am) + divMat.addValue(j * nbComp, j * nbComp, Am * (-1.)) + + dUi[0] = Uk[(j + 1) * nbComp + 0] - Uk[j * nbComp + 0] + dUi[1] = Uk[(j + 1) * nbComp + 1] - Uk[j * nbComp + 1] + dUi[2] = Uk[(j + 1) * nbComp + 2] - Uk[j * nbComp + 2] + temp = Am * dUi + + if(isImplicit): + 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]) + Rhs[j * nbComp + 2] = -temp[2] - (Uk[j * nbComp + 2] - Un[j * nbComp + 2]) + + 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) + divMat.addValue(j * nbComp, j * nbComp, Ap) + divMat.addValue(j * nbComp, (j - 1) * nbComp, Ap * (-1.)) + + dUi[0] = Uk[(j - 1) * nbComp + 0] - Uk[j * nbComp + 0] + dUi[1] = Uk[(j - 1) * nbComp + 1] - Uk[j * nbComp + 1] + dUi[2] = Uk[(j - 1) * nbComp + 2] - Uk[j * nbComp + 2] + + temp = Ap * dUi + + if(isImplicit): + 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]) + Rhs[j * nbComp + 2] = temp[2] - (Uk[j * nbComp + 2] - Un[j * nbComp + 2]) + +def FillInnerCell(j, Uk, nbComp, divMat, Rhs, Un, dt, dx, isImplicit): + + 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) + + 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) + + 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(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 + + if(isImplicit): + 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]) + +def SetPicture(rho_field_Roe, q_field_Roe, h_field_Roe, p_field_Roe, v_field_Roe, T_field_Roe, dx): + max_initial_rho = max(rho_field_Roe) + min_initial_rho = min(rho_field_Roe) + max_initial_q = max(q_field_Roe) + min_initial_q = min(q_field_Roe) + min_initial_h = min(h_field_Roe) + max_initial_h = max(h_field_Roe) + max_initial_p = max(p_field_Roe) + min_initial_p = min(p_field_Roe) + max_initial_v = max(v_field_Roe) + min_initial_v = min(v_field_Roe) + max_initial_T = max(T_field_Roe) + min_initial_T = min(T_field_Roe) + + fig, ([axDensity, axMomentum, axRhoE],[axPressure, axVitesse, axTemperature]) = plt.subplots(2, 3,sharex=True, figsize=(14,10)) + plt.gcf().subplots_adjust(wspace = 0.5) + + lineDensity_Roe, = axDensity.plot([a+0.5*dx + i*dx for i in range(nx)], rho_field_Roe, label='Roe') + axDensity.set(xlabel='x (m)', ylabel='Densité (kg/m3)') + axDensity.set_xlim(a,b) + #axDensity.set_ylim(min_initial_rho - 0.1 * (max_initial_rho - min_initial_rho), 700.) + axDensity.set_ylim(657, 660.5) + axDensity.legend() + + lineMomentum_Roe, = axMomentum.plot([a+0.5*dx + i*dx for i in range(nx)], q_field_Roe, label='Roe') + axMomentum.set(xlabel='x (m)', ylabel='Momentum (kg/(m2.s))') + axMomentum.set_xlim(a,b) + #axMomentum.set_ylim(min_initial_q - 0.1*(max_initial_q-min_initial_q), max_initial_q * 2.5) + axMomentum.set_ylim(-50, 500) + axMomentum.legend() + + lineRhoE_Roe, = axRhoE.plot([a+0.5*dx + i*dx for i in range(nx)], h_field_Roe, label='Roe') + axRhoE.set(xlabel='x (m)', ylabel='h (J/m3)') + axRhoE.set_xlim(a,b) + #axRhoE.set_ylim(min_initial_h - 0.05*abs(min_initial_h), max_initial_h + 0.5*(max_initial_h-min_initial_h)) + axRhoE.set_ylim(1.495 * 10**6, 1.5*10**6) + axRhoE.legend() + + linePressure_Roe, = axPressure.plot([a+0.5*dx + i*dx for i in range(nx)], p_field_Roe, label='Roe') + axPressure.set(xlabel='x (m)', ylabel='Pression (bar)') + axPressure.set_xlim(a,b) + #axPressure.set_ylim(min_initial_p - 0.05*abs(min_initial_p), max_initial_p + 0.5*(max_initial_p-min_initial_p)) + axPressure.set_ylim(149, 156) + axPressure.ticklabel_format(axis='y', style='sci', scilimits=(0,0)) + axPressure.legend() + + lineVitesse_Roe, = axVitesse.plot([a+0.5*dx + i*dx for i in range(nx)], v_field_Roe, label='Roe') + axVitesse.set(xlabel='x (m)', ylabel='Vitesse (m/s)') + axVitesse.set_xlim(a,b) + #axVitesse.set_ylim(min_initial_v - 0.05*abs(min_initial_v), 15) + axVitesse.set_ylim(-0.5, 1) + axVitesse.ticklabel_format(axis='y', style='sci', scilimits=(0,0)) + axVitesse.legend() + + lineTemperature_Roe, = axTemperature.plot([a+0.5*dx + i*dx for i in range(nx)], T_field_Roe, label='Roe') + axTemperature.set(xlabel='x (m)', ylabel='Température (K)') + axTemperature.set_xlim(a,b) + #axTemperature.set_ylim(min_initial_T - 0.005*abs(min_initial_T), 604) + axTemperature.set_ylim(600, 601) + axTemperature.ticklabel_format(axis='y', style='sci', scilimits=(0,0)) + axTemperature.legend() + + return(fig, lineDensity_Roe, lineMomentum_Roe, lineRhoE_Roe, linePressure_Roe, lineVitesse_Roe, lineTemperature_Roe, lineDensity_Roe, lineMomentum_Roe, lineRhoE_Roe, linePressure_Roe, lineVitesse_Roe, lineTemperature_Roe) + + +def EulerSystemRoe(ntmax, tmax, cfl, a, b, nbCells, output_freq, meshName, state_law, isImplicit): + #state_law_parameters(state_law) + dim = 1 + nbComp = 3 + dt = 0. + time = 0. + it = 0 + isStationary = False + dx = (b - a) / nx + dt = cfl * dx / c0 + #dt = 5*10**(-6) + nbVoisinsMax = 2 + + # iteration vectors + Un_Roe = cdmath.Vector(nbCells * (nbComp)) + dUn_Roe = cdmath.Vector(nbCells * (nbComp)) + dUk_Roe = cdmath.Vector(nbCells * (nbComp)) + Rhs_Roe = cdmath.Vector(nbCells * (nbComp)) + + # Initial conditions + print("Construction of the initial condition …") + + rho_field_Roe, q_field_Roe, rho_E_field_Roe, p_field_Roe, v_field_Roe, T_field_Roe = initial_conditions_Riemann_problem(a, b, nx) + h_field_Roe = (rho_E_field_Roe + p_field_Roe) / rho_field_Roe - 0.5 * (q_field_Roe / rho_field_Roe) **2 + p_field_Roe = p_field_Roe * 10 ** (-5) + + + for k in range(nbCells): + Un_Roe[k * nbComp + 0] = rho_field_Roe[k] + Un_Roe[k * nbComp + 1] = q_field_Roe[k] + Un_Roe[k * nbComp + 2] = rho_E_field_Roe[k] + + divMat_Roe = cdmath.SparseMatrixPetsc(nbCells * nbComp, nbCells * nbComp, (nbVoisinsMax + 1) * nbComp) + + # Picture settings + fig, lineDensity, lineMomentum, lineRhoE, linePressure, lineVitesse, lineTemperature, lineDensity_Roe, lineMomentum_Roe, lineRhoE_Roe, linePressure_Roe, lineVitesse_Roe, lineTemperature_Roe = SetPicture( rho_field_Roe, q_field_Roe, h_field_Roe, p_field_Roe, v_field_Roe, T_field_Roe, dx) + if(isImplicit): + ImplicitOrExplicit="Implicit" + else: + ImplicitOrExplicit="Explicit" + + # Video settings + FFMpegWriter = manimation.writers['ffmpeg'] + metadata = dict(title=ImplicitOrExplicit+" Roe scheme for the 1D Euler system", artist="CEA Saclay", comment="Shock tube") + writer = FFMpegWriter(fps=10, metadata=metadata, codec='h264') + with writer.saving(fig, "1DEulerEquations_RiemannProblem_"+ImplicitOrExplicit+"Roe" + ".mp4", ntmax): + writer.grab_frame() + plt.savefig("EulerEquations_RiemannProblem_" + str(dim) + "D_"+ImplicitOrExplicit+"Roe" + meshName + "0" + ".png") + np.savetxt( "EulerEquations_RiemannProblem_" + str(dim) + "D_"+ImplicitOrExplicit+"Roe" + meshName + "_rho" + "0" + ".txt", rho_field_Roe, delimiter="\n") + np.savetxt( "EulerEquations_RiemannProblem_" + str(dim) + "D_"+ImplicitOrExplicit+"Roe" + meshName + "_q" + "0" + ".txt", q_field_Roe, delimiter="\n") + iterGMRESMax = 50 + newton_max = 100 + + print("Starting computation of the non linear Euler non isentropic system with Roe scheme …") + # STARTING TIME LOOP + while (it < ntmax and time <= tmax and not isStationary): + dUn_Roe = Un_Roe.deepCopy() + Uk_Roe = Un_Roe.deepCopy() + residu_Roe = 1. + + k_Roe = 0 + while (k_Roe < newton_max and residu_Roe > precision): + # STARTING NEWTON LOOP + divMat_Roe.zeroEntries() #sets the matrix coefficients to zero + for j in range(nbCells): + + # traitements des bords + if (j == 0): + FillEdges(j, Uk_Roe, nbComp, divMat_Roe, Rhs_Roe, Un_Roe, dt, dx, isImplicit) + elif (j == nbCells - 1): + FillEdges(j, Uk_Roe, nbComp, divMat_Roe, Rhs_Roe, Un_Roe, dt, dx, isImplicit) + + # traitement des cellules internes + else: + FillInnerCell(j, Uk_Roe, nbComp, divMat_Roe, Rhs_Roe, Un_Roe, dt, dx, isImplicit) + + #solving the linear system divMat * dUk = Rhs + + if(isImplicit): + divMat_Roe.diagonalShift(1.) + LS_Roe = cdmath.LinearSolver(divMat_Roe, Rhs_Roe, iterGMRESMax, precision, "GMRES", "LU") + dUk_Roe = LS_Roe.solve() + residu_Roe = dUk_Roe.norm() + else: + dUk_Roe=Rhs_Roe - divMat_Roe*Un_Roe + residu_Roe = 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 : ", residu_Roe) + print("Linear system converged in ", LS_Roe.getNumberOfIter(), " GMRES iterations") + + #updates for Newton loop + Uk_Roe += dUk_Roe + k_Roe = k_Roe + 1 + if (isImplicit and not LS_Roe.getStatus()): + print("Linear system did not converge ", LS_Roe.getNumberOfIter(), " GMRES iterations") + raise ValueError("No convergence of the linear system") + + if k_Roe == newton_max: + raise ValueError("No convergence of Newton Roe Scheme") + + #updating fields + Un_Roe = Uk_Roe.deepCopy() + dUn_Roe -= Un_Roe + if (dUn_Roe.norm()= ntmax or isStationary or time >= tmax): + + print("-- Time step : " + str(it) + ", Time: " + str(time) + ", dt: " + str(dt)) + + plt.savefig("EulerEquations_RiemannProblem_" + str(dim) + "D_"+ImplicitOrExplicit+"Roe" + meshName + str(it) + '_time' + str(time) + ".png") + + print("-- Iter: " + str(it) + ", Time: " + str(time) + ", dt: " + str(dt)) + if (it >= ntmax): + print("Maximum number of time steps ntmax= ", ntmax, " reached") + return + + elif (isStationary): + print("Stationary regime reached at time step ", it, ", t= ", time) + print("------------------------------------------------------------------------------------") + np.savetxt("EulerEquations_RiemannProblem_" + str(dim) + "D_"+ImplicitOrExplicit+"Roe" + meshName + "_rho_Stat.txt", rho_field_Roe, delimiter="\n") + np.savetxt("EulerEquations_RiemannProblem_" + str(dim) + "D_"+ImplicitOrExplicit+"Roe" + meshName + "_q_Stat.txt", q_field_Roe, delimiter="\n") + np.savetxt("EulerEquations_RiemannProblem_" + str(dim) + "D_"+ImplicitOrExplicit+"Roe" + meshName + "_rhoE_Stat.txt", rho_E_field_Roe, delimiter="\n") + np.savetxt("EulerEquations_RiemannProblem_" + str(dim) + "D_"+ImplicitOrExplicit+"Roe" + meshName + "_p_Stat.txt", p_field_Roe, delimiter="\n") + plt.savefig("EulerEquations_RiemannProblem_" + str(dim) + "D_"+ImplicitOrExplicit+"Roe" + meshName + "_Stat.png") + return + else: + print("Maximum time Tmax= ", tmax, " reached") + return + + +def solve(a, b, nx, meshName, meshType, cfl, state_law, isImplicit): + print("Resolution of the Euler System in dimension 1 on " + str(nx) + " cells") + print("State Law Stiffened Gaz, " + state_law) + print("Initial data : ", "Riemann problem") + print("Boundary conditions : ", "Neumann") + print("Mesh name : ", meshName, ", ", nx, " cells") + # Problem data + tmax = 10. + ntmax = 25 + output_freq = 1 + EulerSystemRoe(ntmax, tmax, cfl, a, b, nx, output_freq, meshName, state_law, isImplicit) + return + + +if __name__ == """__main__""": + a = 0. + b = 1. + nx = 50 + cfl = 0.95 + isImplicit = bool(int(sys.argv[1])) + #state_law = "Hermite590K" + #state_law_parameters(state_law) + solve(a, b, nx, "RegularSquares", "", cfl, state_law, isImplicit)