471 INTEGER,
PARAMETER :: select_test = 6
472 LOGICAL,
PARAMETER :: run_all = .true.
473 REAL,
PARAMETER :: gamma = 1.4
475 INTEGER,
PARAMETER :: mgeo = cartesian
478 INTEGER,
PARAMETER :: xres = 100
479 INTEGER,
PARAMETER :: yres = 100
480 INTEGER,
PARAMETER :: zres = 1
481 REAL,
PARAMETER :: rmin = 1.0e-4
483 INTEGER,
PARAMETER :: onum = 1
484 CHARACTER(LEN=256),
PARAMETER &
486 CHARACTER(LEN=256),
PARAMETER &
487 :: ofname =
'riemann2d'
489 INTEGER,
PARAMETER :: num_tests = 19
490 REAL,
PARAMETER,
DIMENSION(NUM_TESTS) :: &
491 test_stoptime = (/ 0.2,0.2,0.3,0.25,0.23,0.3,0.25, &
492 0.25,0.3,0.15,0.3,0.25,0.3,0.1, &
493 0.2,0.2,0.3,0.2,0.3/)
495 CLASS(
fosite),
ALLOCATABLE :: sim
496 LOGICAL :: ok(num_tests), mpifinalize = .false.
507 IF (.NOT.run_all.AND.ic.NE.select_test) cycle
509 IF (
ALLOCATED(sim))
DEALLOCATE(sim)
510 ALLOCATE(sim,stat=err)
512 CALL sim%Error(
"riemann2d",
"cannot allocate memory for Sim")
514 CALL sim%InitFosite()
522 CALL initdata(sim%Mesh, sim%Physics, sim%Timedisc, ic)
524 ok(ic) = .NOT.sim%aborted
526 IF ((run_all.AND.ic.EQ.num_tests).OR.(.NOT.run_all)) mpifinalize=.true.
527 CALL sim%Finalize(mpifinalize)
532 IF (.NOT.run_all.AND.ic.NE.select_test) cycle
533tap_check(ok(ic),
"Simulation finished")
536 IF (
ALLOCATED(sim))
DEALLOCATE(sim)
547 TYPE(dict_typ),
POINTER :: config
552 TYPE(dict_typ),
POINTER :: mesh, physics, boundary, datafile, &
554 REAL :: x1,x2,y1,y2,z1,z2,sc
555 CHARACTER(LEN=16) :: fext
559 CHARACTER(LEN=8) :: geo_str
561 IF (ic.LT.1.OR.ic.GT.num_tests) &
562 CALL sim%Error(
"riemann2d::MakeConfig",
"invalid test number selected")
574 bc(west) = no_gradients
575 bc(east) = no_gradients
576 bc(south) = no_gradients
577 bc(north) = no_gradients
578 bc(bottom)= no_gradients
579 bc(top) = no_gradients
590 bc(east) = no_gradients
593 bc(bottom) = reflecting
599 x2 = log(0.5*sqrt(2.0)/sc)
605 bc(east) = no_gradients
608 bc(bottom) = no_gradients
609 bc(top) = no_gradients
612 CALL sim%Error(
"riemann2d::MakeConfig",
"geometry not supported in this test")
616 mesh => dict(
"meshtype" / midpoint, &
627 "output/commutator" / 1,&
631 boundary => dict(
"western" / bc(west), &
632 "eastern" / bc(east), &
633 "southern" / bc(south), &
634 "northern" / bc(north), &
635 "bottomer" / bc(bottom),&
639 physics => dict(
"problem" / euler, &
643 fluxes => dict(
"order" / linear, &
645 "variables" / conservative, &
646 "limiter" / monocent, &
651 "method" / modified_euler, &
654 "stoptime" / test_stoptime(ic), &
655 "dtlimit" / 1.0e-10, &
656 "maxiter" / 10000000, &
657 "output/geometrical_sources" / 1)
661 WRITE (fext,
'(A,I2.2)') trim(geo_str) //
"_", ic
664 "fileformat" / vtk, &
666 "filename" / (trim(odir) // trim(ofname) // trim(fext)), &
669 config => dict(
"mesh" / mesh, &
670 "physics" / physics, &
671 "boundary" / boundary, &
673 "timedisc" / timedisc, &
675 "datafile" / datafile)
682 CLASS(physics_base) :: Physics
683 CLASS(mesh_base) :: Mesh
684 CLASS(timedisc_base):: Timedisc
688 TYPE(marray_base) :: vcart,vcurv
689 REAL :: xmin,ymin,xmax,ymax,x0,y0
691 REAL :: xmin_all,xmax_all,ymin_all,ymax_all
694 REAL,
DIMENSION(NUM_TESTS,4):: den,pre,vx,vy
695 CHARACTER(LEN=64) :: teststr
697 INTENT(IN) :: mesh,physics,ic
698 INTENT(INOUT) :: timedisc
701 xmin = minval(mesh%bccart(:,:,:,1))
702 ymin = minval(mesh%bccart(:,:,:,2))
703 xmax = maxval(mesh%bccart(:,:,:,1))
704 ymax = maxval(mesh%bccart(:,:,:,2))
706 CALL mpi_allreduce(xmin,xmin_all,1,default_mpi_real,mpi_min,mesh%comm_cart,ierr)
708 CALL mpi_allreduce(ymin,ymin_all,1,default_mpi_real,mpi_min,mesh%comm_cart,ierr)
710 CALL mpi_allreduce(xmax,xmax_all,1,default_mpi_real,mpi_max,mesh%comm_cart,ierr)
712 CALL mpi_allreduce(ymax,ymax_all,1,default_mpi_real,mpi_max,mesh%comm_cart,ierr)
715 x0 = xmin + 0.5*abs(xmax-xmin)
716 y0 = ymin + 0.5*abs(ymax-ymin)
718 vcart = marray_base(3)
719 vcurv = marray_base(3)
727 IF (ic.LT.1.OR.ic.GT.num_tests) &
728 CALL mesh%Error(
"InitData",
"Sorry, this 2D Riemann problem is currently not supported!")
730 WRITE (teststr,
'(A,I2)')
"2D Riemann problem no. ", ic
1136 SELECT TYPE(pvar => timedisc%pvar)
1137 TYPE IS(statevector_euler)
1138 WHERE ( (mesh%bccart(:,:,:,1).GE.x0).AND.(mesh%bccart(:,:,:,2).GE.y0) )
1139 pvar%density%data3d(:,:,:) = den(ic,1)
1140 pvar%pressure%data3d(:,:,:) = pre(ic,1)
1141 vcart%data4d(:,:,:,1) = vx(ic,1)
1142 vcart%data4d(:,:,:,2) = vy(ic,1)
1143 ELSEWHERE ( (mesh%bccart(:,:,:,1).LT.x0).AND.(mesh%bccart(:,:,:,2).GE.y0) )
1144 pvar%density%data3d(:,:,:) = den(ic,2)
1145 pvar%pressure%data3d(:,:,:) = pre(ic,2)
1146 vcart%data4d(:,:,:,1) = vx(ic,2)
1147 vcart%data4d(:,:,:,2) = vy(ic,2)
1148 ELSEWHERE ( (mesh%bccart(:,:,:,1).LT.x0).AND.(mesh%bccart(:,:,:,2).LT.y0) )
1149 pvar%density%data3d(:,:,:) = den(ic,3)
1150 pvar%pressure%data3d(:,:,:) = pre(ic,3)
1151 vcart%data4d(:,:,:,1) = vx(ic,3)
1152 vcart%data4d(:,:,:,2) = vy(ic,3)
1153 ELSEWHERE ( (mesh%bccart(:,:,:,1).GE.x0).AND.(mesh%bccart(:,:,:,2).LT.y0) )
1154 pvar%density%data3d(:,:,:) = den(ic,4)
1155 pvar%pressure%data3d(:,:,:) = pre(ic,4)
1156 vcart%data4d(:,:,:,1) = vx(ic,4)
1157 vcart%data4d(:,:,:,2) = vy(ic,4)
1159 vcart%data2d(:,3) = 0.0
1160 CALL mesh%geometry%Convert2Curvilinear(mesh%bcenter,vcart%data4d,vcurv%data4d)
1161 pvar%velocity%data2d(:,1:2) = vcurv%data2d(:,1:2)
1164 CALL physics%Convert2Conservative(timedisc%pvar,timedisc%cvar)
1165 CALL mesh%Info(
" DATA-----> initial condition: " // trim(teststr))
subroutine initdata(Mesh, Physics, Fluxes, Timedisc)
subroutine makeconfig(Sim, config)
elemental real function, public gamma(x)
Compute the Gamma function for arguments != 0,-1,-2,..
elemental real function, public asinh(x)
inverse hyperbolic sine function