[Commits] [svn:einsteintoolkit] GRHydro/trunk/ (Rev. 635)

User: rhaas Date: 2014/04/15 02:49 PM

Modified: /trunk/ schedule.ccl /trunk/src/ GRHydro_Analysis.F90, GRHydro_Boundaries.F90, GRHydro_CalcUpdate.F90, GRHydro_UpdateMask.F90

Log: GRHydro: clean up handling of Avec evolution and boundaries

Remove almost empty Analysis_Init routine, apply boundary conditions consistently.

File Changes:

Directory: /trunk/src/ ======================

File [modified]: GRHydro_Analysis.F90 Delta lines: +70 -55 =================================================================== --- trunk/src/GRHydro_Analysis.F90 2014-04-15 19:49:52 UTC (rev 634) +++ trunk/src/GRHydro_Analysis.F90 2014-04-15 19:49:55 UTC (rev 635) @@ -5,18 +5,6 @@ #include "cctk_Functions.h" #include "SpaceMask.h"

-subroutine GRHydro_Analysis_Init(CCTK_ARGUMENTS) - - implicit none - - DECLARE_CCTK_ARGUMENTS - DECLARE_CCTK_PARAMETERS - DECLARE_CCTK_FUNCTIONS - - divB = 0.0d0 - -end subroutine GRHydro_Analysis_Init - subroutine GRHydro_CalcDivB(CCTK_ARGUMENTS)

implicit none @@ -27,56 +15,83 @@

CCTK_INT :: i,j,k,itracer CCTK_REAL :: idx, idy, idz, Bcons_l1, Bcons_r1, Bcons_l2, Bcons_r2, Bcons_l3, Bcons_r3 + logical :: do_Avec, do_Bcons

+ idx = 1.d0 / CCTK_DELTA_SPACE(1) idy = 1.d0 / CCTK_DELTA_SPACE(2) idz = 1.d0 / CCTK_DELTA_SPACE(3)

- !$OMP PARALLEL DO PRIVATE(i,j,k,itracer,Bcons_l1,Bcons_r1,Bcons_l2,Bcons_r2,Bcons_l3,Bcons_r3) - 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 + do_Avec = CCTK_EQUALS(Bvec_evolution_method,"GRHydro_Avec") + do_Bcons = CCTK_EQUALS(Bvec_evolution_method,"GRHydro")

- if(evolve_mhd.ne.0) then - if(track_divB.ne.0) then - if(transport_constraints.ne.0) then - ! edge based divergence (see WhiskyMHD & Bruno's thesis, Eq. 7.27) - divB(i,j,k) = & - 0.25d0*(Bcons(i+1,j,k,1)-Bcons(i,j,k,1)+ & - Bcons(i+1,j+1,k,1)-Bcons(i,j+1,k,1)+ & - Bcons(i+1,j,k+1,1)-Bcons(i,j,k+1,1)+ & - Bcons(i+1,j+1,k+1,1)-Bcons(i,j+1,k+1,1))*idx + & - 0.25d0*(Bcons(i,j+1,k,2)-Bcons(i,j,k,2)+ & - Bcons(i+1,j+1,k,2)-Bcons(i+1,j,k,2)+ & - Bcons(i,j+1,k+1,2)-Bcons(i,j,k+1,2)+ & - Bcons(i+1,j+1,k+1,2)-Bcons(i+1,j,k+1,2))*idy + & - 0.25d0*(Bcons(i,j,k+1,3)-Bcons(i,j,k,3)+ & - Bcons(i+1,j,k+1,3)-Bcons(i+1,j,k,3)+ & - Bcons(i,j+1,k+1,3)-Bcons(i,j+1,k,3)+ & - Bcons(i+1,j+1,k+1,3)-Bcons(i+1,j+1,k,3))*idz - else - Bcons_l1 = 0.5d0 * (Bcons(i,j,k,1) + & - Bcons(i-1,j,k,1)) - Bcons_l2 = 0.5d0 * (Bcons(i,j,k,2) + & - Bcons(i,j-1,k,2)) - Bcons_l3 = 0.5d0 * (Bcons(i,j,k,3) + & - Bcons(i,j,k-1,3)) - Bcons_r1 = 0.5d0 * (Bcons(i,j,k,1) + & - Bcons(i+1,j,k,1)) - Bcons_r2 = 0.5d0 * (Bcons(i,j,k,2) + & - Bcons(i,j+1,k,2)) - Bcons_r3 = 0.5d0 * (Bcons(i,j,k,3) + & - Bcons(i,j,k+1,3)) - - divB(i,j,k) = (Bcons_l1 - Bcons_r1) * idx + (Bcons_l2 - Bcons_r2) * idy + (Bcons_l3 - Bcons_r3) * idz + divB = 0d0

- endif - endif - endif - - enddo - enddo + !$OMP PARALLEL DO PRIVATE(i,j,k,itracer,Bcons_l1,Bcons_r1,Bcons_l2,Bcons_r2,Bcons_l3,Bcons_r3) + 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 + + if(evolve_mhd.ne.0) then + if(track_divB.ne.0) then + if(transport_constraints.ne.0) then + ! edge based divergence (see WhiskyMHD & Bruno's thesis, Eq. 7.27) + divB(i,j,k) = & + 0.25d0*(Bcons(i+1,j,k,1)-Bcons(i,j,k,1)+ & + Bcons(i+1,j+1,k,1)-Bcons(i,j+1,k,1)+ & + Bcons(i+1,j,k+1,1)-Bcons(i,j,k+1,1)+ & + Bcons(i+1,j+1,k+1,1)-Bcons(i,j+1,k+1,1))*idx + & + 0.25d0*(Bcons(i,j+1,k,2)-Bcons(i,j,k,2)+ & + Bcons(i+1,j+1,k,2)-Bcons(i+1,j,k,2)+ & + Bcons(i,j+1,k+1,2)-Bcons(i,j,k+1,2)+ & + Bcons(i+1,j+1,k+1,2)-Bcons(i+1,j,k+1,2))*idy + & + 0.25d0*(Bcons(i,j,k+1,3)-Bcons(i,j,k,3)+ & + Bcons(i+1,j,k+1,3)-Bcons(i+1,j,k,3)+ & + Bcons(i,j+1,k+1,3)-Bcons(i,j+1,k,3)+ & + Bcons(i+1,j+1,k+1,3)-Bcons(i+1,j+1,k,3))*idz + else if(do_Avec) then + Bcons_l1 = 0.5d0 * (sdetg(i,j,k)*Bvec(i,j,k,1) + & + sdetg(i-1,j,k)*Bvec(i-1,j,k,1)) + Bcons_l2 = 0.5d0 * (sdetg(i,j,k)*Bvec(i,j,k,2) + & + sdetg(i,j-1,k)*Bvec(i,j-1,k,2)) + Bcons_l3 = 0.5d0 * (sdetg(i,j,k)*Bvec(i,j,k,3) + & + sdetg(i,j,k-1)*Bvec(i,j,k-1,3)) + Bcons_r1 = 0.5d0 * (sdetg(i,j,k)*Bvec(i,j,k,1) + & + sdetg(i+1,j,k)*Bvec(i+1,j,k,1)) + Bcons_r2 = 0.5d0 * (sdetg(i,j,k)*Bvec(i,j,k,2) + & + sdetg(i,j+1,k)*Bvec(i,j+1,k,2)) + Bcons_r3 = 0.5d0 * (sdetg(i,j,k)*Bvec(i,j,k,3) + & + sdetg(i,j,k+1)*Bvec(i,j,k+1,3)) + + divB(i,j,k) = (Bcons_l1 - Bcons_r1) * idx + (Bcons_l2 - Bcons_r2) * idy + (Bcons_l3 - Bcons_r3) * idz + + else if(do_Bcons) then + Bcons_l1 = 0.5d0 * (Bcons(i,j,k,1) + & + Bcons(i-1,j,k,1)) + Bcons_l2 = 0.5d0 * (Bcons(i,j,k,2) + & + Bcons(i,j-1,k,2)) + Bcons_l3 = 0.5d0 * (Bcons(i,j,k,3) + & + Bcons(i,j,k-1,3)) + Bcons_r1 = 0.5d0 * (Bcons(i,j,k,1) + & + Bcons(i+1,j,k,1)) + Bcons_r2 = 0.5d0 * (Bcons(i,j,k,2) + & + Bcons(i,j+1,k,2)) + Bcons_r3 = 0.5d0 * (Bcons(i,j,k,3) + & + Bcons(i,j,k+1,3)) + + divB(i,j,k) = (Bcons_l1 - Bcons_r1) * idx + (Bcons_l2 - Bcons_r2) * idy + (Bcons_l3 - Bcons_r3) * idz + + else + call CCTK_ERROR("Internal error. Do now know how to compute divB") + STOP + + endif + endif + endif + enddo - !$OMP END PARALLEL DO + enddo + enddo + !$OMP END PARALLEL DO

end subroutine GRHydro_CalcDivB

File [modified]: GRHydro_Boundaries.F90 Delta lines: +86 -149 =================================================================== --- trunk/src/GRHydro_Boundaries.F90 2014-04-15 19:49:52 UTC (rev 634) +++ trunk/src/GRHydro_Boundaries.F90 2014-04-15 19:49:55 UTC (rev 635) @@ -165,199 +165,136 @@

@@*/

-subroutine GRHydro_Boundaries(CCTK_ARGUMENTS) - +subroutine GRHydro_RegisterBoundariesWithMethod(CCTK_ARGUMENTS, bcname) + implicit none - + DECLARE_CCTK_ARGUMENTS DECLARE_CCTK_PARAMETERS DECLARE_CCTK_FUNCTIONS - + + character(len=*), INTENT(IN) :: bcname + integer, dimension(3) :: sw integer :: ierr CCTK_INT :: GRHydro_UseGeneralCoordinates, general_coordinates

CCTK_INT, parameter :: faces=CCTK_ALL_FACES CCTK_INT, parameter :: ione=1 + logical :: do_Avec, do_Bcons

ierr = 0

sw = GRHydro_stencil

general_coordinates = GRHydro_UseGeneralCoordinates(cctkGH) + do_Avec = CCTK_EQUALS(Bvec_evolution_method,"GRHydro_Avec") + do_Bcons = CCTK_EQUALS(Bvec_evolution_method,"GRHydro")

- if (verbose.eq.1) call CCTK_INFO("Selecting conserved BC (and primitive BC if selected)") - -!!$Flat boundaries if required - - if (CCTK_EQUALS(bound,"flat")) then + ierr = ierr + Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & + "GRHydro::dens", bcname) + ierr = ierr + Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & + "GRHydro::tau", bcname) + ierr = ierr + Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & + "GRHydro::scon", bcname) + if (sync_conserved_only .eq. 0) then ierr = ierr + Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & - "GRHydro::dens", "Flat") + "HydroBase::w_lorentz", bcname) ierr = ierr + Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & - "GRHydro::tau", "Flat") + "HydroBase::rho", bcname) ierr = ierr + Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & - "GRHydro::scon", "Flat") - if (sync_conserved_only .eq. 0) then - ierr = ierr + Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & - "HydroBase::w_lorentz", "Flat") - ierr = ierr + Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & - "HydroBase::rho", "Flat") - ierr = ierr + Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & - "HydroBase::press", "Flat") - ierr = ierr + Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & - "HydroBase::eps", "Flat") - if (general_coordinates .ne. 0) then - ierr = ierr + Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & - "GRHydro::lvel", "Flat") - else - ierr = ierr + Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & - "HydroBase::vel", "Flat") - endif - endif - if(evolve_mhd.ne.0) then - if (sync_conserved_only .eq. 0) then - if (general_coordinates .ne. 0) then - ierr = ierr + Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & - "GRHydro::lBvec", "Flat") - else - ierr = ierr + Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & - "HydroBase::Bvec", "Flat") - endif - endif + "HydroBase::press", bcname) + ierr = ierr + Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & + "HydroBase::eps", bcname) + if (general_coordinates .ne. 0) then ierr = ierr + Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & - "GRHydro::Bcons", "Flat") - if(clean_divergence.ne.0) then + "GRHydro::lvel", bcname) + else + ierr = ierr + Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & + "HydroBase::vel", bcname) + endif + endif + if(evolve_mhd.ne.0) then + if (sync_conserved_only .eq. 0) then + if (general_coordinates .ne. 0) then ierr = ierr + Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & - "GRHydro::psidc", "Flat") + "GRHydro::lBvec", bcname) + else + ierr = ierr + Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & + "HydroBase::Bvec", bcname) endif - endif + endif + if(do_Bcons) then + ierr = ierr + Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & + "GRHydro::Bcons", bcname) + end if + if (do_Avec) then + ierr = ierr + Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & + "HydroBase::Avec", bcname) + endif + if(clean_divergence.ne.0) then + ierr = ierr + Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & + "GRHydro::psidc", bcname) + endif + endif

- if (CCTK_EQUALS(Bvec_evolution_method, "GRHydro_Avec")) then - ierr = ierr + Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & - "HydroBase::Avec", "Flat") - endif

- if(evolve_tracer.ne.0) then - if (sync_conserved_only .eq. 0) then - ierr = ierr + Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & - "GRHydro::GRHydro_tracers", "Flat") - endif - ierr = ierr + Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & - "GRHydro::GRHydro_cons_tracers", "Flat") + if(evolve_tracer.ne.0) then + if (sync_conserved_only .eq. 0) then + ierr = ierr + Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & + "GRHydro::GRHydro_tracers", bcname) endif + ierr = ierr + Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & + "GRHydro::GRHydro_cons_tracers", bcname) + endif

- if(evolve_y_e.ne.0) then - if (sync_conserved_only .eq. 0) then - ierr = ierr + Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & - "HydroBase::Y_e", "Flat") - endif - ierr = ierr + Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & - "GRHydro::Y_e_con", "Flat") + if(evolve_y_e.ne.0) then + if (sync_conserved_only .eq. 0) then + ierr = ierr + Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & + "HydroBase::Y_e", bcname) endif - - if(evolve_temper.ne.0) then - if (sync_conserved_only .eq. 0) then - ierr = ierr + Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & - "HydroBase::temperature", "Flat") - ierr = ierr + Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & - "HydroBase::entropy", "Flat") - endif - endif - - + ierr = ierr + Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & + "GRHydro::Y_e_con", bcname) endif

- if (CCTK_EQUALS(bound,"none")) then - ierr = ierr + Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & - "GRHydro::dens", "None") - ierr = ierr + Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & - "GRHydro::tau", "None") - ierr = ierr + Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & - "GRHydro::scon", "None") - + if(evolve_temper.ne.0) then if (sync_conserved_only .eq. 0) then - ierr = ierr + Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & - "HydroBase::w_lorentz", "None") - ierr = ierr + Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & - "HydroBase::rho", "None") - ierr = ierr + Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & - "HydroBase::press", "None") - ierr = ierr + Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & - "HydroBase::eps", "None") - if (general_coordinates .ne. 0) then - ierr = ierr + Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & - "GRHydro::lvel", "None") - else - ierr = ierr + Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & - "HydroBase::vel", "None") - endif - endif - - if(evolve_mhd.ne.0) then - if (sync_conserved_only .eq. 0) then - if (general_coordinates .ne. 0) then - ierr = ierr + Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & - "GRHydro::lBvec", "None") - else - ierr = ierr + Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & - "HydroBase::Bvec", "None") - endif - endif ierr = ierr + Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & - "GRHydro::Bcons", "None") - if(clean_divergence.ne.0) then - ierr = ierr + Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & - "GRHydro::psidc", "None") - endif + "HydroBase::temperature", bcname) + ierr = ierr + Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & + "HydroBase::entropy", bcname) endif + endif

- if (CCTK_EQUALS(Bvec_evolution_method, "GRHydro_Avec")) then - ierr = ierr + Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & - "HydroBase::Avec", "None") - if (general_coordinates .ne. 0) then - ierr = ierr + Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & - "GRHydro::lBvec", "Flat") - else - ierr = ierr + Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & - "HydroBase::Bvec", "Flat") - endif - endif + if (ierr < 0) call CCTK_WARN(0, "problems with applying the chosen boundary condition")

- if(evolve_tracer.ne.0) then - if (sync_conserved_only .eq. 0) then - ierr = ierr + Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & - "GRHydro::GRHydro_tracers", "None") - endif - ierr = ierr + Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & - "GRHydro::GRHydro_cons_tracers", "None") - endif +end subroutine

- if(evolve_y_e.ne.0) then - if (sync_conserved_only .eq. 0) then - ierr = ierr + Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & - "HydroBase::Y_e", "None") - endif - ierr = ierr + Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & - "GRHydro::Y_e_con", "None") - endif +subroutine GRHydro_Boundaries(CCTK_ARGUMENTS)

- if(evolve_temper.ne.0) then - if (sync_conserved_only .eq. 0) then - ierr = ierr + Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & - "HydroBase::temperature", "None") - ierr = ierr + Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & - "HydroBase::entropy", "None") - endif - endif + implicit none

+ DECLARE_CCTK_ARGUMENTS + DECLARE_CCTK_PARAMETERS + DECLARE_CCTK_FUNCTIONS + + if (verbose.eq.1) call CCTK_INFO("Selecting conserved BC (and primitive BC if selected)") + +!!$Flat boundaries if required + if (CCTK_EQUALS(bound,"flat")) then + call GRHydro_RegisterBoundariesWithMethod(CCTK_PASS_FTOF, "Flat") endif

+!!$None boundaries if required + if (CCTK_EQUALS(bound,"none")) then + call GRHydro_RegisterBoundariesWithMethod(CCTK_PASS_FTOF, "None") + endif + +!!$Scalar boundaries if required if (CCTK_EQUALS(bound,"scalar")) then call CCTK_WARN(0, "Until somebody uses this I see no reason to support it") end if

- if (ierr < 0) call CCTK_WARN(0, "problems with applying the chosen boundary condition") - end subroutine GRHydro_Boundaries

File [modified]: GRHydro_CalcUpdate.F90 Delta lines: +16 -17 =================================================================== --- trunk/src/GRHydro_CalcUpdate.F90 2014-04-15 19:49:52 UTC (rev 634) +++ trunk/src/GRHydro_CalcUpdate.F90 2014-04-15 19:49:55 UTC (rev 635) @@ -108,24 +108,23 @@ 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.1) then - Bconsrhs(i,j,k,1) = Bconsrhs(i,j,k,1) + & - (alp_l * Bconsxflux(i-xoffset,j-yoffset,k-zoffset) - & - alp_r * Bconsxflux(i,j,k)) * idx - Bconsrhs(i,j,k,2) = Bconsrhs(i,j,k,2) + & - (alp_l * Bconsyflux(i-xoffset,j-yoffset,k-zoffset) - & - alp_r * Bconsyflux(i,j,k)) * idx - Bconsrhs(i,j,k,3) = Bconsrhs(i,j,k,3) + & - (alp_l * Bconszflux(i-xoffset,j-yoffset,k-zoffset) - & - alp_r * Bconszflux(i,j,k)) * idx - endif - 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(transport_constraints.eq.0) then + Bconsrhs(i,j,k,1) = Bconsrhs(i,j,k,1) + & + (alp_l * Bconsxflux(i-xoffset,j-yoffset,k-zoffset) - & + alp_r * Bconsxflux(i,j,k)) * idx + Bconsrhs(i,j,k,2) = Bconsrhs(i,j,k,2) + & + (alp_l * Bconsyflux(i-xoffset,j-yoffset,k-zoffset) - & + alp_r * Bconsyflux(i,j,k)) * idx + Bconsrhs(i,j,k,3) = Bconsrhs(i,j,k,3) + & + (alp_l * Bconszflux(i-xoffset,j-yoffset,k-zoffset) - & + alp_r * Bconszflux(i,j,k)) * idx + endif + 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 endif

if (evolve_tracer .ne. 0) then

File [modified]: GRHydro_UpdateMask.F90 Delta lines: +6 -0 =================================================================== --- trunk/src/GRHydro_UpdateMask.F90 2014-04-15 19:49:52 UTC (rev 634) +++ trunk/src/GRHydro_UpdateMask.F90 2014-04-15 19:49:55 UTC (rev 635) @@ -65,6 +65,9 @@ if (evolve_Bvec) then Bconsrhs(i,j,k,:) = 0.0d0 endif + if (evolve_Avec) then + Avecrhs(i,j,k,:) = 0.0d0 + endif

! TODO: Need to set Avecrhs and Aphirhs to zero for centered / or staggered case if vector potential is evolved

@@ -93,6 +96,9 @@ if (evolve_Bvec) then Bconsrhs(i,j,k,:) = 0.0d0 endif + if (evolve_Avec) then + Avecrhs(i,j,k,:) = 0.0d0 + endif

! TODO: Need to set Avecrhs and Aphirhs to zero for centered / or staggered case if vector potential is evolved

Directory: /trunk/ ==================

File [modified]: schedule.ccl Delta lines: +1 -7 =================================================================== --- trunk/schedule.ccl 2014-04-15 19:49:52 UTC (rev 634) +++ trunk/schedule.ccl 2014-04-15 19:49:55 UTC (rev 635) @@ -80,7 +80,6 @@ } else if(CCTK_Equals(Bvec_evolution_method,"GRHydro_Avec")) { STORAGE: HydroBase::Bvec[timelevels] STORAGE: HydroBase::Avec[timelevels] - STORAGE: GRHydro::Bcons[timelevels] if (CCTK_Equals(Avec_gauge,"Lorenz")) { STORAGE: HydroBase::Aphi[timelevels] @@ -1611,16 +1610,11 @@ { } "Calculate analysis quantities"

- schedule GRHydro_Analysis_Init IN GRHydroAnalysis - { - LANG: Fortran - } "Initialize variables" - schedule GRHydro_CalcDivB IN GRHydroAnalysis AFTER GRHydro_Analysis_Init { LANG: Fortran } "Calculate divB" -} +}

if (constrain_to_1D) {