60 INTEGER,
PARAMETER :: testnum = 8
61 CHARACTER(LEN=32),
PARAMETER :: teststr(testnum) = (/ &
68 "Isothermal shock tube ", &
69 "Isothermal Noh problem " /)
72 CHARACTER(LEN=3),
PARAMETER :: testdir =
"all"
73 INTEGER,
PARAMETER :: res = 400
75 INTEGER,
PARAMETER :: onum = 1
76 CHARACTER(LEN=8),
PARAMETER :: ofname(testnum) = (/ &
77 "sod ",
"toro2 ",
"toro3 ",
"toro4 ",
"toro5 ",
"noh ", &
78 "sodiso ",
"nohiso " /)
79 CHARACTER(LEN=256),
PARAMETER &
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
96 SELECT CASE (trim(testdir))
112 ALLOCATE(sim,sigma(testnum,dir_min:dir_max))
115 CALL sim%InitFosite()
118 IF (sim%GetRank().EQ.0)
THEN
120tap_plan(testnum*(dir_max-dir_min+1))
126 tests:
DO ic=1,testnum
128 DO sd=dir_min,dir_max
139 IF (sim%aborted)
THEN
140 sigma(ic,sd) = huge(1.0)
145 IF (
ASSOCIATED(sim%Timedisc%solution))
THEN
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(:))
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)
157 IF (abs(sum_denom).GT.4*tiny(sum_denom)) &
158 sigma(ic,sd) = sum_numer / sum_denom
163 IF (.NOT.(ic.EQ.testnum.AND.sd.EQ.dir_max))
THEN
164 CALL sim%Finalize(mpifinalize_=.false.)
167 CALL sim%InitFosite()
175 IF (sim%GetRank().EQ.0)
THEN
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)
183tap_check_small(sigma(ic,sd),sigma_tol(ic),trim(verbose_tap_output))
201 TYPE(dict_typ),
POINTER :: config
205 TYPE(dict_typ),
POINTER :: mesh, physics, boundary, datafile, &
207 INTEGER :: xres,yres,zres
208 REAL :: xmax,ymax,zmax
237 CALL sim%Error(
"riemann1d::InitData",
"direction must be one of 1,2,3")
241 mesh => dict(
"meshtype" / midpoint, &
242 "geometry" / cartesian, &
280 CALL sim%Mesh%Error(
"riemann1d::InitData",
"Test problem number should be 1,2,3,4,5,6,7 or 8")
285 IF (csiso.NE.0.0)
THEN
286 physics => dict(
"problem" / euler_isotherm, &
289 physics => dict(
"problem" / euler, &
297 "variables" / conservative, &
298 "limiter" / monocent, &
303 "method" / modified_euler, &
307 "dtlimit" / 1.0e-10, &
313 "output/solution" / 0, &
318 IF (csiso.EQ.0.0)
CALL setattr(timedisc,
"output/solution",1)
322 "western" / no_gradients, &
323 "eastern" / no_gradients, &
324 "southern" / no_gradients, &
325 "northern" / no_gradients, &
326 "bottomer" / no_gradients, &
327 "topper" / no_gradients)
331 "fileformat" / gnuplot, &
332 "filename" / (trim(odir) // trim(ofname(ic)) //
"-" // achar(119+dir)), &
336 config => dict(
"mesh" / mesh, &
337 "physics" / physics, &
338 "boundary" / boundary, &
340 "timedisc" / timedisc, &
341 "datafile" / datafile)
344 SUBROUTINE run(this,ic,dir)
348 CLASS(
fosite),
INTENT(INOUT) :: this
350 INTEGER,
INTENT(IN) :: ic,dir
353 REAL :: x0, rho_l, rho_r, u_l, u_r, p_l, p_r
354 REAL,
DIMENSION(:),
ALLOCATABLE :: x
418 CALL this%Error(
"InitData",
"Test problem should be 1,2,3,4,5,6,7 or 8")
423 SELECT TYPE(pvar => this%Timedisc%pvar)
424 CLASS IS(statevector_eulerisotherm)
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
431 pvar%density%data3d(:,:,:) = rho_r
432 pvar%velocity%data4d(:,:,:,1) = u_r
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
442 pvar%pressure%data3d(:,:,:) = p_r
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)))
450 IF (
ASSOCIATED(this%Timedisc%solution))
THEN
452 this%Timedisc%solution = this%Timedisc%pvar
455 ALLOCATE(x(
SIZE(this%Mesh%cart%data3d(:,2,dir),1)), &
456 source=this%Mesh%cart%data3d(:,2,dir))
458 DO WHILE((this%Timedisc%maxiter.LE.0).OR.(this%iter.LE.this%Timedisc%maxiter))
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(:,:))
subroutine makeconfig(Sim, config)
subroutine, private run(this)
subroutine, public riemann(x0, gamma_, rho_l_, u_l_, p_l_, rho_r_, u_r_, p_r_, t, x, pvar)
program riemann1d
1D Riemann problems