############################################################################### # This file is part of SWIFT. # Coypright (c) 2012 Pedro Gonnet (pedro.gonnet@durham.ac.uk), # Matthieu Schaller (matthieu.schaller@durham.ac.uk) # # This program is free software: you can redistribute it and/or modify # it under the terms of the GNU Lesser General Public License as published # by the Free Software Foundation, either version 3 of the License, or # (at your option) any later version. # # This program is distributed in the hope that it will be useful, # but WITHOUT ANY WARRANTY; without even the implied warranty of # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the # GNU General Public License for more details. # # You should have received a copy of the GNU Lesser General Public License # along with this program. If not, see <http://www.gnu.org/licenses/>. # ############################################################################## import random from numpy import * # Generates the analytical solution for the Sod shock test case # The script works for a given left (x<0) and right (x>0) state and computes the solution at a later time t. # The code writes five files rho.dat, P.dat, v.dat, u.dat and s.dat with the density, pressure, internal energy and # entropic function on N points between x_min and x_max. # This follows the solution given in (Toro, 2009) # Parameters rho_L = 1 P_L = 1 v_L = 0. rho_R = 0.25 P_R = 0.1795 v_R = 0. gamma = 5./3. # Polytropic index t = 0.12 # Time of the evolution N = 1000 # Number of points x_min = -0.25 x_max = 0.25 # --------------------------------------------------------------- # Don't touch anything after this. # --------------------------------------------------------------- c_L = sqrt(gamma * P_L / rho_L) # Speed of the rarefaction wave c_R = sqrt(gamma * P_R / rho_R) # Speed of the shock front # Helpful variable Gama = (gamma - 1.) / (gamma + 1.) beta = (gamma - 1.) / (2. * gamma) # Characteristic function and its derivative, following Toro (2009) def compute_f(P_3, P, c): u = P_3 / P if u > 1: term1 = gamma*((gamma+1.)*u + gamma-1.) term2 = sqrt(2./term1) fp = (u - 1.)*c*term2 dfdp = c*term2/P + (u - 1.)*c/term2*(-1./term1**2)*gamma*(gamma+1.)/P else: fp = (u**beta - 1.)*(2.*c/(gamma-1.)) dfdp = 2.*c/(gamma-1.)*beta*u**(beta-1.)/P return (fp, dfdp) # Solution of the Riemann problem following Toro (2009) def RiemannProblem(rho_L, P_L, v_L, rho_R, P_R, v_R): P_new = ((c_L + c_R + (v_L - v_R)*0.5*(gamma-1.))/(c_L / P_L**beta + c_R / P_R**beta))**(1./beta) P_3 = 0.5*(P_R + P_L) f_L = 1. while fabs(P_3 - P_new) > 1e-6: P_3 = P_new (f_L, dfdp_L) = compute_f(P_3, P_L, c_L) (f_R, dfdp_R) = compute_f(P_3, P_R, c_R) f = f_L + f_R + (v_R - v_L) df = dfdp_L + dfdp_R dp = -f/df prnew = P_3 + dp v_3 = v_L - f_L return (P_new, v_3) # Solve Riemann problem for post-shock region (P_3, v_3) = RiemannProblem(rho_L, P_L, v_L, rho_R, P_R, v_R) # Check direction of shocks and wave shock_R = (P_3 > P_R) shock_L = (P_3 > P_L) # Velocity of shock front and and rarefaction wave if shock_R: v_right = v_R + c_R**2*(P_3/P_R - 1.)/(gamma*(v_3-v_R)) else: v_right = c_R + 0.5*(gamma+1.)*v_3 - 0.5*(gamma-1.)*v_R if shock_L: v_left = v_L + c_L**2*(P_3/p_L - 1.)/(gamma*(v_3-v_L)) else: v_left = c_L - 0.5*(gamma+1.)*v_3 + 0.5*(gamma-1.)*v_L # Compute position of the transitions x_23 = -fabs(v_left) * t if shock_L : x_12 = -fabs(v_left) * t else: x_12 = -(c_L - v_L) * t x_34 = v_3 * t x_45 = fabs(v_right) * t if shock_R: x_56 = fabs(v_right) * t else: x_56 = (c_R + v_R) * t # Prepare arrays delta_x = (x_max - x_min) / N x_s = arange(x_min, x_max, delta_x) rho_s = zeros(N) P_s = zeros(N) v_s = zeros(N) # Compute solution in the different regions for i in range(N): if x_s[i] <= x_12: rho_s[i] = rho_L P_s[i] = P_L v_s[i] = v_L if x_s[i] >= x_12 and x_s[i] < x_23: if shock_L: rho_s[i] = rho_L*(Gama + P_3/P_L)/(1. + Gama * P_3/P_L) P_s[i] = P_3 v_s[i] = v_3 else: rho_s[i] = rho_L*(Gama * (0. - x_s[i])/(c_L * t) + Gama * v_L/c_L + (1.-Gama))**(2./(gamma-1.)) P_s[i] = P_L*(rho_s[i] / rho_L)**gamma v_s[i] = (1.-Gama)*(c_L -(0. - x_s[i]) / t) + Gama*v_L if x_s[i] >= x_23 and x_s[i] < x_34: if shock_L: rho_s[i] = rho_L*(Gama + P_3/P_L)/(1+Gama * P_3/p_L) else: rho_s[i] = rho_L*(P_3 / P_L)**(1./gamma) P_s[i] = P_3 v_s[i] = v_3 if x_s[i] >= x_34 and x_s[i] < x_45: if shock_R: rho_s[i] = rho_R*(Gama + P_3/P_R)/(1. + Gama * P_3/P_R) else: rho_s[i] = rho_R*(P_3 / P_R)**(1./gamma) P_s[i] = P_3 v_s[i] = v_3 if x_s[i] >= x_45 and x_s[i] < x_56: if shock_R: rho_s[i] = rho_R P_s[i] = P_R v_s[i] = v_R else: rho_s[i] = rho_R*(Gama*(x_s[i])/(c_R*t) - Gama*v_R/c_R + (1.-Gama))**(2./(gamma-1.)) P_s[i] = p_R*(rho_s[i]/rho_R)**gamma v_s[i] = (1.-Gama)*(-c_R - (-x_s[i])/t) + Gama*v_R if x_s[i] >= x_56: rho_s[i] = rho_R P_s[i] = P_R v_s[i] = v_R # Additional arrays u_s = P_s / (rho_s * (gamma - 1.)) #internal energy s_s = P_s / rho_s**gamma # entropic function #--------------------------------------------------------------- # Print arrays savetxt("rho.dat", column_stack((x_s, rho_s))) savetxt("P.dat", column_stack((x_s, P_s))) savetxt("v.dat", column_stack((x_s, v_s))) savetxt("u.dat", column_stack((x_s, u_s))) savetxt("s.dat", column_stack((x_s, s_s)))