diff options
author | knarf <knarf@c83d129a-5a75-4d5a-9c4d-ed3a5855bf45> | 2010-04-06 17:50:33 +0000 |
---|---|---|
committer | knarf <knarf@c83d129a-5a75-4d5a-9c4d-ed3a5855bf45> | 2010-04-06 17:50:33 +0000 |
commit | d7d1a02dc4d3a3e1202185b4fffaf87bac2fcf32 (patch) | |
tree | 4be7230b09f7e692343d3a25ff73046295d6f464 | |
parent | fda4ba17d2e54a4b4fdcdd27d025c9980c04a610 (diff) |
add excision
however: be carefull with ENO and excision for the time being
git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinEvolve/GRHydro/trunk@94 c83d129a-5a75-4d5a-9c4d-ed3a5855bf45
-rw-r--r-- | param.ccl | 6 | ||||
-rw-r--r-- | src/Whisky_Con2Prim.F90 | 109 | ||||
-rw-r--r-- | src/Whisky_ENOReconstruct.F90 | 70 | ||||
-rw-r--r-- | src/Whisky_PPM.F90 | 100 | ||||
-rw-r--r-- | src/Whisky_Prim2Con.F90 | 7 | ||||
-rw-r--r-- | src/Whisky_Reconstruct.F90 | 138 | ||||
-rw-r--r-- | src/Whisky_ReconstructPoly.F90 | 114 | ||||
-rw-r--r-- | src/Whisky_TVDReconstruct.F90 | 76 | ||||
-rw-r--r-- | src/Whisky_UpdateMask.F90 | 5 |
9 files changed, 340 insertions, 285 deletions
@@ -25,6 +25,7 @@ shares: HydroBase USES CCTK_INT timelevels USES KEYWORD prolongation_type +USES INT hydro_excision EXTENDS KEYWORD evolution_method "" { "whisky" :: "Use Whisky to evolve the hydro variables" @@ -52,6 +53,11 @@ USES INT spatial_order restricted: +CCTK_INT whisky_hydro_excision "Turns excision automatically on in HydroBase" ACCUMULATOR-BASE=HydroBase::hydro_excision +{ + 1:1 :: "Only '1' allowed" +} 1 + CCTK_INT Whisky_MaxNumEvolvedVars "The maximum number of evolved variables used by Whisky" ACCUMULATOR-BASE=MethodofLines::MoL_Num_Evolved_Vars { 5:8 :: "dens scon[3] tau" diff --git a/src/Whisky_Con2Prim.F90 b/src/Whisky_Con2Prim.F90 index 5902c77..19d9a5c 100644 --- a/src/Whisky_Con2Prim.F90 +++ b/src/Whisky_Con2Prim.F90 @@ -146,14 +146,13 @@ subroutine Conservative2Primitive(CCTK_ARGUMENTS) character(len=100) warnline CCTK_INT :: type_bits, atmosphere - CCTK_INT :: type2_bits, excised + CCTK_INT :: type2_bits CCTK_REAL :: local_min_tracer call SpaceMask_GetTypeBits(type_bits, "Hydro_Atmosphere") call SpaceMask_GetStateBits(atmosphere, "Hydro_Atmosphere", "in_atmosphere") - type2_bits = -1 - excised = -1 + type2_bits = -1 nx = cctk_lsh(1) ny = cctk_lsh(2) @@ -175,8 +174,9 @@ subroutine Conservative2Primitive(CCTK_ARGUMENTS) do j = 1, ny do i = 1, nx - !do not compute if in atmosphere - if (SpaceMask_CheckStateBitsF90(space_mask, i, j, k, type_bits, atmosphere)) cycle + !do not compute if in atmosphere or in excised region + if (SpaceMask_CheckStateBitsF90(space_mask, i, j, k, type_bits, atmosphere) .or. & + (hydro_excision_mask(i,j,k) .ne. 0)) cycle epsnegative = .false. @@ -481,6 +481,14 @@ subroutine Con2Prim_pt(handle, dens, sx, sy, sz, tau, rho, velx, vely, & f = pnew - EOS_Pressure(handle, rho, epsilon) + !check whether rho is NaN + if (rho .ne. rho) then + !$OMP CRITICAL + write(warnline,'(a20,g20.7)') 'rho is NaN: ',rho + call CCTK_WARN(Whisky_NaN_verbose,warnline) + !$OMP END CRITICAL + goto 51 + end if enddo @@ -594,7 +602,7 @@ subroutine Conservative2PrimitiveBounds(CCTK_ARGUMENTS) character(len=100) warnline CCTK_INT :: type_bits, atmosphere - CCTK_INT :: type2_bits, excised + CCTK_INT :: type2_bits CCTK_REAL :: local_min_tracer @@ -608,8 +616,7 @@ subroutine Conservative2PrimitiveBounds(CCTK_ARGUMENTS) call SpaceMask_GetTypeBits(type_bits, "Hydro_Atmosphere") call SpaceMask_GetStateBits(atmosphere, "Hydro_Atmosphere", "in_atmosphere") - type2_bits = -1 - excised = -1 + type2_bits = -1 nx = cctk_lsh(1) ny = cctk_lsh(2) @@ -625,8 +632,9 @@ subroutine Conservative2PrimitiveBounds(CCTK_ARGUMENTS) do j = whisky_stencil, ny - whisky_stencil + 1 do i = whisky_stencil, nx - whisky_stencil + 1 - !do not compute if in atmosphere - if (SpaceMask_CheckStateBitsF90(space_mask, i, j, k, type_bits, atmosphere)) cycle + !do not compute if in atmosphere or in an excised region + if (SpaceMask_CheckStateBitsF90(space_mask, i, j, k, type_bits, atmosphere) .or. & + (hydro_excision_mask(i,j,k) .ne. 0)) cycle gxxl = 0.5d0 * (gxx(i,j,k) + gxx(i-xoffset,j-yoffset,k-zoffset)) gxyl = 0.5d0 * (gxy(i,j,k) + gxy(i-xoffset,j-yoffset,k-zoffset)) @@ -809,7 +817,7 @@ subroutine Conservative2PrimitivePolytype(CCTK_ARGUMENTS) CCTK_REAL :: uxx, uxy, uxz, uyy, uyz, uzz, det, psi4pt CCTK_INT :: type_bits, atmosphere - CCTK_INT :: type2_bits, excised + CCTK_INT :: type2_bits CCTK_REAL :: local_min_tracer ! character(len=400) :: warnline @@ -817,8 +825,7 @@ subroutine Conservative2PrimitivePolytype(CCTK_ARGUMENTS) call SpaceMask_GetTypeBits(type_bits, "Hydro_Atmosphere") call SpaceMask_GetStateBits(atmosphere, "Hydro_Atmosphere", "in_atmosphere") - type2_bits = -1 - excised = -1 + type2_bits = -1 nx = cctk_lsh(1) ny = cctk_lsh(2) @@ -839,9 +846,9 @@ subroutine Conservative2PrimitivePolytype(CCTK_ARGUMENTS) do j = 1, ny do i = 1, nx - !do not compute if in atmosphere - !or in an excised region - if (SpaceMask_CheckStateBitsF90(space_mask, i, j, k, type_bits, atmosphere)) cycle + !do not compute if in atmosphere or in an excised region + if (SpaceMask_CheckStateBitsF90(space_mask, i, j, k, type_bits, atmosphere) .or. & + (hydro_excision_mask(i,j,k) .ne. 0)) cycle if (conformal_state .eq. 0) then call SpatialDeterminant(gxx(i,j,k),gxy(i,j,k),gxz(i,j,k),& @@ -929,7 +936,6 @@ subroutine Con2Prim_ptPolytype(handle, dens, sx, sy, sz, tau, rho, & #undef _EOS_BASE_INC_ #endif #include "EOS_Base.inc" - !!$ Undensitize the variables @@ -954,7 +960,6 @@ subroutine Con2Prim_ptPolytype(handle, dens, sx, sy, sz, tau, rho, & !$OMP CRITICAL write(warnline,'(a70,8g15.6)') 'NaN produced in sqrt(): (dens, det, s2, udens, enthalpy, x, y, z)', dens, det, s2, udens, enthalpy, x, y, z call CCTK_WARN(Whisky_NaN_verbose, warnline) - call CCTK_WARN(0, warnline) !$OMP END CRITICAL endif !!$ END: Check for NaN value (1st check) @@ -1018,6 +1023,12 @@ subroutine Con2Prim_ptPolytype(handle, dens, sx, sy, sz, tau, rho, & rhoold = rhonew rhonew = rhoold - f/df + if (rhonew .lt. 0.d0) then + whisky_C2P_failed = 1 + !$OMP CRITICAL + call CCTK_WARN(Whisky_NaN_verbose, 'rhonew went below 0') + !$OMP END CRITICAL + end if !!$ Recalculate primitive variables and function @@ -1145,12 +1156,11 @@ subroutine Con2PrimBoundsPolytype(CCTK_ARGUMENTS) gxxr,gxyr,gxzr,gyyr,gyzr,gzzr,avg_detr,psi4r,psi4l CCTK_INT :: type_bits, atmosphere - CCTK_INT :: type2_bits, excised + CCTK_INT :: type2_bits call SpaceMask_GetTypeBits(type_bits, "Hydro_Atmosphere") call SpaceMask_GetStateBits(atmosphere, "Hydro_Atmosphere", "in_atmosphere") - type2_bits = -1 - excised = -1 + type2_bits = -1 nx = cctk_lsh(1) ny = cctk_lsh(2) @@ -1160,8 +1170,9 @@ subroutine Con2PrimBoundsPolytype(CCTK_ARGUMENTS) do j = whisky_stencil, ny - whisky_stencil + 1 do i = whisky_stencil, nx - whisky_stencil + 1 - !do not compute if in atmosphere - if (SpaceMask_CheckStateBitsF90(space_mask, i, j, k, type_bits, atmosphere)) cycle + !do not compute if in atmosphere or in an excised region + if (SpaceMask_CheckStateBitsF90(space_mask, i, j, k, type_bits, atmosphere) .or. & + (hydro_excision_mask(i,j,k) .ne. 0)) cycle gxxl = 0.5d0 * (gxx(i,j,k) + gxx(i-xoffset,j-yoffset,k-zoffset)) gxyl = 0.5d0 * (gxy(i,j,k) + gxy(i-xoffset,j-yoffset,k-zoffset)) @@ -1447,7 +1458,7 @@ subroutine Conservative2PrimitiveGeneral(CCTK_ARGUMENTS) logical, dimension(:,:,:), allocatable :: point_loop_condition CCTK_INT :: type_bits, atmosphere - CCTK_INT :: type2_bits, excised + CCTK_INT :: type2_bits CCTK_REAL :: local_min_tracer @@ -1463,7 +1474,6 @@ subroutine Conservative2PrimitiveGeneral(CCTK_ARGUMENTS) call SpaceMask_GetTypeBits(type_bits, "Hydro_Atmosphere") call SpaceMask_GetStateBits(atmosphere, "Hydro_Atmosphere", "in_atmosphere") type2_bits = -1 - excised = -1 nx = cctk_lsh(1) ny = cctk_lsh(2) @@ -1482,11 +1492,6 @@ subroutine Conservative2PrimitiveGeneral(CCTK_ARGUMENTS) !!$ Set up rho and epsilon -!!! write(6,"(A15,i3,1P10E18.9)") "Xatmo 34 4 4: ", atmosphere_mask(34,4,4), rho(34,4,4), & -!!! whisky_rho_min, & -!!! dens(34,4,4), sqrt(Whisky_Det(34,4,4))*whisky_rho_min*(1.0d0+whisky_atmo_tolerance) -!!! call flush(6) - do k = 1, nz do j = 1, ny do i = 1, nx @@ -1518,11 +1523,7 @@ subroutine Conservative2PrimitiveGeneral(CCTK_ARGUMENTS) SpaceMask_SetStateBitsF90(space_mask, i, j, k, type_bits, atmosphere) atmosphere_mask(i,j,k) = 1 -! call Whisky_DumpPointInAtmos(rho(i,j,k), & -! vel(i,j,k,1), vel(i,j,k,2), vel(i,j,k,3), & -! eps(i,j,k), press(i,j,k), w_lorentz(i,j,k), & -! dens(i,j,k), scon(i,j,k,1), scon(i,j,k,2), scon(i,j,k,3), tau(i,j,k), & -! Whisky_det(i,j,k), whisky_rho_min) + end if udens = dens(i,j,k) / sqrt(Whisky_det(i,j,k)) @@ -1652,13 +1653,13 @@ subroutine Conservative2PrimitiveGeneral(CCTK_ARGUMENTS) do j = 1, ny do i = 1, nx - !do not compute if in atmosphere - !or in an excised region - - if (SpaceMask_CheckStateBitsF90(space_mask, i, j, k, type_bits, atmosphere)) & + !do not compute if in atmosphere or in an excised region + if (SpaceMask_CheckStateBitsF90(space_mask, i, j, k, type_bits, atmosphere) .or. & + (hydro_excision_mask(i,j,k) .ne. 0)) & then press_new(i,j,k) = press(i,j,k) press_old(i,j,k) = press(i,j,k) + cycle endif @@ -1756,8 +1757,9 @@ subroutine Conservative2PrimitiveGeneral(CCTK_ARGUMENTS) do j = 1, ny do i = 1, nx - !do not compute if in atmosphere - if (SpaceMask_CheckStateBitsF90(space_mask, i, j, k, type_bits, atmosphere)) cycle + !do not compute if in atmosphere or in an excised region + if (SpaceMask_CheckStateBitsF90(space_mask, i, j, k, type_bits, atmosphere) .or. & + (hydro_excision_mask(i,j,k) .ne. 0)) cycle f = press_new(i,j,k) - press(i,j,k) @@ -1809,8 +1811,9 @@ subroutine Conservative2PrimitiveGeneral(CCTK_ARGUMENTS) do j = 1, ny do i = 1, nx - !do not compute if in atmosphere - if (SpaceMask_CheckStateBitsF90(space_mask, i, j, k, type_bits, atmosphere)) cycle + !do not compute if in atmosphere or in an excised region + if (SpaceMask_CheckStateBitsF90(space_mask, i, j, k, type_bits, atmosphere) .or. & + (hydro_excision_mask(i,j,k) .ne. 0)) cycle !!$ Note that in the following we enforce eps > 0 !!$ This makes the later check redundant @@ -1903,11 +1906,6 @@ subroutine Conservative2PrimitiveGeneral(CCTK_ARGUMENTS) deallocate(f) -!!! write(6,"(A15,i3,1P10E18.9)") "Yatmo 34 4 4: ", atmosphere_mask(34,4,4), rho(34,4,4), & -!!! whisky_rho_min, & -!!! dens(34,4,4), sqrt(Whisky_Det(34,4,4))*whisky_rho_min*(1.0d0+whisky_atmo_tolerance) -!!! call flush(6) - end subroutine Conservative2PrimitiveGeneral @@ -1949,14 +1947,13 @@ subroutine Con2PrimPolytypeGeneral(CCTK_ARGUMENTS) CCTK_INT :: type_bits CCTK_INT :: atmosphere - CCTK_INT :: type2_bits, excised + CCTK_INT :: type2_bits CCTK_REAL :: local_min_tracer call SpaceMask_GetTypeBits(type_bits, "Hydro_Atmosphere") call SpaceMask_GetStateBits(atmosphere, "Hydro_Atmosphere", "in_atmosphere") - type2_bits = -1 - excised = -1 + type2_bits = -1 nx = cctk_lsh(1) ny = cctk_lsh(2) @@ -2044,8 +2041,9 @@ subroutine Con2PrimPolytypeGeneral(CCTK_ARGUMENTS) do j = 1, ny do i = 1, nx - !do not compute if in atmosphere - if (SpaceMask_CheckStateBitsF90(space_mask, i, j, k, type_bits, atmosphere)) cycle + !do not compute if in atmosphere or in an excised region + if (SpaceMask_CheckStateBitsF90(space_mask, i, j, k, type_bits, atmosphere) .or. & + (hydro_excision_mask(i,j,k) .ne. 0)) cycle enthalpy = 1.d0 + eps(i,j,k) + press(i,j,k) / press_new(i,j,k) @@ -2096,8 +2094,9 @@ subroutine Con2PrimPolytypeGeneral(CCTK_ARGUMENTS) do j = 1, ny do i = 1, nx - !do not compute if in atmosphere - if (SpaceMask_CheckStateBitsF90(space_mask, i, j, k, type_bits, atmosphere)) cycle + !do not compute if in atmosphere or in an excised region + if (SpaceMask_CheckStateBitsF90(space_mask, i, j, k, type_bits, atmosphere) .or. & + (hydro_excision_mask(i,j,k) .ne. 0)) cycle udens = dens(i,j,k) / sqrt(Whisky_det(i,j,k)) diff --git a/src/Whisky_ENOReconstruct.F90 b/src/Whisky_ENOReconstruct.F90 index 093e524..42f0661 100644 --- a/src/Whisky_ENOReconstruct.F90 +++ b/src/Whisky_ENOReconstruct.F90 @@ -153,7 +153,7 @@ end subroutine Whisky_ENOShutdown (iand(mask((i)),(type_bits)).eq.(state_bits)) subroutine Whisky_ENOReconstruct1d(order, nx, v, vminus, vplus, trivial_rp, & - space_mask, excision_descriptors) + hydro_excision_mask) USE Whisky_ENOScalars implicit none @@ -164,17 +164,13 @@ subroutine Whisky_ENOReconstruct1d(order, nx, v, vminus, vplus, trivial_rp, & CCTK_REAL, dimension(nx) :: v, vplus, vminus CCTK_REAL, dimension(order, 1-order:nx+order) :: vdiff - CCTK_INT, dimension(nx) :: space_mask + CCTK_INT, dimension(nx) :: hydro_excision_mask logical, dimension(nx) :: trivial_rp logical, dimension(nx) :: excise - CCTK_INT :: excision_bits, excision_mask - CCTK_INT, dimension(3) :: excision_descriptors + logical :: normal_eno CCTK_REAL :: large = 1.d10 - excision_bits=excision_descriptors(1) - excision_mask=excision_descriptors(2) - vminus = 0.d0 vplus = 0.d0 vdiff = 0.d0 @@ -183,39 +179,81 @@ subroutine Whisky_ENOReconstruct1d(order, nx, v, vminus, vplus, trivial_rp, & excise = .false. trivial_rp = .false. - +!!$ Initialize excision do i = 1, nx + if (hydro_excision_mask(i) .ne. 0) then + vdiff(1, i) = large * (-1.d0*(1.d0+mod(i,10)))**i + trivial_rp(i) = .true. + excise(i) = .true. + if (i > 1) then + trivial_rp(i-1) = .true. + end if + end if + end do + do i = 1, nx +!!$ Handle excision + normal_eno = .true. + if (i < nx) then + if (excise(i+1)) then + vminus(i) = v(i) + vplus(i) = v(i) + normal_eno = .false. + end if + end if + if (i > 1) then + if (excise(i-1)) then + vminus(i) = v(i) + vplus(i) = v(i) + normal_eno = .false. + end if + end if + if (normal_eno) then !!$ Calculate the undivided differences - do k = 2, order do j = max(1, i - order), min(nx, i + order) vdiff(k, j) = vdiff(k - 1, j + 1) - vdiff(k - 1, j) end do end do - + !!$ Ensure the stencil stays within the grid - vdiff(:, 1 - order : 0) = 1.d10 vdiff(:, nx + 1 : nx + order) = 1.d10 - + !!$ Find the stencil - r = 0 - do j = 2, order if ( abs(vdiff(j, i-r-1)) < abs(vdiff(j, i-r)) ) r = r + 1 end do !!$ Calculate the reconstruction - do j = 0, order - 1 vminus(i) = vminus(i) + eno_coeffs(r-1, j) * vdiff(1, i-r+j) vplus(i) = vplus(i) + eno_coeffs(r, j) * vdiff(1, i-r+j) end do + end if + end do - + do i = 1, nx + if (excise(i)) then + if (i > 1) then + if (.not. excise(i-1)) then + vminus(i) = vplus(i-1) + end if + end if + vplus(i) = vminus(i) + end if + end do + do i = nx, 1, -1 + if (excise(i)) then + if (i < nx) then + if (.not. excise(i+1)) then + vplus(i) = vminus(i+1) + end if + end if + vminus(i) = vplus(i) + end if end do end subroutine Whisky_ENOReconstruct1d diff --git a/src/Whisky_PPM.F90 b/src/Whisky_PPM.F90 index e249831..e56b910 100644 --- a/src/Whisky_PPM.F90 +++ b/src/Whisky_PPM.F90 @@ -53,7 +53,7 @@ end subroutine PPM_TVD subroutine SimplePPM_1d(handle,poly,& nx,dx,rho,velx,vely,velz,eps,press,rhominus,& velxminus,velyminus,velzminus,epsminus,rhoplus,velxplus,velyplus,& - velzplus,epsplus,trivial_rp,space_mask, excision_descriptors,& + velzplus,epsplus,trivial_rp, hydro_excision_mask,& gxx, gxy, gxz, gyy, gyz, gzz, psi4, beta, alp, w_lorentz, & dir, ni, nj, nrx, nry, nrz, ev_l, ev_r, xw) @@ -89,9 +89,7 @@ subroutine SimplePPM_1d(handle,poly,& logical, dimension(nx) :: trivial_rp - CCTK_INT, dimension(nx) :: space_mask - CCTK_INT :: excision_bits, excision_mask - CCTK_INT, dimension(3) :: excision_descriptors + CCTK_INT, dimension(nx) :: hydro_excision_mask CCTK_REAL, dimension(nx) :: gxx, gxy, gxz, gyy, gyz, gzz, & psi4, beta, alp, w_lorentz @@ -432,7 +430,99 @@ do i = whisky_stencil, nx - whisky_stencil + 1 end if !!$ excision - return + if (hydro_excision) then + do i = 1, nx + if (hydro_excision_mask(i) .ne. 0) then + if (i .gt. 1) then + trivial_rp(i-1)=.true. + end if + trivial_rp(i)=.true. + else + !!$ Do not optimize cond away by combining the 'if's. Fortran does not + !!$ have to follow the order of sub-expressions given here and might + !!$ access outside the array range + cond = i .gt. 1 + if (cond) cond = hydro_excision_mask(i-1) .ne. 0 + if (cond) then + rhominus(i)=rho(i) + rhoplus(i)=rho(i) + velxminus(i)=velx(i) + velxplus(i)=velx(i) + velyminus(i)=vely(i) + velyplus(i)=vely(i) + velzminus(i)=velz(i) + velzplus(i)=velz(i) + rhominus(i-1)=rho(i) + rhoplus(i-1)=rho(i) + velxminus(i-1)=velx(i) + velxplus(i-1)=velx(i) + velyminus(i-1)=vely(i) + velyplus(i-1)=vely(i) + velzminus(i-1)=velz(i) + velzplus(i-1)=velz(i) + if (poly .eq. 0) then + epsminus(i)=eps(i) + epsplus(i)=eps(i) + epsminus(i-1)=eps(i) + epsplus(i-1)=eps(i) + end if + else + cond = (i.gt.2) .and. (i.lt.nx) + if (cond) cond = (ppm_mppm .eq. 0) .and.& + (hydro_excision_mask(i-2) .ne. 0) + if (cond) then + call PPM_TVD(rho(i-1), rho(i), rho(i+1), rhominus(i), rhoplus(i)) + call PPM_TVD(velx(i-1), velx(i), velx(i+1), velxminus(i), velxplus(i)) + call PPM_TVD(vely(i-1), vely(i), vely(i+1), velyminus(i), velyplus(i)) + call PPM_TVD(velz(i-1), velz(i), velz(i+1), velzminus(i), velzplus(i)) + if (poly .eq. 0) then + call PPM_TVD(eps(i-1), eps(i), eps(i+1), epsminus(i), epsplus(i)) + end if + end if + end if + cond = i.lt.nx + if (cond) cond = hydro_excision_mask(i+1) .ne. 0; + if (cond) then + rhominus(i)=rho(i) + rhoplus(i)=rho(i) + velxminus(i)=velx(i) + velxplus(i)=velx(i) + velyminus(i)=vely(i) + velyplus(i)=vely(i) + velzminus(i)=velz(i) + velzplus(i)=velz(i) + rhominus(i+1)=rho(i) + rhoplus(i+1)=rho(i) + velxminus(i+1)=velx(i) + velxplus(i+1)=velx(i) + velyminus(i+1)=vely(i) + velyplus(i+1)=vely(i) + velzminus(i+1)=velz(i) + velzplus(i+1)=velz(i) + if (poly .eq. 0) then + epsminus(i)=eps(i) + epsplus(i)=eps(i) + epsminus(i+1)=eps(i) + epsplus(i+1)=eps(i) + endif + else + cond = (i.lt.nx-1) .and. (i.gt.1) + if (cond) cond = (ppm_mppm .eq. 0) .and.& + (hydro_excision_mask(i+2) .ne. 0) + if (cond) then + call PPM_TVD(rho(i-1), rho(i), rho(i+1), rhominus(i), rhoplus(i)) + call PPM_TVD(velx(i-1), velx(i), velx(i+1), velxminus(i), velxplus(i)) + call PPM_TVD(vely(i-1), vely(i), vely(i+1), velyminus(i), velyplus(i)) + call PPM_TVD(velz(i-1), velz(i), velz(i+1), velzminus(i), velzplus(i)) + if (poly .eq. 0) then + call PPM_TVD(eps(i-1), eps(i), eps(i+1), epsminus(i), epsplus(i)) + end if + end if + end if + end if + end do + end if +return end subroutine SimplePPM_1d diff --git a/src/Whisky_Prim2Con.F90 b/src/Whisky_Prim2Con.F90 index 01f3d9f..2c823f5 100644 --- a/src/Whisky_Prim2Con.F90 +++ b/src/Whisky_Prim2Con.F90 @@ -329,10 +329,13 @@ subroutine prim2conpolytype(handle, gxx, gxy, gxz, gyy, gyz, & implicit none + DECLARE_CCTK_PARAMETERS + CCTK_REAL :: gxx, gxy, gxz, gyy, gyz, gzz, det CCTK_REAL :: ddens, dsx, dsy, dsz, dtau, drho, dvelx, dvely, dvelz,& deps, dpress, w_tmp, w, vlowx, vlowy, vlowz, sqrtdet CCTK_INT :: handle + character(len=256) NaN_WarnLine #ifdef _EOS_BASE_INC_ #undef _EOS_BASE_INC_ @@ -348,6 +351,10 @@ subroutine prim2conpolytype(handle, gxx, gxy, gxz, gyy, gyz, & ! and hard to trace wrong values below. There is no good value to ! choose in this case, but something small is probably the best of ! all bad choices. + !$OMP CRITICAL + write(NaN_WarnLine,'(a100,1g15.6)') 'Infinite Lorentz factor reset. rho: ', drho + call CCTK_WARN(Whisky_NaN_verbose, NaN_WarnLine) + !$OMP END CRITICAL w = 1.d-20 else w = 1.d0 / sqrt(1.d0 - w_tmp) diff --git a/src/Whisky_Reconstruct.F90 b/src/Whisky_Reconstruct.F90 index 89bedb7..7b2365d 100644 --- a/src/Whisky_Reconstruct.F90 +++ b/src/Whisky_Reconstruct.F90 @@ -55,7 +55,6 @@ subroutine Reconstruction(CCTK_ARGUMENTS) &type_bitsy, trivialy, not_trivialy, & &type_bitsz, trivialz, not_trivialz - CCTK_INT, dimension(3) :: excision_descriptors CCTK_REAL, dimension(:,:,:),allocatable :: & &psi4, lbetax, lbetay, lbetaz !!$ CCTK_REAL, dimension(cctk_lsh(1),cctk_lsh(2),cctk_lsh(3)) :: & @@ -78,10 +77,6 @@ subroutine Reconstruction(CCTK_ARGUMENTS) call CCTK_WARN(0, "Allocation problems with lbeta") end if - excision_descriptors(1)=-1 - excision_descriptors(2)=-1 - excision_descriptors(3)=-1 - call SpaceMask_GetTypeBits(type_bitsx, "Hydro_RiemannProblemX") call SpaceMask_GetStateBits(trivialx, "Hydro_RiemannProblemX", & &"trivial") @@ -147,42 +142,32 @@ subroutine Reconstruction(CCTK_ARGUMENTS) call tvdreconstruct(nx, ny, nz, xoffset, yoffset, zoffset, & tracer(:,:,:,itracer), tracerplus(:,:,:,itracer), & tracerminus(:,:,:,itracer), & - trivial_rp, space_mask, excision_descriptors) + trivial_rp, hydro_excision_mask) enddo end if if (CCTK_EQUALS(recon_vars,"primitive")) then call tvdreconstruct(nx, ny, nz, xoffset, yoffset, zoffset, & - rho, rhoplus, rhominus, trivial_rp, space_mask, & - excision_descriptors) + rho, rhoplus, rhominus, trivial_rp, hydro_excision_mask) call tvdreconstruct(nx, ny, nz, xoffset, yoffset, zoffset, & - vel(:,:,:,1), velxplus, velxminus, trivial_rp, space_mask, & - excision_descriptors) + vel(:,:,:,1), velxplus, velxminus, trivial_rp, hydro_excision_mask) call tvdreconstruct(nx, ny, nz, xoffset, yoffset, zoffset, & - vel(:,:,:,2), velyplus, velyminus, trivial_rp, space_mask, & - excision_descriptors) + vel(:,:,:,2), velyplus, velyminus, trivial_rp, hydro_excision_mask) call tvdreconstruct(nx, ny, nz, xoffset, yoffset, zoffset, & - vel(:,:,:,3), velzplus, velzminus, trivial_rp, space_mask, & - excision_descriptors) + vel(:,:,:,3), velzplus, velzminus, trivial_rp, hydro_excision_mask) call tvdreconstruct(nx, ny, nz, xoffset, yoffset, zoffset, & - eps, epsplus, epsminus, trivial_rp, space_mask, & - excision_descriptors) + eps, epsplus, epsminus, trivial_rp, hydro_excision_mask) else if (CCTK_EQUALS(recon_vars,"conservative")) then call tvdreconstruct(nx, ny, nz, xoffset, yoffset, zoffset, & - dens, densplus, densminus, trivial_rp, space_mask, & - excision_descriptors) + dens, densplus, densminus, trivial_rp, hydro_excision_mask) call tvdreconstruct(nx, ny, nz, xoffset, yoffset, zoffset, & - scon(:,:,:,1), sxplus, sxminus, trivial_rp, space_mask, & - excision_descriptors) + scon(:,:,:,1), sxplus, sxminus, trivial_rp, hydro_excision_mask) call tvdreconstruct(nx, ny, nz, xoffset, yoffset, zoffset, & - scon(:,:,:,2), syplus, syminus, trivial_rp, space_mask, & - excision_descriptors) + scon(:,:,:,2), syplus, syminus, trivial_rp, hydro_excision_mask) call tvdreconstruct(nx, ny, nz, xoffset, yoffset, zoffset, & - scon(:,:,:,3), szplus, szminus, trivial_rp, space_mask, & - excision_descriptors) + scon(:,:,:,3), szplus, szminus, trivial_rp, hydro_excision_mask) call tvdreconstruct(nx, ny, nz, xoffset, yoffset, zoffset, & - tau, tauplus, tauminus, trivial_rp, space_mask, & - excision_descriptors) + tau, tauplus, tauminus, trivial_rp, hydro_excision_mask) else call CCTK_WARN(0, "Variable type to reconstruct not recognized.") end if @@ -224,8 +209,7 @@ subroutine Reconstruction(CCTK_ARGUMENTS) press(:,j,k),rhominus(:,j,k),velxminus(:,j,k),velyminus(:,j,k),& velzminus(:,j,k),epsminus(:,j,k),rhoplus(:,j,k),& velxplus(:,j,k),velyplus(:,j,k),velzplus(:,j,k),epsplus(:,j,k),& - trivial_rp(:,j,k),space_mask(:,j,k), & - excision_descriptors,& + trivial_rp(:,j,k), hydro_excision_mask(:,j,k),& gxx(:,j,k), gxy(:,j,k), gxz(:,j,k), gyy(:,j,k), gyz(:,j,k), & gzz(:,j,k), psi4(:,j,k), lbetax(:,j,k), alp(:,j,k),& w_lorentz(:,j,k), & @@ -264,8 +248,7 @@ subroutine Reconstruction(CCTK_ARGUMENTS) press(j,:,k),rhominus(j,:,k),velyminus(j,:,k),velzminus(j,:,k),& velxminus(j,:,k),epsminus(j,:,k),rhoplus(j,:,k),& velyplus(j,:,k),velzplus(j,:,k),velxplus(j,:,k),epsplus(j,:,k),& - trivial_rp(j,:,k),space_mask(j,:,k), & - excision_descriptors,& + trivial_rp(j,:,k), hydro_excision_mask(j,:,k),& gyy(j,:,k), gyz(j,:,k), gxy(j,:,k), gzz(j,:,k), gxz(j,:,k), & gxx(j,:,k), psi4(j,:,k), lbetay(j,:,k), alp(j,:,k),& w_lorentz(j,:,k), & @@ -304,8 +287,7 @@ subroutine Reconstruction(CCTK_ARGUMENTS) press(j,k,:),rhominus(j,k,:),velzminus(j,k,:),velxminus(j,k,:),& velyminus(j,k,:),epsminus(j,k,:),rhoplus(j,k,:),& velzplus(j,k,:),velxplus(j,k,:),velyplus(j,k,:),epsplus(j,k,:),& - trivial_rp(j,k,:),space_mask(j,k,:), & - excision_descriptors,& + trivial_rp(j,k,:), hydro_excision_mask(j,k,:),& gzz(j,k,:), gxz(j,k,:), gyz(j,k,:), gxx(j,k,:), gxy(j,k,:), & gyy(j,k,:), psi4(j,k,:), lbetaz(j,k,:), alp(j,k,:),& w_lorentz(j,k,:), & @@ -346,7 +328,7 @@ subroutine Reconstruction(CCTK_ARGUMENTS) call tvdreconstruct(nx, ny, nz, xoffset, yoffset, zoffset, & tracer(:,:,:,itracer), tracerplus(:,:,:,itracer), & tracerminus(:,:,:,itracer), trivial_rp, & - space_mask, excision_descriptors) + hydro_excision_mask) enddo end if @@ -357,45 +339,35 @@ subroutine Reconstruction(CCTK_ARGUMENTS) if (CCTK_EQUALS(recon_vars,"primitive")) then call Whisky_ENOReconstruct1d(eno_order,cctk_lsh(1),& rho(:,j,k),rhominus(:,j,k),rhoplus(:,j,k),& - trivial_rp(:,j,k), space_mask(:,j,k), & - excision_descriptors) + trivial_rp(:,j,k), hydro_excision_mask(:,j,k)) call Whisky_ENOReconstruct1d(eno_order,cctk_lsh(1),& velx(:,j,k),velxminus(:,j,k),velxplus(:,j,k),& - trivial_rp(:,j,k), space_mask(:,j,k), & - excision_descriptors) + trivial_rp(:,j,k), hydro_excision_mask(:,j,k)) call Whisky_ENOReconstruct1d(eno_order,cctk_lsh(1),& vely(:,j,k),velyminus(:,j,k),velyplus(:,j,k),& - trivial_rp(:,j,k), space_mask(:,j,k), & - excision_descriptors) + trivial_rp(:,j,k), hydro_excision_mask(:,j,k)) call Whisky_ENOReconstruct1d(eno_order,cctk_lsh(1),& velz(:,j,k),velzminus(:,j,k),velzplus(:,j,k),& - trivial_rp(:,j,k), space_mask(:,j,k), & - excision_descriptors) + trivial_rp(:,j,k), hydro_excision_mask(:,j,k)) call Whisky_ENOReconstruct1d(eno_order,cctk_lsh(1),& eps(:,j,k),epsminus(:,j,k),epsplus(:,j,k),& - trivial_rp(:,j,k), space_mask(:,j,k), & - excision_descriptors) + trivial_rp(:,j,k), hydro_excision_mask(:,j,k)) else if (CCTK_EQUALS(recon_vars,"conservative")) then call Whisky_ENOReconstruct1d(eno_order,cctk_lsh(1),& dens(:,j,k),densminus(:,j,k),densplus(:,j,k),& - trivial_rp(:,j,k), space_mask(:,j,k), & - excision_descriptors) + trivial_rp(:,j,k), hydro_excision_mask(:,j,k)) call Whisky_ENOReconstruct1d(eno_order,cctk_lsh(1),& sx(:,j,k),sxminus(:,j,k),sxplus(:,j,k),& - trivial_rp(:,j,k), space_mask(:,j,k), & - excision_descriptors) + trivial_rp(:,j,k), hydro_excision_mask(:,j,k)) call Whisky_ENOReconstruct1d(eno_order,cctk_lsh(1),& sy(:,j,k),syminus(:,j,k),syplus(:,j,k),& - trivial_rp(:,j,k), space_mask(:,j,k), & - excision_descriptors) + trivial_rp(:,j,k), hydro_excision_mask(:,j,k)) call Whisky_ENOReconstruct1d(eno_order,cctk_lsh(1),& sz(:,j,k),szminus(:,j,k),szplus(:,j,k),& - trivial_rp(:,j,k), space_mask(:,j,k), & - excision_descriptors) + trivial_rp(:,j,k), hydro_excision_mask(:,j,k)) call Whisky_ENOReconstruct1d(eno_order,cctk_lsh(1),& tau(:,j,k),tauminus(:,j,k),tauplus(:,j,k),& - trivial_rp(:,j,k), space_mask(:,j,k), & - excision_descriptors) + trivial_rp(:,j,k), hydro_excision_mask(:,j,k)) else !$OMP CRITICAL call CCTK_WARN(0, "Variable type to reconstruct not recognized.") @@ -418,45 +390,35 @@ subroutine Reconstruction(CCTK_ARGUMENTS) if (CCTK_EQUALS(recon_vars,"primitive")) then call Whisky_ENOReconstruct1d(eno_order,cctk_lsh(2),& rho(j,:,k),rhominus(j,:,k),rhoplus(j,:,k),& - trivial_rp(j,:,k), space_mask(j,:,k), & - excision_descriptors) + trivial_rp(j,:,k), hydro_excision_mask(j,:,k)) call Whisky_ENOReconstruct1d(eno_order,cctk_lsh(2),& velx(j,:,k),velxminus(j,:,k),velxplus(j,:,k),& - trivial_rp(j,:,k), space_mask(j,:,k), & - excision_descriptors) + trivial_rp(j,:,k), hydro_excision_mask(j,:,k)) call Whisky_ENOReconstruct1d(eno_order,cctk_lsh(2),& vely(j,:,k),velyminus(j,:,k),velyplus(j,:,k),& - trivial_rp(j,:,k), space_mask(j,:,k), & - excision_descriptors) + trivial_rp(j,:,k), hydro_excision_mask(j,:,k)) call Whisky_ENOReconstruct1d(eno_order,cctk_lsh(2),& velz(j,:,k),velzminus(j,:,k),velzplus(j,:,k),& - trivial_rp(j,:,k), space_mask(j,:,k), & - excision_descriptors) + trivial_rp(j,:,k), hydro_excision_mask(j,:,k)) call Whisky_ENOReconstruct1d(eno_order,cctk_lsh(2),& eps(j,:,k),epsminus(j,:,k),epsplus(j,:,k),& - trivial_rp(j,:,k), space_mask(j,:,k), & - excision_descriptors) + trivial_rp(j,:,k), hydro_excision_mask(j,:,k)) else if (CCTK_EQUALS(recon_vars,"conservative")) then call Whisky_ENOReconstruct1d(eno_order,cctk_lsh(2),& dens(j,:,k),densminus(j,:,k),densplus(j,:,k),& - trivial_rp(j,:,k), space_mask(j,:,k), & - excision_descriptors) + trivial_rp(j,:,k), hydro_excision_mask(j,:,k)) call Whisky_ENOReconstruct1d(eno_order,cctk_lsh(2),& sx(j,:,k),sxminus(j,:,k),sxplus(j,:,k),& - trivial_rp(j,:,k), space_mask(j,:,k), & - excision_descriptors) + trivial_rp(j,:,k), hydro_excision_mask(j,:,k)) call Whisky_ENOReconstruct1d(eno_order,cctk_lsh(2),& sy(j,:,k),syminus(j,:,k),syplus(j,:,k),& - trivial_rp(j,:,k), space_mask(j,:,k), & - excision_descriptors) + trivial_rp(j,:,k), hydro_excision_mask(j,:,k)) call Whisky_ENOReconstruct1d(eno_order,cctk_lsh(2),& sz(j,:,k),szminus(j,:,k),szplus(j,:,k),& - trivial_rp(j,:,k), space_mask(j,:,k), & - excision_descriptors) + trivial_rp(j,:,k), hydro_excision_mask(j,:,k)) call Whisky_ENOReconstruct1d(eno_order,cctk_lsh(2),& tau(j,:,k),tauminus(j,:,k),tauplus(j,:,k),& - trivial_rp(j,:,k), space_mask(j,:,k), & - excision_descriptors) + trivial_rp(j,:,k), hydro_excision_mask(j,:,k)) else !$OMP CRITICAL call CCTK_WARN(0, "Variable type to reconstruct not recognized.") @@ -479,45 +441,35 @@ subroutine Reconstruction(CCTK_ARGUMENTS) if (CCTK_EQUALS(recon_vars,"primitive")) then call Whisky_ENOReconstruct1d(eno_order,cctk_lsh(3),& rho(j,k,:),rhominus(j,k,:),rhoplus(j,k,:),& - trivial_rp(j,k,:), space_mask(j,k,:), & - excision_descriptors) + trivial_rp(j,k,:), hydro_excision_mask(j,k,:)) call Whisky_ENOReconstruct1d(eno_order,cctk_lsh(3),& velx(j,k,:),velxminus(j,k,:),velxplus(j,k,:),& - trivial_rp(j,k,:), space_mask(j,k,:), & - excision_descriptors) + trivial_rp(j,k,:), hydro_excision_mask(j,k,:)) call Whisky_ENOReconstruct1d(eno_order,cctk_lsh(3),& vely(j,k,:),velyminus(j,k,:),velyplus(j,k,:),& - trivial_rp(j,k,:), space_mask(j,k,:), & - excision_descriptors) + trivial_rp(j,k,:), hydro_excision_mask(j,k,:)) call Whisky_ENOReconstruct1d(eno_order,cctk_lsh(3),& velz(j,k,:),velzminus(j,k,:),velzplus(j,k,:),& - trivial_rp(j,k,:), space_mask(j,k,:), & - excision_descriptors) + trivial_rp(j,k,:), hydro_excision_mask(j,k,:)) call Whisky_ENOReconstruct1d(eno_order,cctk_lsh(3),& eps(j,k,:),epsminus(j,k,:),epsplus(j,k,:),& - trivial_rp(j,k,:), space_mask(j,k,:), & - excision_descriptors) + trivial_rp(j,k,:), hydro_excision_mask(j,k,:)) else if (CCTK_EQUALS(recon_vars,"conservative")) then call Whisky_ENOReconstruct1d(eno_order,cctk_lsh(3),& dens(j,k,:),densminus(j,k,:),densplus(j,k,:),& - trivial_rp(j,k,:), space_mask(j,k,:), & - excision_descriptors) + trivial_rp(j,k,:), hydro_excision_mask(j,k,:)) call Whisky_ENOReconstruct1d(eno_order,cctk_lsh(3),& sx(j,k,:),sxminus(j,k,:),sxplus(j,k,:),& - trivial_rp(j,k,:), space_mask(j,k,:), & - excision_descriptors) + trivial_rp(j,k,:), hydro_excision_mask(j,k,:)) call Whisky_ENOReconstruct1d(eno_order,cctk_lsh(3),& sy(j,k,:),syminus(j,k,:),syplus(j,k,:),& - trivial_rp(j,k,:), space_mask(j,k,:), & - excision_descriptors) + trivial_rp(j,k,:), hydro_excision_mask(j,k,:)) call Whisky_ENOReconstruct1d(eno_order,cctk_lsh(3),& sz(j,k,:),szminus(j,k,:),szplus(j,k,:),& - trivial_rp(j,k,:), space_mask(j,k,:), & - excision_descriptors) + trivial_rp(j,k,:), hydro_excision_mask(j,k,:)) call Whisky_ENOReconstruct1d(eno_order,cctk_lsh(3),& tau(j,k,:),tauminus(j,k,:),tauplus(j,k,:),& - trivial_rp(j,k,:), space_mask(j,k,:), & - excision_descriptors) + trivial_rp(j,k,:), hydro_excision_mask(j,k,:)) else !$OMP CRITICAL call CCTK_WARN(0, "Variable type to reconstruct not recognized.") diff --git a/src/Whisky_ReconstructPoly.F90 b/src/Whisky_ReconstructPoly.F90 index dc065ce..c7a5dfa 100644 --- a/src/Whisky_ReconstructPoly.F90 +++ b/src/Whisky_ReconstructPoly.F90 @@ -55,7 +55,6 @@ subroutine ReconstructionPolytype(CCTK_ARGUMENTS) &type_bitsy, trivialy, not_trivialy, & &type_bitsz, trivialz, not_trivialz - CCTK_INT, dimension(3) :: excision_descriptors CCTK_REAL, dimension(:,:,:),allocatable :: & &psi4, lbetax, lbetay, lbetaz !!$ CCTK_REAL, dimension(cctk_lsh(1),cctk_lsh(2),cctk_lsh(3)) :: & @@ -78,10 +77,6 @@ subroutine ReconstructionPolytype(CCTK_ARGUMENTS) call CCTK_WARN(0, "Allocation problems with lbeta") end if - excision_descriptors(1)=-1 - excision_descriptors(2)=-1 - excision_descriptors(3)=-1 - call SpaceMask_GetTypeBits(type_bitsx, "Hydro_RiemannProblemX") call SpaceMask_GetStateBits(trivialx, "Hydro_RiemannProblemX", & &"trivial") @@ -147,38 +142,30 @@ subroutine ReconstructionPolytype(CCTK_ARGUMENTS) call tvdreconstruct(nx, ny, nz, xoffset, yoffset, zoffset, & tracer(:,:,:,itracer), tracerplus(:,:,:,itracer), & tracerminus(:,:,:,itracer), & - trivial_rp, space_mask, excision_descriptors) + trivial_rp, hydro_excision_mask) enddo end if if (CCTK_EQUALS(recon_vars,"primitive")) then call tvdreconstruct(nx, ny, nz, xoffset, yoffset, zoffset, & - rho, rhoplus, rhominus, trivial_rp, space_mask, & - excision_descriptors) + rho, rhoplus, rhominus, trivial_rp, hydro_excision_mask) call tvdreconstruct(nx, ny, nz, xoffset, yoffset, zoffset, & - vel(:,:,:,1), velxplus, velxminus, trivial_rp, space_mask, & - excision_descriptors) + vel(:,:,:,1), velxplus, velxminus, trivial_rp, hydro_excision_mask) call tvdreconstruct(nx, ny, nz, xoffset, yoffset, zoffset, & - vel(:,:,:,2), velyplus, velyminus, trivial_rp, space_mask, & - excision_descriptors) + vel(:,:,:,2), velyplus, velyminus, trivial_rp, hydro_excision_mask) call tvdreconstruct(nx, ny, nz, xoffset, yoffset, zoffset, & - vel(:,:,:,3), velzplus, velzminus, trivial_rp, space_mask, & - excision_descriptors) + vel(:,:,:,3), velzplus, velzminus, trivial_rp, hydro_excision_mask) else if (CCTK_EQUALS(recon_vars,"conservative")) then call tvdreconstruct(nx, ny, nz, xoffset, yoffset, zoffset, & - dens, densplus, densminus, trivial_rp, space_mask, & - excision_descriptors) + dens, densplus, densminus, trivial_rp, hydro_excision_mask) call tvdreconstruct(nx, ny, nz, xoffset, yoffset, zoffset, & - scon(:,:,:,1), sxplus, sxminus, trivial_rp, space_mask, & - excision_descriptors) + scon(:,:,:,1), sxplus, sxminus, trivial_rp, hydro_excision_mask) call tvdreconstruct(nx, ny, nz, xoffset, yoffset, zoffset, & - scon(:,:,:,2), syplus, syminus, trivial_rp, space_mask, & - excision_descriptors) + scon(:,:,:,2), syplus, syminus, trivial_rp, hydro_excision_mask) call tvdreconstruct(nx, ny, nz, xoffset, yoffset, zoffset, & - scon(:,:,:,3), szplus, szminus, trivial_rp, space_mask, & - excision_descriptors) + scon(:,:,:,3), szplus, szminus, trivial_rp, hydro_excision_mask) else call CCTK_WARN(0, "Variable type to reconstruct not recognized.") @@ -221,8 +208,7 @@ subroutine ReconstructionPolytype(CCTK_ARGUMENTS) press(:,j,k),rhominus(:,j,k),velxminus(:,j,k),velyminus(:,j,k),& velzminus(:,j,k),epsminus(:,j,k),rhoplus(:,j,k),& velxplus(:,j,k),velyplus(:,j,k),velzplus(:,j,k),epsplus(:,j,k),& - trivial_rp(:,j,k),space_mask(:,j,k), & - excision_descriptors,& + trivial_rp(:,j,k), hydro_excision_mask(:,j,k),& gxx(:,j,k), gxy(:,j,k), gxz(:,j,k), gyy(:,j,k), gyz(:,j,k), & gzz(:,j,k), psi4(:,j,k), lbetax(:,j,k), alp(:,j,k),& w_lorentz(:,j,k), & @@ -261,8 +247,7 @@ subroutine ReconstructionPolytype(CCTK_ARGUMENTS) press(j,:,k),rhominus(j,:,k),velyminus(j,:,k),velzminus(j,:,k),& velxminus(j,:,k),epsminus(j,:,k),rhoplus(j,:,k),& velyplus(j,:,k),velzplus(j,:,k),velxplus(j,:,k),epsplus(j,:,k),& - trivial_rp(j,:,k),space_mask(j,:,k), & - excision_descriptors,& + trivial_rp(j,:,k), hydro_excision_mask(j,:,k),& gyy(j,:,k), gyz(j,:,k), gxy(j,:,k), gzz(j,:,k), gxz(j,:,k), & gxx(j,:,k), psi4(j,:,k), lbetay(j,:,k), alp(j,:,k),& w_lorentz(j,:,k), & @@ -301,8 +286,7 @@ subroutine ReconstructionPolytype(CCTK_ARGUMENTS) press(j,k,:),rhominus(j,k,:),velzminus(j,k,:),velxminus(j,k,:),& velyminus(j,k,:),epsminus(j,k,:),rhoplus(j,k,:),& velzplus(j,k,:),velxplus(j,k,:),velyplus(j,k,:),epsplus(j,k,:),& - trivial_rp(j,k,:),space_mask(j,k,:), & - excision_descriptors,& + trivial_rp(j,k,:), hydro_excision_mask(j,k,:),& gzz(j,k,:), gxz(j,k,:), gyz(j,k,:), gxx(j,k,:), gxy(j,k,:), & gyy(j,k,:), psi4(j,k,:), lbetaz(j,k,:), alp(j,k,:),& w_lorentz(j,k,:), & @@ -343,7 +327,7 @@ subroutine ReconstructionPolytype(CCTK_ARGUMENTS) call tvdreconstruct(nx, ny, nz, xoffset, yoffset, zoffset, & tracer(:,:,:,itracer), tracerplus(:,:,:,itracer), & tracerminus(:,:,:,itracer), trivial_rp, & - space_mask, excision_descriptors) + hydro_excision_mask) enddo end if @@ -354,37 +338,29 @@ subroutine ReconstructionPolytype(CCTK_ARGUMENTS) if (CCTK_EQUALS(recon_vars,"primitive")) then call Whisky_ENOReconstruct1d(eno_order,cctk_lsh(1),& rho(:,j,k),rhominus(:,j,k),rhoplus(:,j,k),& - trivial_rp(:,j,k), space_mask(:,j,k), & - excision_descriptors) + trivial_rp(:,j,k), hydro_excision_mask(:,j,k)) call Whisky_ENOReconstruct1d(eno_order,cctk_lsh(1),& velx(:,j,k),velxminus(:,j,k),velxplus(:,j,k),& - trivial_rp(:,j,k), space_mask(:,j,k), & - excision_descriptors) + trivial_rp(:,j,k), hydro_excision_mask(:,j,k)) call Whisky_ENOReconstruct1d(eno_order,cctk_lsh(1),& vely(:,j,k),velyminus(:,j,k),velyplus(:,j,k),& - trivial_rp(:,j,k), space_mask(:,j,k), & - excision_descriptors) + trivial_rp(:,j,k), hydro_excision_mask(:,j,k)) call Whisky_ENOReconstruct1d(eno_order,cctk_lsh(1),& velz(:,j,k),velzminus(:,j,k),velzplus(:,j,k),& - trivial_rp(:,j,k), space_mask(:,j,k), & - excision_descriptors) + trivial_rp(:,j,k), hydro_excision_mask(:,j,k)) else if (CCTK_EQUALS(recon_vars,"conservative")) then call Whisky_ENOReconstruct1d(eno_order,cctk_lsh(1),& dens(:,j,k),densminus(:,j,k),densplus(:,j,k),& - trivial_rp(:,j,k), space_mask(:,j,k), & - excision_descriptors) + trivial_rp(:,j,k), hydro_excision_mask(:,j,k)) call Whisky_ENOReconstruct1d(eno_order,cctk_lsh(1),& sx(:,j,k),sxminus(:,j,k),sxplus(:,j,k),& - trivial_rp(:,j,k), space_mask(:,j,k), & - excision_descriptors) + trivial_rp(:,j,k), hydro_excision_mask(:,j,k)) call Whisky_ENOReconstruct1d(eno_order,cctk_lsh(1),& sy(:,j,k),syminus(:,j,k),syplus(:,j,k),& - trivial_rp(:,j,k), space_mask(:,j,k), & - excision_descriptors) + trivial_rp(:,j,k), hydro_excision_mask(:,j,k)) call Whisky_ENOReconstruct1d(eno_order,cctk_lsh(1),& sz(:,j,k),szminus(:,j,k),szplus(:,j,k),& - trivial_rp(:,j,k), space_mask(:,j,k), & - excision_descriptors) + trivial_rp(:,j,k), hydro_excision_mask(:,j,k)) else !$OMP CRITICAL call CCTK_WARN(0, "Variable type to reconstruct not recognized.") @@ -407,37 +383,29 @@ subroutine ReconstructionPolytype(CCTK_ARGUMENTS) if (CCTK_EQUALS(recon_vars,"primitive")) then call Whisky_ENOReconstruct1d(eno_order,cctk_lsh(2),& rho(j,:,k),rhominus(j,:,k),rhoplus(j,:,k),& - trivial_rp(j,:,k), space_mask(j,:,k), & - excision_descriptors) + trivial_rp(j,:,k), hydro_excision_mask(j,:,k)) call Whisky_ENOReconstruct1d(eno_order,cctk_lsh(2),& velx(j,:,k),velxminus(j,:,k),velxplus(j,:,k),& - trivial_rp(j,:,k), space_mask(j,:,k), & - excision_descriptors) + trivial_rp(j,:,k), hydro_excision_mask(j,:,k)) call Whisky_ENOReconstruct1d(eno_order,cctk_lsh(2),& vely(j,:,k),velyminus(j,:,k),velyplus(j,:,k),& - trivial_rp(j,:,k), space_mask(j,:,k), & - excision_descriptors) + trivial_rp(j,:,k), hydro_excision_mask(j,:,k)) call Whisky_ENOReconstruct1d(eno_order,cctk_lsh(2),& velz(j,:,k),velzminus(j,:,k),velzplus(j,:,k),& - trivial_rp(j,:,k), space_mask(j,:,k), & - excision_descriptors) + trivial_rp(j,:,k), hydro_excision_mask(j,:,k)) else if (CCTK_EQUALS(recon_vars,"conservative")) then call Whisky_ENOReconstruct1d(eno_order,cctk_lsh(2),& dens(j,:,k),densminus(j,:,k),densplus(j,:,k),& - trivial_rp(j,:,k), space_mask(j,:,k), & - excision_descriptors) + trivial_rp(j,:,k), hydro_excision_mask(j,:,k)) call Whisky_ENOReconstruct1d(eno_order,cctk_lsh(2),& sx(j,:,k),sxminus(j,:,k),sxplus(j,:,k),& - trivial_rp(j,:,k), space_mask(j,:,k), & - excision_descriptors) + trivial_rp(j,:,k), hydro_excision_mask(j,:,k)) call Whisky_ENOReconstruct1d(eno_order,cctk_lsh(2),& sy(j,:,k),syminus(j,:,k),syplus(j,:,k),& - trivial_rp(j,:,k), space_mask(j,:,k), & - excision_descriptors) + trivial_rp(j,:,k), hydro_excision_mask(j,:,k)) call Whisky_ENOReconstruct1d(eno_order,cctk_lsh(2),& sz(j,:,k),szminus(j,:,k),szplus(j,:,k),& - trivial_rp(j,:,k), space_mask(j,:,k), & - excision_descriptors) + trivial_rp(j,:,k), hydro_excision_mask(j,:,k)) else !$OMP CRITICAL call CCTK_WARN(0, "Variable type to reconstruct not recognized.") @@ -460,37 +428,29 @@ subroutine ReconstructionPolytype(CCTK_ARGUMENTS) if (CCTK_EQUALS(recon_vars,"primitive")) then call Whisky_ENOReconstruct1d(eno_order,cctk_lsh(3),& rho(j,k,:),rhominus(j,k,:),rhoplus(j,k,:),& - trivial_rp(j,k,:), space_mask(j,k,:), & - excision_descriptors) + trivial_rp(j,k,:), hydro_excision_mask(j,k,:)) call Whisky_ENOReconstruct1d(eno_order,cctk_lsh(3),& velx(j,k,:),velxminus(j,k,:),velxplus(j,k,:),& - trivial_rp(j,k,:), space_mask(j,k,:), & - excision_descriptors) + trivial_rp(j,k,:), hydro_excision_mask(j,k,:)) call Whisky_ENOReconstruct1d(eno_order,cctk_lsh(3),& vely(j,k,:),velyminus(j,k,:),velyplus(j,k,:),& - trivial_rp(j,k,:), space_mask(j,k,:), & - excision_descriptors) + trivial_rp(j,k,:), hydro_excision_mask(j,k,:)) call Whisky_ENOReconstruct1d(eno_order,cctk_lsh(3),& velz(j,k,:),velzminus(j,k,:),velzplus(j,k,:),& - trivial_rp(j,k,:), space_mask(j,k,:), & - excision_descriptors) + trivial_rp(j,k,:), hydro_excision_mask(j,k,:)) else if (CCTK_EQUALS(recon_vars,"conservative")) then call Whisky_ENOReconstruct1d(eno_order,cctk_lsh(3),& dens(j,k,:),densminus(j,k,:),densplus(j,k,:),& - trivial_rp(j,k,:), space_mask(j,k,:), & - excision_descriptors) + trivial_rp(j,k,:), hydro_excision_mask(j,k,:)) call Whisky_ENOReconstruct1d(eno_order,cctk_lsh(3),& sx(j,k,:),sxminus(j,k,:),sxplus(j,k,:),& - trivial_rp(j,k,:), space_mask(j,k,:), & - excision_descriptors) + trivial_rp(j,k,:), hydro_excision_mask(j,k,:)) call Whisky_ENOReconstruct1d(eno_order,cctk_lsh(3),& sy(j,k,:),syminus(j,k,:),syplus(j,k,:),& - trivial_rp(j,k,:), space_mask(j,k,:), & - excision_descriptors) + trivial_rp(j,k,:), hydro_excision_mask(j,k,:)) call Whisky_ENOReconstruct1d(eno_order,cctk_lsh(3),& sz(j,k,:),szminus(j,k,:),szplus(j,k,:),& - trivial_rp(j,k,:), space_mask(j,k,:), & - excision_descriptors) + trivial_rp(j,k,:), hydro_excision_mask(j,k,:)) else !$OMP CRITICAL call CCTK_WARN(0, "Variable type to reconstruct not recognized.") diff --git a/src/Whisky_TVDReconstruct.F90 b/src/Whisky_TVDReconstruct.F90 index d6f4b95..9788a05 100644 --- a/src/Whisky_TVDReconstruct.F90 +++ b/src/Whisky_TVDReconstruct.F90 @@ -14,8 +14,6 @@ #include "SpaceMask.h" -/* #define WHISKY_UNSTABLE_TVD */ - /*@@ @routine tvdreconstruct @date Sat Jan 26 02:12:12 2002 @@ -32,7 +30,7 @@ @@*/ subroutine tvdreconstruct(nx, ny, nz, xoffset, yoffset, zoffset, & - orig, bextp, bextm, trivial_rp, space_mask, excision_descriptors) + orig, bextp, bextm, trivial_rp, hydro_excision_mask) USE Whisky_Scalars @@ -45,9 +43,7 @@ subroutine tvdreconstruct(nx, ny, nz, xoffset, yoffset, zoffset, & CCTK_REAL, dimension(nx, ny, nz) :: orig, bextp, bextm CCTK_REAL :: dupw, dloc, delta, ratio, hdelta logical, dimension(nx,ny,nz) :: trivial_rp - CCTK_INT, dimension(nx,ny,nz) :: space_mask - CCTK_INT :: excision_bits, excision_mask - CCTK_INT, dimension(3) :: excision_descriptors + CCTK_INT, dimension(nx,ny,nz) :: hydro_excision_mask bextp = 0.d0 @@ -60,41 +56,45 @@ subroutine tvdreconstruct(nx, ny, nz, xoffset, yoffset, zoffset, & do k = whisky_stencil, nz-whisky_stencil+1 do j = whisky_stencil, ny-whisky_stencil+1 do i = whisky_stencil, nx-whisky_stencil+1 - dupw = orig(i, j, k) - orig(i-xoffset, j-yoffset, k-zoffset) - dloc = orig(i+xoffset, j+yoffset, k+zoffset) - orig(i, j, k) -!!$ For minmod set everything here. -!!$ Otherwise call the limiter function. - - if (MINMOD) then -!!$ if (dupw*dloc < 0.d0) then -!!$ delta=0.d0 -!!$ else if (abs(dupw)<abs(dloc)) then -!!$ delta=dupw -!!$ else -!!$ delta=dloc -!!$ end if - - delta = minmod_func(dupw,dloc) - -!!$ This is an alternative equivalent implementation -!!$ of vanLeeer MC slopelimiter - else if (MC2) then - if (dupw*dloc < 0.d0) then - delta=0.d0 - else - delta=sign(min(2.d0*abs(dupw),2.d0*abs(dloc),& - 0.5d0*(abs(dupw)+abs(dloc))),dupw+dloc) + if (hydro_excision_mask(i,j,k) .ne. 0) then + trivial_rp(i-xoffset, j-yoffset, k-zoffset) = .true. + trivial_rp(i, j, k) = .true. + bextm(i, j, k) = orig(i, j, k) + bextp(i, j, k) = orig(i, j, k) + if (hydro_excision_mask(i+xoffset,j+yoffset,k+zoffset) .eq. 0) then + bextm(i, j, k) = orig(i+xoffset, j+yoffset, k+zoffset) + bextp(i, j, k) = orig(i+xoffset, j+yoffset, k+zoffset) end if + else if ((hydro_excision_mask(i-xoffset,j-yoffset,k-zoffset) .ne. 0) .or. & + (hydro_excision_mask(i+xoffset,j+yoffset,k+zoffset) .ne. 0)) then + bextm(i, j, k) = orig(i, j, k) + bextp(i, j, k) = orig(i, j, k) else - delta = 0.5d0*(dupw + dloc) - if (abs(dupw) < myfloor ) dupw = myfloor*sign(1.d0, dupw) - if (abs(dloc) < myfloor ) dloc = myfloor*sign(1.d0, dloc) - ratio = dupw / dloc - call slopelimiter(ratio, delta) + dupw = orig(i, j, k) - orig(i-xoffset, j-yoffset, k-zoffset) + dloc = orig(i+xoffset, j+yoffset, k+zoffset) - orig(i, j, k) + + if (MINMOD) then + delta = minmod_func(dupw,dloc) + else if (MC2) then +!!$ This is an alternative equivalent implementation +!!$ of vanLeeer MC slopelimiter + if (dupw*dloc < 0.d0) then + delta=0.d0 + else + delta=sign(min(2.d0*abs(dupw),2.d0*abs(dloc),& + 0.5d0*(abs(dupw)+abs(dloc))),dupw+dloc) + end if + else + delta = 0.5d0*(dupw + dloc) + if (abs(dupw) < myfloor ) dupw = myfloor*sign(1.d0, dupw) + if (abs(dloc) < myfloor ) dloc = myfloor*sign(1.d0, dloc) + ratio = dupw / dloc + call slopelimiter(ratio, delta) + end if + hdelta = 0.5d0 * delta + bextm(i, j, k) = orig(i, j, k) - hdelta + bextp(i, j, k) = orig(i, j, k) + hdelta end if - hdelta = 0.5d0 * delta - bextm(i, j, k) = orig(i, j, k) - hdelta - bextp(i, j, k) = orig(i, j, k) + hdelta end do end do end do diff --git a/src/Whisky_UpdateMask.F90 b/src/Whisky_UpdateMask.F90 index 49e28c3..0b590ff 100644 --- a/src/Whisky_UpdateMask.F90 +++ b/src/Whisky_UpdateMask.F90 @@ -34,10 +34,12 @@ subroutine WhiskyUpdateAtmosphereMask(CCTK_ARGUMENTS) frac = CCTK_DELTA_TIME + !$OMP PARALLEL DO PRIVATE(j,i) do k = 1, cctk_lsh(3) do j = 1, cctk_lsh(2) do i = 1, cctk_lsh(1) - if ( (SpaceMask_CheckStateBitsF90(space_mask, i, j, k, \ + if ( (hydro_excision_mask(i,j,k) .ne. 0) .or. & + (SpaceMask_CheckStateBitsF90(space_mask, i, j, k, \ type_bits, atmosphere)) .or. & (tau(i,j,k) + frac * taurhs(i,j,k) .le. 0.d0) .or. & (dens(i,j,k) + frac * densrhs(i,j,k) .le. 0.d0) ) then @@ -51,6 +53,7 @@ subroutine WhiskyUpdateAtmosphereMask(CCTK_ARGUMENTS) end do end do end do + !$OMP END PARALLEL DO end subroutine WhiskyUpdateAtmosphereMask |