[Commits] [svn:einsteintoolkit] GRHydro/trunk/ (Rev. 228)
bcmsma at astro.rit.edu
bcmsma at astro.rit.edu
Tue Apr 12 22:29:13 CDT 2011
User: bmundim
Date: 2011/04/12 10:29 PM
Removed:
/trunk/src/
GRHydro_CalcUpdateM.F90
Modified:
/trunk/
interface.ccl, param.ccl, schedule.ccl
/trunk/src/
GRHydro_CalcUpdate.F90, GRHydro_FluxM.F90, GRHydro_HLLEM.F90, GRHydro_PreLoop.F90, GRHydro_SourceM.F90
Log:
RIT MHD development:
1) Fix a bug in the divergence cleaning implementation:
a psidc term in the induction equation was implemented
as a source term where it was supposed to be coded as
part of the flux calculation.
2) Introduce divB as a diagnostic grid function.
3) Remove an old file: GRHydro_CalcUpdateM.F90
File Changes:
Directory: /trunk/src/
======================
File [modified]: GRHydro_CalcUpdate.F90
Delta lines: +16 -2
===================================================================
--- trunk/src/GRHydro_CalcUpdate.F90 2011-04-06 14:37:10 UTC (rev 227)
+++ trunk/src/GRHydro_CalcUpdate.F90 2011-04-13 03:29:13 UTC (rev 228)
@@ -39,7 +39,7 @@
DECLARE_CCTK_FUNCTIONS
CCTK_INT :: i,j,k,itracer
- CCTK_REAL :: idx, alp_l, alp_r
+ CCTK_REAL :: idx, alp_l, alp_r, Bvec_l, Bvec_r
CCTK_INT :: type_bits, atmosphere, not_atmosphere
@@ -55,7 +55,7 @@
if (use_weighted_fluxes == 0) then
- !$OMP PARALLEL DO PRIVATE(i,j,itracer,alp_l, alp_r)
+ !$OMP PARALLEL DO PRIVATE(i,j,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
@@ -95,6 +95,13 @@
(alp_l * psidcflux(i-xoffset,j-yoffset,k-zoffset) - &
alp_r * psidcflux(i,j,k)) * idx
endif
+ 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))
+ Bvec_r = 0.5d0 * (Bvec(i,j,k,flux_direction) + &
+ Bvec(i+xoffset,j+yoffset,k+zoffset,flux_direction))
+ divB(i,j,k) = divB(i,j,k) + ( alp_l * Bvec_l - alp_r * Bvec_r ) * idx
+ endif
endif
if (evolve_tracer .ne. 0) then
@@ -248,6 +255,13 @@
(psidcflux(i-xoffset,j-yoffset,k-zoffset) - &
psidcflux(i,j,k)) * idx
endif
+ 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))
+ Bvec_r = 0.5d0 * (Bvec(i,j,k,flux_direction) + &
+ Bvec(i+xoffset,j+yoffset,k+zoffset,flux_direction))
+ divB(i,j,k) = divB(i,j,k) + ( Bvec_l - Bvec_r ) * idx
+ endif
endif
enddo
File [removed]: GRHydro_CalcUpdateM.F90
Delta lines: +0 -255
===================================================================
--- trunk/src/GRHydro_CalcUpdateM.F90 2011-04-06 14:37:10 UTC (rev 227)
+++ trunk/src/GRHydro_CalcUpdateM.F90 2011-04-13 03:29:13 UTC (rev 228)
@@ -1,255 +0,0 @@
- /*@@
- @file GRHydro_CalcUpdateM.F90
- @date Oct 10, 2010
- @author Joshua Faber, Scott Noble, Bruno Mundim, Ian Hawke
- @desc
- Calculates the update terms given the fluxes. Moved to here so that
- @enddesc
- @@*/
-
-#include "cctk.h"
-#include "cctk_Arguments.h"
-#include "cctk_Parameters.h"
-#include "cctk_Functions.h"
-#include "SpaceMask.h"
-
- /*@@
- @routine UpdateCalculationM
- @date Oct 10, 2010
- @author Joshua Faber, Scott Noble, Bruno Mundim, Ian Hawke
- @desc
- Calculates the update terms from the fluxes.
- @enddesc
- @calls
- @calledby
- @history
- Moved out of the Riemann solver routines to make the FishEye /
- weighted flux calculation easier.
- @endhistory
-
-@@*/
-
-
-subroutine UpdateCalculationM(CCTK_ARGUMENTS)
-
- implicit none
-
- DECLARE_CCTK_ARGUMENTS
- DECLARE_CCTK_PARAMETERS
- DECLARE_CCTK_FUNCTIONS
-
- CCTK_INT :: i,j,k,itracer
- CCTK_REAL :: idx, alp_l, alp_r
-
- CCTK_INT :: type_bits, atmosphere, not_atmosphere
-
- call SpaceMask_GetTypeBits(type_bits, "Hydro_Atmosphere")
- call SpaceMask_GetStateBits(atmosphere, "Hydro_Atmosphere",&
- "in_atmosphere")
- call SpaceMask_GetStateBits(not_atmosphere, "Hydro_Atmosphere",&
- "not_in_atmosphere")
-
- idx = 1.d0 / CCTK_DELTA_SPACE(flux_direction)
-
- if (CCTK_EQUALS(method_type, "RSA FV")) then
-
- if (use_weighted_fluxes == 0) then
-
- !$OMP PARALLEL DO PRIVATE(i,j,itracer,alp_l, alp_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
-
- alp_l = 0.5d0 * (alp(i,j,k) + &
- alp(i-xoffset,j-yoffset,k-zoffset))
- alp_r = 0.5d0 * (alp(i,j,k) + &
- alp(i+xoffset,j+yoffset,k+zoffset))
-
- densrhs(i,j,k) = densrhs(i,j,k) + &
- (alp_l * densflux(i-xoffset,j-yoffset,k-zoffset) - &
- alp_r * densflux(i,j,k)) * idx
- srhs(i,j,k,1) = srhs(i,j,k,1) + &
- (alp_l * sxflux(i-xoffset,j-yoffset,k-zoffset) - &
- alp_r * sxflux(i,j,k)) * idx
- srhs(i,j,k,2) = srhs(i,j,k,2) + &
- (alp_l * syflux(i-xoffset,j-yoffset,k-zoffset) - &
- alp_r * syflux(i,j,k)) * idx
- srhs(i,j,k,3) = srhs(i,j,k,3) + &
- (alp_l * szflux(i-xoffset,j-yoffset,k-zoffset) - &
- alp_r * szflux(i,j,k)) * idx
- taurhs(i,j,k) = taurhs(i,j,k) + &
- (alp_l * tauflux(i-xoffset,j-yoffset,k-zoffset) - &
- alp_r * tauflux(i,j,k)) * idx
- Bvecrhs(i,j,k,1) = Bvecrhs(i,j,k,1) + &
- (alp_l * Bvecxflux(i-xoffset,j-yoffset,k-zoffset) - &
- alp_r * Bvecxflux(i,j,k)) * idx
- Bvecrhs(i,j,k,2) = Bvecrhs(i,j,k,2) + &
- (alp_l * Bvecyflux(i-xoffset,j-yoffset,k-zoffset) - &
- alp_r * Bvecyflux(i,j,k)) * idx
- Bvecrhs(i,j,k,3) = Bvecrhs(i,j,k,3) + &
- (alp_l * Bveczflux(i-xoffset,j-yoffset,k-zoffset) - &
- alp_r * Bveczflux(i,j,k)) * idx
- if(clean_divergence.ne.0) then
- psidcrhs(i,j,k) = psidcrhs(i,j,k) + &
- (alp_l * psidcflux(i-xoffset,j-yoffset,k-zoffset) - &
- alp_r * psidcflux(i,j,k)) * idx
- endif
-
- if (evolve_tracer .ne. 0) then
- do itracer=1,number_of_tracers
- cons_tracerrhs(i,j,k,itracer) = cons_tracerrhs(i,j,k,itracer) + &
- (alp_l * cons_tracerflux(i-xoffset,j-yoffset,k-zoffset,itracer) - &
- alp_r * cons_tracerflux(i,j,k,itracer)) * idx
- enddo
- end if
-
- if (evolve_Y_e .ne. 0) then
- Y_e_con_rhs(i,j,k) = Y_e_con_rhs(i,j,k) + &
- (alp_l * Y_e_con_flux(i-xoffset,j-yoffset,k-zoffset) - &
- alp_r * Y_e_con_flux(i,j,k)) * idx
- end if
-
- if (wk_atmosphere .eq. 1) then
-
- if ( (atmosphere_mask(i,j,k) .eq. 1) .or. &
- (SpaceMask_CheckStateBitsF90(space_mask,i,j,k,type_bits,atmosphere)) ) then
-
-!!$ We are in the atmosphere so the momentum flux must vanish
-
- srhs(i,j,k,:) = 0.d0
-
- if ( ( (atmosphere_mask(i-1,j ,k ) .eq. 1) .and. &
- (atmosphere_mask(i+1,j ,k ) .eq. 1) .and. &
- (atmosphere_mask(i ,j-1,k ) .eq. 1) .and. &
- (atmosphere_mask(i ,j+1,k ) .eq. 1) .and. &
- (atmosphere_mask(i ,j ,k-1) .eq. 1) .and. &
- (atmosphere_mask(i ,j ,k+1) .eq. 1) &
- ) .or. &
- ( (SpaceMask_CheckStateBitsF90(space_mask,i-1,j ,k ,type_bits,atmosphere)) .and. &
- (SpaceMask_CheckStateBitsF90(space_mask,i+1,j ,k ,type_bits,atmosphere)) .and. &
- (SpaceMask_CheckStateBitsF90(space_mask,i ,j-1,k ,type_bits,atmosphere)) .and. &
- (SpaceMask_CheckStateBitsF90(space_mask,i ,j+1,k ,type_bits,atmosphere)) .and. &
- (SpaceMask_CheckStateBitsF90(space_mask,i ,j ,k-1,type_bits,atmosphere)) .and. &
- (SpaceMask_CheckStateBitsF90(space_mask,i ,j ,k+1,type_bits,atmosphere)) &
- ) &
- ) then
-
-!!$ All neighbours are also atmosphere so all rhs vanish
-
- densrhs(i,j,k) = 0.d0
- taurhs(i,j,k) = 0.d0
-!!$
-!!$ We should still evolve the B-field in the atmosphere
-!!$
-
- end if
- end if
-
- end if
-
- enddo
- enddo
- enddo
- !$OMP END PARALLEL DO
-
- else
-
- call CCTK_WARN(0, "Not supported")
-
-!!$ 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
-!!$
-!!$ alp_l = 0.5d0 * (alp(i,j,k) + &
-!!$ alp(i-xoffset,j-yoffset,k-zoffset))
-!!$ alp_r = 0.5d0 * (alp(i,j,k) + &
-!!$ alp(i+xoffset,j+yoffset,k+zoffset))
-!!$
-!!$ densrhs(i,j,k) = densrhs(i,j,k) + &
-!!$ (alp_l * &
-!!$ &cell_surface(i-xoffset,j-yoffset,k-zoffset,flux_direction) * &
-!!$ &densflux(i-xoffset,j-yoffset,k-zoffset) - &
-!!$ alp_r * &
-!!$ &cell_surface(i,j,k,flux_direction) * &
-!!$ &densflux(i,j,k)) * idx / cell_volume(i,j,k)
-!!$ sxrhs(i,j,k) = sxrhs(i,j,k) + &
-!!$ (alp_l * &
-!!$ &cell_surface(i-xoffset,j-yoffset,k-zoffset,flux_direction) * &
-!!$ &sxflux(i-xoffset,j-yoffset,k-zoffset) - &
-!!$ alp_r * &
-!!$ &cell_surface(i,j,k,flux_direction) * &
-!!$ &sxflux(i,j,k)) * idx / cell_volume(i,j,k)
-!!$ syrhs(i,j,k) = syrhs(i,j,k) + &
-!!$ (alp_l * &
-!!$ &cell_surface(i-xoffset,j-yoffset,k-zoffset,flux_direction) * &
-!!$ &syflux(i-xoffset,j-yoffset,k-zoffset) - &
-!!$ alp_r * &
-!!$ &cell_surface(i,j,k,flux_direction) * &
-!!$ &syflux(i,j,k)) * idx / cell_volume(i,j,k)
-!!$ szrhs(i,j,k) = szrhs(i,j,k) + &
-!!$ (alp_l * &
-!!$ &cell_surface(i-xoffset,j-yoffset,k-zoffset,flux_direction) * &
-!!$ &szflux(i-xoffset,j-yoffset,k-zoffset) - &
-!!$ alp_r * &
-!!$ &cell_surface(i,j,k,flux_direction) * &
-!!$ &szflux(i,j,k)) * idx / cell_volume(i,j,k)
-!!$ taurhs(i,j,k) = taurhs(i,j,k) + &
-!!$ (alp_l * &
-!!$ &cell_surface(i-xoffset,j-yoffset,k-zoffset,flux_direction) * &
-!!$ &tauflux(i-xoffset,j-yoffset,k-zoffset) - &
-!!$ alp_r * &
-!!$ &cell_surface(i,j,k,flux_direction) * &
-!!$ &tauflux(i,j,k)) * idx / cell_volume(i,j,k)
-!!$
-!!$ enddo
-!!$ enddo
-!!$ enddo
-
- end if
-
- else if (CCTK_EQUALS(method_type, "Flux split FD")) then
-
- 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
-
- densrhs(i,j,k) = densrhs(i,j,k) + &
- (densflux(i-xoffset,j-yoffset,k-zoffset) - &
- densflux(i,j,k)) * idx
- srhs(i,j,k,1) = srhs(i,j,k,1) + &
- (sxflux(i-xoffset,j-yoffset,k-zoffset) - &
- sxflux(i,j,k)) * idx
- srhs(i,j,k,2) = srhs(i,j,k,2) + &
- (syflux(i-xoffset,j-yoffset,k-zoffset) - &
- syflux(i,j,k)) * idx
- srhs(i,j,k,3) = srhs(i,j,k,3) + &
- (szflux(i-xoffset,j-yoffset,k-zoffset) - &
- szflux(i,j,k)) * idx
- taurhs(i,j,k) = taurhs(i,j,k) + &
- (tauflux(i-xoffset,j-yoffset,k-zoffset) - &
- tauflux(i,j,k)) * idx
- Bvecrhs(i,j,k,1) = Bvecrhs(i,j,k,1) + &
- (Bvecxflux(i-xoffset,j-yoffset,k-zoffset) - &
- Bvecxflux(i,j,k)) * idx
- Bvecrhs(i,j,k,2) = Bvecrhs(i,j,k,2) + &
- (Bvecyflux(i-xoffset,j-yoffset,k-zoffset) - &
- Bvecyflux(i,j,k)) * idx
- Bvecrhs(i,j,k,3) = Bvecrhs(i,j,k,3) + &
- (Bveczflux(i-xoffset,j-yoffset,k-zoffset) - &
- Bveczflux(i,j,k)) * idx
- if(clean_divergence.ne.0) then
- psidcrhs(i,j,k) = psidcrhs(i,j,k) + &
- (psidcflux(i-xoffset,j-yoffset,k-zoffset) - &
- psidcflux(i,j,k)) * idx
- endif
-
- enddo
- enddo
- enddo
-
- end if
-
- return
-
-end subroutine UpdateCalculationM
-
File [modified]: GRHydro_FluxM.F90
Delta lines: +1 -12
===================================================================
--- trunk/src/GRHydro_FluxM.F90 2011-04-06 14:37:10 UTC (rev 227)
+++ trunk/src/GRHydro_FluxM.F90 2011-04-13 03:29:13 UTC (rev 228)
@@ -47,21 +47,10 @@
!!$ [psi^6 tau] vtilde^i +p* v^i - alp b^0 B^i/w
tauf = tau*vxt + psipstar*velm - ab0*psiBx/w
-!!$ [psi^6 (B^k vtilde^i - B^i vtilde^k]
+!!$ [psi^6 (B^k vtilde^i - B^i vtilde^k)]
bxf = 0.0
byf = psiBy * vxt - psiBx*vyt
bzf = psiBz * vxt - psiBx*vzt
end subroutine num_x_fluxM
-
-
-
-
-
-
-
-
-
-
-
File [modified]: GRHydro_HLLEM.F90
Delta lines: +51 -15
===================================================================
--- trunk/src/GRHydro_HLLEM.F90 2011-04-06 14:37:10 UTC (rev 227)
+++ trunk/src/GRHydro_HLLEM.F90 2011-04-13 03:29:13 UTC (rev 228)
@@ -172,6 +172,11 @@
avg_det = SPATIAL_DETERMINANT(gxxh,gxyh,gxzh,gyyh,gyzh,gzzh)
+ call UpperMetric(uxxh, uxyh, uxzh, uyyh, uyzh, uzzh, &
+ avg_det,gxxh, gxyh, gxzh, &
+ gyyh, gyzh, gzzh)
+
+
vxtp = prim_p(2)-avg_betax/avg_alp
vytp = prim_p(3)-avg_betay/avg_alp
vztp = prim_p(4)-avg_betaz/avg_alp
@@ -218,7 +223,10 @@
vxtp,vytp,vztp,pressstarp,bxlowp,bylowp,bzlowp,ab0p,wp, &
avg_det,avg_alp,avg_beta)
if(clean_divergence.ne.0) then
- psidcf = cons_p(6)
+ f1(6)=f1(6)+uxxh*psidcp
+ f1(7)=f1(7)+uxyh*psidcp
+ f1(8)=f1(8)+uxzh*psidcp
+ psidcf = ch_dc*ch_dc*cons_p(6)
endif
else if (flux_direction == 2) then
@@ -228,7 +236,10 @@
vytp,vztp,vxtp,pressstarp,bylowp,bzlowp,bxlowp,ab0p,wp, &
avg_det,avg_alp,avg_beta)
if(clean_divergence.ne.0) then
- psidcf = cons_p(7)
+ f1(6)=f1(6)+uxyh*psidcp
+ f1(7)=f1(7)+uyyh*psidcp
+ f1(8)=f1(8)+uyzh*psidcp
+ psidcf = ch_dc*ch_dc*cons_p(7)
endif
else if (flux_direction == 3) then
@@ -238,7 +249,10 @@
vztp,vxtp,vytp,pressstarp,bzlowp,bxlowp,bylowp,ab0p,wp, &
avg_det,avg_alp,avg_beta)
if(clean_divergence.ne.0) then
- psidcf = cons_p(8)
+ f1(6)=f1(6)+uxzh*psidcp
+ f1(7)=f1(7)+uyzh*psidcp
+ f1(8)=f1(8)+uzzh*psidcp
+ psidcf = ch_dc*ch_dc*cons_p(8)
endif
else
@@ -247,10 +261,6 @@
else !!! The end of this branch is right at the bottom of the routine
- call UpperMetric(uxxh, uxyh, uxzh, uyyh, uyzh, uzzh, &
- avg_det,gxxh, gxyh, gxzh, &
- gyyh, gyzh, gzzh)
-
if (flux_direction == 1) then
usendh = uxxh
else if (flux_direction == 2) then
@@ -302,8 +312,14 @@
avg_det,avg_alp,avg_beta)
if(clean_divergence.ne.0) then
- psidcfp = cons_p(6)
- psidcfm = cons_m(6)
+ fminus(6)=fminus(6)+uxxh*psidcm
+ fminus(7)=fminus(7)+uxyh*psidcm
+ fminus(8)=fminus(8)+uxzh*psidcm
+ fplus(6)=fplus(6)+uxxh*psidcp
+ fplus(7)=fplus(7)+uxyh*psidcp
+ fplus(8)=fplus(8)+uxzh*psidcp
+ psidcfp = ch_dc*ch_dc*cons_p(6)
+ psidcfm = ch_dc*ch_dc*cons_m(6)
endif
else if (flux_direction == 2) then
@@ -329,8 +345,14 @@
vytm,vztm,vxtm,pressstarm,bylowm,bzlowm,bxlowm,ab0m,wm, &
avg_det,avg_alp,avg_beta)
if(clean_divergence.ne.0) then
- psidcfp = cons_p(7)
- psidcfm = cons_m(7)
+ fminus(6)=fminus(6)+uxyh*psidcm
+ fminus(7)=fminus(7)+uyyh*psidcm
+ fminus(8)=fminus(8)+uyzh*psidcm
+ fplus(6)=fplus(6)+uxyh*psidcp
+ fplus(7)=fplus(7)+uyyh*psidcp
+ fplus(8)=fplus(8)+uyzh*psidcp
+ psidcfp = ch_dc*ch_dc*cons_p(7)
+ psidcfm = ch_dc*ch_dc*cons_m(7)
endif
else if (flux_direction == 3) then
@@ -356,8 +378,14 @@
vztm,vxtm,vytm,pressstarm,bzlowm,bxlowm,bylowm,ab0m,wm, &
avg_det,avg_alp,avg_beta)
if(clean_divergence.ne.0) then
- psidcfp = cons_p(8)
- psidcfm = cons_m(8)
+ fminus(6)=fminus(6)+uxzh*psidcm
+ fminus(7)=fminus(7)+uyzh*psidcm
+ fminus(8)=fminus(8)+uzzh*psidcm
+ fplus(6)=fplus(6)+uxzh*psidcp
+ fplus(7)=fplus(7)+uyzh*psidcp
+ fplus(8)=fplus(8)+uzzh*psidcp
+ psidcfp = ch_dc*ch_dc*cons_p(8)
+ psidcfm = ch_dc*ch_dc*cons_m(8)
endif
else
@@ -389,8 +417,16 @@
if(clean_divergence.ne.0) then
psidcdiff = psidcm - psidcp
- psidcf = (charmax * psidcfp - charmin * psidcfm + &
- charmax * charmin * psidcdiff) / charpm
+
+!!$ psidcf = (charmax * psidcfp - charmin * psidcfm + &
+!!$ charmax * charmin * psidcdiff) / charpm
+
+!!$ Wavespeeds for psidc are +/-c, not Fast Magnetosonic?
+
+ psidcf = (1.d0 * psidcfp - (-1.d0) * psidcfm + &
+ 1.d0 * (-1.d0) * psidcdiff) / 2.d0
+
+
endif
end if !!! The end of the SpaceMask check for a trivial RP.
File [modified]: GRHydro_PreLoop.F90
Delta lines: +29 -0
===================================================================
--- trunk/src/GRHydro_PreLoop.F90 2011-04-06 14:37:10 UTC (rev 227)
+++ trunk/src/GRHydro_PreLoop.F90 2011-04-13 03:29:13 UTC (rev 228)
@@ -94,3 +94,32 @@
end subroutine GRHydro_Scalar_Setup
+
+
+
+ /*@@
+ @routine GRHydro_DivBInit
+ @date Apr 06, 2011
+ @author Bruno Mundim
+ @desc
+ Set divB=0 initially.
+ @enddesc
+ @calls
+ @calledby
+ @history
+
+ @endhistory
+
+@@*/
+
+subroutine GRHydro_DivBInit(CCTK_ARGUMENTS)
+
+ implicit none
+
+ DECLARE_CCTK_ARGUMENTS
+
+ divB=0.0
+
+end subroutine GRHydro_DivBInit
+
+
File [modified]: GRHydro_SourceM.F90
Delta lines: +5 -20
===================================================================
--- trunk/src/GRHydro_SourceM.F90 2011-04-06 14:37:10 UTC (rev 227)
+++ trunk/src/GRHydro_SourceM.F90 2011-04-13 03:29:13 UTC (rev 228)
@@ -116,6 +116,10 @@
if (clean_divergence .ne. 0) then
psidcrhs=0.d0
endif
+
+ if (track_divB .ne. 0) then
+ divB=0.d0
+ endif
!!$ Set up the array for checking the order. We switch to second order
@@ -295,22 +299,12 @@
dy_alp = DIFF_Y_2(alp)
dz_alp = DIFF_Z_2(alp)
- if(clean_divergence.ne.0) then
- dx_psidc = DIFF_X_2(psidc)
- dy_psidc = DIFF_Y_2(psidc)
- dz_psidc = DIFF_Z_2(psidc)
- endif
-
else
dx_alp = DIFF_X_4(alp)
dy_alp = DIFF_Y_4(alp)
dz_alp = DIFF_Z_4(alp)
- if(clean_divergence.ne.0) then
- dx_psidc = DIFF_X_4(psidc)
- dy_psidc = DIFF_Y_4(psidc)
- dz_psidc = DIFF_Z_4(psidc)
- endif
+
end if
if (local_spatial_order .eq. 2) then
@@ -436,21 +430,12 @@
rhohstarW2*invalp*(vlowx*dz_betax + vlowy*dz_betay + vlowz*dz_betaz) -&
bt*(bxlow*dz_betax + bylow*dz_betay + bzlow*dz_betaz)
- if(clean_divergence.ne.0) then
- bvcx_source = uxx*dx_psidc+uxy*dy_psidc+uxz*dz_psidc
- bvcy_source = uxy*dx_psidc+uyy*dy_psidc+uyz*dz_psidc
- bvcz_source = uxz*dx_psidc+uyz*dy_psidc+uzz*dz_psidc
- endif
-
densrhs(i,j,k) = 0.d0
srhs(i,j,k,1) = alp(i,j,k)*sqrtdet*sx_source
srhs(i,j,k,2) = alp(i,j,k)*sqrtdet*sy_source
srhs(i,j,k,3) = alp(i,j,k)*sqrtdet*sz_source
taurhs(i,j,k) = alp(i,j,k)*sqrtdet*tau_source
if(clean_divergence.ne.0) then
- Bvecrhs(i,j,k,1) = alp(i,j,k)*sqrtdet*bvcx_source
- Bvecrhs(i,j,k,2) = alp(i,j,k)*sqrtdet*bvcy_source
- Bvecrhs(i,j,k,3) = alp(i,j,k)*sqrtdet*bvcz_source
psidcrhs(i,j,k) = -1.d0*ch_dc*ch_dc/cp_dc/cp_dc*psidc(i,j,k)
endif
Directory: /trunk/
==================
File [modified]: interface.ccl
Delta lines: +2 -0
===================================================================
--- trunk/interface.ccl 2011-04-06 14:37:10 UTC (rev 227)
+++ trunk/interface.ccl 2011-04-13 03:29:13 UTC (rev 228)
@@ -332,6 +332,8 @@
real psidcrhs type = GF Timelevels = 1 tags='Prolongation="None" checkpoint="no"' "Update term for psidc"
+real divB type = GF Timelevels = 1 tags='Prolongation="None" checkpoint="no"' "Magnetic field constraint"
+
##################################################
### These variables are only protected so that ###
### the tests in init_data work. Should fix. ###
File [modified]: param.ccl
Delta lines: +5 -1
===================================================================
--- trunk/param.ccl 2011-04-06 14:37:10 UTC (rev 227)
+++ trunk/param.ccl 2011-04-13 03:29:13 UTC (rev 228)
@@ -492,7 +492,7 @@
0:* :: "Some big value"
} 1.e100
-boolean clean_divergence "Should we advect tracers?"
+boolean clean_divergence "Use hyperbolic divergence cleaning"
{
} "no"
@@ -506,6 +506,10 @@
0:* :: "Any value, but one to 12 is preferred"
} 1.0
+boolean track_divB "Track the magnetic field constraint violations"
+{
+} "no"
+
############### temporary parameters to be removed once connected issues are fixed.
boolean con2prim_oct_hack "Disregard c2p failures in oct/rotsym90 boundary cells with xyz < 0" STEERABLE=ALWAYS
File [modified]: schedule.ccl
Delta lines: +24 -7
===================================================================
--- trunk/schedule.ccl 2011-04-06 14:37:10 UTC (rev 227)
+++ trunk/schedule.ccl 2011-04-13 03:29:13 UTC (rev 228)
@@ -68,10 +68,17 @@
STORAGE:densrhs
STORAGE:taurhs
STORAGE:srhs
-STORAGE:Bvecrhs
-if (clean_divergence)
+if(CCTK_Equals(Bvec_evolution_method,"GRHydro"))
{
- STORAGE:psidcrhs
+ STORAGE:Bvecrhs
+ if (clean_divergence)
+ {
+ STORAGE:psidcrhs
+ }
+ if (track_divB)
+ {
+ STORAGE:divB
+ }
}
STORAGE:GRHydro_eos_scalars
STORAGE:GRHydro_minima
@@ -97,12 +104,22 @@
{
} "GRHydro initial data group"
-if(CCTK_Equals(Bvec_evolution_method,"GRHydro") && clean_divergence)
+if(CCTK_Equals(Bvec_evolution_method,"GRHydro"))
{
- schedule GRHydro_InitDivergenceClean IN HydroBase_Initial
+ if (clean_divergence)
{
- LANG: Fortran
- } "Set psi for divergence cleaning initially to zero"
+ schedule GRHydro_InitDivergenceClean IN HydroBase_Initial
+ {
+ LANG: Fortran
+ } "Set psi for divergence cleaning initially to zero"
+ }
+ if (track_divB)
+ {
+ schedule GRHydro_DivBInit IN HydroBase_Initial
+ {
+ LANG: Fortran
+ } "Set divB initially to zero"
+ }
}
#################################################
More information about the Commits
mailing list