GRHydro: Changes to make GRHydro MHD work with nuclear/hot equation of state:

1) Switch to EOSOmni pointwise C2P routine and modify where necessary.
2) Modify Con2PrimM.F90 to allow for the evolution of temperature and adjust the wrapper routine.
3) Create EigenProblemM_hot pointwise routine and call that from HLLEM.F90 when temperature is evolved.
   Additionally adjust HLLEM where necessary.
4) Adjust InterfacesM.h to incorporate the newly created functions.
5) Fix a loop problem (not taking into account constraint transport) in PPM reconstruction of Y_e
6) Introduce Prim2ConM_hot and call this pointwise routine from Prim2ConM.F90 when temperature is evolved.
   Additionally also make this routine available to initial data routine in GRHydro_InitData
7) Adjust loops in GRHydro_PoloidalMagFieldM.F90 to not set boundary points it cannot set but instead call boundary group afterwards! Pay attention as this will not work with boundary conditions set to "none" in MHD case anymore but is the correct thing to do.
8) Allow StarMapper to extend HydroBase::initial_hydro = "starmapper".
9) Smaller fixes.

From: Philipp Moesta <pmoesta@tapir.caltech.edu>

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

1 files changed, 270 insertions, 1 deletions

diff --git a/src/GRHydro_Prim2ConM.F90 b/src/GRHydro_Prim2ConM.F90
index 59acf0d..a6aef7d 100644
--- a/src/GRHydro_Prim2ConM.F90
+++ b/src/GRHydro_Prim2ConM.F90
@@ -56,9 +56,29 @@ subroutine primitive2conservativeM(CCTK_ARGUMENTS)
   integer :: i, j, k
   CCTK_REAL :: gxxl,gxyl,gxzl,gyyl,gyzl,gzzl,avg_detl,&
        gxxr,gxyr,gxzr,gyyr,gyzr,gzzr,avg_detr
+  !CCTK_REAL :: g11l,g12l,g13l,g22l,g23l,g33l,&
+  !     g11r,g12r,g13r,g22r,g23r,g33r
+  CCTK_REAL :: xtemp(1)
+  character(len=256) NaN_WarnLine
+  !TARGET gxx, gxy, gxz, gyy, gyz, gzz
+
+  !CCTK_REAL, DIMENSION(:,:,:), POINTER :: g11, g12, g13, g22, g23, g33
   
+  !g11 => gxx
+  !g12 => gxy
+  !g13 => gxz
+  !g22 => gyy
+  !g23 => gyz
+  !g33 => gzz
+
   ! constraint transport needs to be able to average fluxes in the directions
   ! other that flux_direction, which in turn need the primitives on interfaces
+
+  if(evolve_temper.ne.1) then
+
+  !$OMP PARALLEL DO PRIVATE(k,j,i,avg_detl,avg_detr,&
+  !$OMP                      gxxl,gxyl,gxzl,gyyl,gyzl,gzzl, &
+  !$OMP                      gxxr,gxyr,gxzr,gyyr,gyzr,gzzr)
   do k = GRHydro_stencil, cctk_lsh(3)-GRHydro_stencil+1 + transport_constraints*(1-zoffset)
     do j = GRHydro_stencil, cctk_lsh(2)-GRHydro_stencil+1 + transport_constraints*(1-yoffset)
       do i = GRHydro_stencil, cctk_lsh(1)-GRHydro_stencil+1 + transport_constraints*(1-xoffset)
@@ -99,6 +119,79 @@ subroutine primitive2conservativeM(CCTK_ARGUMENTS)
       end do
     end do
   end do
+  !$OMP END PARALLEL DO
+
+  else
+
+     !$OMP PARALLEL DO PRIVATE(i, j, k, avg_detl, avg_detr, xtemp,&
+     !$OMP                      gxxl,gxyl,gxzl,gyyl,gyzl,gzzl, &
+     !$OMP                      gxxr,gxyr,gxzr,gyyr,gyzr,gzzr)
+     do k = GRHydro_stencil, cctk_lsh(3)-GRHydro_stencil+1 + transport_constraints*(1-zoffset)
+       do j = GRHydro_stencil, cctk_lsh(2)-GRHydro_stencil+1 + transport_constraints*(1-yoffset)
+         do i = GRHydro_stencil, cctk_lsh(1)-GRHydro_stencil+1 + transport_constraints*(1-xoffset)
+              
+              gxxl = 0.5d0 * (gxx(i,j,k) + gxx(i-xoffset,j-yoffset,k-zoffset))
+              gxyl = 0.5d0 * (gxy(i,j,k) + gxy(i-xoffset,j-yoffset,k-zoffset))
+              gxzl = 0.5d0 * (gxz(i,j,k) + gxz(i-xoffset,j-yoffset,k-zoffset))
+              gyyl = 0.5d0 * (gyy(i,j,k) + gyy(i-xoffset,j-yoffset,k-zoffset))
+              gyzl = 0.5d0 * (gyz(i,j,k) + gyz(i-xoffset,j-yoffset,k-zoffset))
+              gzzl = 0.5d0 * (gzz(i,j,k) + gzz(i-xoffset,j-yoffset,k-zoffset))
+              gxxr = 0.5d0 * (gxx(i,j,k) + gxx(i+xoffset,j+yoffset,k+zoffset))
+              gxyr = 0.5d0 * (gxy(i,j,k) + gxy(i+xoffset,j+yoffset,k+zoffset))
+              gxzr = 0.5d0 * (gxz(i,j,k) + gxz(i+xoffset,j+yoffset,k+zoffset))
+              gyyr = 0.5d0 * (gyy(i,j,k) + gyy(i+xoffset,j+yoffset,k+zoffset))
+              gyzr = 0.5d0 * (gyz(i,j,k) + gyz(i+xoffset,j+yoffset,k+zoffset))
+              gzzr = 0.5d0 * (gzz(i,j,k) + gzz(i+xoffset,j+yoffset,k+zoffset))
+              
+              avg_detl = SPATIAL_DETERMINANT(gxxl,gxyl,gxzl,gyyl,gyzl,gzzl)
+              avg_detr = SPATIAL_DETERMINANT(gxxr,gxyr,gxzr,gyyr,gyzr,gzzr)
+
+              ! we do not have plus/minus vars for temperature since
+              ! eps is reconstructed. Hence, we do not update the
+              ! variable 'temperature' in prim2con at the interfaces 
+              ! We will instead use an average temperature as an initial
+              ! guess.
+              xtemp(1) = 0.5d0*(temperature(i,j,k) + &
+                   temperature(i-xoffset,j-yoffset,k-zoffset))
+
+    !$OMP CRITICAL
+              if (y_e_minus(i,j,k) .le. 0.0d0 .or. y_e_plus(i,j,k) .le. 0.0d0) then
+                  write(NaN_WarnLine,'(a100,7g15.6)') '(y_e_minus,y_e_plus,x,y,z,rho)', y_e(i,j,k), y_e_minus(i,j,k), y_e_plus(i,j,k), x(i,j,k),y(i,j,k),z(i,j,k),rho(i,j,k)
+                  call CCTK_WARN(1, NaN_WarnLine)
+              endif 
+    !$OMP END CRITICAL
+              call prim2conM_hot(GRHydro_eos_handle, GRHydro_reflevel,&
+                   i,j,k,x(i,j,k),y(i,j,k),z(i,j,k), gxxl,gxyl,gxzl,gyyl,gyzl,gzzl, &
+                   avg_detl,densminus(i,j,k),sxminus(i,j,k),&
+                   syminus(i,j,k),szminus(i,j,k),tauminus(i,j,k),&
+                   Bconsxminus(i,j,k),Bconsyminus(i,j,k),Bconszminus(i,j,k), rhominus(i,j,k), &
+                   velxminus(i,j,k),velyminus(i,j,k),velzminus(i,j,k),&
+                   epsminus(i,j,k),pressminus(i,j,k),Bvecxminus(i,j,k), &
+                   Bvecyminus(i,j,k), Bveczminus(i,j,k), w_lorentzminus(i, j, k), xtemp, y_e_minus(i,j,k))
+              ! we do not have plus/minus vars for temperature since
+              ! eps is reconstructed. Hence, we do not update the
+              ! variable 'temperature' in prim2con at the interfaces
+              ! We will instead use an average temperature as an initial
+              ! guess.
+              xtemp(1) = 0.5d0*(temperature(i,j,k) + &
+                   temperature(i+xoffset,j+yoffset,k+zoffset))
+  
+              call prim2conM_hot(GRHydro_eos_handle, GRHydro_reflevel,&
+                   i,j,k,x(i,j,k),y(i,j,k),z(i,j,k), gxxr,gxyr,gxzr,gyyr,gyzr,gzzr, &
+                   avg_detr, densplus(i,j,k),sxplus(i,j,k),&
+                   syplus(i,j,k),szplus(i,j,k),tauplus(i,j,k),&
+                   Bconsxplus(i,j,k),Bconsyplus(i,j,k),Bconszplus(i,j,k),&
+                   rhoplus(i,j,k),velxplus(i,j,k),velyplus(i,j,k),&
+                   velzplus(i,j,k),epsplus(i,j,k),pressplus(i,j,k),&
+                   Bvecxplus(i,j,k), Bvecyplus(i,j,k), Bveczplus(i,j,k), &
+                   w_lorentzplus(i,j,k), xtemp, y_e_plus(i,j,k)) 
+             
+           end do
+        end do
+     end do
+     !$OMP END PARALLEL DO
+  endif
+
 
 end subroutine primitive2conservativeM
 
@@ -117,6 +210,152 @@ end subroutine primitive2conservativeM
 
 @@*/
 
+subroutine prim2conM_hot(handle, GRHydro_reflevel, ii, jj, kk, x, y, z, gxx, gxy, gxz, gyy, gyz, gzz, det, ddens, &
+     dsx, dsy, dsz, dtau , dBconsx, dBconsy, dBconsz, drho, dvelx, dvely, dvelz, deps, dpress, dBvcx, dBvcy, dBvcz, w, temp, ye) 
+  
+  implicit none
+
+  DECLARE_CCTK_PARAMETERS
+  DECLARE_CCTK_FUNCTIONS
+
+  CCTK_REAL :: gxx, gxy, gxz, gyy, gyz, gzz, det
+  CCTK_REAL, dimension(1) :: ddens, dsx, dsy, dsz, dtau, &
+       dBconsx, dBconsy, dBconsz, &
+       drho, dvelx, dvely, dvelz,&
+       deps, dpress, dBvcx, dBvcy, dBvcz, vlowx, vlowy, vlowz   
+  CCTK_REAL :: temp(1),ye(1), yein(1), x, y, z
+  CCTK_INT :: handle, GRHydro_reflevel, ii, jj, kk
+  CCTK_REAL :: w
+  CCTK_REAL, dimension(1) :: Bdotv,ab0,b2,blowx,blowy,blowz
+  character(len=256) NaN_WarnLine
+  character(len=512) warnline
+
+! begin EOS Omni vars
+  integer :: n, keytemp, anyerr, keyerr(1)
+ ! real*8  :: xpress(1),xeps(1),xtemp(1),xye(1)
+  real*8  :: temp0(1)
+  n = 1; keytemp = 0; anyerr = 0; keyerr(1) = 0
+  !xpress = 0.0d0; xeps = 0.0d0; xtemp = 0.0d0; xye = 0.0d0
+! end EOS Omni vars
+  
+  temp0 = temp
+  w = 1.d0 / sqrt(1.d0 - DOT2(dvelx(1),dvely(1),dvelz(1)))
+
+!!$ BEGIN: Check for NaN value
+  if (w .ne. w) then
+    !$OMP CRITICAL
+    write(NaN_WarnLine,'(a100,6g15.6)') 'NaN produced in sqrt(): (gxx,gxy,gxz,gyy,gyz,gzz)', gxx, gxy, gxz, gyy, gyz, gzz
+    call CCTK_WARN(1, NaN_WarnLine)
+    write(NaN_WarnLine,'(a100,3g15.6)') 'NaN produced in sqrt(): (dvelx,dvely,dvelz)', dvelx, dvely, dvelz
+    call CCTK_WARN(1, NaN_WarnLine)
+    write(NaN_WarnLine,'(a100,3g15.6)') 'NaN produced in sqrt(): (x,y,z)', x, y, z
+    call CCTK_WARN(1, NaN_WarnLine)
+    !write(NaN_WarnLine,'(a100,g15.6)') 'NaN produced in sqrt(): w', w
+    !call CCTK_WARN(GRHydro_NaN_verbose, NaN_WarnLine)
+    !$OMP END CRITICAL
+  endif
+!!$ END: Check for NaN value
+
+  ye = max(min(ye,GRHydro_Y_e_max),GRHydro_Y_e_min)
+ 
+  yein = ye 
+ 
+  call EOS_Omni_press(handle,keytemp,GRHydro_eos_rf_prec,n,&
+       drho,deps,temp,ye,dpress,keyerr,anyerr)  
+  ! error handling
+  if(anyerr.ne.0) then
+     if(GRHydro_reflevel.lt.GRHydro_c2p_warn_from_reflevel) then
+        ! in this case (coarse grid error that is hopefully restricted
+        ! away), we use the average temperature between cells and call
+        ! the EOS with keytemp=1
+        keytemp=1
+        call EOS_Omni_press(handle,keytemp,GRHydro_eos_rf_prec,n,&
+             drho,deps,temp,ye,dpress,keyerr,anyerr)  
+        keytemp=0
+        if(anyerr.ne.0) then
+           !$OMP CRITICAL
+           call CCTK_WARN(1,"EOS error in prim2con_hot: lev 2")
+           write(warnline,"(3i5,1P10E15.6)") ii,jj,kk,x,y,z
+           call CCTK_WARN(1,warnline)
+           write(warnline,"(1P10E15.6)") drho,deps,temp,ye,yein
+           call CCTK_WARN(1,warnline)
+           write(warnline,"(A7,i8)") "code: ",keyerr(1)
+           call CCTK_WARN(1,warnline)
+           write(warnline,"(A10,i5)") "reflevel: ", GRHydro_reflevel
+           call CCTK_WARN(1,warnline)
+           !$OMP END CRITICAL
+        endif
+     else
+        ! This is a way of recovering even on finer refinement levels:
+        ! Use the average temperature at the interface instead of the
+        ! reconstructed specific internal energy.
+        !$OMP CRITICAL
+        call CCTK_WARN(1,"EOS error in prim2con_hot: NOW using averaged temp!")
+        write(warnline,"(3i5,1P10E15.6)") ii,jj,kk,x,y,z
+        call CCTK_WARN(1,warnline)
+        write(warnline,"(1P10E15.6)") drho,deps,temp0,ye
+        call CCTK_WARN(1,warnline)
+        write(warnline,"(A7,i8)") "code: ",keyerr(1)
+        call CCTK_WARN(1,warnline)
+        write(warnline,"(A10,i5)") "reflevel: ", GRHydro_reflevel
+        call CCTK_WARN(1,warnline)
+        !$OMP END CRITICAL
+        keytemp=1
+        temp = temp0
+        call EOS_Omni_press(handle,keytemp,GRHydro_eos_rf_prec,n,&
+             drho,deps,temp,ye,dpress,keyerr,anyerr)  
+        keytemp=0
+        if(anyerr.ne.0) then
+           !$OMP CRITICAL
+           call CCTK_WARN(1,"EOS error in prim2con_hot")
+           write(warnline,"(3i5,1P10E15.6)") ii,jj,kk,x,y,z
+           call CCTK_WARN(1,warnline)
+           write(warnline,"(1P10E15.6)") drho,deps,temp,ye
+           call CCTK_WARN(1,warnline)
+           write(warnline,"(A7,i8)") "code: ",keyerr(1)
+           call CCTK_WARN(1,warnline)
+           write(warnline,"(A10,i5)") "reflevel: ", GRHydro_reflevel
+           call CCTK_WARN(1,warnline)
+           call CCTK_WARN(0,"Aborting!!!")
+           !$OMP END CRITICAL
+        endif
+     endif
+  endif
+
+  vlowx = gxx*dvelx + gxy*dvely + gxz*dvelz
+  vlowy = gxy*dvelx + gyy*dvely + gyz*dvelz
+  vlowz = gxz*dvelx + gyz*dvely + gzz*dvelz
+
+  Bdotv=DOT(dvelx,dvely,dvelz,dBvcx,dBvcy,dBvcz)
+  ab0=w*Bdotv
+  b2 = DOT2(dBvcx,dBvcy,dBvcz)/w**2+Bdotv**2
+  blowx = (gxx*dBvcx + gxy*dBvcy + gxz*dBvcz)/w + &
+       w*Bdotv*vlowx
+  blowy = (gxy*dBvcx + gyy*dBvcy + gyz*dBvcz)/w + &
+       w*Bdotv*vlowy
+  blowz = (gxz*dBvcx + gyz*dBvcy + gzz*dBvcz)/w + &
+       w*Bdotv*vlowz
+
+  ddens = sqrt(det) * drho * w 
+  dsx = sqrt(det) * ((drho*(1+deps)+dpress+b2)*w*w * vlowx - &
+       ab0*blowx)
+  dsy = sqrt(det) * ((drho*(1+deps)+dpress+b2)*w*w * vlowy - &
+       ab0*blowy)
+  dsz = sqrt(det) * ((drho*(1+deps)+dpress+b2)*w*w * vlowz - &
+       ab0*blowz)
+  dtau = sqrt(det) * ((drho*(1+deps)+dpress+b2)*w*w - dpress-b2/2.0-ab0**2) - ddens 
+
+  dBconsx = sqrt(det)*dBvcx
+  dBconsy = sqrt(det)*dBvcy
+  dBconsz = sqrt(det)*dBvcz
+
+  !!$OMP CRITICAL
+  !write(NaN_WarnLine,'(a100,3g15.6)') '(dens out, sqrt(det), w:)', ddens,sqrt(det),w
+  !call CCTK_WARN(1, NaN_WarnLine)
+  !!$OMP END CRITICAL
+end subroutine prim2conM_hot
+
+
 subroutine prim2conM(handle, gxx, gxy, gxz, gyy, gyz, gzz, det, ddens, &
      dsx, dsy, dsz, dtau , dBconsx, dBconsy, dBconsz, drho, dvelx, dvely, dvelz, deps, dpress, dBvcx, dBvcy, dBvcz, w) 
   
@@ -202,7 +441,6 @@ end subroutine prim2conM
 
 @@*/
 
-
 subroutine Primitive2ConservativeCellsM(CCTK_ARGUMENTS)
 
   implicit none
@@ -210,9 +448,13 @@ subroutine Primitive2ConservativeCellsM(CCTK_ARGUMENTS)
   DECLARE_CCTK_ARGUMENTS
   DECLARE_CCTK_PARAMETERS
   
+  CCTK_REAL :: xtemp(1)
   CCTK_INT :: i, j, k
   CCTK_REAL :: det
   
+  if(evolve_temper.ne.1) then
+
+  !$OMP PARALLEL DO PRIVATE(k,j,i,det)
   do k = GRHydro_stencil,cctk_lsh(3)-GRHydro_stencil+1
     do j = GRHydro_stencil,cctk_lsh(2)-GRHydro_stencil+1
       do i = GRHydro_stencil,cctk_lsh(1)-GRHydro_stencil+1
@@ -232,6 +474,31 @@ subroutine Primitive2ConservativeCellsM(CCTK_ARGUMENTS)
       end do
     end do
   end do
+  !$OMP END PARALLEL DO
+  else 
+  !$OMP PARALLEL DO PRIVATE(k,j,i,det)
+  do k = GRHydro_stencil,cctk_lsh(3)-GRHydro_stencil+1
+    do j = GRHydro_stencil,cctk_lsh(2)-GRHydro_stencil+1
+      do i = GRHydro_stencil,cctk_lsh(1)-GRHydro_stencil+1
+        
+         det = SPATIAL_DETERMINANT(gxx(i,j,k),gxy(i,j,k),gxz(i,j,k), \
+              gyy(i,j,k),gyz(i,j,k),gzz(i,j,k))
+
+         call prim2conM_hot(GRHydro_eos_handle,GRHydro_reflevel,&
+              i,j,k,x(i,j,k),y(i,j,k),z(i,j,k),&
+              gxx(i,j,k),gxy(i,j,k),gxz(i,j,k),&
+              gyy(i,j,k),gyz(i,j,k),gzz(i,j,k),&
+              det, dens(i,j,k),sx(i,j,k),sy(i,j,k),sz(i,j,k),&
+              tau(i,j,k),Bconsx(i,j,k),Bconsy(i,j,k),Bconsz(i,j,k),&
+              rho(i,j,k),velx(i,j,k),vely(i,j,k),velz(i,j,k),&
+              eps(i,j,k),press(i,j,k),Bvecx(i,j,k), &
+              Bvecy(i,j,k), Bvecz(i,j,k), w_lorentz(i,j,k), temperature(i,j,k), Y_e(i,j,k))
+
+      end do
+    end do
+  end do
+  !$OMP END PARALLEL DO
+  endif
 
   if(evolve_Y_e.ne.0) then
      !$OMP PARALLEL DO PRIVATE(i, j, k)
@@ -451,6 +718,7 @@ subroutine Primitive2ConservativePolyCellsM(CCTK_ARGUMENTS)
   CCTK_INT :: i, j, k
   CCTK_REAL :: det
   
+  !$OMP PARALLEL DO PRIVATE(k,j,i,det)
   do k = GRHydro_stencil,cctk_lsh(3)-GRHydro_stencil+1
     do j = GRHydro_stencil,cctk_lsh(2)-GRHydro_stencil+1
       do i = GRHydro_stencil,cctk_lsh(1)-GRHydro_stencil+1
@@ -469,6 +737,7 @@ subroutine Primitive2ConservativePolyCellsM(CCTK_ARGUMENTS)
       end do
     end do
   end do
+  !$OMP END PARALLEL DO
 
   if(evolve_Y_e.ne.0) then
      !$OMP PARALLEL DO PRIVATE(i, j, k)