]> SALOME platform Git repositories - tools/solverlab.git/commitdiff
Salome HOME
Added validation tests for the Reference solutions
authormichael <michael@localhost.localdomain>
Sat, 9 Oct 2021 22:02:53 +0000 (00:02 +0200)
committermichael <michael@localhost.localdomain>
Sat, 9 Oct 2021 22:02:53 +0000 (00:02 +0200)
CDMATH/tests/validation/ReferenceSolutions/CMakeLists.txt [new file with mode: 0644]
CDMATH/tests/validation/ReferenceSolutions/ExactRiemannProblemSolutions.py [new file with mode: 0644]
CDMATH/tests/validation/ReferenceSolutions/plot_output.sh [new file with mode: 0755]
CDMATH/tests/validation/scripts/ReferenceSolutions/exact_rs_stiffenedgas.py [new file with mode: 0644]

diff --git a/CDMATH/tests/validation/ReferenceSolutions/CMakeLists.txt b/CDMATH/tests/validation/ReferenceSolutions/CMakeLists.txt
new file mode 100644 (file)
index 0000000..37ce57a
--- /dev/null
@@ -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 (file)
index 0000000..b1ad296
--- /dev/null
@@ -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 (executable)
index 0000000..4877bc0
--- /dev/null
@@ -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 (file)
index 0000000..0a60dce
--- /dev/null
@@ -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);
+
+