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" + } }

#################################################