sedov3d.f90
Go to the documentation of this file.
1!#############################################################################
2!# #
3!# fosite - 3D hydrodynamical simulation program #
4!# module: sedov3d.f90 #
5!# #
6!# Copyright (C) 2006-2014 #
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!----------------------------------------------------------------------------!
34!----------------------------------------------------------------------------!
35PROGRAM sedov3d
36 USE fosite_mod
37 USE solutions
38#include "tap.h"
39 IMPLICIT NONE
40 !--------------------------------------------------------------------------!
41 ! simulation parameters
42 REAL, PARAMETER :: tsim = 1.0 ! simulation stop time
43 REAL, PARAMETER :: gamma = 1.4 ! ratio of specific heats
44 ! initial condition (dimensionless units)
45 REAL, PARAMETER :: rho0 = 1.0 ! ambient density
46 REAL, PARAMETER :: p0 = 1.0e-05 ! ambient pressure
47 REAL, PARAMETER :: e1 = 1.0 ! initial energy input
48 ! Spatial with of the initial pulse should be at least 5 cells;
49 ! if you wish to compare the results on different grids
50 ! R0 should be of the same order
51 REAL, PARAMETER :: r0 = 5.0e-2
52 ! mesh settings
53! INTEGER, PARAMETER :: MGEO = LOGSPHERICAL ! geometry
54 INTEGER, PARAMETER :: mgeo = spherical ! geometry
55! INTEGER, PARAMETER :: MGEO = CARTESIAN ! geometry
56 INTEGER, PARAMETER :: xres = 200 ! x-resolution
57 INTEGER, PARAMETER :: yres = 6 ! y-resolution
58 INTEGER, PARAMETER :: zres = 1 ! z-resolution
59 REAL, PARAMETER :: rmax = 2.0 ! outer radius of comput. domain
60 REAL, PARAMETER :: gpar = 0.2 ! geometry scaling parameter
61 ! output parameters
62 INTEGER, PARAMETER :: onum = 10 ! number of output data sets
63 CHARACTER(LEN=256), PARAMETER & ! output data dir
64 :: odir = './'
65 CHARACTER(LEN=256), PARAMETER & ! output data file name
66 :: ofname = 'sedov3d'
67 !-------------------------------------------------------------------------!
68 CLASS(fosite), ALLOCATABLE :: sim
69 REAL, DIMENSION(:), ALLOCATABLE :: sigma, sum_numer, sum_denom
70 INTEGER :: n,den,vel,pre
71 LOGICAL :: ok
72 !-------------------------------------------------------------------------!
73 ALLOCATE(sim)
74 CALL sim%InitFosite()
75
76#ifdef PARALLEL
77 IF (sim%GetRank().EQ.0) THEN
78#endif
79tap_plan(4)
80#ifdef PARALLEL
81 END IF
82#endif
83
84 CALL makeconfig(sim, sim%config)
85 CALL sim%Setup()
86 ALLOCATE(sum_numer(sim%Physics%VNUM),sum_denom(sim%Physics%VNUM),sigma(sim%Physics%VNUM))
87 sigma(:) = 0.0
88
89 CALL initdata(sim%Mesh, sim%Physics, sim%Timedisc)
90 CALL run(sim%Mesh, sim%Physics, sim%Timedisc)
91 ok = .NOT.sim%aborted
92
93 ! compare with exact solution if requested
94 IF (ASSOCIATED(sim%Timedisc%solution)) THEN
95 DO n=1,sim%Physics%VNUM
96 ! use L1 norm to estimate the deviation from the exact solution:
97 ! Σ |pvar - pvar_exact| / Σ |pvar_exact|
98 sum_numer(n) = sum(abs(sim%Timedisc%pvar%data2d(:,n)-sim%Timedisc%solution%data2d(:,n)), &
99 mask=sim%Mesh%without_ghost_zones%mask1d)
100 sum_denom(n) = sum(abs(sim%Timedisc%solution%data2d(:,n)),mask=sim%Mesh%without_ghost_zones%mask1d)
101 END DO
102#ifdef PARALLEL
103 IF (sim%GetRank().GT.0) THEN
104 CALL mpi_reduce(sum_numer,sum_numer,sim%Physics%VNUM,default_mpi_real,mpi_sum,0,mpi_comm_world,sim%ierror)
105 CALL mpi_reduce(sum_denom,sum_denom,sim%Physics%VNUM,default_mpi_real,mpi_sum,0,mpi_comm_world,sim%ierror)
106 ELSE
107 CALL mpi_reduce(mpi_in_place,sum_numer,sim%Physics%VNUM,default_mpi_real,mpi_sum,0,mpi_comm_world,sim%ierror)
108 CALL mpi_reduce(mpi_in_place,sum_denom,sim%Physics%VNUM,default_mpi_real,mpi_sum,0,mpi_comm_world,sim%ierror)
109 END IF
110#endif
111 WHERE (abs(sum_denom(:)).GT.2*tiny(sum_denom))
112 sigma(:) = sum_numer(:) / sum_denom(:)
113 END WHERE
114 END IF
115
116 den = sim%Physics%DENSITY
117 vel = sim%Physics%XVELOCITY
118 pre = sim%Physics%PRESSURE
119
120#ifdef PARALLEL
121 IF (sim%GetRank().EQ.0) THEN
122#endif
123tap_check(ok,"stoptime reached")
124! These lines are very long if expanded. So we can't indent it or it will be cropped.
125tap_check_small(sigma(den),5.0e-02,"density deviation < 5%")
126tap_check_small(sigma(vel),4.0e-02,"radial velocity deviation < 4%")
127tap_check_small(sigma(pre),5.0e-02,"pressure deviation < 5%")
128tap_done
129! PRINT *,sigma(:)
130#ifdef PARALLEL
131 END IF
132#endif
133
134 CALL sim%Finalize()
135 DEALLOCATE(sim,sigma,sum_numer,sum_denom)
136
137CONTAINS
138
139 SUBROUTINE makeconfig(Sim, config)
140 USE functions, ONLY : asinh,acosh
141 IMPLICIT NONE
142 !------------------------------------------------------------------------!
143 CLASS(fosite) :: Sim
144 TYPE(dict_typ),POINTER :: config
145 !------------------------------------------------------------------------!
146 ! Local variable declaration
147 INTEGER :: bc(6)
148 TYPE(dict_typ), POINTER :: mesh, physics, boundary, datafile, &
149 timedisc, fluxes
150 REAL :: x1,x2,y1,y2,z1,z2
151 !------------------------------------------------------------------------!
152 INTENT(INOUT) :: sim
153 !------------------------------------------------------------------------!
154 ! mesh settings and boundary conditions
155 SELECT CASE(mgeo)
156 CASE(cartesian)
157 x1 = -0.5
158 x2 = 0.5
159 y1 = -0.5
160 y2 = 0.5
161 z1 = -0.5
162 z2 = 0.5
163 bc(west) = no_gradients
164 bc(east) = no_gradients
165 bc(south) = no_gradients
166 bc(north) = no_gradients
167 bc(bottom)= no_gradients
168 bc(top) = no_gradients
169 CASE(spherical)
170 x1 = 0.0
171 x2 = rmax
172 y1 = 0.0
173 y2 = pi
174 z1 = 0.0
175 z2 = 0.0 !2*PI
176 bc(west) = reflecting
177 bc(east) = no_gradients
178 bc(south) = axis
179 bc(north) = axis
180 bc(bottom)= periodic
181 bc(top) = periodic
182 CASE(logspherical)
183 x1 = log(1e-3/gpar)
184 x2 = log(rmax/gpar)
185 y1 = 0.0
186 y2 = pi
187 z1 = 0.0
188 z2 = 2*pi
189 bc(west) = reflecting
190 bc(east) = no_gradients
191 bc(south) = axis
192 bc(north) = axis
193 bc(bottom)= periodic
194 bc(top) = periodic
195 CASE(cylindrical)
196 x1 = -rmax
197 x2 = rmax
198 y1 = 0.0
199 y2 = 2.0*pi
200 z1 = 0.0
201 z2 = rmax
202 bc(west) = no_gradients !default: ABSORBING
203 bc(east) = no_gradients !default: ABSORBING
204 bc(south) = no_gradients !default: AXIS
205 bc(north) = no_gradients !default: ABSORBING
206 bc(bottom)= no_gradients
207 bc(top) = no_gradients
208! CASE(OBLATE_SPHEROIDAL) !TODO anderen Meshs zu 3D konvertieren
209! x1 = 0.0
210! x2 = Acosh(RMAX/GPAR)
211! y1 = -0.5*PI
212! y2 = 0.5*PI
213! bc(WEST) = FOLDED
214! bc(EAST) = NO_GRADIENTS
215! bc(SOUTH) = AXIS
216! bc(NORTH) = AXIS
217! CASE(TANCYLINDRICAL)
218! x1 = ATAN(-RMAX/GPAR)
219! x2 = ATAN(RMAX/GPAR)
220! y1 = 0.0
221! y2 = RMAX
222! bc(WEST) = NO_GRADIENTS
223! bc(EAST) = NO_GRADIENTS
224! bc(SOUTH) = AXIS
225! bc(NORTH) = NO_GRADIENTS
226! CASE(SINHSPHERICAL)
227! x1 = 0.0
228! x2 = Asinh(RMAX/GPAR)
229! y1 = 0.0
230! y2 = PI
231! bc(WEST) = REFLECTING
232! bc(EAST) = NO_GRADIENTS
233! bc(SOUTH) = AXIS
234! bc(NORTH) = AXIS
235 CASE DEFAULT
236 CALL sim%Physics%Error("InitProgram","geometry not supported for 3D Sedov explosion")
237 END SELECT
238
239 ! mesh settings
240 mesh => dict( &
241 "meshtype" / midpoint, &
242 "geometry" / mgeo, &
243 "inum" / xres, &
244 "jnum" / yres, &
245 "knum" / zres, &
246 "xmin" / x1, &
247 "xmax" / x2, &
248 "ymin" / y1, &
249 "ymax" / y2, &
250 "zmin" / z1, &
251 "zmax" / z2, &
252 "gparam" / gpar )
253
254 ! boundary conditions
255 boundary => dict( &
256 "western" / bc(west), &
257 "eastern" / bc(east), &
258 "southern" / bc(south), &
259 "northern" / bc(north), &
260 "bottomer" / bc(bottom), &
261 "topper" / bc(top))
262
263 ! physics settings
264 physics => dict( &
265 "problem" / euler, &
266 "gamma" / gamma)
267
268 ! flux calculation and reconstruction method
269 fluxes => dict( &
270 "order" / linear, &
271 "fluxtype" / kt, &
272 "variables" / primitive, &
273 "limiter" / vanleer)
274
275 ! time discretization settings
276 timedisc => dict( &
277 "method" / dormand_prince, &
278 "order" / 5, &
279 "cfl" / 0.4, &
280 "stoptime" / tsim, &
281 "dtlimit" / 1.0e-13, &
282 "tol_rel" / 0.01, &
283! "tol_abs" / (/1e-5,1e-5,1e-5,1e-5,1e-5/), &
284 "output/solution" / 1, &
285 "maxiter" / 1000000 )
286
287 ! initialize data input/output
288 datafile => dict( &
289 "fileformat" / vtk, &
290 "filename" / (trim(odir) // trim(ofname)), &
291 "count" / onum)
292
293 config => dict( &
294 "mesh" / mesh, &
295 "physics" / physics, &
296 "boundary" / boundary, &
297 "fluxes" / fluxes, &
298 "timedisc" / timedisc, &
299 "datafile" / datafile )
300 END SUBROUTINE makeconfig
301
302 SUBROUTINE initdata(Mesh,Physics,Timedisc)
304 IMPLICIT NONE
305 !------------------------------------------------------------------------!
306 CLASS(physics_base), INTENT(IN) :: Physics
307 CLASS(mesh_base), INTENT(IN) :: Mesh
308 CLASS(timedisc_base), INTENT(INOUT) :: Timedisc
309 !------------------------------------------------------------------------!
310 ! Local variable declaration
311 INTEGER :: n
312 REAL :: P1
313 !------------------------------------------------------------------------!
314 ! isothermal modules are excluded
315 SELECT TYPE (phys => physics)
316 CLASS IS(physics_euler)
317 SELECT TYPE (pvar => timedisc%pvar)
318 CLASS IS(statevector_euler)
319 ! uniform density everywhere
320 pvar%density%data1d(:) = rho0
321 ! vanishing initial velocities
322 pvar%velocity%data1d(:) = 0.0
323 ! set initial peak pressure P1 inside sphere with radius R0 centered on the origin
324 n = 3 ! 3 for 3D
325 p1 = 3.*(phys%gamma - 1.0)*e1 / ((n + 1)*pi*r0**n)
326 WHERE (mesh%radius%bcenter(:,:,:).LE.r0)
327 ! behind the shock front
328 pvar%pressure%data3d(:,:,:) = p1
329 ELSEWHERE
330 ! in front of the shock front (ambient medium)
331 pvar%pressure%data3d(:,:,:) = p0
332 END WHERE
333 CLASS DEFAULT
334 ! abort
335 CALL phys%Error("sedov3d:InitData","statevector must be of class euler")
336 END SELECT
337 CLASS DEFAULT
338 ! abort
339 CALL phys%Error("sedov3d:InitData","physics not supported")
340 END SELECT
341
342 CALL physics%Convert2Conservative(timedisc%pvar,timedisc%cvar)
343 CALL mesh%Info(" DATA-----> initial condition: 3D Sedov explosion")
344
345 END SUBROUTINE initdata
346
347 SUBROUTINE run(Mesh,Physics,Timedisc)
348 USE sedov_mod
349 IMPLICIT NONE
350 !------------------------------------------------------------------------!
351 CLASS(mesh_base), INTENT(IN) :: Mesh
352 CLASS(physics_base), INTENT(IN) :: Physics
353 CLASS(timedisc_base), INTENT(INOUT) :: Timedisc
354 !------------------------------------------------------------------------!
355 TYPE(sedov_typ), ALLOCATABLE :: sedov
356 REAL :: T0,vr
357 TYPE(marray_base) :: er
358 INTEGER :: n,m,i,j,k
359 !------------------------------------------------------------------------!
360 IF (ASSOCIATED(timedisc%solution)) THEN
361 ! initialize sedov class
362 ! parameters: ndim, gamma, w, rho0, E1, P0
363 sedov = sedov_typ(3,gamma,0.0,rho0,e1,p0)
364#ifdef PARALLEL
365 IF (sim%GetRank().EQ.0) &
366#endif
367 CALL sedov%PrintConfiguration()
368 ! store initial condition in exact solution array
369 timedisc%solution = timedisc%pvar
370
371 ! compute initial time using R0 as the initial position of the shock
372 t0 = sqrt((e1/rho0)*(r0/sedov%alpha)**5)
373
374 ! set radial unit vector to compute full 3D radial velocities in any geometry
375 er = marray_base(3)
376 DO n=1,3
377 er%data4d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,n) = &
378 mesh%posvec%bcenter(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,n) &
379 / mesh%radius%bcenter(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX)
380 END DO
381 DO WHILE((timedisc%maxiter.LE.0).OR.(sim%iter.LE.timedisc%maxiter))
382 SELECT TYPE(solution => timedisc%solution)
383 TYPE IS (statevector_euler)
384 ! compute exact solution of the spherical Sedov explosion problem
385 DO k=mesh%KMIN,mesh%KMAX
386 DO j=mesh%JMIN,mesh%JMAX
387 DO i=mesh%IMIN,mesh%IMAX
388 CALL sedov%ComputeSolution(timedisc%time+t0,mesh%radius%bcenter(i,j,k), &
389 solution%density%data3d(i,j,k),vr,solution%pressure%data3d(i,j,k))
390 ! set the velocities
391 m = 1
392 DO n=1,3
393 IF (btest(mesh%VECTOR_COMPONENTS,n-1)) THEN
394 solution%velocity%data4d(i,j,k,m) = vr * er%data4d(i,j,k,n)
395 m = m + 1
396 END IF
397 END DO
398 END DO
399 END DO
400 END DO
401 END SELECT
402 IF(sim%Step()) EXIT
403 END DO
404 CALL sedov%Destroy()
405 ELSE
406 CALL sim%Run()
407 END IF
408 END SUBROUTINE run
409
410END PROGRAM sedov3d
subroutine initdata(Mesh, Physics, Fluxes, Timedisc)
Definition: bondi2d.f90:274
subroutine makeconfig(Sim, config)
Definition: bondi2d.f90:165
subroutine, private run(this)
Definition: fosite.f90:352
Mathematical functions.
Definition: functions.f90:48
elemental real function, public acosh(x)
inverse hyperbolic cosine function
Definition: functions.f90:110
real, parameter pi
Definition: functions.f90:52
elemental real function, public asinh(x)
inverse hyperbolic sine function
Definition: functions.f90:100
physics module for 1D,2D and 3D non-isothermal Euler equations
module providing sedov class with basic init and solve subroutines
Definition: sedov.f90:32
real gamma
Definition: solutions.f90:38
program sedov3d
3D Sedov explosion
Definition: sedov3d.f90:35
main fosite class
Definition: fosite.f90:71