From: michael Date: Sat, 9 Oct 2021 22:02:53 +0000 (+0200) Subject: Added validation tests for the Reference solutions X-Git-Tag: V9_8_0~85 X-Git-Url: http://git.salome-platform.org/gitweb/?a=commitdiff_plain;h=a44cffcfd58ffc199e1db5effdeb2358005db6b6;p=tools%2Fsolverlab.git Added validation tests for the Reference solutions --- diff --git a/CDMATH/tests/validation/ReferenceSolutions/CMakeLists.txt b/CDMATH/tests/validation/ReferenceSolutions/CMakeLists.txt new file mode 100644 index 0000000..37ce57a --- /dev/null +++ b/CDMATH/tests/validation/ReferenceSolutions/CMakeLists.txt @@ -0,0 +1,21 @@ + +SET(SCRIPT + ../scripts/ReferenceSolutions/exact_rs_stiffenedgas.py + ExactRiemannProblemSolutions.py + ) + +file(COPY ${SCRIPT} DESTINATION ${CMAKE_CURRENT_BINARY_DIR}) + +if (CDMATH_WITH_PYTHON AND CDMATH_WITH_PETSC AND CDMATH_WITH_POSTPRO) + + ADD_TEST(validationReferenceSolutions_RiemannProblem ${PYTHON_EXECUTABLE} ${CMAKE_CURRENT_SOURCE_DIR}/ExactRiemannProblemSolutions.py ) + +endif (CDMATH_WITH_PYTHON AND CDMATH_WITH_PETSC AND CDMATH_WITH_POSTPRO) + +install( DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}/${dir}/ DESTINATION share/validation/ReferenceSolutions + FILES_MATCHING PATTERN "*.dat" + PATTERN "*.sh" + PATTERN "*.py" + PATTERN "*.png" +) + diff --git a/CDMATH/tests/validation/ReferenceSolutions/ExactRiemannProblemSolutions.py b/CDMATH/tests/validation/ReferenceSolutions/ExactRiemannProblemSolutions.py new file mode 100644 index 0000000..b1ad296 --- /dev/null +++ b/CDMATH/tests/validation/ReferenceSolutions/ExactRiemannProblemSolutions.py @@ -0,0 +1,442 @@ +#!/usr/bin/env python3 +# -*-coding:utf-8 -* + +###################################################################################################################### +# This file runs a set of Riemann problems for the one dimensional Euler equations with stiffened gas equation of state. +# It uses the class exact_rs_stiffenedgas to compute the exact solution of the Riemann Problems. +# +# Author: Michael Ndjinga +# Date: 19/02/2021 +# Description : Translated from C++ package developped by Murray Cutforth +####################################################################################################################### + +import exact_rs_stiffenedgas +from math import fabs +import matplotlib +matplotlib.use("Agg") +import matplotlib.pyplot as plt + +def run_Riemann_problems(numsamples = 100): + # Output test solution for many different Riemann problems + + print("") + print("Determination of the exact solutions of some Riemann problems for the Euler equations on " + str(numsamples) + " sample points.") + + for TC in range(1, 23): + WL = [0.]*3 + WR = [0.]*3 + + gammaL = 1.4; + gammaR = 1.4; + pinf_L = 0.0; + pinf_R = 0.0; + offset = 0.5;#position of the initial discontinuity + xmin = 0.0; + xmax = 1.0; + + + # TC1 to TC5 are the 7 shock tube problems from Toro + + if (TC == 1): + WL[0] = 1.0; + WL[1] = 0.75; + WL[2] = 1.0; + WR[0] = 0.125; + WR[1] = 0.0; + WR[2] = 0.1; + t = 0.2; + offset = 0.3; + filename = "TTC1"; + elif (TC ==2): + WL[0] = 1.0; + WL[1] = -2.0; + WL[2] = 0.4; + WR[0] = 1.0; + WR[1] = 2.0; + WR[2] = 0.4; + t = 0.15; + filename = "TTC2"; + elif (TC ==3): + WL[0] = 1.0; + WL[1] = 0.0; + WL[2] = 1000.0; + WR[0] = 1.0; + WR[1] = 0.0; + WR[2] = 0.01; + t = 0.012; + filename = "TTC3"; + elif (TC ==4): + WL[0] = 5.99924; + WL[1] = 19.5975; + WL[2] = 460.894; + WR[0] = 5.99242; + WR[1] = -6.19633; + WR[2] = 46.0950; + t = 0.035; + offset = 0.4; + filename = "TTC4"; + elif (TC ==5): + WL[0] = 1.0; + WL[1] = -19.59745; + WL[2] = 1000.; + WR[0] = 1.0; + WR[1] = -19.59745; + WR[2] = 0.01; + t = 0.012; + offset = 0.8; + filename = "TTC5"; + elif (TC == 21): + WL[0] = 1.4; + WL[1] = 0.0; + WL[2] = 1.0; + WR[0] = 1.; + WR[1] = 0.0; + WR[2] = 1.; + t = 2.; + filename = "TTC6"; + elif (TC == 22): + WL[0] = 1.4; + WL[1] = 0.1; + WL[2] = 1.0; + WR[0] = 1.; + WR[1] = 0.1; + WR[2] = 1.; + t = 2.; + filename = "TTC7"; + elif (TC == 6): + # Air - Helium shock tube from Sambasivan 2009 + + WL[0] = 1.0; + WL[1] = 0.0; + WL[2] = 1.0; + WR[0] = 0.125; + WR[1] = 0.0; + WR[2] = 0.1; + t = 0.25; + filename = "Samb1"; + gammaR = 1.667; + elif (TC == 7): + # Gaseous shock tube from rGFM + + WL[0] = 1.0; + WL[1] = 0.0; + WL[2] = 100000.0; + WR[0] = 0.125; + WR[1] = 0.0; + WR[2] = 10000.0; + t = 0.0007; + filename = "NE1"; + gammaR = 1.2; + elif (TC == 8): + # Air - water shock from rGFM + + WL[0] = 0.00596521; + WL[1] = 911.8821; + WL[2] = 1000.0; + WR[0] = 1.0; + WR[1] = 0.0; + WR[2] = 1.0; + gammaR = 7.15; + pinf_R = 3309.0; + t = 0.0007; + filename = "rGFM2"; + elif (TC == 9): + # Air - water jet from rGFM + + WL[0] = 1.0; + WL[1] = 90.0; + WL[2] = 1.0; + WR[0] = 1000.0; + WR[1] = 0.0; + WR[2] = 1.0; + gammaR = 7.15; + pinf_R = 3309.0; + t = 0.015; + filename = "rGFM4"; + offset = 0.6; + elif (TC == 10): + # Reversed version of TC 9 + + WR[0] = 1.0; + WR[1] = -90.0; + WR[2] = 1.0; + WL[0] = 1000.0; + WL[1] = 0.0; + WL[2] = 1.0; + gammaL = 7.15; + pinf_L = 3309.0; + t = 0.015; + filename = "rGFM_reversed"; + offset = 0.4; + elif (TC == 11): + # Water - air shock from Saurel 1999 + + WL[0] = 1.0; + WL[1] = 0.0; + WL[2] = 4.0; + WR[0] = 0.05; + WR[1] = 0.0; + WR[2] = 0.0004; + gammaL = 4.4; + pinf_L = 2.4; + t = 0.12; + filename = "Saurel1"; + offset = 0.7; + elif (TC == 12): + # Reversed version of TC 11 + + WR[0] = 1.0; + WR[1] = 0.0; + WR[2] = 4.0; + WL[0] = 0.05; + WL[1] = 0.0; + WL[2] = 0.0004; + gammaR = 4.4; + pinf_R = 2.4; + t = 0.12; + filename = "Saurel1_reversed"; + offset = 0.3; + elif (TC == 13): + # Numerical experiment 2 + + WL[0] = 3.175962; + WL[1] = 9.434992; + WL[2] = 100.0; + WR[0] = 1.0; + WR[1] = 0.0; + WR[2] = 1.0; + gammaL = 1.667; + gammaR = 1.2; + t = 0.045; + filename = "NE2"; + offset = 0.2; + elif (TC == 14): + # Numerical experiment 3 + + WL[0] = 0.00596521; + WL[1] = 911.8821; + WL[2] = 100.0; + WR[0] = 1.0; + WR[1] = 0.0; + WR[2] = 1.0; + gammaL = 1.4; + gammaR = 7.15; + pinf_L = 0.0; + pinf_R = 3309.0; + t = 0.0007; + filename = "NE3"; + offset = 0.5; + elif (TC == 15): + # Gaseous shock tube from So/Hu/Adams 2012 + + WL[0] = 1.0; + WL[1] = 0.0; + WL[2] = 1.0; + WR[0] = 0.125; + WR[1] = 0.0; + WR[2] = 0.1; + t = 0.15; + filename = "ST3"; + gammaR = 1.667; + elif (TC == 16): + # Gaseous shock tube from Garrick/Owkes/Regele 2016 + + WL[0] = 1.0; + WL[1] = 0.0; + WL[2] = 1.0; + WR[0] = 0.125; + WR[1] = 0.0; + WR[2] = 0.1; + t = 0.14; + filename = "ST1"; + gammaR = 2.4; + elif (TC == 17): + # Water - air shock tube from Murrone/Guillard 2004 + + WL[0] = 1000.0; + WL[1] = 0.0; + WL[2] = 1000000000.0; + WR[0] = 50.0; + WR[1] = 0.0; + WR[2] = 100000.0; + t = 0.0009; + filename = "ST2"; + offset = 0.7; + gammaL = 4.4; + gammaR = 1.4; + pinf_L = 600000000.0; + pinf_R = 0.0; + xmin = -2.0; + xmax = 2.0; + elif (TC == 18): + # Water shock tube from PWR - version 1 + + WL[0] = 700.0; + WL[1] = 0.0; + WL[2] = 15500000.0; + WR[0] = 700.0; + WR[1] = 0.0; + WR[2] = 100000.0; + t = 3.e-4; + filename = "PWR-ShockTube1"; + gammaL = 1.58; + gammaR = 1.58; + pinf_L = 353637173.0; + pinf_R = 353637173.0; + elif (TC == 19): + # Water shock tube from PWR - version 2 + + WL[0] = 700.0; + WL[1] = 0.0; + WL[2] = 15500000.0; + WR[0] = 700.0; + WR[1] = 20.0; + WR[2] = 15500000.0; + t = 3.e-4; + filename = "PWR-ShockTube2"; + gammaL = 1.58; + gammaR = 1.58; + pinf_L = 353637173.0; + pinf_R = 353637173.0; + elif (TC == 20): + # Water shock tube from PWR - version 3 + + WL[0] = 700.0; + WL[1] = 0.0; + WL[2] = 15500000.0; + WR[0] = 650.0; + WR[1] = 0.0; + WR[2] = 15500000.0; + t = 3.e-4; + filename = "PWR-ShockTube3"; + gammaL = 1.58; + gammaR = 1.58; + pinf_L = 353637173.0; + pinf_R = 353637173.0; + else: + print( "Unknown test case" ) + return 1; + + + RS = exact_rs_stiffenedgas.exact_rs_stiffenedgas(gammaL, gammaR, pinf_L, pinf_R); + RS.solve_RP(WL,WR); + + print( "") + print( "Solved Riemann problem for TC = " , TC ) + print( "Star state pressure calculated as " , RS.P_STAR ) + print( "Star state velocity calculated as " , RS.S_STAR ) + print( "Left star state density calculated as " , RS.rho_star_L ) + print( "Right state state density calculated as " , RS.rho_star_R ) + print( "Left shock speed calculated as " , RS.S_L ) + print( "Right shock speed calculated as " , RS.S_R ) + print( "Left rarefaction head speed calculated as " , RS.S_HL ) + print( "Left rarefaction tail speed calculated as " , RS.S_TL ) + print( "Right rarefaction head speed calculated as " , RS.S_HR ) + print( "Right rarefaction tail speed calculated as " , RS.S_TR ) + + + dx = (xmax - xmin)/numsamples; + + outfile=open(filename+'.dat', 'w') + + rho_field = [0]*(numsamples+1) + u_field = [0]*(numsamples+1) + p_field = [0]*(numsamples+1) + q_field = [0]*(numsamples+1) + h_field = [0]*(numsamples+1) + e_field = [0]*(numsamples+1) + + x = xmin; + for i in range(numsamples+1): + + S = x/t; + soln = RS.sample_solution(WL, WR, S - offset/t); + + thisgamma = gammaL if S - offset/t < RS.S_STAR else gammaR; + thispinf = pinf_L if S - offset/t < RS.S_STAR else pinf_R; + thisz = 1.0 if S - offset/t < RS.S_STAR else 0.0; + + rho_field[i]=soln[0] + u_field[i] =soln[1] + p_field[i] =soln[2] + q_field[i] =rho_field[i]*u_field[i] + h_field[i] =exact_rs_stiffenedgas.stiffenedgas_h(rho_field[i], p_field[i], thisgamma, thispinf) + e_field[i] =exact_rs_stiffenedgas.stiffenedgas_e(rho_field[i], p_field[i], thisgamma, thispinf) + + outfile.write( str( x ) + " " + str( soln[0]) + " " + str( soln[1]) + " " + str( soln[2]) + " " + str(exact_rs_stiffenedgas.stiffenedgas_e(soln[0], soln[2], thisgamma, thispinf)) + " " + str(exact_rs_stiffenedgas.stiffenedgas_h(soln[0], soln[2], thisgamma, thispinf)) + " " + str(fabs(soln[1])/RS.a(soln[0], soln[2], thisgamma, thispinf)) + " " + str(thisz) + "\n") + + x += dx; + + outfile.close() + + #Set picture parameters + max_initial_rho = max(WL[0],WR[0]) + min_initial_rho = min(WL[0],WR[0]) + max_initial_p = max(WL[2],WR[2]) + min_initial_p = min(WL[2],WR[2]) + max_initial_v = max(WL[2],WR[2]) + min_initial_v = min(WL[2],WR[2]) + max_initial_q = max_initial_rho*max_initial_v + min_initial_q = min_initial_rho*min_initial_v + + e_L=exact_rs_stiffenedgas.p_to_e_StiffenedGaz(WL[2], WL[0], gammaL, pinf_L) + e_R=exact_rs_stiffenedgas.p_to_e_StiffenedGaz(WR[2], WR[0], gammaR, pinf_R) + h_L=e_L+WL[2]/WL[0] + h_R=e_R+WR[2]/WR[0] + 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) + + # Create plots + fig, ([axDensity, axMomentum, axEnthalpie],[axPressure, axVitesse, axEinterne]) = plt.subplots(2, 3,sharex=True, figsize=(14,10)) + plt.gcf().subplots_adjust(wspace = 0.5) + + axDensity.plot( [xmin + i*dx for i in range(numsamples+1)], rho_field, label='Density, '+str(numsamples)+' cells') #new picture for video # Returns a tuple of line objects, thus the comma + axDensity.set(xlabel='x (m)', ylabel='Densité (Kg/m3)') + axDensity.set_xlim(xmin,xmax) + axDensity.set_ylim(0.9*min_initial_rho , 6*max_initial_rho ) + axDensity.legend() + + axMomentum.plot( [xmin + i*dx for i in range(numsamples+1)], q_field, label='Momentum, '+str(numsamples)+' cells') + axMomentum.set(xlabel='x (m)', ylabel='Momentum (Kg/m2/s)') + axMomentum.set_xlim(xmin,xmax) + axMomentum.set_ylim(0.9*min_initial_q , 2.5*max_initial_q ) + axMomentum.legend() + + axEnthalpie.plot([xmin + i*dx for i in range(numsamples+1)], h_field, label='Enthalpy, '+str(numsamples)+' cells') + axEnthalpie.set(xlabel='x (m)', ylabel='Enthalpy (J/Kg)') + axEnthalpie.set_xlim(xmin,xmax) + axEnthalpie.set_ylim(0.9*min_initial_h , 1.75*max_initial_h ) + axEnthalpie.legend() + + axPressure.plot( [xmin + i*dx for i in range(numsamples+1)], p_field, label='Pressure, '+str(numsamples)+' cells') + axPressure.set(xlabel='x (m)', ylabel='Pressure (Pa)') + axPressure.set_xlim(xmin,xmax) + axPressure.set_ylim(0.9*min_initial_p , 4*max_initial_p) + axPressure.ticklabel_format(axis='y', style='sci', scilimits=(0,0)) + axPressure.legend() + + axVitesse.plot( [xmin + i*dx for i in range(numsamples+1)], u_field, label=' Velocity, '+str(numsamples)+' cells') + axVitesse.set(xlabel='x (m)', ylabel='Velocity (m/s)') + axVitesse.set_xlim(xmin,xmax) + axVitesse.set_ylim(0.9*min_initial_v , 1.1*max_initial_v) + axVitesse.ticklabel_format(axis='y', style='sci', scilimits=(0,0)) + axVitesse.legend() + + axEinterne.plot( [xmin + i*dx for i in range(numsamples+1)], e_field, label='Internal energy, '+str(numsamples)+' cells') + axEinterne.set(xlabel='x (m)', ylabel='Internal energy (J/Kg)') + axEinterne.set_xlim(xmin,xmax) + axEinterne.set_ylim(0.9*min_initial_e , 1.75*max_initial_e) + axEinterne.ticklabel_format(axis='y', style='sci', scilimits=(0,0)) + axEinterne.legend() + + plt.suptitle('Euler equations : exact solution of Riemann problem ' + filename) + plt.savefig("EulerRiemannProblemExactSolution_" + filename + ".png") + + + return 0.0; + + +if __name__ == """__main__""": + run_Riemann_problems() diff --git a/CDMATH/tests/validation/ReferenceSolutions/plot_output.sh b/CDMATH/tests/validation/ReferenceSolutions/plot_output.sh new file mode 100755 index 0000000..4877bc0 --- /dev/null +++ b/CDMATH/tests/validation/ReferenceSolutions/plot_output.sh @@ -0,0 +1,7 @@ +#!/bin/bash + +for filepath in ./*.dat +do + filename=$(basename $filepath | cut -d "." -f1) + gnuplot -c ./epsplotsoln.gnu $filename +done diff --git a/CDMATH/tests/validation/scripts/ReferenceSolutions/exact_rs_stiffenedgas.py b/CDMATH/tests/validation/scripts/ReferenceSolutions/exact_rs_stiffenedgas.py new file mode 100644 index 0000000..0a60dce --- /dev/null +++ b/CDMATH/tests/validation/scripts/ReferenceSolutions/exact_rs_stiffenedgas.py @@ -0,0 +1,263 @@ +#!/usr/bin/env python3 +# -*-coding:utf-8 -* + +###################################################################################################################### +# This file contains a class to solve for the exact solution of the Riemann Problem for the one dimensional Euler +# equations with stiffened gas equation of state +# +# Author: Michael Ndjinga +# Date: 18/02/2021 +# Description : Translated from C++ package developped by Murray Cutforth +####################################################################################################################### + +from math import pow, fabs, sqrt + +def stiffenedgas_e (rho, p, gamma, pinf): + return (p+gamma*pinf)/(rho*(gamma-1)); + +def stiffenedgas_h (rho, p, gamma, pinf): + return gamma*(p+pinf)/(rho*(gamma-1)); + +def p_to_e_StiffenedGaz(p, rho, gamma, pinf): + e_field = (p + gamma*pinf) / (gamma - 1.) / rho + return e_field + +class exact_rs_stiffenedgas : + + def __init__(self, gamma_L, gamma_R, pinf_L, pinf_R, tol=1.e-6, max_iter=100): + self.TOL = tol + self.MAX_NB_ITER = max_iter + + self.gamma_L = gamma_L + self.gamma_R = gamma_R + self.pinf_L = pinf_L + self.pinf_R = pinf_R + + self.S_STAR = 0. + self.P_STAR = 0. + self.rho_star_L = 0. + self.rho_star_R = 0. + + self.S_L = 0. + self.S_R = 0. + self.S_HL = 0. + self.S_TL = 0. + self.S_HR = 0. + self.S_TR = 0. + + + + # Functions used to generate exact solutions to Riemann problems + + def solve_RP (self, W_L, W_R): + assert len(W_L) == 3, "Left state should have three components (rho, u p)" + assert len(W_R) == 3, "Right state should have three components (rho, u p)" + assert W_L[0] >= 0.0, "Left density should be positive" + assert W_R[0] >= 0.0, "Right density should be positive" + # assert W_L[2] >= 0.0 # Since stiffened gases will often exhibit p<0.. + # assert W_R[2] >= 0.0 # + + print("") + print("Solving Riemann problem for left state W_L=", W_L, ", and right state W_R=",W_R) + + # Calculate p_star + + self.P_STAR = self.find_p_star_newtonraphson(W_L[0], W_L[1], W_L[2], W_R[0], W_R[1], W_R[2]) + + + # Calculate u_star + + self.S_STAR = 0.5*(W_L[1]+W_R[1]) + 0.5*(self.f(self.P_STAR,W_R[0],W_R[2],self.gamma_R,self.pinf_R) - self.f(self.P_STAR,W_L[0],W_L[2],self.gamma_L,self.pinf_L)) + + + # Solution now depends on character of 1st and 3rd waves + + if (self.P_STAR > W_L[2]): + # Left shock + + self.rho_star_L = W_L[0]*((2.0*self.gamma_L*self.pinf_L + (self.gamma_L+1.0)*self.P_STAR + (self.gamma_L-1.0)*W_L[2])/(2.0*(W_L[2] + self.gamma_L*self.pinf_L) + (self.gamma_L-1.0)*self.P_STAR + (self.gamma_L-1.0)*W_L[2])) + self.S_L = W_L[1] - (self.Q_K(self.P_STAR,W_L[0],W_L[2],self.gamma_L,self.pinf_L)/W_L[0]) + else: + # Left rarefaction + + self.rho_star_L = W_L[0]*pow((self.P_STAR + self.pinf_L)/(W_L[2] + self.pinf_L), 1.0/self.gamma_L) + + a_L = self.a(W_L[0], W_L[2], self.gamma_L, self.pinf_L) + a_star_L = a_L*pow((self.P_STAR + self.pinf_L)/(W_L[2] + self.pinf_L), (self.gamma_L-1.0)/(2.0*self.gamma_L)) + + self.S_HL = W_L[1] - a_L + self.S_TL = self.S_STAR - a_star_L + + if (self.P_STAR > W_R[2]): + # Right shock + + self.rho_star_R = W_R[0]*((2.0*self.gamma_R*self.pinf_R + (self.gamma_R+1.0)*self.P_STAR + (self.gamma_R-1.0)*W_R[2])/(2.0*(W_R[2] + self.gamma_R*self.pinf_R) + (self.gamma_R-1.0)*self.P_STAR + (self.gamma_R-1.0)*W_R[2])) + + self.S_R = W_R[1] + (self.Q_K(self.P_STAR,W_R[0],W_R[2],self.gamma_R,self.pinf_R)/W_R[0]) + else: + # Right rarefaction + + self.rho_star_R = W_R[0]*pow((self.P_STAR + self.pinf_R)/(W_R[2] + self.pinf_R), 1.0/self.gamma_R) + + a_R = self.a(W_R[0],W_R[2],self.gamma_R, self.pinf_R) + a_star_R = a_R*pow((self.P_STAR + self.pinf_R)/(W_R[2] + self.pinf_R), (self.gamma_R-1.0)/(2.0*self.gamma_R)) + + self.S_HR = W_R[1] + a_R + self.S_TR = self.S_STAR + a_star_R + + def sample_solution (self, W_L, W_R, S): + W = [0.]*3 + + # Find appropriate part of solution and return primitives + + if (S < self.S_STAR): + # To the left of the contact + + if (self.P_STAR > W_L[2]): + # Left shock + + if (S < self.S_L): + W = W_L + else: + W[0] = self.rho_star_L + W[1] = self.S_STAR + W[2] = self.P_STAR + else: + # Left rarefaction + + if (S < self.S_HL): + W = W_L + else: + if (S > self.S_TL): + W[0] = self.rho_star_L + W[1] = self.S_STAR + W[2] = self.P_STAR + else: + self.set_left_rarefaction_fan_state(W_L, S, W) + else: + # To the right of the contact + + if (self.P_STAR > W_R[2]): + # Right shock + + if (S > self.S_R): + W = W_R + else: + W[0] = self.rho_star_R + W[1] = self.S_STAR + W[2] = self.P_STAR + else: + # Right rarefaction + + if (S > self.S_HR): + W = W_R + else: + if (S < self.S_TR): + W[0] = self.rho_star_R + W[1] = self.S_STAR + W[2] = self.P_STAR + else: + self.set_right_rarefaction_fan_state(W_R, S, W) + + return W + + # Functions used to solve for p_star iteratively + + def find_p_star_newtonraphson (self, rho_L, u_L, p_L, rho_R, u_R, p_R ): + + # First we set the initial guess for p_star using a simple mean-value approximation + + p_star_next = 0.5*(p_L+p_R) + n = 0 + + + # Now use the Newton-Raphson algorithm + + while True:#conversion of do ... while by while True... if (...) break + p_star = p_star_next + + p_star_next = p_star - self.total_pressure_function(p_star,rho_L,u_L,p_L,rho_R,u_R,p_R)/self.total_pressure_function_deriv(p_star,rho_L,p_L,rho_R,p_R) + + p_star_next = max(p_star_next, self.TOL) + + n+=1 + + if not ((fabs(p_star_next - p_star)/(0.5*(p_star+p_star_next)) > self.TOL) and n < self.MAX_NB_ITER): + break + + if (n == self.MAX_NB_ITER): + raise ValueError("!!!!!!!!!!Newton algorithm did not converge. Increase tolerance or maximum number of time steps. Current values : tol=" + str(self.TOL) + ", max_iter=" + str(self.MAX_NB_ITER) ) + #p_star_next = 0.5*(p_L+p_R) + + return p_star_next + + def total_pressure_function (self, p_star, rho_L, u_L, p_L, rho_R, u_R, p_R ): + + return self.f(p_star, rho_L, p_L, self.gamma_L, self.pinf_L) + self.f(p_star, rho_R, p_R, self.gamma_R, self.pinf_R) + u_R - u_L + + def total_pressure_function_deriv (self, p_star, rho_L, p_L, rho_R, p_R ): + + return self.f_deriv (p_star, rho_L, p_L, self.gamma_L, self.pinf_L) + self.f_deriv (p_star, rho_R, p_R, self.gamma_R, self.pinf_R) + + + def f (self, p_star, rho, p, gamma, pinf): + if (p_star > p): + + return (p_star - p)/self.Q_K(p_star, rho, p, gamma, pinf) + + else: + + return (2.0*self.a(rho,p,gamma,pinf)/(gamma-1.0))*(pow((p_star + pinf)/(p + pinf), (gamma-1.0)/(2.0*gamma)) - 1.0) + + + def f_deriv (self, p_star, rho, p, gamma, pinf): + A = 2.0/((gamma+1.0)*rho) + B = (p+pinf)*(gamma-1.0)/(gamma+1.0) + + if (p_star > p): + + return sqrt(A/(B+p_star+pinf))*(1.0 - ((p_star-p)/(2.0*(B+p_star+pinf)))) + + else: + + return (1.0/(rho*self.a(rho,p,gamma,pinf)))*pow((p_star+pinf)/(p+pinf), -(gamma+1.0)/(2.0*gamma)) + + + + # Functions to find the state inside a rarefaction fan + + def set_left_rarefaction_fan_state (self, W_L, S, W): + a_L = self.a(W_L[0],W_L[2],self.gamma_L,self.pinf_L) + W[0] = W_L[0]*pow((2.0/(self.gamma_L+1.0)) + ((self.gamma_L-1.0)/(a_L*(self.gamma_L+1.0)))*(W_L[1] - S), 2.0/(self.gamma_L - 1.0)) + W[1] = (2.0/(self.gamma_L+1.0))*(a_L + S + ((self.gamma_L-1.0)/2.0)*W_L[1]) + W[2] = (W_L[2] + self.pinf_L)*pow((2.0/(self.gamma_L+1.0)) + ((self.gamma_L-1.0)/(a_L*(self.gamma_L+1.0)))*(W_L[1] - S), (2.0*self.gamma_L)/(self.gamma_L-1.0)) - self.pinf_L + + def set_right_rarefaction_fan_state (self, W_R, S, W): + a_R = self.a(W_R[0],W_R[2],self.gamma_R,self.pinf_R) + W[0] = W_R[0]*pow((2.0/(self.gamma_R+1.0)) - ((self.gamma_R-1.0)/(a_R*(self.gamma_R+1.0)))*(W_R[1] - S), 2.0/(self.gamma_R - 1.0)) + W[1] = (2.0/(self.gamma_R+1.0))*(- a_R + S + ((self.gamma_R-1.0)/2.0)*W_R[1]) + W[2] = (W_R[2] + self.pinf_R)*pow((2.0/(self.gamma_R+1.0)) - ((self.gamma_R-1.0)/(a_R*(self.gamma_R+1.0)))*(W_R[1] - S), (2.0*self.gamma_R)/(self.gamma_R-1.0)) - self.pinf_R + + + + # Misc functions + + def Q_K (self, p_star, rho, p, gamma, pinf): + A = 2.0/((gamma+1.0)*rho) + B = (p+pinf)*(gamma-1.0)/(gamma+1.0) + return sqrt((p_star+pinf+B)/A) + + + + # Equation of state functions + + def a (self, rho, p, gamma, pinf):#sound speed + return sqrt(gamma*((p+pinf)/rho)) + + #Determine the solution value at position x and time t + def rho_u_p_solution (initialLeftState, initialRightState, x, t, gamma, pinf, offset=0): + RS = exact_rs_stiffenedgas(gamma, gamma, pinf, pinf); + RS.solve_RP(initialLeftState, initialRightState); + return RS.sample_solution(initialLeftState, initialRightState, (x - offset)/t); + +