diff options
| -rw-r--r-- | param.ccl | 17 | ||||
| -rw-r--r-- | schedule.ccl | 24 | ||||
| -rw-r--r-- | src/GRHydro_Reconstruct.F90 | 5 | ||||
| -rw-r--r-- | src/GRHydro_WENOReconstruct.F90 | 276 | ||||
| -rw-r--r-- | src/GRHydro_WENOReconstruct_drv.F90 | 471 | ||||
| -rw-r--r-- | src/GRHydro_WENOScalars.F90 | 23 | ||||
| -rw-r--r-- | src/make.code.defn | 6 |
7 files changed, 814 insertions, 8 deletions
@@ -133,9 +133,10 @@ CCTK_INT GRHydro_MaxNumSandRVars "The maximum number of save and restore variabl keyword recon_method "Which reconstruction method to use" { - "tvd" :: "Slope limited TVD" - "ppm" :: "PPM reconstruction" - "eno" :: "ENO reconstruction" + "tvd" :: "Slope limited TVD" + "ppm" :: "PPM reconstruction" + "eno" :: "ENO reconstruction" + "weno" :: "WENO reconstruction" } "tvd" keyword method_type "Which type of method to use" @@ -256,6 +257,16 @@ int eno_order "The order of accuracy of the ENO reconstruction" 1:* :: "Anything above 1, but above 5 is pointless" } 2 +int WENO_order "The order of accuracy of the WENO reconstruction" +{ + 5 :: "Fifth-order" +} 5 + +real weno_eps "WENO epsilon parameter" +{ + 0:* :: "small and positive" +} 1e-6 + keyword riemann_solver "Which Riemann solver to use" { "Roe" :: "Standard Roe solver" diff --git a/schedule.ccl b/schedule.ccl index fa71107..fb017da 100644 --- a/schedule.ccl +++ b/schedule.ccl @@ -427,9 +427,9 @@ if(CCTK_Equals(Bvec_evolution_method,"GRHydro")) { } "Evaluate ADM constraints, and perform symmetry boundary conditions" -###################################################################### -### Setup and remove the coefficient arrays for ENO reconstruction ### -###################################################################### +########################################################################### +### Setup and remove the coefficient arrays for ENO/WENO reconstruction ### +########################################################################### if (CCTK_Equals(recon_method,"eno")) { @@ -448,6 +448,24 @@ if (CCTK_Equals(recon_method,"eno")) } +if (CCTK_Equals(recon_method,"weno")) +{ + + schedule GRHydro_WENOSetup AT CCTK_Basegrid + { + OPTIONS: global + LANG:Fortran + } "Coefficients for WENO reconstruction" + + schedule GRHydro_WENOShutdown AT CCTK_Terminate BEFORE Driver_Terminate + { + OPTIONS: global + LANG:Fortran + } "Deallocate WENO coefficients" + +} + + ###################################################################### ### This routine is a bit strangely situated, as it's really to do ### ### with initial data. There may be a better place for it later. ### diff --git a/src/GRHydro_Reconstruct.F90 b/src/GRHydro_Reconstruct.F90 index 8ac033a..6a2fc27 100644 --- a/src/GRHydro_Reconstruct.F90 +++ b/src/GRHydro_Reconstruct.F90 @@ -89,8 +89,13 @@ subroutine Reconstruction(CCTK_ARGUMENTS) ! this handles MHD and non-MHD call GRHydro_ENOReconstruct_drv(CCTK_PASS_FTOF) + else if (CCTK_EQUALS(recon_method,"weno")) then + ! this handles MHD and non-MHD + call GRHydro_WENOReconstruct_drv(CCTK_PASS_FTOF) + else call CCTK_WARN(0, "Reconstruction method not recognized!") + end if if (evolve_tracer .ne. 0) then diff --git a/src/GRHydro_WENOReconstruct.F90 b/src/GRHydro_WENOReconstruct.F90 new file mode 100644 index 0000000..e7273a0 --- /dev/null +++ b/src/GRHydro_WENOReconstruct.F90 @@ -0,0 +1,276 @@ + /*@@ + @file GRHydro_WENOReconstruct.F90 + @date Fri Jan 3 2013 + @author Ian Hawke, Christian Reisswig + @desc + Routines to set up the coefficient array and to perform one dimensional + ENO reconstruction of arbitrary order. + @enddesc + @@*/ + +#include "cctk.h" +#include "cctk_Parameters.h" +#include "cctk_Arguments.h" + + /*@@ + @routine GRHydro_WENOSetup + @date Fri Jan 3 2013 + @author Christian Reisswig + @desc + Sets up the coefficient array for WENO reconstruction. + Uses the notation of Shu, equation (2.21), in + ''High Order ENO and WENO Schemes for CFD'' + (see ThornGuide for full reference). + One exception: (Shu) r -> (Here) i: avoiding name clash. + @enddesc + @calls + @calledby + @history + + @endhistory + +@@*/ + +subroutine GRHydro_WENOSetup(CCTK_ARGUMENTS) + + USE GRHydro_WENOScalars + + implicit none + + DECLARE_CCTK_ARGUMENTS + DECLARE_CCTK_PARAMETERS + + CCTK_INT :: i, j, k, l, m, q, allocstat + CCTK_REAL :: denominator, numerator, numerator_product + + if(.not.coeffs_allocated) then + ! Right now we hardcode to 5th order + allocate(weno_coeffs(3, 5), STAT=allocstat) + if (allocstat .ne. 0) call CCTK_WARN(0, "Failed to allocate WENO coefficient arrays!") + allocate(beta_shu(3, 6), STAT=allocstat) + if (allocstat .ne. 0) call CCTK_WARN(0, "Failed to allocate smoothness indicator stencil coefficient arrays!") + coeffs_allocated = .true. + endif + + ! Set stencils + weno_coeffs(1,1) = 3.0d0/8.0d0 + weno_coeffs(1,2) = -5.0d0/4.0d0 + weno_coeffs(1,3) = 15.0d0/8.0d0 + weno_coeffs(1,4) = 0.0d0 + weno_coeffs(1,5) = 0.0d0 + + weno_coeffs(2,1) = 0.0d0 + weno_coeffs(2,2) = -1.0d0/8.0d0 + weno_coeffs(2,3) = 3.0d0/4.0d0 + weno_coeffs(2,4) = 3.0d0/8.0d0 + weno_coeffs(2,5) = 0.0d0 + + weno_coeffs(3,1) = 0.0d0 + weno_coeffs(3,2) = 0.0d0 + weno_coeffs(3,3) = 3.0d0/8.0d0 + weno_coeffs(3,4) = 3.0d0/4.0d0 + weno_coeffs(3,5) = -1.0d0/8.0d0 + + ! Shu smoothness indicator stencil coefficients + beta_shu(1,1) = 4.0d0/3.0d0 + beta_shu(1,2) = -19.0d0/3.0d0 + beta_shu(1,3) = 25.0d0/3.0d0 + beta_shu(1,4) = 11.0d0/3.0d0 + beta_shu(1,5) = -31.0d0/3.0d0 + beta_shu(1,6) = 10.0d0/3.0d0 + + beta_shu(2,1) = 4.0d0/3.0d0 + beta_shu(2,2) = -13.0d0/3.0d0 + beta_shu(2,3) = 13.0d0/3.0d0 + beta_shu(2,4) = 5.0d0/3.0d0 + beta_shu(2,5) = -13.0d0/3.0d0 + beta_shu(2,6) = 4.0d0/3.0d0 + + beta_shu(3,1) = 10.0d0/3.0d0 + beta_shu(3,2) = -31.0d0/3.0d0 + beta_shu(3,3) = 25.0d0/3.0d0 + beta_shu(3,4) = 11.0d0/3.0d0 + beta_shu(3,5) = -19.0d0/3.0d0 + beta_shu(3,6) = 4.0d0/3.0d0 + +end subroutine GRHydro_WENOSetup + + /*@@ + @routine GRHydro_ENOShutdown + @date Fri Jan 3 2013 + @author Christian Reisswig + @desc + Deallocates the coefficient arrays + @enddesc + @calls + @calledby + @history + + @endhistory + +@@*/ + +subroutine GRHydro_WENOShutdown(CCTK_ARGUMENTS) + + USE GRHydro_WENOScalars + + implicit none + + DECLARE_CCTK_ARGUMENTS + + CCTK_INT :: deallocstat + + if(coeffs_allocated) then + deallocate(weno_coeffs, STAT = deallocstat) + if (deallocstat .ne. 0) call CCTK_WARN(0, "Failed to deallocate WENO coefficients.") + deallocate(beta_shu, STAT = deallocstat) + if (deallocstat .ne. 0) call CCTK_WARN(0, "Failed to deallocate shu smoothness indicator coefficients.") + coeffs_allocated = .false. + endif + +end subroutine GRHydro_WENOShutdown + + /*@@ + @routine GRHydro_ENOReconstruct1d + @date Fri Jan 3 2013 + @author Christian Reisswig + @desc + Perform a one dimensional reconstruction of a given array using WENO. + @enddesc + @calls + @calledby + @history + + @endhistory + +@@*/ + +#define SpaceMask_CheckStateBitsF90_1D(mask,i,type_bits,state_bits) \ + (iand(mask((i)),(type_bits)).eq.(state_bits)) + +subroutine GRHydro_WENOReconstruct1d(order, nx, v, vminus, vplus, trivial_rp, & + hydro_excision_mask) + USE GRHydro_WENOScalars + + implicit none + + DECLARE_CCTK_PARAMETERS + + CCTK_INT :: order, nx, i, j, k, r + CCTK_REAL, dimension(nx) :: v, vplus, vminus + CCTK_REAL, dimension(order, 1-order:nx+order) :: vdiff + + CCTK_INT, dimension(nx) :: hydro_excision_mask + logical, dimension(nx) :: trivial_rp + logical, dimension(nx) :: excise + logical :: normal_weno + + CCTK_REAL :: large = 1.d10, gamma1, gamma2, gamma3, beta1, beta2, beta3, w1, w2, w3, wbar1, wbar2, wbar3 + + vminus = 0.d0 + vplus = 0.d0 + + excise = .false. + trivial_rp = .false. + + + +!!$ Initialize excision + do i = 1, nx + if (GRHydro_enable_internal_excision /= 0 .and. (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 = 3, nx-2 +!!$ Handle excision + normal_weno = .true. + if (i < nx) then + if (excise(i+1)) then + vminus(i) = v(i) + vplus(i) = v(i) + normal_weno = .false. + end if + end if + if (i > 1) then + if (excise(i-1)) then + vminus(i) = v(i) + vplus(i) = v(i) + normal_weno = .false. + end if + end if + + if (normal_weno) then + +!!$ Compute smoothness indicators + beta1 = beta_shu(1,1)*v(i-2)**2 & + + beta_shu(1,2)*v(i-2)*v(i-1) & + + beta_shu(1,3)*v(i-1)**2 & + + beta_shu(1,4)*v(i-2)*v(i) & + + beta_shu(1,5)*v(i-1)*v(i) & + + beta_shu(1,6)*v(i)**2 + + beta2 = beta_shu(2,1)*v(i-1)**2 & + + beta_shu(2,2)*v(i-1)*v(i) & + + beta_shu(2,3)*v(i)**2 & + + beta_shu(2,4)*v(i-1)*v(i+1) & + + beta_shu(2,5)*v(i)*v(i+1) & + + beta_shu(2,6)*v(i+1)**2 + + beta3 = beta_shu(3,1)*v(i)**2 & + + beta_shu(3,2)*v(i)*v(i+1) & + + beta_shu(3,3)*v(i+1)**2 & + + beta_shu(3,4)*v(i)*v(i+2) & + + beta_shu(3,5)*v(i+1)*v(i+2) & + + beta_shu(3,6)*v(i+2)**2 + + + wbar1 = 1.0d0/16.d0 / (weno_eps + beta1)**2 + wbar2 = 5.0d0/8.d0 / (weno_eps + beta2)**2 + wbar3 = 5.0d0/16.d0 / (weno_eps + beta3)**2 + + w1 = wbar1 / (wbar1 + wbar2 + wbar3) + w2 = wbar2 / (wbar1 + wbar2 + wbar3) + w3 = wbar3 / (wbar1 + wbar2 + wbar3) + +!!$ Calculate the reconstruction + do j = 1, 5 + vplus(i) = vplus(i) + w1 * weno_coeffs(1,j)*v(i-3+j) & + + w2 * weno_coeffs(2,j)*v(i-3+j) & + + w3 * weno_coeffs(3,j)*v(i-3+j) + vminus(i) = vminus(i) + w1 * weno_coeffs(3,6-j)*v(i-3+j) & + + w2 * weno_coeffs(2,6-j)*v(i-3+j) & + + w3 * weno_coeffs(1,6-j)*v(i-3+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 GRHydro_WENOReconstruct1d diff --git a/src/GRHydro_WENOReconstruct_drv.F90 b/src/GRHydro_WENOReconstruct_drv.F90 new file mode 100644 index 0000000..2c748c5 --- /dev/null +++ b/src/GRHydro_WENOReconstruct_drv.F90 @@ -0,0 +1,471 @@ + /*@@ + @file GRHydro_WENOReconstruct_drv.F90 + @date Fri Jan 3 2013 + @author Christian Reisswig, Bruno C. Mundim, Joshua Faber, Christian D. Ott + @desc + Driver routine to perform the WENO reconstruction. + @enddesc + @@*/ + +#include "cctk.h" +#include "cctk_Arguments.h" +#include "cctk_Parameters.h" +#include "cctk_Functions.h" + +#include "SpaceMask.h" + +#define velx(i,j,k) vup(i,j,k,1) +#define vely(i,j,k) vup(i,j,k,2) +#define velz(i,j,k) vup(i,j,k,3) +#define sx(i,j,k) scon(i,j,k,1) +#define sy(i,j,k) scon(i,j,k,2) +#define sz(i,j,k) scon(i,j,k,3) +#define Bvecx(i,j,k) Bprim(i,j,k,1) +#define Bvecy(i,j,k) Bprim(i,j,k,2) +#define Bvecz(i,j,k) Bprim(i,j,k,3) +#define Bconsx(i,j,k) Bcons(i,j,k,1) +#define Bconsy(i,j,k) Bcons(i,j,k,2) +#define Bconsz(i,j,k) Bcons(i,j,k,3) + + + /*@@ + @routine GRHydro_WENOReconstruct_drv + @date Fri Jan 3 2013 + @author Christian Reisswig, Luca Baiotti, Ian Hawke, Bruno C. Mundim, Joshua Faber, Christian D. Ott + @desc + A driver routine to do WENO reconstruction. + @enddesc + @calls + @calledby + @history + + @endhistory + +@@*/ + +subroutine GRHydro_WENOReconstruct_drv(CCTK_ARGUMENTS) + + implicit none + + ! save memory when MP is not used + ! TARGET as to be before DECLARE_CCTK_ARGUMENTS for gcc 4.1 + TARGET vel, lvel + TARGET Bvec, lBvec + + DECLARE_CCTK_ARGUMENTS + DECLARE_CCTK_PARAMETERS + DECLARE_CCTK_FUNCTIONS + + integer :: nx, ny, nz, i, j, k, itracer + + logical, dimension(:,:,:), allocatable :: trivial_rp + + CCTK_INT :: type_bitsx, trivialx, not_trivialx, & + &type_bitsy, trivialy, not_trivialy, & + &type_bitsz, trivialz, not_trivialz + + CCTK_REAL, dimension(:,:,:),allocatable :: & + &dum, dump, dumm + + CCTK_INT :: ierr + + ! save memory when MP is not used + CCTK_INT :: GRHydro_UseGeneralCoordinates + CCTK_REAL, DIMENSION(:,:,:,:), POINTER :: vup, Bprim + + if (GRHydro_UseGeneralCoordinates(cctkGH).ne.0) then + vup => lvel + Bprim => lBvec + else + vup => vel + Bprim => Bvec + end if + + allocate(trivial_rp(cctk_lsh(1),cctk_lsh(2),cctk_lsh(3)),STAT=ierr) + if (ierr .ne. 0) then + call CCTK_WARN(0, "Allocation problems with trivial_rp") + end if + +!!$ The dum variables are used as dummies if MHD on but divergence cleaning off + allocate(dum(cctk_lsh(1),cctk_lsh(2),cctk_lsh(3)),& + dump(cctk_lsh(1),cctk_lsh(2),cctk_lsh(3)),& + dumm(cctk_lsh(1),cctk_lsh(2),cctk_lsh(3)),STAT=ierr) + + if (ierr .ne. 0) then + call CCTK_WARN(0, "Allocation problems with dum dump dumm") + end if + + call SpaceMask_GetTypeBits(type_bitsx, "Hydro_RiemannProblemX") + call SpaceMask_GetStateBits(trivialx, "Hydro_RiemannProblemX", & + &"trivial") + call SpaceMask_GetStateBits(not_trivialx, "Hydro_RiemannProblemX", & + &"not_trivial") + call SpaceMask_GetTypeBits(type_bitsy, "Hydro_RiemannProblemY") + call SpaceMask_GetStateBits(trivialy, "Hydro_RiemannProblemY", & + &"trivial") + call SpaceMask_GetStateBits(not_trivialy, "Hydro_RiemannProblemY", & + &"not_trivial") + call SpaceMask_GetTypeBits(type_bitsz, "Hydro_RiemannProblemZ") + call SpaceMask_GetStateBits(trivialz, "Hydro_RiemannProblemZ", & + &"trivial") + call SpaceMask_GetStateBits(not_trivialz, "Hydro_RiemannProblemZ", & + &"not_trivial") + + nx = cctk_lsh(1) + ny = cctk_lsh(2) + nz = cctk_lsh(3) + +!!$ Initialize variables that store reconstructed quantities: + + !$OMP PARALLEL DO PRIVATE(i,j,k) + do k=1,cctk_lsh(3) + do j=1,cctk_lsh(2) + do i=1,cctk_lsh(1) + trivial_rp(i,j,k) = .false. + rhoplus(i,j,k) = 0.0d0 + rhominus(i,j,k)= 0.0d0 + epsplus(i,j,k) = 0.0d0 + epsminus(i,j,k) = 0.0d0 + velxplus(i,j,k) = 0.0d0 + velxminus(i,j,k) = 0.0d0 + velyplus(i,j,k) = 0.0d0 + velyminus(i,j,k) = 0.0d0 + velzplus(i,j,k) = 0.0d0 + velzminus(i,j,k) = 0.0d0 + + if(evolve_mhd.ne.0) then + Bvecxplus(i,j,k) = 0.0d0 + Bvecxminus(i,j,k) = 0.0d0 + Bvecyplus(i,j,k) = 0.0d0 + Bvecyminus(i,j,k) = 0.0d0 + Bveczplus(i,j,k) = 0.0d0 + Bveczminus(i,j,k) = 0.0d0 + if(clean_divergence.ne.0) then + psidcplus(i,j,k) = 0.0d0 + psidcminus(i,j,k) = 0.0d0 + endif + endif + + if (evolve_entropy .ne. 0) then + entropyplus(i,j,k) = 0.0d0 + entropyminus(i,j,k) = 0.0d0 + endif + + if (evolve_tracer .ne. 0) then + tracerplus(i,j,k,:) = 0.0d0 + tracerminus(i,j,k,:) = 0.0d0 + endif + + if (evolve_Y_e .ne. 0) then + Y_e_plus(i,j,k) = 0.0d0 + Y_e_minus(i,j,k) = 0.0d0 + endif + + enddo + enddo + enddo + !$OMP END PARALLEL DO + +!!$ WENO starts: + if (evolve_tracer .ne. 0) then + do itracer=1,number_of_tracers + call tvdreconstruct(nx, ny, nz, xoffset, yoffset, zoffset, & + tracer(:,:,:,itracer), tracerplus(:,:,:,itracer), & + tracerminus(:,:,:,itracer), trivial_rp, & + hydro_excision_mask) + enddo + end if + + if (evolve_Y_e .ne. 0) then + call tvdreconstruct(nx, ny, nz, xoffset, yoffset, zoffset, & + Y_e(:,:,:), Y_e_plus(:,:,:), & + Y_e_minus(:,:,:), & + trivial_rp, hydro_excision_mask) + endif + + if (flux_direction == 1) then + ! constraint transport needs to be able to average fluxes in the directions + ! other that flux_direction, which in turn need the primitives on interfaces + !$OMP PARALLEL DO PRIVATE(i, j, k) + do k = GRHydro_stencil, cctk_lsh(3) - GRHydro_stencil + 1 + transport_constraints + do j = GRHydro_stencil, cctk_lsh(2) - GRHydro_stencil + 1 + transport_constraints + if (CCTK_EQUALS(recon_vars,"primitive")) then + call GRHydro_WENOReconstruct1d(WENO_order,cctk_lsh(1),& + rho(:,j,k),rhominus(:,j,k),rhoplus(:,j,k),& + trivial_rp(:,j,k), hydro_excision_mask(:,j,k)) + call GRHydro_WENOReconstruct1d(WENO_order,cctk_lsh(1),& + velx(:,j,k),velxminus(:,j,k),velxplus(:,j,k),& + trivial_rp(:,j,k), hydro_excision_mask(:,j,k)) + call GRHydro_WENOReconstruct1d(WENO_order,cctk_lsh(1),& + vely(:,j,k),velyminus(:,j,k),velyplus(:,j,k),& + trivial_rp(:,j,k), hydro_excision_mask(:,j,k)) + call GRHydro_WENOReconstruct1d(WENO_order,cctk_lsh(1),& + velz(:,j,k),velzminus(:,j,k),velzplus(:,j,k),& + trivial_rp(:,j,k), hydro_excision_mask(:,j,k)) + call GRHydro_WENOReconstruct1d(WENO_order,cctk_lsh(1),& + eps(:,j,k),epsminus(:,j,k),epsplus(:,j,k),& + trivial_rp(:,j,k), hydro_excision_mask(:,j,k)) + if(evolve_mhd.ne.0) then + call GRHydro_WENOReconstruct1d(WENO_order,cctk_lsh(1),& + Bvecx(:,j,k),Bvecxminus(:,j,k),Bvecxplus(:,j,k),& + trivial_rp(:,j,k), hydro_excision_mask(:,j,k)) + call GRHydro_WENOReconstruct1d(WENO_order,cctk_lsh(1),& + Bvecy(:,j,k),Bvecyminus(:,j,k),Bvecyplus(:,j,k),& + trivial_rp(:,j,k), hydro_excision_mask(:,j,k)) + call GRHydro_WENOReconstruct1d(WENO_order,cctk_lsh(1),& + Bvecz(:,j,k),Bveczminus(:,j,k),Bveczplus(:,j,k),& + trivial_rp(:,j,k), hydro_excision_mask(:,j,k)) + endif + if (evolve_entropy .ne. 0) then + call GRHydro_WENOReconstruct1d(WENO_order,cctk_lsh(1),& + entropy(:,j,k),entropyminus(:,j,k),entropyplus(:,j,k),& + trivial_rp(:,j,k), hydro_excision_mask(:,j,k)) + endif + else if (CCTK_EQUALS(recon_vars,"conservative")) then + call GRHydro_WENOReconstruct1d(WENO_order,cctk_lsh(1),& + dens(:,j,k),densminus(:,j,k),densplus(:,j,k),& + trivial_rp(:,j,k), hydro_excision_mask(:,j,k)) + call GRHydro_WENOReconstruct1d(WENO_order,cctk_lsh(1),& + sx(:,j,k),sxminus(:,j,k),sxplus(:,j,k),& + trivial_rp(:,j,k), hydro_excision_mask(:,j,k)) + call GRHydro_WENOReconstruct1d(WENO_order,cctk_lsh(1),& + sy(:,j,k),syminus(:,j,k),syplus(:,j,k),& + trivial_rp(:,j,k), hydro_excision_mask(:,j,k)) + call GRHydro_WENOReconstruct1d(WENO_order,cctk_lsh(1),& + sz(:,j,k),szminus(:,j,k),szplus(:,j,k),& + trivial_rp(:,j,k), hydro_excision_mask(:,j,k)) + call GRHydro_WENOReconstruct1d(WENO_order,cctk_lsh(1),& + tau(:,j,k),tauminus(:,j,k),tauplus(:,j,k),& + trivial_rp(:,j,k), hydro_excision_mask(:,j,k)) + if(evolve_mhd.ne.0) then + call GRHydro_WENOReconstruct1d(WENO_order,cctk_lsh(1),& + Bconsx(:,j,k),Bconsxminus(:,j,k),Bconsxplus(:,j,k),& + trivial_rp(:,j,k), hydro_excision_mask(:,j,k)) + call GRHydro_WENOReconstruct1d(WENO_order,cctk_lsh(1),& + Bconsy(:,j,k),Bconsyminus(:,j,k),Bconsyplus(:,j,k),& + trivial_rp(:,j,k), hydro_excision_mask(:,j,k)) + call GRHydro_WENOReconstruct1d(WENO_order,cctk_lsh(1),& + Bconsz(:,j,k),Bconszminus(:,j,k),Bconszplus(:,j,k),& + trivial_rp(:,j,k), hydro_excision_mask(:,j,k)) + endif + if (evolve_entropy .ne. 0) then + call GRHydro_WENOReconstruct1d(WENO_order,cctk_lsh(1),& + entropycons(:,j,k),entropyconsminus(:,j,k),entropyconsplus(:,j,k),& + trivial_rp(:,j,k), hydro_excision_mask(:,j,k)) + endif + else + !$OMP CRITICAL + call CCTK_WARN(0, "Variable type to reconstruct not recognized.") + !$OMP END CRITICAL + end if + + if(evolve_mhd.ne.0.and.clean_divergence.ne.0) then + call GRHydro_WENOReconstruct1d(WENO_order,cctk_lsh(1),& + psidc(:,j,k),psidcminus(:,j,k),psidcplus(:,j,k),& + trivial_rp(:,j,k), hydro_excision_mask(:,j,k)) + endif + + do i = 1, cctk_lsh(1) + if (trivial_rp(i,j,k)) then + SpaceMask_SetStateBitsF90(space_mask, i, j, k, type_bitsx, trivialx) + else + SpaceMask_SetStateBitsF90(space_mask, i, j, k, type_bitsx, not_trivialx) + end if + end do + end do + end do + !$OMP END PARALLEL DO + else if (flux_direction == 2) then + ! constraint transport needs to be able to average fluxes in the directions + ! other that flux_direction, which in turn need the primitives on interfaces + !$OMP PARALLEL DO PRIVATE(i, j, k) + do k = GRHydro_stencil, cctk_lsh(3) - GRHydro_stencil + 1 + transport_constraints + do j = GRHydro_stencil, cctk_lsh(1) - GRHydro_stencil + 1 + transport_constraints + if (CCTK_EQUALS(recon_vars,"primitive")) then + call GRHydro_WENOReconstruct1d(WENO_order,cctk_lsh(2),& + rho(j,:,k),rhominus(j,:,k),rhoplus(j,:,k),& + trivial_rp(j,:,k), hydro_excision_mask(j,:,k)) + call GRHydro_WENOReconstruct1d(WENO_order,cctk_lsh(2),& + velx(j,:,k),velxminus(j,:,k),velxplus(j,:,k),& + trivial_rp(j,:,k), hydro_excision_mask(j,:,k)) + call GRHydro_WENOReconstruct1d(WENO_order,cctk_lsh(2),& + vely(j,:,k),velyminus(j,:,k),velyplus(j,:,k),& + trivial_rp(j,:,k), hydro_excision_mask(j,:,k)) + call GRHydro_WENOReconstruct1d(WENO_order,cctk_lsh(2),& + velz(j,:,k),velzminus(j,:,k),velzplus(j,:,k),& + trivial_rp(j,:,k), hydro_excision_mask(j,:,k)) + call GRHydro_WENOReconstruct1d(WENO_order,cctk_lsh(2),& + eps(j,:,k),epsminus(j,:,k),epsplus(j,:,k),& + trivial_rp(j,:,k), hydro_excision_mask(j,:,k)) + if(evolve_mhd.ne.0) then + call GRHydro_WENOReconstruct1d(WENO_order,cctk_lsh(2),& + Bvecx(j,:,k),Bvecxminus(j,:,k),Bvecxplus(j,:,k),& + trivial_rp(j,:,k), hydro_excision_mask(j,:,k)) + call GRHydro_WENOReconstruct1d(WENO_order,cctk_lsh(2),& + Bvecy(j,:,k),Bvecyminus(j,:,k),Bvecyplus(j,:,k),& + trivial_rp(j,:,k), hydro_excision_mask(j,:,k)) + call GRHydro_WENOReconstruct1d(WENO_order,cctk_lsh(2),& + Bvecz(j,:,k),Bveczminus(j,:,k),Bveczplus(j,:,k),& + trivial_rp(j,:,k), hydro_excision_mask(j,:,k)) + endif + if (evolve_entropy .ne. 0) then + call GRHydro_WENOReconstruct1d(WENO_order,cctk_lsh(2),& + entropy(j,:,k),entropyminus(j,:,k),entropyplus(j,:,k),& + trivial_rp(j,:,k), hydro_excision_mask(j,:,k)) + endif + else if (CCTK_EQUALS(recon_vars,"conservative")) then + call GRHydro_WENOReconstruct1d(WENO_order,cctk_lsh(2),& + dens(j,:,k),densminus(j,:,k),densplus(j,:,k),& + trivial_rp(j,:,k), hydro_excision_mask(j,:,k)) + call GRHydro_WENOReconstruct1d(WENO_order,cctk_lsh(2),& + sx(j,:,k),sxminus(j,:,k),sxplus(j,:,k),& + trivial_rp(j,:,k), hydro_excision_mask(j,:,k)) + call GRHydro_WENOReconstruct1d(WENO_order,cctk_lsh(2),& + sy(j,:,k),syminus(j,:,k),syplus(j,:,k),& + trivial_rp(j,:,k), hydro_excision_mask(j,:,k)) + call GRHydro_WENOReconstruct1d(WENO_order,cctk_lsh(2),& + sz(j,:,k),szminus(j,:,k),szplus(j,:,k),& + trivial_rp(j,:,k), hydro_excision_mask(j,:,k)) + call GRHydro_WENOReconstruct1d(WENO_order,cctk_lsh(2),& + tau(j,:,k),tauminus(j,:,k),tauplus(j,:,k),& + trivial_rp(j,:,k), hydro_excision_mask(j,:,k)) + if(evolve_mhd.ne.0) then + call GRHydro_WENOReconstruct1d(WENO_order,cctk_lsh(2),& + Bconsx(j,:,k),Bconsxminus(j,:,k),Bconsxplus(j,:,k),& + trivial_rp(j,:,k), hydro_excision_mask(j,:,k)) + call GRHydro_WENOReconstruct1d(WENO_order,cctk_lsh(2),& + Bconsy(j,:,k),Bconsyminus(j,:,k),Bconsyplus(j,:,k),& + trivial_rp(j,:,k), hydro_excision_mask(j,:,k)) + call GRHydro_WENOReconstruct1d(WENO_order,cctk_lsh(2),& + Bconsz(j,:,k),Bconszminus(j,:,k),Bconszplus(j,:,k),& + trivial_rp(j,:,k), hydro_excision_mask(j,:,k)) + endif + if (evolve_entropy .ne. 0) then + call GRHydro_WENOReconstruct1d(WENO_order,cctk_lsh(2),& + entropycons(j,:,k),entropyconsminus(j,:,k),entropyconsplus(j,:,k),& + trivial_rp(j,:,k), hydro_excision_mask(j,:,k)) + endif + else + !$OMP CRITICAL + call CCTK_WARN(0, "Variable type to reconstruct not recognized.") + !$OMP END CRITICAL + end if + + if(evolve_mhd.ne.0.and.clean_divergence.ne.0) then + call GRHydro_WENOReconstruct1d(WENO_order,cctk_lsh(2),& + psidc(j,:,k),psidcminus(j,:,k),psidcplus(j,:,k),& + trivial_rp(j,:,k), hydro_excision_mask(j,:,k)) + endif + + do i = 1, cctk_lsh(2) + if (trivial_rp(j,i,k)) then + SpaceMask_SetStateBitsF90(space_mask, j, i, k, type_bitsy, trivialy) + else + SpaceMask_SetStateBitsF90(space_mask, j, i, k, type_bitsy, not_trivialy) + end if + end do + end do + end do + !$OMP END PARALLEL DO + else if (flux_direction == 3) then + ! constraint transport needs to be able to average fluxes in the directions + ! other that flux_direction, which in turn need the primitives on interfaces + !$OMP PARALLEL DO PRIVATE(i, j, k) + do k = GRHydro_stencil, cctk_lsh(2) - GRHydro_stencil + 1 + transport_constraints + do j = GRHydro_stencil, cctk_lsh(1) - GRHydro_stencil + 1 + transport_constraints + if (CCTK_EQUALS(recon_vars,"primitive")) then + call GRHydro_WENOReconstruct1d(WENO_order,cctk_lsh(3),& + rho(j,k,:),rhominus(j,k,:),rhoplus(j,k,:),& + trivial_rp(j,k,:), hydro_excision_mask(j,k,:)) + call GRHydro_WENOReconstruct1d(WENO_order,cctk_lsh(3),& + velx(j,k,:),velxminus(j,k,:),velxplus(j,k,:),& + trivial_rp(j,k,:), hydro_excision_mask(j,k,:)) + call GRHydro_WENOReconstruct1d(WENO_order,cctk_lsh(3),& + vely(j,k,:),velyminus(j,k,:),velyplus(j,k,:),& + trivial_rp(j,k,:), hydro_excision_mask(j,k,:)) + call GRHydro_WENOReconstruct1d(WENO_order,cctk_lsh(3),& + velz(j,k,:),velzminus(j,k,:),velzplus(j,k,:),& + trivial_rp(j,k,:), hydro_excision_mask(j,k,:)) + call GRHydro_WENOReconstruct1d(WENO_order,cctk_lsh(3),& + eps(j,k,:),epsminus(j,k,:),epsplus(j,k,:),& + trivial_rp(j,k,:), hydro_excision_mask(j,k,:)) + if(evolve_mhd.ne.0) then + call GRHydro_WENOReconstruct1d(WENO_order,cctk_lsh(3),& + Bvecx(j,k,:),Bvecxminus(j,k,:),Bvecxplus(j,k,:),& + trivial_rp(j,k,:), hydro_excision_mask(j,k,:)) + call GRHydro_WENOReconstruct1d(WENO_order,cctk_lsh(3),& + Bvecy(j,k,:),Bvecyminus(j,k,:),Bvecyplus(j,k,:),& + trivial_rp(j,k,:), hydro_excision_mask(j,k,:)) + call GRHydro_WENOReconstruct1d(WENO_order,cctk_lsh(3),& + Bvecz(j,k,:),Bveczminus(j,k,:),Bveczplus(j,k,:),& + trivial_rp(j,k,:), hydro_excision_mask(j,k,:)) + endif + if (evolve_entropy .ne. 0) then + call GRHydro_WENOReconstruct1d(WENO_order,cctk_lsh(2),& + entropy(j,k,:),entropyminus(j,k,:),entropyplus(j,k,:),& + trivial_rp(j,k,:), hydro_excision_mask(j,k,:)) + endif + else if (CCTK_EQUALS(recon_vars,"conservative")) then + call GRHydro_WENOReconstruct1d(WENO_order,cctk_lsh(3),& + dens(j,k,:),densminus(j,k,:),densplus(j,k,:),& + trivial_rp(j,k,:), hydro_excision_mask(j,k,:)) + call GRHydro_WENOReconstruct1d(WENO_order,cctk_lsh(3),& + sx(j,k,:),sxminus(j,k,:),sxplus(j,k,:),& + trivial_rp(j,k,:), hydro_excision_mask(j,k,:)) + call GRHydro_WENOReconstruct1d(WENO_order,cctk_lsh(3),& + sy(j,k,:),syminus(j,k,:),syplus(j,k,:),& + trivial_rp(j,k,:), hydro_excision_mask(j,k,:)) + call GRHydro_WENOReconstruct1d(WENO_order,cctk_lsh(3),& + sz(j,k,:),szminus(j,k,:),szplus(j,k,:),& + trivial_rp(j,k,:), hydro_excision_mask(j,k,:)) + call GRHydro_WENOReconstruct1d(WENO_order,cctk_lsh(3),& + tau(j,k,:),tauminus(j,k,:),tauplus(j,k,:),& + trivial_rp(j,k,:), hydro_excision_mask(j,k,:)) + if(evolve_mhd.ne.0) then + call GRHydro_WENOReconstruct1d(WENO_order,cctk_lsh(3),& + Bconsx(j,k,:),Bconsxminus(j,k,:),Bconsxplus(j,k,:),& + trivial_rp(j,k,:), hydro_excision_mask(j,k,:)) + call GRHydro_WENOReconstruct1d(WENO_order,cctk_lsh(3),& + Bconsy(j,k,:),Bconsyminus(j,k,:),Bconsyplus(j,k,:),& + trivial_rp(j,k,:), hydro_excision_mask(j,k,:)) + call GRHydro_WENOReconstruct1d(WENO_order,cctk_lsh(3),& + Bconsz(j,k,:),Bconszminus(j,k,:),Bconszplus(j,k,:),& + trivial_rp(j,k,:), hydro_excision_mask(j,k,:)) + endif + if (evolve_entropy .ne. 0) then + call GRHydro_WENOReconstruct1d(WENO_order,cctk_lsh(3),& + entropycons(j,k,:),entropyconsminus(j,k,:),entropyconsplus(j,k,:),& + trivial_rp(j,k,:), hydro_excision_mask(j,k,:)) + endif + else + !$OMP CRITICAL + call CCTK_WARN(0, "Variable type to reconstruct not recognized.") + !$OMP END CRITICAL + end if + + if(evolve_mhd.ne.0.and.clean_divergence.ne.0) then + call GRHydro_WENOReconstruct1d(WENO_order,cctk_lsh(3),& + psidc(j,k,:),psidcminus(j,k,:),psidcplus(j,k,:),& + trivial_rp(j,k,:), hydro_excision_mask(j,k,:)) + endif + + do i = 1, cctk_lsh(3) + if (trivial_rp(j,k,i)) then + SpaceMask_SetStateBitsF90(space_mask, j, k, i, type_bitsz, trivialz) + else + SpaceMask_SetStateBitsF90(space_mask, j, k, i, type_bitsz, not_trivialz) + end if + end do + end do + end do + !$OMP END PARALLEL DO + else + call CCTK_WARN(0, "Flux direction not x,y,z") + end if +!!$ WENO ends. + + deallocate(trivial_rp) + deallocate(dum,dump,dumm) + +end subroutine GRHydro_WENOReconstruct_drv + diff --git a/src/GRHydro_WENOScalars.F90 b/src/GRHydro_WENOScalars.F90 new file mode 100644 index 0000000..2aee7f1 --- /dev/null +++ b/src/GRHydro_WENOScalars.F90 @@ -0,0 +1,23 @@ + /*@@ + @file GRHydro_WENOScalars.F90 + @date Fri Jan 3 2013 + @author Ian Hawke, Christian Reisswig + @desc + Module containing the coefficient array for WENO reconstruction. + @enddesc + @@*/ + +#include "cctk.h" +#include "cctk_Parameters.h" +#include "cctk_Arguments.h" + + module GRHydro_WENOScalars + + implicit none + + CCTK_REAL, allocatable :: weno_coeffs(:, :) + CCTK_REAL, allocatable :: beta_shu(:, :) + + logical, save :: coeffs_allocated = .false. + + end module GRHydro_WENOScalars diff --git a/src/make.code.defn b/src/make.code.defn index 5a1b96e..7244b0d 100644 --- a/src/make.code.defn +++ b/src/make.code.defn @@ -68,10 +68,12 @@ SRCS = Utils.F90 \ GRHydro_PPMMReconstruct_drv.F90\ GRHydro_ENOReconstruct_drv.F90\ GRHydro_PPMReconstruct_drv.F90\ + GRHydro_WENOReconstruct_drv.F90\ GRHydro_LastMoLPostStep.c\ GRHydro_TVDReconstruct_drv.F90\ - - + GRHydro_EvolutionMask.F90\ + GRHydro_WENOReconstruct.F90\ + GRHydro_WENOScalars.F90 # Subdirectories containing source files SUBDIRS = |