From 04b230f0831053faf2f6a1e798613f39ea378a37 Mon Sep 17 00:00:00 2001 From: rhaas Date: Wed, 6 Mar 2013 00:13:26 +0000 Subject: GRHydro: Reconstruct Wv for WENO5. From: Christian Reisswig git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinEvolve/GRHydro/trunk@486 c83d129a-5a75-4d5a-9c4d-ed3a5855bf45 --- src/GRHydro_WENOReconstruct_drv.F90 | 176 ++++++++++++++++++++++++++++++------ 1 file changed, 149 insertions(+), 27 deletions(-) diff --git a/src/GRHydro_WENOReconstruct_drv.F90 b/src/GRHydro_WENOReconstruct_drv.F90 index 0a6374c..bf3ff13 100644 --- a/src/GRHydro_WENOReconstruct_drv.F90 +++ b/src/GRHydro_WENOReconstruct_drv.F90 @@ -28,6 +28,63 @@ #define Bconsz(i,j,k) Bcons(i,j,k,3) +!!$ transform back to vel if Wv was used +subroutine undo_Wv(nx,velxminus,velyminus,velzminus,& + velxplus,velyplus,velzplus,& + gxx,gxy,gxz,gyy,gyz,gzz) + + implicit none + + DECLARE_CCTK_PARAMETERS + + CCTK_INT :: i,nx + CCTK_REAL :: w, agxx, agxy, agxz, agyy, agyz, agzz + CCTK_REAL, dimension(nx) :: gxx, gxy, gxz, gyy, gyz, gzz + CCTK_REAL, dimension(nx) :: velxminus,velyminus,velzminus + CCTK_REAL, dimension(nx) :: velxplus,velyplus,velzplus + + do i = grhydro_stencil, nx - grhydro_stencil + 1 + ! divide out the Loretnz factor obtained from w_lorentz = + ! sqrt(1+g_{ij} w^i w^j) for both the + ! plus and minus quantities this should by construction ensure + ! that any Lorentz factor calculated from them later on is + ! physical (ie. > 1.d0) + agxx = 0.5d0*( gxx(i) + gxx(i-1) ) + agxy = 0.5d0*( gxy(i) + gxy(i-1) ) + agxz = 0.5d0*( gxz(i) + gxz(i-1) ) + agyy = 0.5d0*( gyy(i) + gyy(i-1) ) + agyz = 0.5d0*( gyz(i) + gyz(i-1) ) + agzz = 0.5d0*( gzz(i) + gzz(i-1) ) + w = sqrt( 1.d0 + agxx*velxminus(i)*velxminus(i) + & + agyy*velyminus(i)*velyminus(i) & + + agzz*velzminus(i)*velzminus(i) + & + 2.d0*agxy*velxminus(i)*velyminus(i) & + + 2.d0*agxz*velxminus(i)*velzminus(i) + & + 2.d0*agyz*velyminus(i)*velzminus(i) ) + velxminus(i) = velxminus(i)/w + velyminus(i) = velyminus(i)/w + velzminus(i) = velzminus(i)/w + + agxx = 0.5d0*( gxx(i) + gxx(i+1) ) + agxy = 0.5d0*( gxy(i) + gxy(i+1) ) + agxz = 0.5d0*( gxz(i) + gxz(i+1) ) + agyy = 0.5d0*( gyy(i) + gyy(i+1) ) + agyz = 0.5d0*( gyz(i) + gyz(i+1) ) + agzz = 0.5d0*( gzz(i) + gzz(i+1) ) + w = sqrt( 1.d0 + agxx*velxplus(i)*velxplus(i) + & + agyy*velyplus(i)*velyplus(i) & + + agzz*velzplus(i)*velzplus(i) + & + 2.d0*agxy*velxplus(i)*velyplus(i) & + + 2.d0*agxz*velxplus(i)*velzplus(i) + & + 2.d0*agyz*velyplus(i)*velzplus(i) ) + velxplus(i) = velxplus(i)/w + velyplus(i) = velyplus(i)/w + velzplus(i) = velzplus(i)/w + end do +end subroutine undo_Wv + + + /*@@ @routine GRHydro_WENOReconstruct_drv @date Fri Jan 3 2013 @@ -49,6 +106,8 @@ subroutine GRHydro_WENOReconstruct_drv(CCTK_ARGUMENTS) ! save memory when MP is not used ! TARGET as to be before DECLARE_CCTK_ARGUMENTS for gcc 4.1 + TARGET gaa, gab, gac, gbb, gbc, gcc + TARGET gxx, gxy, gxz, gyy, gyz, gzz TARGET vel, lvel TARGET Bvec, lBvec @@ -69,14 +128,32 @@ subroutine GRHydro_WENOReconstruct_drv(CCTK_ARGUMENTS) CCTK_INT :: ierr + ! 1D arrays in x, y and z direction for reconstructing Wv instead of v +! CCTK_REAL, dimension(nx) :: vxX,vyX,vzX +! CCTK_REAL, dimension(ny) :: vxY,vyY,vzY +! CCTK_REAL, dimension(nz) :: vxZ,vyZ,vzZ + ! save memory when MP is not used CCTK_INT :: GRHydro_UseGeneralCoordinates + CCTK_REAL, DIMENSION(:,:,:), POINTER :: g11, g12, g13, g22, g23, g33 CCTK_REAL, DIMENSION(:,:,:,:), POINTER :: vup, Bprim if (GRHydro_UseGeneralCoordinates(cctkGH).ne.0) then + g11 => gaa + g12 => gab + g13 => gac + g22 => gbb + g23 => gbc + g33 => gcc vup => lvel Bprim => lBvec else + g11 => gxx + g12 => gxy + g13 => gxz + g22 => gyy + g23 => gyz + g33 => gzz vup => vel Bprim => Bvec end if @@ -194,15 +271,30 @@ subroutine GRHydro_WENOReconstruct_drv(CCTK_ARGUMENTS) 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)) + if (reconstruct_Wv.eq.0) then + 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)) + else + call GRHydro_WENOReconstruct1d(WENO_order,cctk_lsh(1),& + w_lorentz(:,j,k)*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),& + w_lorentz(:,j,k)*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),& + w_lorentz(:,j,k)*velz(:,j,k),velzminus(:,j,k),velzplus(:,j,k),& + trivial_rp(:,j,k), hydro_excision_mask(:,j,k)) + call undo_Wv(nx, velxminus(:,j,k), velyminus(:,j,k), velzminus(:,j,k),& + velxplus(:,j,k),velyplus(:,j,k),velzplus(:,j,k),& + g11(:,j,k),g12(:,j,k),g13(:,j,k),g22(:,j,k),g23(:,j,k),g33(:,j,k)) + endif 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)) @@ -296,15 +388,30 @@ subroutine GRHydro_WENOReconstruct_drv(CCTK_ARGUMENTS) 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)) + if (reconstruct_Wv.eq.0) then + 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)) + else + call GRHydro_WENOReconstruct1d(WENO_order,cctk_lsh(2),& + w_lorentz(j,:,k)*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),& + w_lorentz(j,:,k)*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),& + w_lorentz(j,:,k)*velz(j,:,k),velzminus(j,:,k),velzplus(j,:,k),& + trivial_rp(j,:,k), hydro_excision_mask(j,:,k)) + call undo_Wv(ny, velyminus(j,:,k), velzminus(j,:,k), velxminus(j,:,k),& + velyplus(j,:,k),velzplus(j,:,k),velxplus(j,:,k),& + g22(j,:,k), g23(j,:,k), g12(j,:,k), g33(j,:,k), g13(j,:,k), g11(j,:,k)) + endif 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)) @@ -398,15 +505,30 @@ subroutine GRHydro_WENOReconstruct_drv(CCTK_ARGUMENTS) 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,:)) + if (reconstruct_Wv.eq.0) then + 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,:)) + else + call GRHydro_WENOReconstruct1d(WENO_order,cctk_lsh(3),& + w_lorentz(j,k,:)*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),& + w_lorentz(j,k,:)*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),& + w_lorentz(j,k,:)*velz(j,k,:),velzminus(j,k,:),velzplus(j,k,:),& + trivial_rp(j,k,:), hydro_excision_mask(j,k,:)) + call undo_Wv(nz, velzminus(j,k,:), velxminus(j,k,:), velyminus(j,k,:),& + velzplus(j,k,:),velxplus(j,k,:),velyplus(j,k,:),& + g33(j,k,:), g13(j,k,:), g23(j,k,:), g11(j,k,:), g12(j,k,:),g22(j,k,:)) + endif 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,:)) -- cgit v1.2.3