From 386965127ce8969481a837335e9fe38e21351ab7 Mon Sep 17 00:00:00 2001 From: michael Date: Sat, 9 Oct 2021 23:59:25 +0200 Subject: [PATCH] Added examples for the Full Euler system --- .../FullEulerSystem/CMakeLists.txt | 8 +- .../CMakeLists.txt | 2 +- .../Euler_complet_HeatedChanel_MAC.py | 893 ++++++++++++++++++ .../CMakeLists.txt | 2 +- .../Euler_complet_HeatedChanel_Roe.py} | 4 +- .../CMakeLists.txt | 9 + .../Euler_complet_Toro4_convergence.py | 548 +++++++++++ .../EulerSystem1D_RiemannProblem_MAC/TTC4.dat | 100 ++ .../CMakeLists.txt | 12 + .../EulerEquations_RiemannProblem_Roe.py} | 46 +- 10 files changed, 1594 insertions(+), 30 deletions(-) rename CDMATH/tests/examples/EulerEquations/FullEulerSystem/{EulerSystem1D_RiemannProblem => EulerSystem1D_HeatedChannel_MAC}/CMakeLists.txt (51%) create mode 100644 CDMATH/tests/examples/EulerEquations/FullEulerSystem/EulerSystem1D_HeatedChannel_MAC/Euler_complet_HeatedChanel_MAC.py rename CDMATH/tests/examples/EulerEquations/FullEulerSystem/{EulerSystem1D_HeatedChannel => EulerSystem1D_HeatedChannel_Roe}/CMakeLists.txt (95%) rename CDMATH/tests/examples/EulerEquations/FullEulerSystem/{EulerSystem1D_HeatedChannel/Euler_complet_HeatedChanel.py => EulerSystem1D_HeatedChannel_Roe/Euler_complet_HeatedChanel_Roe.py} (97%) create mode 100755 CDMATH/tests/examples/EulerEquations/FullEulerSystem/EulerSystem1D_RiemannProblem_MAC/CMakeLists.txt create mode 100644 CDMATH/tests/examples/EulerEquations/FullEulerSystem/EulerSystem1D_RiemannProblem_MAC/Euler_complet_Toro4_convergence.py create mode 100644 CDMATH/tests/examples/EulerEquations/FullEulerSystem/EulerSystem1D_RiemannProblem_MAC/TTC4.dat create mode 100755 CDMATH/tests/examples/EulerEquations/FullEulerSystem/EulerSystem1D_RiemannProblem_Roe/CMakeLists.txt rename CDMATH/tests/examples/EulerEquations/FullEulerSystem/{EulerSystem1D_RiemannProblem/Euler_complet_RP.py => EulerSystem1D_RiemannProblem_Roe/EulerEquations_RiemannProblem_Roe.py} (90%) 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_RiemannProblem/CMakeLists.txt b/CDMATH/tests/examples/EulerEquations/FullEulerSystem/EulerSystem1D_HeatedChannel_MAC/CMakeLists.txt similarity index 51% rename from CDMATH/tests/examples/EulerEquations/FullEulerSystem/EulerSystem1D_RiemannProblem/CMakeLists.txt rename to CDMATH/tests/examples/EulerEquations/FullEulerSystem/EulerSystem1D_HeatedChannel_MAC/CMakeLists.txt index 895c803..6cc4002 100755 --- a/CDMATH/tests/examples/EulerEquations/FullEulerSystem/EulerSystem1D_RiemannProblem/CMakeLists.txt +++ b/CDMATH/tests/examples/EulerEquations/FullEulerSystem/EulerSystem1D_HeatedChannel_MAC/CMakeLists.txt @@ -1,7 +1,7 @@ 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 ) + 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/CMakeLists.txt b/CDMATH/tests/examples/EulerEquations/FullEulerSystem/EulerSystem1D_HeatedChannel_Roe/CMakeLists.txt similarity index 95% rename from CDMATH/tests/examples/EulerEquations/FullEulerSystem/EulerSystem1D_HeatedChannel/CMakeLists.txt rename to CDMATH/tests/examples/EulerEquations/FullEulerSystem/EulerSystem1D_HeatedChannel_Roe/CMakeLists.txt index 8415053..bf750c8 100755 --- a/CDMATH/tests/examples/EulerEquations/FullEulerSystem/EulerSystem1D_HeatedChannel/CMakeLists.txt +++ b/CDMATH/tests/examples/EulerEquations/FullEulerSystem/EulerSystem1D_HeatedChannel_Roe/CMakeLists.txt @@ -1,7 +1,7 @@ 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 ) + 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/Euler_complet_HeatedChanel.py b/CDMATH/tests/examples/EulerEquations/FullEulerSystem/EulerSystem1D_HeatedChannel_Roe/Euler_complet_HeatedChanel_Roe.py similarity index 97% rename from CDMATH/tests/examples/EulerEquations/FullEulerSystem/EulerSystem1D_HeatedChannel/Euler_complet_HeatedChanel.py rename to CDMATH/tests/examples/EulerEquations/FullEulerSystem/EulerSystem1D_HeatedChannel_Roe/Euler_complet_HeatedChanel_Roe.py index e63c38c..859cde8 100644 --- a/CDMATH/tests/examples/EulerEquations/FullEulerSystem/EulerSystem1D_HeatedChannel/Euler_complet_HeatedChanel.py +++ b/CDMATH/tests/examples/EulerEquations/FullEulerSystem/EulerSystem1D_HeatedChannel_Roe/Euler_complet_HeatedChanel_Roe.py @@ -31,9 +31,7 @@ import matplotlib matplotlib.use("Agg") import matplotlib.pyplot as plt -import matplotlib.animation as manimation -import sys -from math import sqrt, atan, pi +from math import sqrt from numpy import sign 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/Euler_complet_RP.py b/CDMATH/tests/examples/EulerEquations/FullEulerSystem/EulerSystem1D_RiemannProblem_Roe/EulerEquations_RiemannProblem_Roe.py similarity index 90% rename from CDMATH/tests/examples/EulerEquations/FullEulerSystem/EulerSystem1D_RiemannProblem/Euler_complet_RP.py rename to CDMATH/tests/examples/EulerEquations/FullEulerSystem/EulerSystem1D_RiemannProblem_Roe/EulerEquations_RiemannProblem_Roe.py index 8785225..2fd3d69 100644 --- a/CDMATH/tests/examples/EulerEquations/FullEulerSystem/EulerSystem1D_RiemannProblem/Euler_complet_RP.py +++ b/CDMATH/tests/examples/EulerEquations/FullEulerSystem/EulerSystem1D_RiemannProblem_Roe/EulerEquations_RiemannProblem_Roe.py @@ -27,14 +27,12 @@ To do correct the computation of the time step : lambda_max (maximum eigenvalue) import cdmath import numpy as np import matplotlib - +import matplotlib.animation as manimation matplotlib.use("Agg") import matplotlib.pyplot as plt -import matplotlib.animation as manimation -import sys -from math import sqrt, atan, pi +from math import sqrt from numpy import sign - +import sys ## Numerical parameter precision = 1e-5 @@ -461,7 +459,7 @@ def SetPicture(rho_field_Roe, q_field_Roe, h_field_Roe, p_field_Roe, v_field_Roe 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): +def EulerSystemRoe(ntmax, tmax, cfl, a, b, nbCells, output_freq, meshName, state_law, isImplicit): #state_law_parameters(state_law) dim = 1 nbComp = 3 @@ -469,7 +467,6 @@ def EulerSystemRoe(ntmax, tmax, cfl, a, b, nbCells, output_freq, meshName, state time = 0. it = 0 isStationary = False - isImplicit = False dx = (b - a) / nx dt = cfl * dx / c0 #dt = 5*10**(-6) @@ -498,16 +495,20 @@ def EulerSystemRoe(ntmax, tmax, cfl, a, b, nbCells, output_freq, meshName, state # 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="Roe scheme for the 1D Euler system", artist="CEA Saclay", comment="Shock tube") + 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, "1DEuler_complet_RP_Roe" + ".mp4", ntmax): + with writer.saving(fig, "1DEulerEquations_RiemannProblem_"+ImplicitOrExplicit+"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") + 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 @@ -597,7 +598,7 @@ def EulerSystemRoe(ntmax, tmax, cfl, a, b, nbCells, output_freq, meshName, state 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") + 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): @@ -607,18 +608,18 @@ def EulerSystemRoe(ntmax, tmax, cfl, a, b, nbCells, output_freq, meshName, state 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") + 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): +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") @@ -628,7 +629,7 @@ def solve(a, b, nx, meshName, meshType, cfl, state_law): tmax = 10. ntmax = 25 output_freq = 1 - EulerSystemRoe(ntmax, tmax, cfl, a, b, nx, output_freq, meshName, state_law) + EulerSystemRoe(ntmax, tmax, cfl, a, b, nx, output_freq, meshName, state_law, isImplicit) return @@ -636,7 +637,8 @@ if __name__ == """__main__""": a = 0. b = 1. nx = 50 - cfl = 0.5 + 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) + solve(a, b, nx, "RegularSquares", "", cfl, state_law, isImplicit) -- 2.39.2