RIT's divergence cleaning development.

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

1 files changed, 44 insertions, 3 deletions

diff --git a/src/GRHydro_RiemannSolveM.F90 b/src/GRHydro_RiemannSolveM.F90
index 46ab310..d0eb991 100644
--- a/src/GRHydro_RiemannSolveM.F90
+++ b/src/GRHydro_RiemannSolveM.F90
@@ -181,6 +181,8 @@ subroutine RiemannSolveGeneralM(CCTK_ARGUMENTS)
   CCTK_REAL :: pressstarp,pressstarm,rhoenth_p,rhoenth_m
   CCTK_REAL :: velxlowp,velxlowm,velylowp,velylowm,velzlowp,velzlowm
 
+  CCTK_REAL :: psidcp, psidcm, psidcf, psidcdiff, psidcfp, psidcfm
+
   densflux = 0.d0
   sxflux   = 0.d0
   syflux   = 0.d0
@@ -189,6 +191,9 @@ subroutine RiemannSolveGeneralM(CCTK_ARGUMENTS)
   Bvecxflux = 0.d0
   Bvecyflux = 0.d0
   Bveczflux = 0.d0
+  if(clean_divergence.ne.0) then
+     psidcflux = 0.d0
+  endif
 
 !!$  Do the EOS call to set the pressure, derivative and cs2
 
@@ -248,7 +253,12 @@ subroutine RiemannSolveGeneralM(CCTK_ARGUMENTS)
               prim_m(5) = epsminus(i+xoffset,j+yoffset,k+zoffset) 
               prim_m(6) = pressminus(i+xoffset,j+yoffset,k+zoffset) 
               rhoenth_m    = prim_m(1)*(1.0d0+prim_m(5))+prim_m(6)
-              
+
+              if(clean_divergence.ne.0) then
+                 psidcp = psidcplus(i,j,k)
+                 psidcm = psidcminus(i+xoffset,j+yoffset,k+zoffset)
+              endif
+             
 !!$        Set metric terms at interface
               
               if (shift_state .ne. 0) then
@@ -322,7 +332,7 @@ subroutine RiemannSolveGeneralM(CCTK_ARGUMENTS)
                       tmp_flux(6),tmp_flux(7),tmp_flux(8), &
                       vxtp,vytp,vztp,pressstarp,bxlowp,bylowp,bzlowp,ab0p,wp, &
                       avg_det,avg_alp,avg_beta)
-                 
+
                  densflux(i,j,k) = 0.5d0 * tmp_flux(1)
                  sxflux(i,j,k)   = 0.5d0 * tmp_flux(2)
                  syflux(i,j,k)   = 0.5d0 * tmp_flux(3)
@@ -332,13 +342,18 @@ subroutine RiemannSolveGeneralM(CCTK_ARGUMENTS)
                  Bvecyflux(i,j,k)   = 0.5d0 * tmp_flux(7)
                  Bveczflux(i,j,k)   = 0.5d0 * tmp_flux(8)
 
+                 if(clean_divergence.ne.0) then
+                    psidcf = cons_p(6)
+                    psidcflux(i,j,k) = 0.5d0 * psidcf
+                 endif
+                 
                  call num_x_fluxM(cons_m(1),cons_m(2),cons_m(3),cons_m(4),cons_m(5),&
                       cons_m(6),cons_m(7),cons_m(8),&
                       tmp_flux(1),tmp_flux(2),tmp_flux(3),tmp_flux(4),tmp_flux(5),&
                       tmp_flux(6),tmp_flux(7),tmp_flux(8),&
                       vxtm,vytm,vztm,pressstarm,bxlowm,bylowm,bzlowm,ab0m,wm, &
                       avg_det,avg_alp,avg_beta)
-                 
+
                  densflux(i,j,k) = densflux(i,j,k) + 0.5d0 * tmp_flux(1)
                  sxflux(i,j,k)   = sxflux(i,j,k)   + 0.5d0 * tmp_flux(2)
                  syflux(i,j,k)   = syflux(i,j,k)   + 0.5d0 * tmp_flux(3)
@@ -347,7 +362,13 @@ subroutine RiemannSolveGeneralM(CCTK_ARGUMENTS)
                  Bvecxflux(i,j,k)= Bvecxflux(i,j,k)+ 0.5d0 * tmp_flux(6)
                  Bvecyflux(i,j,k)= Bvecyflux(i,j,k)+ 0.5d0 * tmp_flux(7)
                  Bveczflux(i,j,k)= Bveczflux(i,j,k)+ 0.5d0 * tmp_flux(8)
+                 psidcflux(i,j,k)= psidcflux(i,j,k)+ 0.5d0 * psidcf
                  
+                 if(clean_divergence.ne.0) then
+                    psidcf = cons_m(6)
+                    psidcflux(i,j,k) = psidcflux(i,j,k) + 0.5d0 * psidcf
+                 endif
+
               else if (flux_direction == 2) then
                  
                  call num_x_fluxM(cons_p(1),cons_p(3),cons_p(4),cons_p(2),cons_p(5),&
@@ -366,6 +387,11 @@ subroutine RiemannSolveGeneralM(CCTK_ARGUMENTS)
                  Bvecyflux(i,j,k)= Bvecyflux(i,j,k)+ 0.5d0 * tmp_flux(7)
                  Bveczflux(i,j,k)= Bveczflux(i,j,k)+ 0.5d0 * tmp_flux(8)
                  
+                 if(clean_divergence.ne.0) then
+                    psidcf = cons_p(7)
+                    psidcflux(i,j,k) = psidcflux(i,j,k) + 0.5d0 * psidcf
+                 endif
+
                  call num_x_fluxM(cons_m(1),cons_m(3),cons_m(4),cons_m(2),cons_m(5),&
                       cons_m(7),cons_m(8),cons_m(6),&
                       tmp_flux(1),tmp_flux(3),tmp_flux(4),tmp_flux(2),tmp_flux(5),&
@@ -382,6 +408,11 @@ subroutine RiemannSolveGeneralM(CCTK_ARGUMENTS)
                  Bvecyflux(i,j,k)= Bvecyflux(i,j,k)+ 0.5d0 * tmp_flux(7)
                  Bveczflux(i,j,k)= Bveczflux(i,j,k)+ 0.5d0 * tmp_flux(8)
                  
+                 if(clean_divergence.ne.0) then
+                    psidcf = cons_m(7)
+                    psidcflux(i,j,k) = psidcflux(i,j,k) + 0.5d0 * psidcf
+                 endif
+
               else if (flux_direction == 3) then
                  
                  call num_x_fluxM(cons_p(1),cons_p(4),cons_p(2),cons_p(3),cons_p(5),&
@@ -400,6 +431,11 @@ subroutine RiemannSolveGeneralM(CCTK_ARGUMENTS)
                  Bvecyflux(i,j,k)= Bvecyflux(i,j,k)+ 0.5d0 * tmp_flux(7)
                  Bveczflux(i,j,k)= Bveczflux(i,j,k)+ 0.5d0 * tmp_flux(8)
                  
+                 if(clean_divergence.ne.0) then
+                    psidcf = cons_p(8)
+                    psidcflux(i,j,k) = psidcflux(i,j,k) + 0.5d0 * psidcf
+                 endif
+
                  call num_x_fluxM(cons_m(1),cons_m(4),cons_m(2),cons_m(3),cons_m(5),&
                       cons_m(8),cons_m(6),cons_m(7),&
                       tmp_flux(1),tmp_flux(4),tmp_flux(2),tmp_flux(3),tmp_flux(5), &
@@ -416,6 +452,11 @@ subroutine RiemannSolveGeneralM(CCTK_ARGUMENTS)
                  Bvecyflux(i,j,k)= Bvecyflux(i,j,k)+ 0.5d0 * tmp_flux(7)
                  Bveczflux(i,j,k)= Bveczflux(i,j,k)+ 0.5d0 * tmp_flux(8)
                  
+                 if(clean_divergence.ne.0) then
+                    psidcf = cons_m(8)
+                    psidcflux(i,j,k) = psidcflux(i,j,k) + 0.5d0 * psidcf
+                 endif
+
               else
                  
                  call CCTK_WARN(0, "Flux direction not x,y,z")