76 INTEGER,
PARAMETER ::
pp = 6
81 "minmod ",
"mc ",
"sweby ",
"superbee ",
"ospre ", &
82 "pp ",
"van leer ",
"no limit "/)
104 CHARACTER(LEN=60) :: key
105 INTEGER :: err,limiter,valwrite,i,j,k,l,m
109 ALLOCATE(this%slopes,this%dx, &
112 CALL this%Error(
"InitReconstruction_linear",
"Unable to allocate memory.")
120 CALL getattr(config,
"limiter", limiter,
minmod)
121 CALL getattr(config,
"theta", theta, 1.0)
129 this%limiter_param = theta
132 CALL this%Error(
"InitReconstruction_linear",
"Unknown limiter")
136 IF (limiter.EQ.
pp)
THEN
137 CALL this%Error(
"InitReconstruction_linear",
"PP limiter currently not supported.")
153 this%slopes%data1d(:) = 0.0
157 IF (mesh%INUM.GT.1)
THEN
159 DO concurrent(i=mesh%IGMIN:mesh%IGMAX,j=mesh%JGMIN:mesh%JGMAX,k=mesh%KGMIN:mesh%KGMAX)
160 this%dx%data4d(i,j,k,m) = mesh%curv%faces(i,j,k,1,1) - mesh%curv%bcenter(i,j,k,1)
161 this%dx%data4d(i,j,k,m+1) = mesh%curv%faces(i,j,k,2,1) - mesh%curv%bcenter(i,j,k,1)
165 IF (mesh%JNUM.GT.1)
THEN
166 DO concurrent(i=mesh%IGMIN:mesh%IGMAX,j=mesh%JGMIN:mesh%JGMAX,k=mesh%KGMIN:mesh%KGMAX)
167 this%dx%data4d(i,j,k,m) = mesh%curv%faces(i,j,k,3,2) - mesh%curv%bcenter(i,j,k,2)
168 this%dx%data4d(i,j,k,m+1) = mesh%curv%faces(i,j,k,4,2) - mesh%curv%bcenter(i,j,k,2)
172 IF (mesh%KNUM.GT.1)
THEN
173 DO concurrent(i=mesh%IGMIN:mesh%IGMAX,j=mesh%JGMIN:mesh%JGMAX,k=mesh%KGMIN:mesh%KGMAX)
174 this%dx%data4d(i,j,k,m) = mesh%curv%faces(i,j,k,5,3) - mesh%curv%bcenter(i,j,k,3)
175 this%dx%data4d(i,j,k,m+1) = mesh%curv%faces(i,j,k,6,3) - mesh%curv%bcenter(i,j,k,3)
179 CALL getattr(config,
"output/slopes", valwrite, 0)
180 IF(valwrite.EQ.1)
THEN
182 IF (mesh%INUM.GT.1)
THEN
185 key = trim(physics%pvarname(l)) //
"_xslope"
186 CALL setattr(io, trim(key), &
187 this%slopes%data5d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,m,l))
191 IF (mesh%JNUM.GT.1)
THEN
194 key = trim(physics%pvarname(l)) //
"_yslope"
195 CALL setattr(io, trim(key), &
196 this%slopes%data5d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,m,l))
200 IF (mesh%KNUM.GT.1)
THEN
203 key = trim(physics%pvarname(l)) //
"_zslope"
204 CALL setattr(io, trim(key), &
205 this%slopes%data5d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,m,l))
210 CALL this%Info(
" limiter: " // trim(this%GetName()))
227 CALL this%CalculateSlopes(mesh,physics,rvar)
234 DO concurrent(i=1:
SIZE(rvar%data2d,dim=1))
235 rstates%data3d(i,n,l) = rvar%data2d(i,l) + &
236 this%slopes%data3d(i,((n+1)/2),l)*this%dx%data2d(i,n)
255 SELECT CASE(this%GetType())
258 IF (mesh%INUM.GT.1)
THEN
261 DO k=mesh%KGMIN,mesh%KGMAX
262 DO j=mesh%JGMIN,mesh%JGMAX
264 DO concurrent(i=mesh%IGMIN+mesh%IP1:mesh%IGMAX+mesh%IM1)
266 rvar%data4d(i,j,k,l) - rvar%data4d(i-1,j,k,l), &
267 rvar%data4d(i+1,j,k,l) - rvar%data4d(i,j,k,l))
275 IF (mesh%JNUM.GT.1)
THEN
276 i = mesh%IGMAX-mesh%IGMIN+1
280 DO k=mesh%KGMIN,mesh%KGMAX
283 DO j=i+1,
SIZE(rvar%data3d,1)-i
285 rvar%data3d(j,k,l) - rvar%data3d(j-i,k,l), &
286 rvar%data3d(j+i,k,l) - rvar%data3d(j,k,l))
293 IF (mesh%KNUM.GT.1)
THEN
294 i = (mesh%IGMAX-mesh%IGMIN+1)*(mesh%JGMAX-mesh%JGMIN+1)
299 DO j=i+1,
SIZE(rvar%data2d,1)-i
301 rvar%data2d(j,l) - rvar%data2d(j-i,l), &
302 rvar%data2d(j+i,l) - rvar%data2d(j,l))
308 IF (mesh%INUM.GT.1)
THEN
311 DO k=mesh%KGMIN,mesh%KGMAX
312 DO j=mesh%JGMIN,mesh%JGMAX
314 DO concurrent(i=mesh%IGMIN+mesh%IP1:mesh%IGMAX+mesh%IM1)
316 this%limiter_param*(rvar%data4d(i,j,k,l) - rvar%data4d(i-1,j,k,l)),&
317 this%limiter_param*(rvar%data4d(i+1,j,k,l) - rvar%data4d(i,j,k,l)),&
318 0.5*(rvar%data4d(i+1,j,k,l) - rvar%data4d(i-1,j,k,l)))
326 IF (mesh%JNUM.GT.1)
THEN
327 i = mesh%IGMAX-mesh%IGMIN+1
330 DO k=mesh%KGMIN,mesh%KGMAX
333 DO j=i+1,
SIZE(rvar%data3d,1)-i
335 this%limiter_param*(rvar%data3d(j,k,l) - rvar%data3d(j-i,k,l)),&
336 this%limiter_param*(rvar%data3d(j+i,k,l) - rvar%data3d(j,k,l)),&
337 0.5*(rvar%data3d(j+i,k,l) - rvar%data3d(j-i,k,l)))
344 IF (mesh%KNUM.GT.1)
THEN
345 i = (mesh%IGMAX-mesh%IGMIN+1)*(mesh%JGMAX-mesh%JGMIN+1)
350 DO j=i+1,
SIZE(rvar%data2d,1)-i
352 this%limiter_param*(rvar%data2d(j,l) - rvar%data2d(j-i,l)),&
353 this%limiter_param*(rvar%data2d(j+i,l) - rvar%data2d(j,l)),&
354 0.5*(rvar%data2d(j+i,l) - rvar%data2d(j-i,l)))
360 IF (mesh%INUM.GT.1)
THEN
363 DO k=mesh%KGMIN,mesh%KGMAX
364 DO j=mesh%JGMIN,mesh%JGMAX
366 DO concurrent(i=mesh%IGMIN+mesh%IP1:mesh%IGMAX+mesh%IM1)
368 rvar%data4d(i,j,k,l) - rvar%data4d(i-1,j,k,l), &
369 rvar%data4d(i+1,j,k,l) - rvar%data4d(i,j,k,l),&
378 IF (mesh%JNUM.GT.1)
THEN
379 i = mesh%IGMAX-mesh%IGMIN+1
382 DO k=mesh%KGMIN,mesh%KGMAX
385 DO j=i+1,
SIZE(rvar%data3d,1)-i
387 rvar%data3d(j,k,l) - rvar%data3d(j-i,k,l),&
388 rvar%data3d(j+i,k,l) - rvar%data3d(j,k,l),&
396 IF (mesh%KNUM.GT.1)
THEN
397 i = (mesh%IGMAX-mesh%IGMIN+1)*(mesh%JGMAX-mesh%JGMIN+1)
402 DO j=i+1,
SIZE(rvar%data2d,1)-i
404 rvar%data2d(j,l) - rvar%data2d(j-i,l),&
405 rvar%data2d(j+i,l) - rvar%data2d(j,l),&
412 IF (mesh%INUM.GT.1)
THEN
415 DO k=mesh%KGMIN,mesh%KGMAX
416 DO j=mesh%JGMIN,mesh%JGMAX
418 DO concurrent(i=mesh%IGMIN+mesh%IP1:mesh%IGMAX+mesh%IM1)
420 rvar%data4d(i,j,k,l) - rvar%data4d(i-1,j,k,l), &
421 rvar%data4d(i+1,j,k,l) - rvar%data4d(i,j,k,l), 2.0)
429 IF (mesh%JNUM.GT.1)
THEN
430 i = mesh%IGMAX-mesh%IGMIN+1
433 DO k=mesh%KGMIN,mesh%KGMAX
436 DO j=i+1,
SIZE(rvar%data3d,1)-i
438 rvar%data3d(j,k,l) - rvar%data3d(j-i,k,l),&
439 rvar%data3d(j+i,k,l) - rvar%data3d(j,k,l), 2.0)
446 IF (mesh%KNUM.GT.1)
THEN
447 i = (mesh%IGMAX-mesh%IGMIN+1)*(mesh%JGMAX-mesh%JGMIN+1)
452 DO j=i+1,
SIZE(rvar%data2d,1)-i
454 rvar%data2d(j,l) - rvar%data2d(j-i,l),&
455 rvar%data2d(j+i,l) - rvar%data2d(j,l), 2.0)
461 IF (mesh%INUM.GT.1)
THEN
464 DO k=mesh%KGMIN,mesh%KGMAX
465 DO j=mesh%JGMIN,mesh%JGMAX
467 DO concurrent(i=mesh%IGMIN+mesh%IP1:mesh%IGMAX+mesh%IM1)
469 rvar%data4d(i,j,k,l) - rvar%data4d(i-1,j,k,l), &
470 rvar%data4d(i+1,j,k,l) - rvar%data4d(i,j,k,l))
478 IF (mesh%JNUM.GT.1)
THEN
479 i = mesh%IGMAX-mesh%IGMIN+1
482 DO k=mesh%KGMIN,mesh%KGMAX
485 DO j=i+1,
SIZE(rvar%data3d,1)-i
487 rvar%data3d(j,k,l) - rvar%data3d(j-i,k,l),&
488 rvar%data3d(j+i,k,l) - rvar%data3d(j,k,l))
495 IF (mesh%KNUM.GT.1)
THEN
496 i = (mesh%IGMAX-mesh%IGMIN+1)*(mesh%JGMAX-mesh%JGMIN+1)
501 DO j=i+1,
SIZE(rvar%data2d,1)-i
503 rvar%data2d(j,l) - rvar%data2d(j-i,l),&
504 rvar%data2d(j+i,l) - rvar%data2d(j,l))
533 IF (mesh%INUM.GT.1)
THEN
536 DO k=mesh%KGMIN,mesh%KGMAX
537 DO j=mesh%JGMIN,mesh%JGMAX
539 DO concurrent(i=mesh%IGMIN+mesh%IP1:mesh%IGMAX+mesh%IM1)
541 rvar%data4d(i,j,k,l) - rvar%data4d(i-1,j,k,l), &
542 rvar%data4d(i+1,j,k,l) - rvar%data4d(i,j,k,l))
550 IF (mesh%JNUM.GT.1)
THEN
551 i = mesh%IGMAX-mesh%IGMIN+1
554 DO k=mesh%KGMIN,mesh%KGMAX
557 DO j=i+1,
SIZE(rvar%data3d,1)-i
559 rvar%data3d(j,k,l) - rvar%data3d(j-i,k,l),&
560 rvar%data3d(j+i,k,l) - rvar%data3d(j,k,l))
567 IF (mesh%KNUM.GT.1)
THEN
568 i = (mesh%IGMAX-mesh%IGMIN+1)*(mesh%JGMAX-mesh%JGMIN+1)
573 DO j=i+1,
SIZE(rvar%data2d,1)-i
575 rvar%data2d(j,l) - rvar%data2d(j-i,l),&
576 rvar%data2d(j+i,l) - rvar%data2d(j,l))
582 IF (mesh%INUM.GT.1)
THEN
585 DO k=mesh%KGMIN,mesh%KGMAX
586 DO j=mesh%JGMIN,mesh%JGMAX
588 DO concurrent(i=mesh%IGMIN+mesh%IP1:mesh%IGMAX+mesh%IM1)
590 rvar%data4d(i,j,k,l) - rvar%data4d(i-1,j,k,l), &
591 rvar%data4d(i+1,j,k,l) - rvar%data4d(i,j,k,l))
599 IF (mesh%JNUM.GT.1)
THEN
600 i = mesh%IGMAX-mesh%IGMIN+1
603 DO k=mesh%KGMIN,mesh%KGMAX
606 DO j=i+1,
SIZE(rvar%data3d,1)-i
608 rvar%data3d(j,k,l) - rvar%data3d(j-i,k,l),&
609 rvar%data3d(j+i,k,l) - rvar%data3d(j,k,l))
616 IF (mesh%KNUM.GT.1)
THEN
617 i = (mesh%IGMAX-mesh%IGMIN+1)*(mesh%JGMAX-mesh%JGMIN+1)
622 DO j=i+1,
SIZE(rvar%data2d,1)-i
624 rvar%data2d(j,l) - rvar%data2d(j-i,l),&
625 rvar%data2d(j+i,l) - rvar%data2d(j,l))
637 DEALLOCATE(this%slopes,this%dx)
639 CALL this%Finalize_base()
649 REAL,
INTENT(IN) :: arg1, arg2
651 IF (sign(1.0,arg1)*sign(1.0,arg2).GT.0)
THEN
652 limarg = sign(min(abs(arg1),abs(arg2)),arg1)
662 REAL,
INTENT(IN) :: arg1, arg2, arg3
664 IF (((sign(1.0,arg1)*sign(1.0,arg2)).GT.0).AND.&
665 ((sign(1.0,arg2)*sign(1.0,arg3)).GT.0))
THEN
666 limarg = sign(min(abs(arg1),abs(arg2),abs(arg3)),arg1)
676 REAL,
INTENT(IN) :: arg1, arg2, param
678 IF (sign(1.0,arg1)*sign(1.0,arg2).GT.0)
THEN
679 limarg = sign(max(min(param*abs(arg1),abs(arg2)), &
680 min(abs(arg1),param*abs(arg2))),arg1)
690 REAL,
INTENT(IN) :: arg1, arg2
692 limarg = 1.5*arg1*arg2*(arg1 + arg2) / (arg1*(arg1 + 0.5*arg2) &
693 + arg2*(arg2 + 0.5*arg1) + tiny(1.0))
717 REAL,
INTENT(IN) :: arg1, arg2
723 limarg = (arg1 * a2 + arg2 * a1)/(a1 + a2 + tiny(a1))
730 REAL,
INTENT(IN) :: arg1, arg2
732 limarg = 0.5*(arg1+arg2)
Dictionary for generic data types.
subroutine finalize(this)
Destructor of logging_base class.
base class for mesh arrays
derived class for compound of mesh arrays
base module for reconstruction process
integer, parameter, public linear
module for linear (first order) TVD reconstruction using slope limiters
subroutine calculateslopes(this, Mesh, Physics, rvar)
computes the limited slopes
integer, parameter, public nolimit
integer, parameter, public sweby
integer, parameter, public pp
elemental real function nolimit_limiter(arg1, arg2)
integer, parameter, public vanleer
integer, parameter, public superbee
integer, parameter, public minmod
subroutine initreconstruction_linear(this, Mesh, Physics, config, IO)
Constructor of linear reconstruction module.
character(len=32), parameter recontype_name
elemental real function sweby_limiter(arg1, arg2, param)
elemental real function minmod3_limiter(arg1, arg2, arg3)
elemental real function ospre_limiter(arg1, arg2)
character(len=32), dimension(8), parameter limitertype_name
subroutine calculatestates(this, Mesh, Physics, rvar, rstates)
Reconstructes states at cell boundaries.
elemental real function minmod2_limiter(arg1, arg2)
integer, parameter, public ospre
elemental real function vanleer_limiter(arg1, arg2)
integer, parameter, public monocent