diff options
Diffstat (limited to 'src')
| -rw-r--r-- | src/GRHydro_CalcUpdate.F90 | 64 | ||||
| -rw-r--r-- | src/GRHydro_ParamCheck.F90 | 4 |
2 files changed, 68 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 + diff --git a/src/GRHydro_ParamCheck.F90 b/src/GRHydro_ParamCheck.F90 index 9310448..132bb7d 100644 --- a/src/GRHydro_ParamCheck.F90 +++ b/src/GRHydro_ParamCheck.F90 @@ -112,6 +112,10 @@ subroutine GRHydro_ParamCheck(CCTK_ARGUMENTS) if (CCTK_Equals(riemann_solver,"HLLE").eq.0.and.coordinates_is_active.ne.0) then call CCTK_WARN(0, "There is currently no Riemann solver other than HLLE compatible with multipatch!") end if + + if (constrain_to_1D .eq. 1 .and. coordinates_is_active.eq.1) then + call CCTK_WARN(0, "Constrain to 1D option does not work with other than Cartesian coordinates. Ask Christian R. to implement it if required!") + endif if (CCTK_EQUALS(riemann_solver,"Roe").and.CCTK_EQUALS(Bvec_evolution_method,"GRHydro")) then call CCTK_PARAMWARN("Roe solver is not implemented yet for MHD") |