Blog: Meeting 3/24: solvers.py

File solvers.py, 2.8 KB (added by smurugan, 11 years ago)
Line 
1import numpy as np
2from problem import *
3import math
4import ERS
5
6
7def find_timestep(U, t, tf):
8 smax = 0
9 for j in range(npts):
10 wave_vel = math.fabs(U[j,1]) + (gamma*U[j,2]/U[j,0])**(.5)
11 if wave_vel > smax:
12 smax = wave_vel
13 dt = CFL*dx/(smax)
14 if dt > (tf -t):
15 dt = tf - t
16 return dt
17
18def FindStar(U,k, flag):
19 if flag == '-1' or flag == '-2':
20 return AdaptiveSolver(U,k, flag)
21 else:
22 return ERS.ExactSolve(U,k)
23
24
25def getUpdate(U, t, tf, flag):
26 u = np.zeros(npts)
27 p = np.zeros(npts)
28 rho = np.zeros(npts)
29 while t < tf:
30 for k in range(npts-1):
31 p[k], rho[k], u[k] = FindStar(U,k, flag)
32 dt = find_timestep(U, t, tf)
33 for i in range(npts-1):
34 Cons = np.array([U[i,0], U[i,0]*U[i,1], U[i,2]/(gamma -1) + .5*U[i,0]*U[i,1]**2])
35 if i > 0:
36 FluxL = ERS.getFlux(rho[i -1], u[i-1], p[i-1])
37 FluxR = ERS.getFlux(rho[i], u[i], p[i])
38 for j in range(3):
39 Cons[j] += (dt/dx)*(FluxL[j] - FluxR[j])
40 U[i,0] = Cons[0]
41 U[i,1] = Cons[1]/Cons[0]
42 U[i,2] = (Cons[2] -(.5*Cons[1]**2)/Cons[0])*(gamma -1)
43 t += dt
44 print t
45 return U, t
46
47def g(U,k,p):
48 A = 2/((gamma+1)*U[k,0])
49 B = (gamma - 1)/(gamma + 1) *U[k,2]
50 return (A/(p +B))**(.5)
51
52def TSrho(U,k,p):
53 return U[k,1]* ( (p/U[k,2]) + (gamma-1)/(gamma +1))/(((gamma -1)*p/((gamma+1)*U[k,2])) +1)
54
55def TSRS(U,k, p0):
56 p_num = g(U, k, p0)*U[k, 2] + g(U, k+ 1, p0)*U[k+1, 2] + U[k, 1] - U[k+1, 1]
57 p_den = g(U, k ,p0) + g(U, k+1, p0)
58 p = p_num/p_den
59 u = .5*(U[k,1] + U[k+1,1]) + .5*((p-U[k+1,2])*g(U,k+1,p0) -(p-U[k,2])*g(U,k,p0))
60 if u > 0:
61 rho = TSrho(U, k, p)
62 elif u < 0:
63 rho = TSrho(U,k+1, p)
64 else:
65 rho = U[k,0]
66 return p, rho, u
67
68def TRRS(U,k):
69 z = (gamma -1)/(2*gamma)
70 aR = gamma*(U[k,2]/U[k,0])**(.5)
71 aL = gamma*(U[k+1,2]/U[k+1,0])**(.5)
72 p_num = aL + aR - .5*(gamma -1)*(U[k+1,1] - U[k,1])
73 p_den = aL/(U[k,2]**z) + aR/(U[k+1,2]**z)
74 p = (p_num/p_den)**(1/z)
75 u = U[k, 1] - (2*aL/(gamma -1))* ( (p/U[k, 2])**z -1)
76 if u > 0:
77 rho = U[k,0] * (p/U[k,2])**(1/gamma)
78 elif u < 0:
79 rho = U[k+1 ,0] * (p/U[k +1,2])**(1/gamma)
80 else:
81 rho = U[k, 0]
82 return p, rho, u
83
84def AdaptiveSolver(U,k, flag):
85 pmin = min(U[k,2],U[k+1,2])
86 pmax = max(U[k,2],U[k+1,2])
87 if pmin == pmax:
88 Q = 1
89 else:
90 Q = pmin/pmax
91 if Q < Q_TOL:
92 rhobar = .5*(U[k,0] + U[k+1, 0])
93 abar = .5*((U[k,2]/U[k,0])**(.5) + (U[k+1, 2]/U[k+1,0])**(.5))*gamma**(.5)
94 p = .5*(U[k, 2] + U[k+1, 2]) + .5*(U[k,1] - U[k+1,1])*(rhobar*abar)
95 if p < pmax and p > pmin:
96 u = .5*(U[k, 1] + U[k+1, 1]) + .5*(U[k,2] - U[k+1,2])/(rhobar*abar)
97 if u > 0:
98 rho = U[k, 0] + (U[k, 1] - u)*(rhobar/abar)
99 elif u < 0:
100 rho = U[k+1, 0] + (u - U[k+1, 1])*(rhobar/abar)
101 else:
102 rho = U[k,0]
103 return p, rho, u
104 else:
105 if flag == '-1':
106 return ERS.ExactSolve(U,k)
107 else:
108 if p <= pmin:
109 return TRRS(U,k)
110 else:
111 p0 = max(0, p)
112 return TSRS(U,k, p0)
113 else:
114 return ERS.ExactSolve(U,k)
115
116