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