solutions.f90
Go to the documentation of this file.
1!#############################################################################
2!# #
3!# fosite - 3D hydrodynamical simulation program #
4!# module: solutions.f90 #
5!# #
6!# Copyright (C) 2013 Manuel Jung <mjung@astrophysik.uni-kiel.de> #
7!# #
8!# This program is free software; you can redistribute it and/or modify #
9!# it under the terms of the GNU General Public License as published by #
10!# the Free Software Foundation; either version 2 of the License, or (at #
11!# your option) any later version. #
12!# #
13!# This program is distributed in the hope that it will be useful, but #
14!# WITHOUT ANY WARRANTY; without even the implied warranty of #
15!# MERCHANTABILITY OR FITNESS FOR A PARTICULAR PURPOSE, GOOD TITLE or #
16!# NON INFRINGEMENT. See the GNU General Public License for more #
17!# details. #
18!# #
19!# You should have received a copy of the GNU General Public License #
20!# along with this program; if not, write to the Free Software #
21!# Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA. #
22!# #
23!#############################################################################
24
25!----------------------------------------------------------------------------!
31!----------------------------------------------------------------------------!
32
33MODULE solutions
34 USE roots, ONLY: getroot, getroot_newton
35 IMPLICIT NONE
36 !--------------------------------------------------------------------------!
37 PRIVATE
38 REAL :: rho_l,u_l,p_l,a_l,b_l,c_l,&
39 rho_r,u_r,p_r,a_r,b_r,c_r,gamma
40 !--------------------------------------------------------------------------!
41 PUBLIC :: &
42 riemann
43 !--------------------------------------------------------------------------!
44CONTAINS
45 SUBROUTINE riemann(x0,gamma_,rho_l_,u_l_,p_l_,rho_r_,u_r_,p_r_,t,x,pvar)
46 IMPLICIT NONE
47 !------------------------------------------------------------------------!
48 REAL :: gamma_,t,x0,rho_l_,u_l_,p_l_,rho_r_,u_r_,p_r_
49 REAL,DIMENSION(:) :: x
50 REAL,DIMENSION(:,:) :: pvar
51 !------------------------------------------------------------------------!
52 INTEGER :: i
53 Real :: rho_lstar, rho_rstar, p_star, u_star, &
54 c_lstar, c_rstar, &
55 s_l, s_hl, s_tl, s_r, s_hr, s_tr
56 INTEGER :: error
57 !------------------------------------------------------------------------!
58 gamma = gamma_ !1.4 !heat capacity ratio
59
60 rho_l = rho_l_
61 rho_r = rho_r_
62
63 u_l = u_l_
64 u_r = u_r_
65
66 p_l = p_l_
67 p_r = p_r_
68
69 IF (t.LE.0) RETURN
70
71 ! sound speeds of far regions
72 c_l = sqrt(gamma*p_l/rho_l)
73 c_r = sqrt(gamma*p_r/rho_r)
74
75 a_l = 2.0/((gamma+1.0)*rho_l)
76 a_r = 2.0/((gamma+1.0)*rho_r)
77
78 b_l = p_l*(gamma-1.0)/(gamma+1.0)
79 b_r = p_r*(gamma-1.0)/(gamma+1.0)
80
81 ! calculate p_star by findung the root of function f
82 ! upper limit is just a somewhat big number. Is this really always
83 ! sufficient?
84 CALL getroot(f, 0., 2.e+3, p_star, error)
85 !p_star = GetRoot_andbjo(f, 0., 1.E+4,1.E-6)
86
87 u_star = 0.5*((u_l + u_r) + ( f_x(p_star, rho_r, u_r, p_r, a_r, b_r, c_r) &
88 - f_x(p_star, rho_l, u_l, p_l, a_l, b_l, c_l)))
89
90 IF (p_star > p_l) THEN
91 !left shock
92 rho_lstar = rho_l * ( ((p_star/p_l) + (gamma-1.0)/(gamma+1.0)) &
93 / (((gamma-1.0)/(gamma+1.0))*(p_star/p_l) + 1.0))
94 ELSE
95 !left rarefaction
96 rho_lstar = rho_l * (p_star/p_l)**(1.0/gamma)
97 END IF
98
99 IF (p_star > p_r) THEN
100 !right shock
101 rho_rstar = rho_r * ( ((p_star/p_r) + (gamma-1.0)/(gamma+1.0)) &
102 /(((gamma-1.0)/(gamma+1.0))*(p_star/p_r) + 1.0))
103 ELSE
104 !right rarefaction
105 rho_rstar = rho_r * (p_star/p_r)**(1.0/gamma)
106 END IF
107
108 ! sound speeds of star regions
109 c_lstar = c_l * (p_star/p_l)**((gamma-1.0)/(2.0*gamma))
110 c_rstar = c_r * (p_star/p_r)**((gamma-1.0)/(2.0*gamma))
111
112 ! shock speeds
113 s_l = u_l - c_l * sqrt( ((gamma+1.0)/(2.0*gamma))*(p_star/p_l) + (gamma-1.0)/(2.0*gamma))
114 s_r = u_r + c_r * sqrt( ((gamma+1.0)/(2.0*gamma))*(p_star/p_r) + (gamma-1.0)/(2.0*gamma))
115
116 ! head speed of rarefraction waves
117 s_hl = u_l - c_l
118 s_hr = u_r + c_r
119 ! tail speed of rarefraction waves
120 s_tl = u_star - c_lstar
121 s_tr = u_star + c_rstar
122
123 DO i=1,SIZE(x)
124 CALL sample(x(i),pvar(i,1),pvar(i,2),pvar(i,3))
125 END DO
126
127 CONTAINS
128 SUBROUTINE sample(x,rho,u,p)
129 IMPLICIT NONE
130 !----------------------------------------------------------------------!
131 REAL :: x, rho, u, p
132 !----------------------------------------------------------------------!
133 REAL :: S
134 !----------------------------------------------------------------------!
135 INTENT(IN) :: x
136 INTENT(OUT) :: rho, u, p
137 !----------------------------------------------------------------------!
138
139 !Sampling the solution
140 s = (x-x0)/t
141 IF (s <= u_star) THEN
142 !left side of contact
143 IF (p_star > p_l) THEN
144 !left shock
145 IF (s < s_l) THEN
146 !far left region
147 rho = rho_l
148 u = u_l
149 p = p_l
150 ELSE
151 !left star region
152 rho = rho_lstar
153 u = u_star
154 p = p_star
155 END IF
156 ELSE
157 !left fan
158 IF (s < s_hl) THEN
159 !far left region
160 rho = rho_l
161 u = u_l
162 p = p_l
163 ELSE IF (s > s_tl) THEN
164 !left star region
165 rho = rho_lstar
166 u = u_star
167 p = p_star
168 ELSE
169 !left fan region
170 rho = rho_l * ( 2.0/(gamma+1.0) + ((gamma-1.0)/((gamma+1.0)*c_l)) &
171 *(u_l - (x-x0)/t) ) ** (2.0/(gamma-1.0))
172 u = 2.0/(gamma+1.0) * (c_l + (gamma-1.0)*u_l/2.0 + (x-x0)/t)
173 p = p_l * ( 2.0/(gamma+1.0) + ((gamma-1.0)/((gamma+1.0)*c_l)) &
174 * (u_l - (x-x0)/t) ) ** (2.0*gamma/(gamma-1.0))
175 END IF
176 END IF
177 ELSE
178 !right side of contact
179 IF (p_star > p_r) THEN
180 !right shock
181 IF (s > s_r) THEN
182 rho = rho_r
183 u = u_r
184 p = p_r
185 ELSE
186 rho = rho_rstar
187 u = u_star
188 p = p_star
189 END IF
190 ELSE
191 !right fan
192 IF (s > s_hr) THEN
193 rho = rho_r
194 u = u_r
195 p = p_r
196 ELSE IF (s < s_tr) THEN
197 rho = rho_rstar
198 u = u_star
199 p = p_star
200 ELSE
201 rho = rho_r * ( 2.0/(gamma+1.0) - ((gamma-1.0)/((gamma+1.0)*c_r)) &
202 * (u_r - (x-x0)/t) ) ** (2.0/(gamma-1.0))
203 u = 2.0/(gamma+1.0) * ( -c_r + (gamma-1.0)*u_r/2.0 + (x-x0)/t)
204 p = p_r * ( 2.0/(gamma+1.0) - ((gamma-1.0)/((gamma+1.0)*c_r)) &
205 * (u_r - (x-x0)/t) ) ** (2.0*gamma/(gamma-1.0))
206 END IF
207 END IF
208 END IF
209 END SUBROUTINE sample
210 END SUBROUTINE riemann
211
212 ! the root of this function gives p_star
213 !REAL FUNCTION f(p)
214 PURE SUBROUTINE f(p,fx,plist)
215 IMPLICIT NONE
216 !----------------------------------------------------------------------!
217 REAL, INTENT(IN) :: p
218 REAL, INTENT(OUT) :: fx
219 REAL, INTENT(IN), DIMENSION(:), OPTIONAL :: plist
220 !----------------------------------------------------------------------!
221 fx = f_x(p, rho_l, u_l, p_l, a_l, b_l, c_l) &
222 + f_x(p, rho_r, u_r, p_r, a_r, b_r, c_r) &
223 + (u_r - u_l)
224 END SUBROUTINE f
225
226 !left side or right side function
227 REAL PURE function f_x(p, rho_x, u_x, p_x, a_x, b_x, c_x)
228 IMPLICIT NONE
229 !----------------------------------------------------------------------!
230 REAL, INTENT(IN) :: p, rho_x, u_x, p_x, a_x, b_x, c_x
231 !----------------------------------------------------------------------!
232 IF (p > p_x) THEN
233 f_x = (p-p_x) * sqrt(a_x/(p+b_x))
234 ELSE
235 f_x = (2.0*c_x/(gamma-1.0)) * ((p/p_x)**((gamma-1.0)/(2.0*gamma)) - 1.0)
236 END IF
237 END FUNCTION f_x
238
239END MODULE
240
roots
root finding subroutines
Definition: roots.f90:45
roots::getroot_newton
subroutine, public getroot_newton(funcd, x1, x2, root, error, plist, xm, iterations)
Definition: roots.f90:247
solutions::p_r
real p_r
Definition: solutions.f90:38
solutions::riemann
subroutine, public riemann(x0, gamma_, rho_l_, u_l_, p_l_, rho_r_, u_r_, p_r_, t, x, pvar)
Definition: solutions.f90:46
solutions::rho_r
real rho_r
Definition: solutions.f90:38
solutions::b_r
real b_r
Definition: solutions.f90:38
solutions::f_x
real pure function f_x(p, rho_x, u_x, p_x, A_x, B_x, c_x)
Definition: solutions.f90:228
solutions::u_r
real u_r
Definition: solutions.f90:38
solutions::f
pure subroutine f(p, fx, plist)
Definition: solutions.f90:215
solutions::a_r
real a_r
Definition: solutions.f90:38
solutions::u_l
real u_l
Definition: solutions.f90:38
solutions::rho_l
real rho_l
Definition: solutions.f90:38
solutions::b_l
real b_l
Definition: solutions.f90:38
solutions::a_l
real a_l
Definition: solutions.f90:38
solutions::c_r
real c_r
Definition: solutions.f90:38
solutions::gamma
real gamma
Definition: solutions.f90:38
solutions::c_l
real c_l
Definition: solutions.f90:38
solutions::p_l
real p_l
Definition: solutions.f90:38
sample
subroutine sample(x, rho, u, p)
Definition: solutions.f90:129
Generated by   doxygen 1.9.4