riemann1d.f90
Go to the documentation of this file.
1!#############################################################################
2!# #
3!# fosite - 3D hydrodynamical simulation program #
4!# module: riemann1d.f90 #
5!# #
6!# Copyright (C) 2006-2024 #
7!# Tobias Illenseer <tillense@astrophysik.uni-kiel.de> #
8!# #
9!# This program is free software; you can redistribute it and/or modify #
10!# it under the terms of the GNU General Public License as published by #
11!# the Free Software Foundation; either version 2 of the License, or (at #
12!# your option) any later version. #
13!# #
14!# This program is distributed in the hope that it will be useful, but #
15!# WITHOUT ANY WARRANTY; without even the implied warranty of #
16!# MERCHANTABILITY OR FITNESS FOR A PARTICULAR PURPOSE, GOOD TITLE or #
17!# NON INFRINGEMENT. See the GNU General Public License for more #
18!# details. #
19!# #
20!# You should have received a copy of the GNU General Public License #
21!# along with this program; if not, write to the Free Software #
22!# Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA. #
23!# #
24!#############################################################################
25
26!----------------------------------------------------------------------------!
53!----------------------------------------------------------------------------!
54PROGRAM riemann1d
55 USE fosite_mod
56 USE solutions
57#include "tap.h"
58 IMPLICIT NONE
59 !--------------------------------------------------------------------------!
60 INTEGER, PARAMETER :: testnum = 8 ! number of available tests (do not modify)
61 CHARACTER(LEN=32), PARAMETER :: teststr(testnum) = (/ &
62 "Sod shock tube ", &
63 "Toro test no. 2 ", &
64 "Toro test no. 3 ", &
65 "Toro test no. 4 ", &
66 "Toro test no. 5 ", &
67 "Noh problem ", &
68 "Isothermal shock tube ", &
69 "Isothermal Noh problem " /)
70 !--------------------------------------------------------------------------!
71 ! mesh settings
72 CHARACTER(LEN=3), PARAMETER :: testdir = "all" !"x" ! direction: x,y,z or all
73 INTEGER, PARAMETER :: res = 400 ! resolution
74 ! output parameters
75 INTEGER, PARAMETER :: onum = 1 ! number of output data sets
76 CHARACTER(LEN=8), PARAMETER :: ofname(testnum) = (/ & ! output file name
77 "sod ","toro2 ","toro3 ","toro4 ","toro5 ","noh ", &
78 "sodiso ","nohiso " /)
79 CHARACTER(LEN=256), PARAMETER & ! output data dir
80 :: odir = './'
81 ! some global constants
82 REAL :: gamma ! ratio of specific heats
83 REAL :: tsim ! simulation time
84 REAL :: csiso = 0.0 ! isothermal sound speed (test no. 6)
85 !--------------------------------------------------------------------------!
86 CLASS(fosite), ALLOCATABLE :: sim
87 CHARACTER(LEN=128) :: verbose_tap_output
88 INTEGER :: ic,sd,dir_min,dir_max,n
89 REAL, DIMENSION(:,:), ALLOCATABLE :: sigma
90 REAL, PARAMETER :: sigma_tol(testnum) &
91 = (/ 3.0e-3, 1.5e-2, 5.0e-3, 7.0e-3, 4.0e-3, 5.0e-3, 1.0e+0, 1.0e+0 /)
92 REAL :: sum_numer, sum_denom
93 !--------------------------------------------------------------------------!
94
95 ! check whether we perform tests in all directions
96 SELECT CASE (trim(testdir))
97 CASE("all")
98 dir_min = 1
99 dir_max = 3
100 CASE("y")
101 dir_min = 2
102 dir_max = 2
103 CASE("z")
104 dir_min = 3
105 dir_max = 3
106 CASE DEFAULT
107 ! x-direction is also the default
108 dir_min = 1
109 dir_max = 1
110 END SELECT
111
112 ALLOCATE(sim,sigma(testnum,dir_min:dir_max))
113
114 ! initialize Fosite
115 CALL sim%InitFosite()
116
117#ifdef PARALLEL
118 IF (sim%GetRank().EQ.0) THEN
119#endif
120tap_plan(testnum*(dir_max-dir_min+1))
121#ifdef PARALLEL
122 END IF
123#endif
124
125 ! loop over all tests
126 tests: DO ic=1,testnum
127 ! loop over selected directions
128 DO sd=dir_min,dir_max
129 ! get configuration for the particular test and direction
130 CALL makeconfig(sim, sim%config,ic,sd)
131! CALL PrintDict(config)
132
133 ! setup the simulation
134 CALL sim%Setup()
135
136 ! set initial conditions and run the simulation
137 CALL run(sim,ic,sd)
138
139 IF (sim%aborted) THEN
140 sigma(ic,sd) = huge(1.0)
141 ELSE
142 sigma(ic,sd) = 0.0 ! default
143 sum_numer = 0.0
144 sum_denom = 0.0
145 IF (ASSOCIATED(sim%Timedisc%solution)) THEN
146 ! use L1 norm to estimate the deviation from the exact solution:
147 ! Σ |pvar - pvar_exact| / Σ |pvar_exact|
148 DO n=1,sim%Physics%VNUM
149 sum_numer = sum(abs(sim%Timedisc%pvar%data2d(:,n)-sim%Timedisc%solution%data2d(:,n)), &
150 mask=sim%Mesh%without_ghost_zones%mask1d(:))
151 sum_denom = sum(abs(sim%Timedisc%solution%data2d(:,n)),mask=sim%Mesh%without_ghost_zones%mask1d(:))
152 END DO
153#ifdef PARALLEL
154 CALL mpi_allreduce(mpi_in_place,sum_numer,1,default_mpi_real,mpi_sum,mpi_comm_world,sim%ierror)
155 CALL mpi_allreduce(mpi_in_place,sum_denom,1,default_mpi_real,mpi_sum,mpi_comm_world,sim%ierror)
156#endif
157 IF (abs(sum_denom).GT.4*tiny(sum_denom)) & ! avoid division by zero error
158 sigma(ic,sd) = sum_numer / sum_denom
159 END IF
160 END IF
161
162 ! reset fosite
163 IF (.NOT.(ic.EQ.testnum.AND.sd.EQ.dir_max)) THEN
164 CALL sim%Finalize(mpifinalize_=.false.)
165 DEALLOCATE(sim)
166 ALLOCATE(sim)
167 CALL sim%InitFosite()
168 END IF
169 END DO
170 END DO tests
171
172 ! check results
173 ! loop over all tests
174#ifdef PARALLEL
175 IF (sim%GetRank().EQ.0) THEN
176#endif
177 DO ic=1,testnum
178 ! loop over selected directions
179 DO sd=dir_min,dir_max
180 WRITE (verbose_tap_output,'(A,ES10.2)') teststr(ic) // " " // achar(119+sd) &
181 // "-direction: " // " sigma = ", sigma(ic,sd)
182! This line is long if expanded. So we can't indent it or it will be cropped.
183tap_check_small(sigma(ic,sd),sigma_tol(ic),trim(verbose_tap_output))
184 END DO
185 END DO
186
187tap_done
188#ifdef PARALLEL
189 END IF
190#endif
191
192CALL sim%Finalize()
193DEALLOCATE(sim,sigma)
194
195CONTAINS
196
197 SUBROUTINE makeconfig(Sim,config,ic,dir)
198 IMPLICIT NONE
199 !------------------------------------------------------------------------!
200 CLASS(fosite) :: Sim
201 TYPE(dict_typ),POINTER :: config
202 INTEGER :: ic,dir
203 !------------------------------------------------------------------------!
204 ! Local variable declaration
205 TYPE(dict_typ), POINTER :: mesh, physics, boundary, datafile, &
206 timedisc, fluxes
207 INTEGER :: xres,yres,zres
208 REAL :: xmax,ymax,zmax
209 !------------------------------------------------------------------------!
210 INTENT(INOUT) :: sim
211 INTENT(IN) :: ic,dir
212 !------------------------------------------------------------------------!
213 ! mesh settings depends on the direction selected above
214 SELECT CASE(dir)
215 CASE(1)
216 xmax = 1.0
217 ymax = 0.0
218 zmax = 0.0
219 xres = res
220 yres = 1
221 zres = 1
222 CASE(2)
223 xmax = 0.0
224 ymax = 1.0
225 zmax = 0.0
226 xres = 1
227 yres = res
228 zres = 1
229 CASE(3)
230 xmax = 0.0
231 ymax = 0.0
232 zmax = 1.0
233 xres = 1
234 yres = 1
235 zres = res
236 CASE DEFAULT
237 CALL sim%Error("riemann1d::InitData","direction must be one of 1,2,3")
238 END SELECT
239
240 ! mesh settings
241 mesh => dict("meshtype" / midpoint, &
242 "geometry" / cartesian, &
243 "inum" / xres, & ! resolution in x and !
244 "jnum" / yres, & ! y direction !
245 "knum" / zres, & ! z direction !
246 "xmin" / (0.0), &
247 "xmax" / xmax, &
248 "ymin" / (-0.0), &
249 "ymax" / ymax, &
250 "zmin" / (0.0), &
251 "zmax" / zmax &
252 )
253
254 SELECT CASE(ic)
255 CASE(1)
256 gamma = 1.4
257 tsim = 0.25
258 CASE(2)
259 gamma = 1.4
260 tsim = 0.15
261 CASE(3)
262 gamma = 1.4
263 tsim = 0.012
264 CASE(4)
265 gamma = 1.4
266 tsim = 0.035
267 CASE(5)
268 gamma = 1.4
269 tsim = 0.012
270 CASE(6)
271 gamma = 5./3.
272 tsim = 1.0
273 CASE(7)
274 csiso = 1.0
275 tsim = 0.25
276 CASE(8)
277 csiso = 1.0
278 tsim = 0.5
279 CASE DEFAULT
280 CALL sim%Mesh%Error("riemann1d::InitData", "Test problem number should be 1,2,3,4,5,6,7 or 8")
281 END SELECT
282
283
284 ! physics settings
285 IF (csiso.NE.0.0) THEN
286 physics => dict("problem" / euler_isotherm, &
287 "cs" / csiso)
288 ELSE
289 physics => dict("problem" / euler, &
290 "gamma" / gamma) ! ratio of specific heats !
291 END IF
292
293 ! flux calculation and reconstruction method
294 fluxes => dict(&
295 "fluxtype" / kt, &
296 "order" / linear, &
297 "variables" / conservative, & ! vars. to use for reconstruction!
298 "limiter" / monocent, & ! one of: minmod, monocent,... !
299 "theta" / 1.2 & ! optional parameter for limiter !
300 )
301 ! time discretization settings
302 timedisc => dict( &
303 "method" / modified_euler, &
304 "order" / 3, &
305 "cfl" / 0.4, &
306 "stoptime" / tsim, &
307 "dtlimit" / 1.0e-10, &
308! "output/xmomentum" / 1, &
309! "output/ymomentum" / 1, &
310! "output/zmomentum" / 1, &
311! "output/energy" / 1, &
312! "output/rhs" / 1 &
313 "output/solution" / 0, &
314 "maxiter" / 100000 &
315 )
316
317 ! compute and output exact solution in non-isothermal simulations
318 IF (csiso.EQ.0.0) CALL setattr(timedisc,"output/solution",1)
319
320 ! boundary conditions
321 boundary => dict( &
322 "western" / no_gradients, &
323 "eastern" / no_gradients, &
324 "southern" / no_gradients, &
325 "northern" / no_gradients, &
326 "bottomer" / no_gradients, &
327 "topper" / no_gradients)
328
329 ! initialize data input/output
330 datafile => dict( &
331 "fileformat" / gnuplot, &
332 "filename" / (trim(odir) // trim(ofname(ic)) // "-" // achar(119+dir)), &
333 "filecycles" / 0, &
334 "count" / onum)
335
336 config => dict("mesh" / mesh, &
337 "physics" / physics, &
338 "boundary" / boundary, &
339 "fluxes" / fluxes, &
340 "timedisc" / timedisc, &
341 "datafile" / datafile)
342 END SUBROUTINE makeconfig
343
344 SUBROUTINE run(this,ic,dir)
345 USE solutions
346 IMPLICIT NONE
347 !------------------------------------------------------------------------!
348 CLASS(fosite), INTENT(INOUT) :: this
349! CLASS(marray_compound), POINTER :: pvar
350 INTEGER, INTENT(IN) :: ic,dir
351 !------------------------------------------------------------------------!
352 ! Local variable declaration
353 REAL :: x0, rho_l, rho_r, u_l, u_r, p_l, p_r
354 REAL, DIMENSION(:), ALLOCATABLE :: x
355 !------------------------------------------------------------------------!
356 SELECT CASE(ic)
357 CASE(1) ! Sod shock tube (Toro's test no. 1)
358 x0 = 0.5
359 rho_l = 1.0
360 rho_r = 0.125
361 u_l = 0.0
362 u_r = 0.0
363 p_l = 1.0
364 p_r = 0.1
365 CASE(2) ! Toro's test no. 2
366 x0 = 0.5
367 rho_l = 1.0
368 rho_r = 1.0
369 u_l = -2.0
370 u_r = 2.0
371 p_l = 0.4
372 p_r = 0.4
373 CASE(3) ! Toro's test no. 3
374 x0 = 0.5
375 rho_l = 1.0
376 rho_r = 1.0
377 u_l = 0.0
378 u_r = 0.0
379 p_l = 1000.
380 p_r = 0.01
381 CASE(4) ! Toro's test no. 4
382 x0 = 0.4
383 rho_l = 5.99924
384 rho_r = 5.99924
385 u_l = 19.5975
386 u_r = -6.19633
387 p_l = 460.894
388 p_r = 46.0950
389 CASE(5) ! Toro's test no. 5
390 x0 = 0.8
391 rho_l = 1.0
392 rho_r = 1.0
393 u_l = -19.5975
394 u_r = -19.5975
395 p_l = 1000.0
396 p_r = 0.01
397 CASE(6) ! Noh problem
398 x0 = 0.5
399 rho_l = 1.0
400 rho_r = 1.0
401 u_l = 1.0
402 u_r = -1.0
403 p_l = 1.0e-05
404 p_r = 1.0e-05
405 CASE(7) ! isothermal shock tube
406 x0 = 0.5
407 rho_l = 1.0
408 rho_r = 0.125
409 u_l = 0.0
410 u_r = 0.0
411 CASE(8) ! isothermal Noh problem
412 x0 = 0.5
413 rho_l = 1.0
414 rho_r = 1.0
415 u_l = 1.0
416 u_r = -1.0
417 CASE DEFAULT
418 CALL this%Error("InitData", "Test problem should be 1,2,3,4,5,6,7 or 8")
419 END SELECT
420
421 ! initial condition
422 ! always set density and first velocity
423 SELECT TYPE(pvar => this%Timedisc%pvar)
424 CLASS IS(statevector_eulerisotherm)
425 ! set all velocities to 0
426 pvar%velocity%data1d(:) = 0.0
427 WHERE (this%Mesh%bcenter(:,:,:,dir).LT.x0)
428 pvar%density%data3d(:,:,:) = rho_l
429 pvar%velocity%data4d(:,:,:,1) = u_l
430 ELSEWHERE
431 pvar%density%data3d(:,:,:) = rho_r
432 pvar%velocity%data4d(:,:,:,1) = u_r
433 END WHERE
434 END SELECT
435
436 ! set pressure for non-isothermal HD
437 SELECT TYPE(pvar => this%Timedisc%pvar)
438 TYPE IS(statevector_euler)
439 WHERE (this%Mesh%bcenter(:,:,:,dir).LT.x0)
440 pvar%pressure%data3d(:,:,:) = p_l
441 ELSEWHERE
442 pvar%pressure%data3d(:,:,:) = p_r
443 END WHERE
444 END SELECT
445
446 CALL this%Physics%Convert2Conservative(this%Timedisc%pvar,this%Timedisc%cvar)
447 CALL this%Info(" DATA-----> initial condition: " &
448 // "1D Riemann problem: " // trim(teststr(ic)))
449
450 IF (ASSOCIATED(this%Timedisc%solution)) THEN
451 ! store initial condition in exact solution array
452 this%Timedisc%solution = this%Timedisc%pvar
453
454 ! set 1D coordinate array for exact solution
455 ALLOCATE(x(SIZE(this%Mesh%cart%data3d(:,2,dir),1)), &
456 source=this%Mesh%cart%data3d(:,2,dir))
457
458 DO WHILE((this%Timedisc%maxiter.LE.0).OR.(this%iter.LE.this%Timedisc%maxiter))
459 ! compute exact solution of the 1D Riemann problem (only non-isothermal HD)
460 ! if output is written during next integration step
461 IF (this%Datafile%time-this%Timedisc%time.LT.this%Timedisc%dt) THEN
462 SELECT TYPE(phys => this%Physics)
463 TYPE IS (physics_euler)
464 CALL riemann(x0,phys%gamma,rho_l,u_l,p_l,rho_r,u_r,p_r, &
465 this%Datafile%time,x,this%Timedisc%solution%data2d(:,:))
466 END SELECT
467 END IF
468 IF(this%Step()) EXIT
469 END DO
470 ELSE
471 CALL this%Run()
472 END IF
473 END SUBROUTINE run
474
475END PROGRAM riemann1d
subroutine makeconfig(Sim, config)
Definition: bondi2d.f90:165
subroutine, private run(this)
Definition: fosite.f90:352
subroutine, public riemann(x0, gamma_, rho_l_, u_l_, p_l_, rho_r_, u_r_, p_r_, t, x, pvar)
Definition: solutions.f90:46
real gamma
Definition: solutions.f90:38
program riemann1d
1D Riemann problems
Definition: riemann1d.f90:54
main fosite class
Definition: fosite.f90:71