diff options
author | rhaas <rhaas@c83d129a-5a75-4d5a-9c4d-ed3a5855bf45> | 2011-11-02 15:04:18 +0000 |
---|---|---|
committer | rhaas <rhaas@c83d129a-5a75-4d5a-9c4d-ed3a5855bf45> | 2011-11-02 15:04:18 +0000 |
commit | 34aeeb9e2d7e870f8d698aad89ea7eecf649e72d (patch) | |
tree | e5646b6d09d47ea7c8a7b647a965ec2b88db5caa /src | |
parent | c43f5e7261a9413d87c3369e15775be19eb8979e (diff) |
First try of a constraint transport scheme for unigrid
git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinEvolve/GRHydro/trunk@296 c83d129a-5a75-4d5a-9c4d-ed3a5855bf45
Diffstat (limited to 'src')
-rw-r--r-- | src/GRHydro_CalcUpdate.F90 | 51 | ||||
-rw-r--r-- | src/GRHydro_Source.F90 | 100 | ||||
-rw-r--r-- | src/GRHydro_SourceM.F90 | 7 |
3 files changed, 155 insertions, 3 deletions
diff --git a/src/GRHydro_CalcUpdate.F90 b/src/GRHydro_CalcUpdate.F90 index fffb42d..cc7b136 100644 --- a/src/GRHydro_CalcUpdate.F90 +++ b/src/GRHydro_CalcUpdate.F90 @@ -48,9 +48,9 @@ subroutine UpdateCalculation(CCTK_ARGUMENTS) if (use_weighted_fluxes == 0) then !$OMP PARALLEL DO PRIVATE(i,j,k,itracer,alp_l,alp_r,Bvec_l,Bvec_r) - do k = GRHydro_stencil + 1, cctk_lsh(3) - GRHydro_stencil - do j = GRHydro_stencil + 1, cctk_lsh(2) - GRHydro_stencil - do i = GRHydro_stencil + 1, cctk_lsh(1) - GRHydro_stencil + do k = GRHydro_stencil + 1 - transport_constraints, cctk_lsh(3) - GRHydro_stencil ! mean we need to compute fluxes for one more point the GRHydro_stencil would indicate + do j = GRHydro_stencil + 1 - transport_constraints, cctk_lsh(2) - GRHydro_stencil + do i = GRHydro_stencil + 1 - transport_constraints, cctk_lsh(1) - GRHydro_stencil alp_l = 0.5d0 * (alp(i,j,k) + & alp(i-xoffset,j-yoffset,k-zoffset)) @@ -87,6 +87,19 @@ subroutine UpdateCalculation(CCTK_ARGUMENTS) (alp_l * psidcflux(i-xoffset,j-yoffset,k-zoffset) - & alp_r * psidcflux(i,j,k)) * idx endif + if(transport_constraints.ne.0) then + ! Evec lives on edges of cell: Evec(i,j,k,1) is at edge i,j+1/2,k+1/2 ie. the lower-front edge of cell (i,j,k) + if(flux_direction.eq.1) then + Evec(i,j,k,2) = Evec(i,j,k,2) + 0.25d0 * alp_r*(Bconszflux(i,j,k) + Bconszflux(i,j ,k+1)) + Evec(i,j,k,3) = Evec(i,j,k,3) - 0.25d0 * alp_r*(Bconsxflux(i,j,k) + Bconsyflux(i,j+1,k )) + elseif(flux_direction.eq.2) then + Evec(i,j,k,1) = Evec(i,j,k,1) - 0.25d0 * alp_r*(Bconszflux(i,j,k) + Bconszflux(i ,j,k+1)) + Evec(i,j,k,3) = Evec(i,j,k,3) + 0.25d0 * alp_r*(Bconsxflux(i,j,k) + Bconsxflux(i+1,j,k )) + elseif(flux_direction.eq.3) then + Evec(i,j,k,1) = Evec(i,j,k,1) + 0.25d0 * alp_r*(Bconsyflux(i,j,k) + Bconsyflux(i ,j+1,k)) + Evec(i,j,k,2) = Evec(i,j,k,2) - 0.25d0 * alp_r*(Bconsxflux(i,j,k) + Bconsxflux(i+1,j ,k)) + end if + end if if(track_divB.ne.0) then Bvec_l = 0.5d0 * (Bvec(i,j,k,flux_direction) + & Bvec(i-xoffset,j-yoffset,k-zoffset,flux_direction)) @@ -205,6 +218,10 @@ subroutine UpdateCalculation(CCTK_ARGUMENTS) else if (CCTK_EQUALS(method_type, "Flux split FD")) then + if (transport_constraints .ne. 0) then + call CCTK_WARN(0, "Not supported") + end if + do k = GRHydro_stencil + 1, cctk_lsh(3) - GRHydro_stencil do j = GRHydro_stencil + 1, cctk_lsh(2) - GRHydro_stencil do i = GRHydro_stencil + 1, cctk_lsh(1) - GRHydro_stencil @@ -253,6 +270,34 @@ subroutine UpdateCalculation(CCTK_ARGUMENTS) enddo end if + + if (transport_constraints.ne.0 .and. flux_direction.eq.1) then ! HACK: x direction is last + ! RH: I think one could wrap all of this into a single do loop and remove the + ! Evec storage + do k = GRHydro_stencil + 1, cctk_lsh(3) - GRHydro_stencil + do j = GRHydro_stencil + 1, cctk_lsh(2) - GRHydro_stencil + do i = GRHydro_stencil + 1, cctk_lsh(1) - GRHydro_stencil + Bconsrhs(i,j,k,1) = & !Bconsrhs(i,j,k,1) & + - 0.5d0 * (Evec(i-1,j-1,k-1,2)-Evec(i-1,j-1,k ,2)) / CCTK_DELTA_SPACE(3) & + - 0.5d0 * (Evec(i-1,j ,k-1,3)-Evec(i-1,j-1,k-1,3)) / CCTK_DELTA_SPACE(2) & + - 0.5d0 * (Evec(i ,j-1,k-1,2)-Evec(i ,j-1,k ,2)) / CCTK_DELTA_SPACE(3) & + - 0.5d0 * (Evec(i ,j ,k-1,3)-Evec(i ,j-1,k-1,3)) / CCTK_DELTA_SPACE(2) + Bconsrhs(i,j,k,2) = & !Bconsrhs(i,j,k,2) & + - 0.5d0 * (Evec(i-1,j-1,k-1,3)-Evec(i-1,j ,k-1,3)) / CCTK_DELTA_SPACE(1) & + - 0.5d0 * (Evec(i ,j-1,k-1,1)-Evec(i-1,j-1,k-1,1)) / CCTK_DELTA_SPACE(3) & + - 0.5d0 * (Evec(i-1,j-1,k ,3)-Evec(i-1,j ,k ,3)) / CCTK_DELTA_SPACE(1) & + - 0.5d0 * (Evec(i ,j-1,k ,1)-Evec(i-1,j-1,k ,1)) / CCTK_DELTA_SPACE(3) + Bconsrhs(i,j,k,3) = & !Bconsrhs(i,j,k,3) & + - 0.5d0 * (Evec(i-1,j-1,k-1,1)-Evec(i ,j-1,k-1,1)) / CCTK_DELTA_SPACE(1) & + - 0.5d0 * (Evec(i-1,j-1,k ,2)-Evec(i-1,j-1,k-1,2)) / CCTK_DELTA_SPACE(2) & + - 0.5d0 * (Evec(i-1,j ,k-1,1)-Evec(i ,j ,k-1,1)) / CCTK_DELTA_SPACE(1) & + - 0.5d0 * (Evec(i-1,j ,k ,2)-Evec(i-1,j ,k-1,2)) / CCTK_DELTA_SPACE(2) + enddo + enddo + enddo + end if + + return diff --git a/src/GRHydro_Source.F90 b/src/GRHydro_Source.F90 index 6758c90..7b04b8c 100644 --- a/src/GRHydro_Source.F90 +++ b/src/GRHydro_Source.F90 @@ -93,6 +93,11 @@ subroutine SourceTerms(CCTK_ARGUMENTS) CCTK_REAL, DIMENSION(:,:,:), POINTER :: beta1, beta2, beta3 CCTK_REAL, DIMENSION(:,:,:,:), POINTER :: vup + ! poison hack + CCTK_INT, SAVE :: reflevel = -1 + CCTK_INT, SAVE :: last_iteration_seen = -1 + CCTK_INT, SAVE :: mol_substep = -1 + if (GRHydro_UseGeneralCoordinates(cctkGH).ne.0) then g11 => gaa g12 => gab @@ -449,6 +454,101 @@ subroutine SourceTerms(CCTK_ARGUMENTS) deallocate(force_spatial_second_order) +#if(0) /* poison edges of domain */ + if(last_iteration_seen .ne. cctk_iteration .or. reflevel .ne. grhydro_reflevel) then + last_iteration_seen = cctk_iteration + reflevel = grhydro_reflevel + mol_substep = 0 + else + mol_substep = mol_substep + 1 + end if + do k = 1, GRHydro_stencil*mol_substep + do j = 1, cctk_lsh(2) + do i = 1, cctk_lsh(1) + dens(i,j,k) = -1d100 + Scon(i,j,k,1) = -1d100 + Scon(i,j,k,2) = -1d100 + Scon(i,j,k,3) = -1d100 + tau(i,j,k) = -1d100 + Bcons(i,j,k,1) = -1d100 + Bcons(i,j,k,2) = -1d100 + Bcons(i,j,k,3) = -1d100 + end do + end do + end do + do k = cctk_lsh(3)-GRHydro_stencil*mol_substep+1, cctk_lsh(3) + do j = 1, cctk_lsh(2) + do i = 1, cctk_lsh(1) + dens(i,j,k) = -1d100 + Scon(i,j,k,1) = -1d100 + Scon(i,j,k,2) = -1d100 + Scon(i,j,k,3) = -1d100 + tau(i,j,k) = -1d100 + Bcons(i,j,k,1) = -1d100 + Bcons(i,j,k,2) = -1d100 + Bcons(i,j,k,3) = -1d100 + end do + end do + end do + do i = 1, GRHydro_stencil*mol_substep + do k = 1, cctk_lsh(3) + do j = 1, cctk_lsh(2) + dens(i,j,k) = -1d100 + Scon(i,j,k,1) = -1d100 + Scon(i,j,k,2) = -1d100 + Scon(i,j,k,3) = -1d100 + tau(i,j,k) = -1d100 + Bcons(i,j,k,1) = -1d100 + Bcons(i,j,k,2) = -1d100 + Bcons(i,j,k,3) = -1d100 + end do + end do + end do + do i = cctk_lsh(1)-GRHydro_stencil*mol_substep+1, cctk_lsh(1) + do k = 1, cctk_lsh(3) + do j = 1, cctk_lsh(2) + dens(i,j,k) = -1d100 + Scon(i,j,k,1) = -1d100 + Scon(i,j,k,2) = -1d100 + Scon(i,j,k,3) = -1d100 + tau(i,j,k) = -1d100 + Bcons(i,j,k,1) = -1d100 + Bcons(i,j,k,2) = -1d100 + Bcons(i,j,k,3) = -1d100 + end do + end do + end do + do j = 1, GRHydro_stencil*mol_substep + do i = 1, cctk_lsh(1) + do k = 1, cctk_lsh(3) + dens(i,j,k) = -1d100 + Scon(i,j,k,1) = -1d100 + Scon(i,j,k,2) = -1d100 + Scon(i,j,k,3) = -1d100 + tau(i,j,k) = -1d100 + Bcons(i,j,k,1) = -1d100 + Bcons(i,j,k,2) = -1d100 + Bcons(i,j,k,3) = -1d100 + end do + end do + end do + do j = cctk_lsh(2)-GRHydro_stencil*mol_substep+1, cctk_lsh(2) + do i = 1, cctk_lsh(1) + do k = 1, cctk_lsh(3) + dens(i,j,k) = -1d100 + Scon(i,j,k,1) = -1d100 + Scon(i,j,k,2) = -1d100 + Scon(i,j,k,3) = -1d100 + tau(i,j,k) = -1d100 + Bcons(i,j,k,1) = -1d100 + Bcons(i,j,k,2) = -1d100 + Bcons(i,j,k,3) = -1d100 + end do + end do + end do +#endif + + end subroutine SourceTerms diff --git a/src/GRHydro_SourceM.F90 b/src/GRHydro_SourceM.F90 index 73334e3..1456500 100644 --- a/src/GRHydro_SourceM.F90 +++ b/src/GRHydro_SourceM.F90 @@ -127,6 +127,13 @@ subroutine SourceTermsM(CCTK_ARGUMENTS) if (track_divB .ne. 0) then divB=0.d0 endif + + if (transport_constraints .ne. 0) then + Evec = 0.d0 + Bconsrhsx(i,j,k) = 0.d0 + Bconsrhsy(i,j,k) = 0.d0 + Bconsrhsz(i,j,k) = 0.d0 + endif !!$ Set up the array for checking the order. We switch to second order |