diff options
Diffstat (limited to 'Carpet/CarpetLib/src/prolongate_3d_real8_minmod.F77')
| -rw-r--r-- | Carpet/CarpetLib/src/prolongate_3d_real8_minmod.F77 | 253 |
1 files changed, 0 insertions, 253 deletions
diff --git a/Carpet/CarpetLib/src/prolongate_3d_real8_minmod.F77 b/Carpet/CarpetLib/src/prolongate_3d_real8_minmod.F77 deleted file mode 100644 index 32c3e6227..000000000 --- a/Carpet/CarpetLib/src/prolongate_3d_real8_minmod.F77 +++ /dev/null @@ -1,253 +0,0 @@ -c -*-Fortran-*- - -#include "cctk.h" -#include "cctk_Parameters.h" - - -c$$$ This routine performs "TVD" prolongation. It is intended to be used -c$$$ with GFs that are not expected to be smooth, particularly those -c$$$ that must also obey certain constraints. The obvious example is the -c$$$ density in hydrodynamics, which may be discontinuous yet must be -c$$$ strictly positive. -c$$$ -c$$$ To ensure that this prolongation method is used you should add the -c$$$ tag -c$$$ -c$$$ tags='Prolongation="TVD"' -c$$$ -c$$$ to the interface.ccl on the appropriate group. -c$$$ -c$$$ This applies minmod type limiting to the slope, checking over the -c$$$ entire coarse grid cell for the minimum modulus in each direction. -c$$$ -c$$$ The actual minmod function is defined in the routine -c$$$ -c$$$ prolongate_3d_real8_minmod.F77 - - - function minmod(a, b) - - implicit none - - CCTK_REAL8 minmod - CCTK_REAL8 a, b - CCTK_REAL8 zero - parameter (zero = 0) - - if (a * b .lt. zero) then - minmod = zero - else if (abs(a) .lt. abs(b)) then - minmod = a - else - minmod = b - end if - - end - - subroutine prolongate_3d_real8_minmod ( - $ src, srciext, srcjext, srckext, - $ dst, dstiext, dstjext, dstkext, - $ srcbbox, dstbbox, regbbox) - - implicit none - - DECLARE_CCTK_PARAMETERS - - CCTK_REAL8 one - parameter (one = 1) - - integer srciext, srcjext, srckext - CCTK_REAL8 src(srciext,srcjext,srckext) - integer dstiext, dstjext, dstkext - CCTK_REAL8 dst(dstiext,dstjext,dstkext) -c bbox(:,1) is lower boundary (inclusive) -c bbox(:,2) is upper boundary (inclusive) -c bbox(:,3) is stride - integer srcbbox(3,3), dstbbox(3,3), regbbox(3,3) - - integer offsetlo, offsethi - - integer regiext, regjext, regkext - - integer dstifac, dstjfac, dstkfac - - integer srcioff, srcjoff, srckoff - integer dstioff, dstjoff, dstkoff - - integer i, j, k - integer i0, j0, k0 - integer fi, fj, fk - integer ii, jj, kk - integer d - - external minmod - CCTK_REAL8 minmod - - CCTK_REAL8 half, zero - parameter (half = 0.5) - parameter (zero = 0) - CCTK_REAL8 dupw, dloc, slopex, slopey, slopez - logical firstloop - - - do d=1,3 - if (srcbbox(d,3).eq.0 .or. dstbbox(d,3).eq.0 - $ .or. regbbox(d,3).eq.0) then - call CCTK_WARN (0, "Internal error: stride is zero") - end if - if (srcbbox(d,3).le.regbbox(d,3) - $ .or. dstbbox(d,3).ne.regbbox(d,3)) then - call CCTK_WARN (0, "Internal error: strides disagree") - end if - if (mod(srcbbox(d,3), dstbbox(d,3)).ne.0) then - call CCTK_WARN (0, "Internal error: destination strides are not integer multiples of the source strides") - end if - if (mod(srcbbox(d,1), srcbbox(d,3)).ne.0 - $ .or. mod(dstbbox(d,1), dstbbox(d,3)).ne.0 - $ .or. mod(regbbox(d,1), regbbox(d,3)).ne.0) then - call CCTK_WARN (0, "Internal error: array origins are not integer multiples of the strides") - end if - if (regbbox(d,1).gt.regbbox(d,2)) then -c This could be handled, but is likely to point to an error elsewhere - call CCTK_WARN (0, "Internal error: region extent is empty") - end if - if (mod(srcbbox(d,2) - srcbbox(d,1), srcbbox(d,3)).ne.0 - $ .or. mod(dstbbox(d,2) - dstbbox(d,1), dstbbox(d,3)).ne.0 - $ .or. mod(regbbox(d,2) - regbbox(d,1), regbbox(d,3)).ne.0) then - call CCTK_WARN (0, "Internal error: array extents are not integer multiples of the strides") - end if - regkext = (regbbox(d,2) - regbbox(d,1)) / regbbox(d,3) + 1 - dstkfac = srcbbox(d,3) / dstbbox(d,3) - srckoff = (regbbox(d,1) - srcbbox(d,1)) / dstbbox(d,3) - offsetlo = regbbox(d,3) - if (mod(srckoff + 0, dstkfac).eq.0) then - offsetlo = 0 - if (regkext.gt.1) then - offsetlo = regbbox(d,3) - end if - end if - offsethi = regbbox(d,3) - if (mod(srckoff + regkext-1, dstkfac).eq.0) then - offsethi = 0 - if (regkext.gt.1) then - offsethi = regbbox(d,3) - end if - end if - if (regbbox(d,1)-offsetlo.lt.srcbbox(d,1) - $ .or. regbbox(d,2)+offsethi.gt.srcbbox(d,2) - $ .or. regbbox(d,1).lt.dstbbox(d,1) - $ .or. regbbox(d,2).gt.dstbbox(d,2)) then - call CCTK_WARN (0, "Internal error: region extent is not contained in array extent") - end if - end do - - if (srciext.ne.(srcbbox(1,2)-srcbbox(1,1))/srcbbox(1,3)+1 - $ .or. srcjext.ne.(srcbbox(2,2)-srcbbox(2,1))/srcbbox(2,3)+1 - $ .or. srckext.ne.(srcbbox(3,2)-srcbbox(3,1))/srcbbox(3,3)+1 - $ .or. dstiext.ne.(dstbbox(1,2)-dstbbox(1,1))/dstbbox(1,3)+1 - $ .or. dstjext.ne.(dstbbox(2,2)-dstbbox(2,1))/dstbbox(2,3)+1 - $ .or. dstkext.ne.(dstbbox(3,2)-dstbbox(3,1))/dstbbox(3,3)+1) then - call CCTK_WARN (0, "Internal error: array sizes don't agree with bounding boxes") - end if - - - - regiext = (regbbox(1,2) - regbbox(1,1)) / regbbox(1,3) + 1 - regjext = (regbbox(2,2) - regbbox(2,1)) / regbbox(2,3) + 1 - regkext = (regbbox(3,2) - regbbox(3,1)) / regbbox(3,3) + 1 - - dstifac = srcbbox(1,3) / dstbbox(1,3) - dstjfac = srcbbox(2,3) / dstbbox(2,3) - dstkfac = srcbbox(3,3) / dstbbox(3,3) - - srcioff = (regbbox(1,1) - srcbbox(1,1)) / dstbbox(1,3) - srcjoff = (regbbox(2,1) - srcbbox(2,1)) / dstbbox(2,3) - srckoff = (regbbox(3,1) - srcbbox(3,1)) / dstbbox(3,3) - - dstioff = (regbbox(1,1) - dstbbox(1,1)) / dstbbox(1,3) - dstjoff = (regbbox(2,1) - dstbbox(2,1)) / dstbbox(2,3) - dstkoff = (regbbox(3,1) - dstbbox(3,1)) / dstbbox(3,3) - - - -c Loop over fine region - - do k = 0, regkext-1 - k0 = (srckoff + k) / dstkfac - fk = mod(srckoff + k, dstkfac) - - do j = 0, regjext-1 - j0 = (srcjoff + j) / dstjfac - fj = mod(srcjoff + j, dstjfac) - - do i = 0, regiext-1 - i0 = (srcioff + i) / dstifac - fi = mod(srcioff + i, dstifac) - - slopex = zero - slopey = zero - slopez = zero - - firstloop = .true. - - do kk = 1, 2 - do jj = 1, 2 - - dupw = src(i0+1 ,j0+jj,k0+kk) - src(i0+0 ,j0+jj,k0+kk) - dloc = src(i0+2 ,j0+jj,k0+kk) - src(i0+1 ,j0+kk,k0+kk) - if (firstloop) then - slopex = half * dble(fi) * minmod(dupw,dloc) - firstloop = .false. - else - slopex = - $ minmod(slopex, half * dble(fi) * minmod(dupw,dloc)) - end if - end do - end do - - firstloop = .true. - - do kk = 1, 2 - do ii = 1, 2 - - dupw = src(i0+ii,j0+1 ,k0+kk) - src(i0+ii,j0+0 ,k0+kk) - dloc = src(i0+ii,j0+2 ,k0+kk) - src(i0+ii,j0+1 ,k0+kk) - if (firstloop) then - slopey = half * dble(fj) * minmod(dupw,dloc) - firstloop = .false. - else - slopey = - $ minmod(slopey, half * dble(fj) * minmod(dupw,dloc)) - end if - end do - end do - - firstloop = .true. - - do jj = 1, 2 - do ii = 1, 2 - - dupw = src(i0+ii,j0+jj,k0+1 ) - src(i0+ii,j0+jj,k0+0 ) - dloc = src(i0+ii,j0+jj,k0+2 ) - src(i0+ii,j0+jj,k0+1 ) - if (firstloop) then - slopez = half * dble(fk) * minmod(dupw,dloc) - firstloop = .false. - else - slopez = - $ minmod(slopez, half * dble(fk) * minmod(dupw,dloc)) - end if - end do - end do - - if (check_array_accesses.ne.0) then - call checkinde (dstioff+i+1, dstjoff+j+1, dstkoff+k+1, 1,1,1, dstiext,dstjext,dstkext, "destination") - end if - dst (dstioff+i+1, dstjoff+j+1, dstkoff+k+1) = - . src(i0+1,j0+1,k0+1) + slopex + slopey + slopez - - - end do - end do - end do - - end |