42 REAL,
PARAMETER :: tsim = 1.0
43 REAL,
PARAMETER ::
gamma = 1.4
45 REAL,
PARAMETER :: rho0 = 1.0
46 REAL,
PARAMETER :: p0 = 1.0e-05
47 REAL,
PARAMETER :: e1 = 0.311357
51 REAL,
PARAMETER :: r0 = 3.0e-2
54 INTEGER,
PARAMETER :: mgeo = cylindrical
56 INTEGER,
PARAMETER :: xres = 200
57 INTEGER,
PARAMETER :: yres = 6
58 INTEGER,
PARAMETER :: zres = 1
59 REAL,
PARAMETER :: gpar = 0.2
60 REAL,
PARAMETER :: rmax = 2.0
61 REAL,
PARAMETER :: rmin = 1.0e-03
63 INTEGER,
PARAMETER :: onum = 10
64 CHARACTER(LEN=256),
PARAMETER &
66 CHARACTER(LEN=256),
PARAMETER &
69 CLASS(
fosite),
ALLOCATABLE :: sim
70 REAL,
DIMENSION(:),
ALLOCATABLE :: sigma, sum_numer, sum_denom
71 INTEGER :: n,den,vel,pre
78 IF (sim%GetRank().EQ.0)
THEN
87 ALLOCATE(sum_numer(sim%Physics%VNUM),sum_denom(sim%Physics%VNUM),sigma(sim%Physics%VNUM))
90 CALL initdata(sim%Mesh, sim%Physics, sim%Timedisc)
91 CALL run(sim%Mesh, sim%Physics, sim%Timedisc)
95 IF (
ASSOCIATED(sim%Timedisc%solution))
THEN
96 DO n=1,sim%Physics%VNUM
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)
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)
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)
112 WHERE (abs(sum_denom(:)).GT.2*tiny(sum_denom))
113 sigma(:) = sum_numer(:) / sum_denom(:)
117 den = sim%Physics%DENSITY
118 vel = sim%Physics%XVELOCITY
119 pre = sim%Physics%PRESSURE
122 IF (sim%GetRank().EQ.0)
THEN
124tap_check(ok,
"stoptime reached")
126tap_check_small(sigma(den),4.0e-02,
"density deviation < 4%")
127tap_check_small(sigma(vel),4.0e-02,
"radial velocity deviation < 4%")
129tap_check_small(sigma(pre),4.0e-02,
"pressure deviation < 4%")
136 DEALLOCATE(sim,sigma,sum_numer,sum_denom)
145 TYPE(dict_typ),
POINTER :: config
149 TYPE(dict_typ),
POINTER :: mesh, physics, boundary, datafile, &
151 REAL :: x1,x2,y1,y2,z1,z2
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
177 bc(west) = reflecting
190 bc(west) = reflecting
191 bc(east) = no_gradients
194 bc(bottom)= no_gradients
195 bc(top) = no_gradients
203 bc(west) = no_gradients
204 bc(east) = no_gradients
207 bc(bottom)= no_gradients
208 bc(top) = no_gradients
237 CALL sim%Physics%Error(
"InitProgram",
"geometry not supported for 3D Sedov explosion")
242 "meshtype" / midpoint, &
257 "western" / bc(west), &
258 "eastern" / bc(east), &
259 "southern" / bc(south), &
260 "northern" / bc(north), &
261 "bottomer" / bc(bottom), &
273 "variables" / primitive, &
278 "method" / dormand_prince, &
283 "dtlimit" / 1.0e-13, &
285 "tol_abs" / (/0.0,1e-3,1e-3,0.0/), &
286 "maxiter" / 1000000, &
287 "output/solution" / 1, &
289 "output/geometrical_sources" / 1 )
293 "fileformat" / vtk, &
294 "filename" / (trim(odir) // trim(ofname)), &
299 "physics" / physics, &
300 "boundary" / boundary, &
302 "timedisc" / timedisc, &
303 "datafile" / datafile )
310 CLASS(physics_base),
INTENT(IN) :: Physics
311 CLASS(mesh_base),
INTENT(IN) :: Mesh
312 CLASS(timedisc_base),
INTENT(INOUT) :: Timedisc
319 SELECT TYPE (phys => physics)
321 SELECT TYPE (pvar => timedisc%pvar)
324 pvar%density%data1d(:) = rho0
326 pvar%velocity%data1d(:) = 0.0
329 p1 = 3.*(phys%gamma - 1.0)*e1 / ((n + 1)*pi*r0**n)
330 WHERE (mesh%radius%bcenter(:,:,:).LE.r0)
332 pvar%pressure%data3d(:,:,:) = p1
335 pvar%pressure%data3d(:,:,:) = p0
339 CALL phys%Error(
"sedov2d:InitData",
"statevector must be of class euler")
343 CALL phys%Error(
"sedov2d:InitData",
"physics not supported")
346 CALL physics%Convert2Conservative(timedisc%pvar,timedisc%cvar)
347 CALL mesh%Info(
" DATA-----> initial condition: 2D Sedov explosion")
351 SUBROUTINE run(Mesh,Physics,Timedisc)
355 CLASS(mesh_base),
INTENT(IN) :: Mesh
356 CLASS(physics_base),
INTENT(IN) :: Physics
357 CLASS(timedisc_base),
INTENT(INOUT) :: Timedisc
361 REAL,
DIMENSION(:,:,:,:),
ALLOCATABLE :: er
364 IF (
ASSOCIATED(timedisc%solution))
THEN
369 IF (sim%GetRank().EQ.0) &
371 CALL sedov%PrintConfiguration()
373 timedisc%solution = timedisc%pvar
376 t0 = sqrt((e1/rho0)*(r0/sedov%alpha)**5)
379 ALLOCATE(er(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,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)
385 DO WHILE((timedisc%maxiter.LE.0).OR.(sim%iter.LE.timedisc%maxiter))
386 SELECT TYPE(solution => timedisc%solution)
387 TYPE IS (statevector_euler)
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))
397 IF (btest(mesh%VECTOR_COMPONENTS,n-1))
THEN
398 solution%velocity%data4d(i,j,k,m) = vr * er(i,j,k,n)
subroutine initdata(Mesh, Physics, Fluxes, Timedisc)
subroutine makeconfig(Sim, config)
subroutine, private run(this)
elemental real function, public acosh(x)
inverse hyperbolic cosine function
elemental real function, public asinh(x)
inverse hyperbolic sine function
physics module for 1D,2D and 3D non-isothermal Euler equations
module providing sedov class with basic init and solve subroutines
program sedov2d
2D Sedov explosion