CDMATH/tests/validation/scripts/ReferenceSolutions/exact_rs_stiffenedgas.py

   1 #!/usr/bin/env python3
   2 # -*-coding:utf-8 -*
   3
   4 ######################################################################################################################
   5 #       This file contains a class to solve for the exact solution of the Riemann Problem for the one dimensional Euler
   6 #       equations with stiffened gas equation of state
   7 #
   8 #       Author: Michael Ndjinga
   9 #       Date:   18/02/2021
  10 #   Description : Translated from C++ package developped by Murray Cutforth
  11 #######################################################################################################################
  12
  13 from math import pow, fabs, sqrt
  14
  15 def stiffenedgas_e (rho, p, gamma, pinf):
  16         return (p+gamma*pinf)/(rho*(gamma-1));
  17
  18 def stiffenedgas_h (rho, p, gamma, pinf):
  19         return gamma*(p+pinf)/(rho*(gamma-1));
  20
  21 def p_to_e_StiffenedGaz(p, rho, gamma, pinf):
  22         e_field = (p + gamma*pinf) / (gamma - 1.) / rho
  23         return e_field
  24
  25 class exact_rs_stiffenedgas :
  26
  27         def __init__(self, gamma_L, gamma_R, pinf_L, pinf_R, tol=1.e-6, max_iter=100):
  28                 self.TOL = tol
  29                 self.MAX_NB_ITER = max_iter
  30
  31                 self.gamma_L = gamma_L
  32                 self.gamma_R = gamma_R
  33                 self.pinf_L  = pinf_L
  34                 self.pinf_R  = pinf_R
  35
  36                 self.S_STAR = 0.
  37                 self.P_STAR = 0.
  38                 self.rho_star_L = 0.
  39                 self.rho_star_R = 0.
  40
  41                 self.S_L = 0.
  42                 self.S_R = 0.
  43                 self.S_HL = 0.
  44                 self.S_TL = 0.
  45                 self.S_HR = 0.
  46                 self.S_TR = 0.
  47
  48
  49
  50         # Functions used to generate exact solutions to Riemann problems
  51
  52         def solve_RP (self, W_L, W_R):
  53                 assert len(W_L) == 3, "Left state should have three components (rho, u p)"
  54                 assert len(W_R) == 3, "Right state should have three components (rho, u p)"
  55                 assert W_L[0] >= 0.0, "Left density should be positive"
  56                 assert W_R[0] >= 0.0, "Right density should be positive"
  57                 # assert W_L[2] >= 0.0 # Since stiffened gases will often exhibit p<0..
  58                 # assert W_R[2] >= 0.0 #
  59
  60                 print("")
  61                 print("Solving Riemann problem for left state W_L=", W_L, ", and right state W_R=",W_R)
  62
  63                 # Calculate p_star
  64
  65                 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])
  66
  67
  68                 # Calculate u_star
  69
  70                 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))
  71
  72
  73                 # Solution now depends on character of 1st and 3rd waves
  74
  75                 if (self.P_STAR > W_L[2]):
  76                         # Left shock
  77
  78                         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]))
  79                         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])
  80                 else:
  81                         # Left rarefaction
  82
  83                         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)
  84
  85                         a_L = self.a(W_L[0], W_L[2], self.gamma_L, self.pinf_L)
  86                         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))
  87
  88                         self.S_HL = W_L[1] - a_L
  89                         self.S_TL = self.S_STAR - a_star_L
  90
  91                 if (self.P_STAR > W_R[2]):
  92                         # Right shock
  93
  94                         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]))
  95
  96                         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])
  97                 else:
  98                         # Right rarefaction
  99
 100                         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)
 101
 102                         a_R = self.a(W_R[0],W_R[2],self.gamma_R, self.pinf_R)
 103                         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))
 104
 105                         self.S_HR = W_R[1] + a_R
 106                         self.S_TR = self.S_STAR + a_star_R
 107
 108         def sample_solution (self, W_L, W_R, S):
 109                 W = [0.]*3
 110
 111                 # Find appropriate part of solution and return primitives
 112
 113                 if (S < self.S_STAR):
 114                         # To the left of the contact
 115
 116                         if (self.P_STAR > W_L[2]):
 117                                 # Left shock
 118
 119                                 if (S < self.S_L):
 120                                         W = W_L
 121                                 else:
 122                                         W[0] = self.rho_star_L
 123                                         W[1] = self.S_STAR
 124                                         W[2] = self.P_STAR
 125                         else:
 126                                 # Left rarefaction
 127
 128                                 if (S < self.S_HL):
 129                                         W = W_L
 130                                 else:
 131                                         if (S > self.S_TL):
 132                                                 W[0] = self.rho_star_L
 133                                                 W[1] = self.S_STAR
 134                                                 W[2] = self.P_STAR
 135                                         else:
 136                                                 self.set_left_rarefaction_fan_state(W_L, S, W)
 137                 else:
 138                         # To the right of the contact
 139
 140                         if (self.P_STAR > W_R[2]):
 141                                 # Right shock
 142
 143                                 if (S > self.S_R):
 144                                         W = W_R
 145                                 else:
 146                                         W[0] = self.rho_star_R
 147                                         W[1] = self.S_STAR
 148                                         W[2] = self.P_STAR
 149                         else:
 150                                 # Right rarefaction
 151
 152                                 if (S > self.S_HR):
 153                                         W = W_R
 154                                 else:
 155                                         if (S < self.S_TR):
 156                                                 W[0] = self.rho_star_R
 157                                                 W[1] = self.S_STAR
 158                                                 W[2] = self.P_STAR
 159                                         else:
 160                                                 self.set_right_rarefaction_fan_state(W_R, S, W)
 161
 162                 return W
 163
 164         # Functions used to solve for p_star iteratively
 165
 166         def find_p_star_newtonraphson (self, rho_L, u_L, p_L, rho_R, u_R, p_R ):
 167
 168                 # First we set the initial guess for p_star using a simple mean-value approximation
 169
 170                 p_star_next = 0.5*(p_L+p_R)
 171                 n = 0
 172
 173
 174                 # Now use the Newton-Raphson algorithm
 175
 176                 while True:#conversion of do ... while by while True... if (...) break
 177                         p_star = p_star_next
 178
 179                         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)
 180
 181                         p_star_next = max(p_star_next, self.TOL)
 182
 183                         n+=1
 184
 185                         if not ((fabs(p_star_next - p_star)/(0.5*(p_star+p_star_next)) > self.TOL) and n < self.MAX_NB_ITER):
 186                                 break
 187
 188                 if (n == self.MAX_NB_ITER):
 189                         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) )
 190                         #p_star_next = 0.5*(p_L+p_R)
 191
 192                 return p_star_next
 193
 194         def total_pressure_function (self, p_star, rho_L, u_L, p_L, rho_R, u_R, p_R ):
 195
 196                 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
 197
 198         def total_pressure_function_deriv (self, p_star, rho_L, p_L, rho_R, p_R ):
 199
 200                 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)
 201
 202
 203         def f (self, p_star, rho, p, gamma, pinf):
 204                 if (p_star > p):
 205
 206                         return (p_star - p)/self.Q_K(p_star, rho, p, gamma, pinf)
 207
 208                 else:
 209
 210                         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)
 211
 212
 213         def f_deriv (self, p_star, rho, p, gamma, pinf):
 214                 A = 2.0/((gamma+1.0)*rho)
 215                 B = (p+pinf)*(gamma-1.0)/(gamma+1.0)
 216
 217                 if (p_star > p):
 218
 219                         return sqrt(A/(B+p_star+pinf))*(1.0 - ((p_star-p)/(2.0*(B+p_star+pinf))))
 220
 221                 else:
 222
 223                         return (1.0/(rho*self.a(rho,p,gamma,pinf)))*pow((p_star+pinf)/(p+pinf), -(gamma+1.0)/(2.0*gamma))
 224
 225
 226
 227         # Functions to find the state inside a rarefaction fan
 228
 229         def set_left_rarefaction_fan_state (self, W_L, S, W):
 230                 a_L = self.a(W_L[0],W_L[2],self.gamma_L,self.pinf_L)
 231                 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))
 232                 W[1] = (2.0/(self.gamma_L+1.0))*(a_L + S + ((self.gamma_L-1.0)/2.0)*W_L[1])
 233                 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
 234
 235         def set_right_rarefaction_fan_state (self, W_R, S, W):
 236                 a_R = self.a(W_R[0],W_R[2],self.gamma_R,self.pinf_R)
 237                 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))
 238                 W[1] = (2.0/(self.gamma_R+1.0))*(- a_R + S + ((self.gamma_R-1.0)/2.0)*W_R[1])
 239                 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
 240
 241
 242
 243         # Misc functions
 244
 245         def Q_K (self, p_star, rho, p, gamma, pinf):
 246                 A = 2.0/((gamma+1.0)*rho)
 247                 B = (p+pinf)*(gamma-1.0)/(gamma+1.0)
 248                 return sqrt((p_star+pinf+B)/A)
 249
 250
 251
 252         # Equation of state functions
 253
 254         def a (self, rho, p, gamma, pinf):#sound speed
 255                 return sqrt(gamma*((p+pinf)/rho))
 256
 257         #Determine the solution value at position x and time t
 258         def rho_u_p_solution (initialLeftState, initialRightState, x, t, gamma, pinf, offset=0):
 259                 RS = exact_rs_stiffenedgas(gamma, gamma, pinf, pinf);
 260                 RS.solve_RP(initialLeftState, initialRightState);
 261                 return  RS.sample_solution(initialLeftState, initialRightState, (x - offset)/t);
 262
 263