merge branch hot_and_MHD_temp_dev into branch at revision r185

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

1 files changed, 59 insertions, 3 deletions

diff --git a/src/GRHydro_Reconstruct.F90 b/src/GRHydro_Reconstruct.F90
index 3d4b843..9168147 100644
--- a/src/GRHydro_Reconstruct.F90
+++ b/src/GRHydro_Reconstruct.F90
@@ -131,6 +131,12 @@ subroutine Reconstruction(CCTK_ARGUMENTS)
      tracerminus = 0.0d0
   endif
 
+  if (evolve_Y_e .ne. 0) then
+     Y_e_plus = 0.0d0
+     Y_e_minus = 0.0d0
+  endif
+
+
   if (CCTK_EQUALS(recon_method,"tvd")) then
 
      if (evolve_tracer .ne. 0) then
@@ -142,9 +148,17 @@ subroutine Reconstruction(CCTK_ARGUMENTS)
         enddo
      end if
 
+     if (evolve_Y_e .ne. 0) then
+        call tvdreconstruct(nx, ny, nz, xoffset, yoffset, zoffset, &
+             Y_e(:,:,:), Y_e_plus(:,:,:), &
+             Y_e_minus(:,:,:), &
+             trivial_rp, hydro_excision_mask)
+     endif
+
     if (CCTK_EQUALS(recon_vars,"primitive")) then
       call tvdreconstruct(nx, ny, nz, xoffset, yoffset, zoffset, &
            rho, rhoplus, rhominus, trivial_rp, hydro_excision_mask)
+
       call tvdreconstruct(nx, ny, nz, xoffset, yoffset, zoffset, &
            vel(:,:,:,1), velxplus, velxminus, trivial_rp, hydro_excision_mask)
       call tvdreconstruct(nx, ny, nz, xoffset, yoffset, zoffset, &
@@ -234,6 +248,18 @@ subroutine Reconstruction(CCTK_ARGUMENTS)
          end do
          !$OMP END PARALLEL DO
       end if
+      if(evolve_Y_e.ne.0) then
+         !$OMP PARALLEL DO PRIVATE(j)
+         do k = GRHydro_stencil, nz - GRHydro_stencil + 1
+            do j = GRHydro_stencil, ny - GRHydro_stencil + 1
+               call SimplePPM_ye_1d(nx,CCTK_DELTA_SPACE(1),rho(:,j,k), &
+                    velx(:,j,k),vely(:,j,k),velz(:,j,k), &
+                    Y_e(:,j,k),Y_e_minus(:,j,k),Y_e_plus(:,j,k), &
+                    press(:,j,k))
+            end do
+         end do
+         !$OMP END PARALLEL DO
+      end if
 
     else if (flux_direction == 2) then
       !$OMP PARALLEL DO PRIVATE(i, j)
@@ -273,6 +299,18 @@ subroutine Reconstruction(CCTK_ARGUMENTS)
          end do
          !$OMP END PARALLEL DO
       end if
+      if(evolve_Y_e.ne.0) then
+         !$OMP PARALLEL DO PRIVATE(j)
+         do k = GRHydro_stencil, nz - GRHydro_stencil + 1
+            do j = GRHydro_stencil, nx - GRHydro_stencil + 1
+               call SimplePPM_ye_1d(ny,CCTK_DELTA_SPACE(2),rho(j,:,k), &
+                    velx(j,:,k),vely(j,:,k),velz(j,:,k), &
+                    Y_e(j,:,k),Y_e_minus(j,:,k),Y_e_plus(j,:,k), &
+                    press(j,:,k))
+            end do
+         end do
+         !$OMP END PARALLEL DO
+      end if
 
     else if (flux_direction == 3) then
       !$OMP PARALLEL DO PRIVATE(i, j)
@@ -313,6 +351,18 @@ subroutine Reconstruction(CCTK_ARGUMENTS)
          end do
          !$OMP END PARALLEL DO
       end if
+      if(evolve_Y_e.ne.0) then
+         !$OMP PARALLEL DO PRIVATE(j)
+         do k = GRHydro_stencil, ny - GRHydro_stencil + 1
+            do j = GRHydro_stencil, nx - GRHydro_stencil + 1
+               call SimplePPM_ye_1d(nz,CCTK_DELTA_SPACE(3),rho(j,k,:), &
+                    velx(j,k,:),vely(j,k,:),velz(j,k,:), &
+                    Y_e(j,k,:),Y_e_minus(j,k,:),Y_e_plus(j,k,:), &
+                    press(j,k,:))
+            end do
+         end do
+         !$OMP END PARALLEL DO
+      end if
 
     else
       call CCTK_WARN(0, "Flux direction not x,y,z")
@@ -328,6 +378,13 @@ subroutine Reconstruction(CCTK_ARGUMENTS)
         enddo
      end if
 
+     if (evolve_Y_e .ne. 0) then
+        call tvdreconstruct(nx, ny, nz, xoffset, yoffset, zoffset, &
+             Y_e(:,:,:), Y_e_plus(:,:,:), &
+             Y_e_minus(:,:,:), &
+             trivial_rp, hydro_excision_mask)
+     endif
+
     if (flux_direction == 1) then
       !$OMP PARALLEL DO PRIVATE(i, j)
       do k = GRHydro_stencil, cctk_lsh(3) - GRHydro_stencil + 1
@@ -492,8 +549,7 @@ subroutine Reconstruction(CCTK_ARGUMENTS)
   deallocate(psi4, lbetax, lbetay, lbetaz)
 
   !$OMP WORKSHARE
-  where ( (rhoplus < GRHydro_rho_min).or.(rhominus < GRHydro_rho_min).or.&
-          (epsplus < 0.d0).or.(epsminus < 0.d0) )
+  where ( (rhoplus < GRHydro_rho_min).or.(rhominus < GRHydro_rho_min) )
     rhoplus = rho
     rhominus = rho
     velxplus = vel(:,:,:,1)
@@ -511,7 +567,7 @@ subroutine Reconstruction(CCTK_ARGUMENTS)
     if (use_min_tracer .ne. 0) then
       local_min_tracer = min_tracer
     else
-      local_min_tracer = 0d0
+      local_min_tracer = 0.0d0
     end if
 
     !$OMP WORKSHARE