diff options
| author | rhaas <rhaas@c83d129a-5a75-4d5a-9c4d-ed3a5855bf45> | 2012-06-05 20:40:27 +0000 |
|---|---|---|
| committer | rhaas <rhaas@c83d129a-5a75-4d5a-9c4d-ed3a5855bf45> | 2012-06-05 20:40:27 +0000 |
| commit | 82e74d30096e82cd06da78ae4621b3d41583b439 (patch) | |
| tree | c06f258f302978cd9577089a19f7799a9ac4a094 /src/GRHydro_CalcUpdate.F90 | |
| parent | 94b2aa0ebe60a7ee10dd4804e6b6e6cd4b541f9a (diff) | |
GRHydro: Added option to constrain evolution to 1D by
projecting out non-radial parts of the conserved
velocity RHSs.
Patch by Christian Reisswig.
git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinEvolve/GRHydro/trunk@346 c83d129a-5a75-4d5a-9c4d-ed3a5855bf45
Diffstat (limited to 'src/GRHydro_CalcUpdate.F90')
| -rw-r--r-- | src/GRHydro_CalcUpdate.F90 | 64 |
1 files changed, 64 insertions, 0 deletions
diff --git a/src/GRHydro_CalcUpdate.F90 b/src/GRHydro_CalcUpdate.F90 index 62653e9..718d618 100644 --- a/src/GRHydro_CalcUpdate.F90 +++ b/src/GRHydro_CalcUpdate.F90 @@ -72,6 +72,7 @@ subroutine UpdateCalculation(CCTK_ARGUMENTS) taurhs(i,j,k) = taurhs(i,j,k) + & (alp_l * tauflux(i-xoffset,j-yoffset,k-zoffset) - & alp_r * tauflux(i,j,k)) * idx + if(evolve_mhd.ne.0) then if(transport_constraints.ne.0) then ! we have to first compute all components of v\crossB = E and @@ -260,6 +261,7 @@ subroutine UpdateCalculation(CCTK_ARGUMENTS) taurhs(i,j,k) = taurhs(i,j,k) + & (tauflux(i-xoffset,j-yoffset,k-zoffset) - & tauflux(i,j,k)) * idx + if(evolve_mhd.ne.0) then Bconsrhs(i,j,k,1) = Bconsrhs(i,j,k,1) + & (Bconsxflux(i-xoffset,j-yoffset,k-zoffset) - & @@ -319,3 +321,65 @@ subroutine UpdateCalculation(CCTK_ARGUMENTS) end subroutine UpdateCalculation + + + + + /*@@ + @routine ConstrainSconTo1D + @date Tue 24 14:12 2012 + @author Christian Reisswig + @desc + Constrains the conserved fluid velocity to radial direction + @enddesc + @calls + @calledby + @history + @endhistory + +@@*/ + + +subroutine ConstrainSconTo1D(CCTK_ARGUMENTS) + implicit none + + DECLARE_CCTK_ARGUMENTS + DECLARE_CCTK_PARAMETERS + DECLARE_CCTK_FUNCTIONS + + CCTK_INT :: i,j,k + CCTK_REAL :: rnorm, rnormI, scon_tmp1, scon_tmp2, scon_tmp3 + + !$OMP PARALLEL DO PRIVATE(i,j,k,rnorm,rnormI,scon_tmp1,scon_tmp2,scon_tmp3) + do k = GRHydro_stencil + 1 - transport_constraints, cctk_lsh(3) - GRHydro_stencil ! we need to compute Evec on all faces/edges where the fluxes are defined + do j = GRHydro_stencil + 1 - transport_constraints, cctk_lsh(2) - GRHydro_stencil + do i = GRHydro_stencil + 1 - transport_constraints, cctk_lsh(1) - GRHydro_stencil + + + ! Eliminate non-radial fluid velocities to obtain pseudo 1D scheme + rnorm = (x(i,j,k)**2 + y(i,j,k)**2 + z(i,j,k)**2) + if (rnorm.eq.0) then + rnormI = 0 + else + rnormI = 1.0d0/rnorm + endif + + scon_tmp1 = (x(i,j,k)*scon(i,j,k,1) & + + y(i,j,k)*scon(i,j,k,2) & + + z(i,j,k)*scon(i,j,k,3)) * rnormI * x(i,j,k) + scon_tmp2 = (x(i,j,k)*scon(i,j,k,1) & + + y(i,j,k)*scon(i,j,k,2) & + + z(i,j,k)*scon(i,j,k,3)) * rnormI * y(i,j,k) + scon_tmp3 = (x(i,j,k)*scon(i,j,k,1) & + + y(i,j,k)*scon(i,j,k,2) & + + z(i,j,k)*scon(i,j,k,3)) * rnormI * z(i,j,k) + + scon(i,j,k,1) = scon_tmp1 + scon(i,j,k,2) = scon_tmp2 + scon(i,j,k,3) = scon_tmp3 + + end do + end do + end do +end subroutine ConstrainSconTo1D + |