MERGE divergence cleaning feature into trunk.

git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinEvolve/GRHydro/trunk@244 c83d129a-5a75-4d5a-9c4d-ed3a5855bf45

1 files changed, 103 insertions, 54 deletions

diff --git a/src/GRHydro_Reconstruct.F90 b/src/GRHydro_Reconstruct.F90
index 065d2bd..cf84024 100644
--- a/src/GRHydro_Reconstruct.F90
+++ b/src/GRHydro_Reconstruct.F90
@@ -23,6 +23,9 @@
 #define Bvecx(i,j,k) Bvec(i,j,k,1)
 #define Bvecy(i,j,k) Bvec(i,j,k,2)
 #define Bvecz(i,j,k) Bvec(i,j,k,3)
+#define Bconsx(i,j,k) Bcons(i,j,k,1)
+#define Bconsy(i,j,k) Bcons(i,j,k,2)
+#define Bconsz(i,j,k) Bcons(i,j,k,3)
 
 
  /*@@
@@ -208,6 +211,17 @@ subroutine Reconstruction(CCTK_ARGUMENTS)
            vel(:,:,:,3), velzplus, velzminus, trivial_rp, hydro_excision_mask)
       call tvdreconstruct(nx, ny, nz, xoffset, yoffset, zoffset, &
            eps, epsplus, epsminus, trivial_rp, hydro_excision_mask)
+      if(evolve_mhd.ne.0) then
+         call tvdreconstruct(nx, ny, nz, xoffset, yoffset, zoffset, &
+              Bvec(:,:,:,1), Bvecxplus, Bvecxminus, trivial_rp, hydro_excision_mask)
+         call tvdreconstruct(nx, ny, nz, xoffset, yoffset, zoffset, &
+              Bvec(:,:,:,2), Bvecyplus, Bvecyminus, trivial_rp, hydro_excision_mask)
+         call tvdreconstruct(nx, ny, nz, xoffset, yoffset, zoffset, &
+              Bvec(:,:,:,3), Bveczplus, Bveczminus, trivial_rp, hydro_excision_mask)
+      endif
+
+!!$      write(6,*)'recon0:',Bvecxplus(75,5,5),Bvecxminus(75,5,5)
+
     else if (CCTK_EQUALS(recon_vars,"conservative")) then
       call tvdreconstruct(nx, ny, nz, xoffset, yoffset, zoffset, &
            dens, densplus, densminus, trivial_rp, hydro_excision_mask)
@@ -219,22 +233,22 @@ subroutine Reconstruction(CCTK_ARGUMENTS)
            scon(:,:,:,3), szplus, szminus, trivial_rp, hydro_excision_mask)
       call tvdreconstruct(nx, ny, nz, xoffset, yoffset, zoffset, &
            tau, tauplus, tauminus, trivial_rp, hydro_excision_mask)
+      if(evolve_mhd.ne.0) then
+         call tvdreconstruct(nx, ny, nz, xoffset, yoffset, zoffset, &
+              Bcons(:,:,:,1), Bconsxplus, Bconsxminus, trivial_rp, hydro_excision_mask)
+         call tvdreconstruct(nx, ny, nz, xoffset, yoffset, zoffset, &
+              Bcons(:,:,:,2), Bconsyplus, Bconsyminus, trivial_rp, hydro_excision_mask)
+         call tvdreconstruct(nx, ny, nz, xoffset, yoffset, zoffset, &
+              Bcons(:,:,:,3), Bconszplus, Bconszminus, trivial_rp, hydro_excision_mask)
+      endif
+
     else
       call CCTK_WARN(0, "Variable type to reconstruct not recognized.")
     end if
 
-!!$ B-field is both prim and con
-    if(evolve_mhd.ne.0) then
+      if(evolve_mhd.ne.0.and.clean_divergence.ne.0) then
        call tvdreconstruct(nx, ny, nz, xoffset, yoffset, zoffset, &
-            Bvec(:,:,:,1), Bvecxplus, Bvecxminus, trivial_rp, hydro_excision_mask)
-       call tvdreconstruct(nx, ny, nz, xoffset, yoffset, zoffset, &
-            Bvec(:,:,:,2), Bvecyplus, Bvecyminus, trivial_rp, hydro_excision_mask)
-       call tvdreconstruct(nx, ny, nz, xoffset, yoffset, zoffset, &
-            Bvec(:,:,:,3), Bveczplus, Bveczminus, trivial_rp, hydro_excision_mask)
-       if(clean_divergence.ne.0) then
-          call tvdreconstruct(nx, ny, nz, xoffset, yoffset, zoffset, &
-               psidc, psidcplus, psidcminus, trivial_rp, hydro_excision_mask)
-       endif
+            psidc, psidcplus, psidcminus, trivial_rp, hydro_excision_mask)
     endif
 
     !$OMP PARALLEL DO PRIVATE(i, j)
@@ -640,6 +654,17 @@ subroutine Reconstruction(CCTK_ARGUMENTS)
             call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(1),&
                  eps(:,j,k),epsminus(:,j,k),epsplus(:,j,k),&
                  trivial_rp(:,j,k), hydro_excision_mask(:,j,k))
+            if(evolve_mhd.ne.0) then
+               call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(1),&
+                    Bvecx(:,j,k),Bvecxminus(:,j,k),Bvecxplus(:,j,k),&
+                    trivial_rp(:,j,k), hydro_excision_mask(:,j,k))
+               call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(1),&
+                    Bvecy(:,j,k),Bvecyminus(:,j,k),Bvecyplus(:,j,k),&
+                    trivial_rp(:,j,k), hydro_excision_mask(:,j,k))
+               call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(1),&
+                    Bvecz(:,j,k),Bveczminus(:,j,k),Bveczplus(:,j,k),&
+                    trivial_rp(:,j,k), hydro_excision_mask(:,j,k))
+            endif
           else if (CCTK_EQUALS(recon_vars,"conservative")) then
             call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(1),&
                  dens(:,j,k),densminus(:,j,k),densplus(:,j,k),&
@@ -656,28 +681,27 @@ subroutine Reconstruction(CCTK_ARGUMENTS)
             call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(1),&
                  tau(:,j,k),tauminus(:,j,k),tauplus(:,j,k),&
                  trivial_rp(:,j,k), hydro_excision_mask(:,j,k))
+            if(evolve_mhd.ne.0) then
+               call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(1),&
+                    Bconsx(:,j,k),Bconsxminus(:,j,k),Bconsxplus(:,j,k),&
+                    trivial_rp(:,j,k), hydro_excision_mask(:,j,k))
+               call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(1),&
+                    Bconsy(:,j,k),Bconsyminus(:,j,k),Bconsyplus(:,j,k),&
+                    trivial_rp(:,j,k), hydro_excision_mask(:,j,k))
+               call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(1),&
+                    Bconsz(:,j,k),Bconszminus(:,j,k),Bconszplus(:,j,k),&
+                    trivial_rp(:,j,k), hydro_excision_mask(:,j,k))
+            endif
           else
             !$OMP CRITICAL
             call CCTK_WARN(0, "Variable type to reconstruct not recognized.")
             !$OMP END CRITICAL
           end if
 
-!!$ B-fields are both prim and con
-          if(evolve_mhd.ne.0) then
-             call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(1),&
-                  Bvecx(:,j,k),Bvecxminus(:,j,k),Bvecxplus(:,j,k),&
-                  trivial_rp(:,j,k), hydro_excision_mask(:,j,k))
-             call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(1),&
-                  Bvecy(:,j,k),Bvecyminus(:,j,k),Bvecyplus(:,j,k),&
-                  trivial_rp(:,j,k), hydro_excision_mask(:,j,k))
+          if(evolve_mhd.ne.0.and.clean_divergence.ne.0) then
              call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(1),&
-                  Bvecz(:,j,k),Bveczminus(:,j,k),Bveczplus(:,j,k),&
+                  psidc(:,j,k),psidcminus(:,j,k),psidcplus(:,j,k),&
                   trivial_rp(:,j,k), hydro_excision_mask(:,j,k))
-             if(clean_divergence.ne.0) then
-                call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(1),&
-                     psidc(:,j,k),psidcminus(:,j,k),psidcplus(:,j,k),&
-                     trivial_rp(:,j,k), hydro_excision_mask(:,j,k))
-             endif
           endif
 
           do i = 1, cctk_lsh(1)
@@ -710,7 +734,18 @@ subroutine Reconstruction(CCTK_ARGUMENTS)
             call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(2),&
                  eps(j,:,k),epsminus(j,:,k),epsplus(j,:,k),&
                  trivial_rp(j,:,k), hydro_excision_mask(j,:,k))
-          else if (CCTK_EQUALS(recon_vars,"conservative")) then
+            if(evolve_mhd.ne.0) then
+               call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(2),&
+                    Bvecx(j,:,k),Bvecxminus(j,:,k),Bvecxplus(j,:,k),&
+                    trivial_rp(j,:,k), hydro_excision_mask(j,:,k))
+               call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(2),&
+                    Bvecy(j,:,k),Bvecyminus(j,:,k),Bvecyplus(j,:,k),&
+                    trivial_rp(j,:,k), hydro_excision_mask(j,:,k))
+               call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(2),&
+                    Bvecz(j,:,k),Bveczminus(j,:,k),Bveczplus(j,:,k),&
+                    trivial_rp(j,:,k), hydro_excision_mask(j,:,k))
+            endif
+         else if (CCTK_EQUALS(recon_vars,"conservative")) then
             call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(2),&
                  dens(j,:,k),densminus(j,:,k),densplus(j,:,k),&
                  trivial_rp(j,:,k), hydro_excision_mask(j,:,k))
@@ -726,27 +761,27 @@ subroutine Reconstruction(CCTK_ARGUMENTS)
             call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(2),&
                  tau(j,:,k),tauminus(j,:,k),tauplus(j,:,k),&
                  trivial_rp(j,:,k), hydro_excision_mask(j,:,k))
+            if(evolve_mhd.ne.0) then
+               call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(2),&
+                    Bconsx(j,:,k),Bconsxminus(j,:,k),Bconsxplus(j,:,k),&
+                    trivial_rp(j,:,k), hydro_excision_mask(j,:,k))
+               call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(2),&
+                    Bconsy(j,:,k),Bconsyminus(j,:,k),Bconsyplus(j,:,k),&
+                    trivial_rp(j,:,k), hydro_excision_mask(j,:,k))
+               call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(2),&
+                    Bconsz(j,:,k),Bconszminus(j,:,k),Bconszplus(j,:,k),&
+                    trivial_rp(j,:,k), hydro_excision_mask(j,:,k))
+            endif
           else
             !$OMP CRITICAL
             call CCTK_WARN(0, "Variable type to reconstruct not recognized.")
             !$OMP END CRITICAL
           end if
 
-          if(evolve_mhd.ne.0) then
-             call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(2),&
-                  Bvecx(j,:,k),Bvecxminus(j,:,k),Bvecxplus(j,:,k),&
-                  trivial_rp(j,:,k), hydro_excision_mask(j,:,k))
-             call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(2),&
-                  Bvecy(j,:,k),Bvecyminus(j,:,k),Bvecyplus(j,:,k),&
-                  trivial_rp(j,:,k), hydro_excision_mask(j,:,k))
-             call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(2),&
-                  Bvecz(j,:,k),Bveczminus(j,:,k),Bveczplus(j,:,k),&
-                  trivial_rp(j,:,k), hydro_excision_mask(j,:,k))
-             if(clean_divergence.ne.0) then
-                call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(2),&
-                     psidc(j,:,k),psidcminus(j,:,k),psidcplus(j,:,k),&
-                     trivial_rp(j,:,k), hydro_excision_mask(j,:,k))
-             endif
+          if(evolve_mhd.ne.0.and.clean_divergence.ne.0) then
+            call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(2),&
+                 psidc(j,:,k),psidcminus(j,:,k),psidcplus(j,:,k),&
+                 trivial_rp(j,:,k), hydro_excision_mask(j,:,k))
           endif
 
           do i = 1, cctk_lsh(2)
@@ -779,6 +814,17 @@ subroutine Reconstruction(CCTK_ARGUMENTS)
             call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(3),&
                  eps(j,k,:),epsminus(j,k,:),epsplus(j,k,:),&
                  trivial_rp(j,k,:), hydro_excision_mask(j,k,:))
+            if(evolve_mhd.ne.0) then
+               call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(3),&
+                    Bvecx(j,k,:),Bvecxminus(j,k,:),Bvecxplus(j,k,:),&
+                    trivial_rp(j,k,:), hydro_excision_mask(j,k,:))
+               call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(3),&
+                    Bvecy(j,k,:),Bvecyminus(j,k,:),Bvecyplus(j,k,:),&
+                    trivial_rp(j,k,:), hydro_excision_mask(j,k,:))
+               call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(3),&
+                    Bvecz(j,k,:),Bveczminus(j,k,:),Bveczplus(j,k,:),&
+                    trivial_rp(j,k,:), hydro_excision_mask(j,k,:))
+            endif
           else if (CCTK_EQUALS(recon_vars,"conservative")) then
             call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(3),&
                  dens(j,k,:),densminus(j,k,:),densplus(j,k,:),&
@@ -795,27 +841,27 @@ subroutine Reconstruction(CCTK_ARGUMENTS)
             call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(3),&
                  tau(j,k,:),tauminus(j,k,:),tauplus(j,k,:),&
                  trivial_rp(j,k,:), hydro_excision_mask(j,k,:))
+            if(evolve_mhd.ne.0) then
+               call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(3),&
+                    Bconsx(j,k,:),Bconsxminus(j,k,:),Bconsxplus(j,k,:),&
+                    trivial_rp(j,k,:), hydro_excision_mask(j,k,:))
+               call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(3),&
+                    Bconsy(j,k,:),Bconsyminus(j,k,:),Bconsyplus(j,k,:),&
+                    trivial_rp(j,k,:), hydro_excision_mask(j,k,:))
+               call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(3),&
+                    Bconsz(j,k,:),Bconszminus(j,k,:),Bconszplus(j,k,:),&
+                    trivial_rp(j,k,:), hydro_excision_mask(j,k,:))
+            endif
           else
             !$OMP CRITICAL
             call CCTK_WARN(0, "Variable type to reconstruct not recognized.")
             !$OMP END CRITICAL
           end if
 
-          if(evolve_mhd.ne.0) then
-             call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(3),&
-                  Bvecx(j,k,:),Bvecxminus(j,k,:),Bvecxplus(j,k,:),&
-                  trivial_rp(j,k,:), hydro_excision_mask(j,k,:))
-             call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(3),&
-                  Bvecy(j,k,:),Bvecyminus(j,k,:),Bvecyplus(j,k,:),&
-                  trivial_rp(j,k,:), hydro_excision_mask(j,k,:))
+          if(evolve_mhd.ne.0.and.clean_divergence.ne.0) then
              call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(3),&
-                  Bvecz(j,k,:),Bveczminus(j,k,:),Bveczplus(j,k,:),&
+                  psidc(j,k,:),psidcminus(j,k,:),psidcplus(j,k,:),&
                   trivial_rp(j,k,:), hydro_excision_mask(j,k,:))
-             if(clean_divergence.ne.0) then
-                call GRHydro_ENOReconstruct1d(eno_order,cctk_lsh(3),&
-                     psidc(j,k,:),psidcminus(j,k,:),psidcplus(j,k,:),&
-                     trivial_rp(j,k,:), hydro_excision_mask(j,k,:))
-             endif
           endif
 
           do i = 1, cctk_lsh(3)
@@ -877,6 +923,9 @@ subroutine Reconstruction(CCTK_ARGUMENTS)
        CCTK_EQUALS(recon_method,"ppm")) then
      if(evolve_mhd.ne.0) then
         call primitive2conservativeM(CCTK_PASS_FTOF)
+
+!!$      write(6,*)'recon1:',Bvecxplus(75,5,5),Bvecxminus(75,5,5),Bconsxplus(75,5,5),Bconsxminus(75,5,5)
+
      else
         call primitive2conservative(CCTK_PASS_FTOF)
      endif