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 = 1.0
51 REAL,
PARAMETER :: r0 = 5.0e-2
54 INTEGER,
PARAMETER :: mgeo = spherical
56 INTEGER,
PARAMETER :: xres = 200
57 INTEGER,
PARAMETER :: yres = 6
58 INTEGER,
PARAMETER :: zres = 1
59 REAL,
PARAMETER :: rmax = 2.0
60 REAL,
PARAMETER :: gpar = 0.2
62 INTEGER,
PARAMETER :: onum = 10
63 CHARACTER(LEN=256),
PARAMETER &
65 CHARACTER(LEN=256),
PARAMETER &
68 CLASS(
fosite),
ALLOCATABLE :: sim
69 REAL,
DIMENSION(:),
ALLOCATABLE :: sigma, sum_numer, sum_denom
70 INTEGER :: n,den,vel,pre
77 IF (sim%GetRank().EQ.0)
THEN
86 ALLOCATE(sum_numer(sim%Physics%VNUM),sum_denom(sim%Physics%VNUM),sigma(sim%Physics%VNUM))
89 CALL initdata(sim%Mesh, sim%Physics, sim%Timedisc)
90 CALL run(sim%Mesh, sim%Physics, sim%Timedisc)
94 IF (
ASSOCIATED(sim%Timedisc%solution))
THEN
95 DO n=1,sim%Physics%VNUM
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)
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)
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)
111 WHERE (abs(sum_denom(:)).GT.2*tiny(sum_denom))
112 sigma(:) = sum_numer(:) / sum_denom(:)
116 den = sim%Physics%DENSITY
117 vel = sim%Physics%XVELOCITY
118 pre = sim%Physics%PRESSURE
121 IF (sim%GetRank().EQ.0)
THEN
123tap_check(ok,
"stoptime reached")
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%")
135 DEALLOCATE(sim,sigma,sum_numer,sum_denom)
144 TYPE(dict_typ),
POINTER :: config
148 TYPE(dict_typ),
POINTER :: mesh, physics, boundary, datafile, &
150 REAL :: x1,x2,y1,y2,z1,z2
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
176 bc(west) = reflecting
177 bc(east) = no_gradients
189 bc(west) = reflecting
190 bc(east) = no_gradients
202 bc(west) = no_gradients
203 bc(east) = no_gradients
204 bc(south) = no_gradients
205 bc(north) = no_gradients
206 bc(bottom)= no_gradients
207 bc(top) = no_gradients
236 CALL sim%Physics%Error(
"InitProgram",
"geometry not supported for 3D Sedov explosion")
241 "meshtype" / midpoint, &
256 "western" / bc(west), &
257 "eastern" / bc(east), &
258 "southern" / bc(south), &
259 "northern" / bc(north), &
260 "bottomer" / bc(bottom), &
272 "variables" / primitive, &
277 "method" / dormand_prince, &
281 "dtlimit" / 1.0e-13, &
284 "output/solution" / 1, &
285 "maxiter" / 1000000 )
289 "fileformat" / vtk, &
290 "filename" / (trim(odir) // trim(ofname)), &
295 "physics" / physics, &
296 "boundary" / boundary, &
298 "timedisc" / timedisc, &
299 "datafile" / datafile )
306 CLASS(physics_base),
INTENT(IN) :: Physics
307 CLASS(mesh_base),
INTENT(IN) :: Mesh
308 CLASS(timedisc_base),
INTENT(INOUT) :: Timedisc
315 SELECT TYPE (phys => physics)
317 SELECT TYPE (pvar => timedisc%pvar)
320 pvar%density%data1d(:) = rho0
322 pvar%velocity%data1d(:) = 0.0
325 p1 = 3.*(phys%gamma - 1.0)*e1 / ((n + 1)*pi*r0**n)
326 WHERE (mesh%radius%bcenter(:,:,:).LE.r0)
328 pvar%pressure%data3d(:,:,:) = p1
331 pvar%pressure%data3d(:,:,:) = p0
335 CALL phys%Error(
"sedov3d:InitData",
"statevector must be of class euler")
339 CALL phys%Error(
"sedov3d:InitData",
"physics not supported")
342 CALL physics%Convert2Conservative(timedisc%pvar,timedisc%cvar)
343 CALL mesh%Info(
" DATA-----> initial condition: 3D Sedov explosion")
347 SUBROUTINE run(Mesh,Physics,Timedisc)
351 CLASS(mesh_base),
INTENT(IN) :: Mesh
352 CLASS(physics_base),
INTENT(IN) :: Physics
353 CLASS(timedisc_base),
INTENT(INOUT) :: Timedisc
357 TYPE(marray_base) :: er
360 IF (
ASSOCIATED(timedisc%solution))
THEN
365 IF (sim%GetRank().EQ.0) &
367 CALL sedov%PrintConfiguration()
369 timedisc%solution = timedisc%pvar
372 t0 = sqrt((e1/rho0)*(r0/sedov%alpha)**5)
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)
381 DO WHILE((timedisc%maxiter.LE.0).OR.(sim%iter.LE.timedisc%maxiter))
382 SELECT TYPE(solution => timedisc%solution)
383 TYPE IS (statevector_euler)
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))
393 IF (btest(mesh%VECTOR_COMPONENTS,n-1))
THEN
394 solution%velocity%data4d(i,j,k,m) = vr * er%data4d(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 sedov3d
3D Sedov explosion