[svn:einsteintoolkit] GRHydro/trunk/src/ (Rev. 249) - Commits

11 Jun 2011

User: cott
Date: 2011/06/11 04:51 PM
Modified:
 /trunk/src/
  GRHydro_Con2Prim.F90, GRHydro_Prim2Con.F90, GRHydro_UpdateMask.F90
Log:
 * improve atmosphere handling so that it is now possible to use
   hot EOS with an atmosphere
File Changes:
Directory: /trunk/src/
======================
File [modified]: GRHydro_Con2Prim.F90
Delta lines: +26 -11
===================================================================

--- trunk/src/GRHydro_Con2Prim.F90	2011-05-23 22:31:18 UTC (rev 248)
+++ trunk/src/GRHydro_Con2Prim.F90	2011-06-11 21:51:15 UTC (rev 249)
@@ -88,7 +88,7 @@
        GRHydro_rho_min,xeps,xtemp,xye,pmin,epsmin,keyerr,anyerr)
!$OMP PARALLEL DO PRIVATE(i,j,itracer,&
-  !$OMP uxx, uxy, uxz, uyy, uyz, uzz, det, epsnegative)
+  !$OMP uxx, uxy, uxz, uyy, uyz, uzz, det, epsnegative, anyerr, keyerr, keytemp)
   do k = 1, nz 
     do j = 1, ny 
       do i = 1, nx
@@ -128,7 +128,6 @@
            Y_e(i,j,k) = Y_e_con(i,j,k) / dens(i,j,k)
         endif
-
          if ( dens(i,j,k) .le. sqrt(det)*GRHydro_rho_min*(1.d0+GRHydro_atmo_tolerance) ) then
dens(i,j,k) = sqrt(det)*GRHydro_rho_min !/(1.d0+GRHydro_atmo_tolerance)
@@ -137,15 +136,26 @@
            vel(i,j,k,:) = 0.d0
            w_lorentz(i,j,k) = 1.d0
-           call EOS_Omni_press(GRHydro_polytrope_handle,keytemp,GRHydro_eos_rf_prec,n,&
-                rho(i,j,k),eps(i,j,k),xtemp,xye,press(i,j,k),keyerr,anyerr)
+           if(evolve_temper.ne.0) then
+              ! set hot atmosphere values
+              temperature(i,j,k) = grhydro_hot_atmo_temp
+              y_e(i,j,k) = grhydro_hot_atmo_Y_e
+              y_e_con(i,j,k) = y_e(i,j,k) * dens(i,j,k)
+              keytemp = 1
+              call EOS_Omni_press(GRHydro_eos_handle,keytemp,GRHydro_eos_rf_prec,n,&
+                   rho(i,j,k),eps(i,j,k),temperature(i,j,k),y_e(i,j,k),&
+                   press(i,j,k),keyerr,anyerr)
+           else
+              ! use polytropic EOS
+              call EOS_Omni_press(GRHydro_polytrope_handle,keytemp,GRHydro_eos_rf_prec,n,&
+                   rho(i,j,k),eps(i,j,k),xtemp,xye,press(i,j,k),keyerr,anyerr)
-           call EOS_Omni_EpsFromPress(GRHydro_polytrope_handle,keytemp,GRHydro_eos_rf_prec,n,&
-                rho(i,j,k),eps(i,j,k),xtemp,xye,press(i,j,k),eps(i,j,k),keyerr,anyerr)
+              call EOS_Omni_EpsFromPress(GRHydro_polytrope_handle,keytemp,GRHydro_eos_rf_prec,n,&
+                   rho(i,j,k),eps(i,j,k),xtemp,xye,press(i,j,k),eps(i,j,k),keyerr,anyerr)
+           endif
-           ! w_lorentz=1, so the expression for tau reduces to:
            tau(i,j,k)  = sqrt(det) * (rho(i,j,k)+rho(i,j,k)*eps(i,j,k)) - dens(i,j,k)
-
+           
            cycle
end if
@@ -196,6 +206,7 @@
             endif
          endif
+
         if (epsnegative) then
#if 0
@@ -848,10 +859,14 @@
     rho = GRHydro_rho_min
     udens = rho
     dens = sqrt(det) * rho
-!    epsilon = (sqrt( (utau + pnew + udens)**2) - pnew -  udens)/udens
-    epsilon = epsminl
+    temp = GRHydro_hot_atmo_temp
+    ye = GRHydro_hot_atmo_Y_e
+    keytemp=1
+    call EOS_Omni_press(handle,keytemp,GRHydro_eos_rf_prec,n, &
+         rho,epsilon,temp,ye,xpress,keyerr,anyerr)
+    keytemp=0
     ! w_lorentz=1, so the expression for utau reduces to:
-    utau  = rho + rho*epsminl - udens
+    utau  = rho + rho*epsilon - udens
     sx = 0.d0
     sy = 0.d0
     sz = 0.d0
File [modified]: GRHydro_Prim2Con.F90
Delta lines: +0 -1
===================================================================
--- trunk/src/GRHydro_Prim2Con.F90	2011-05-23 22:31:18 UTC (rev 248)
+++ trunk/src/GRHydro_Prim2Con.F90	2011-06-11 21:51:15 UTC (rev 249)
@@ -303,7 +303,6 @@
      endif
   endif
-
   vlowx = gxx*dvelx + gxy*dvely + gxz*dvelz
   vlowy = gxy*dvelx + gyy*dvely + gyz*dvelz
   vlowz = gxz*dvelx + gyz*dvely + gzz*dvelz
File [modified]: GRHydro_UpdateMask.F90
Delta lines: +50 -23
===================================================================
--- trunk/src/GRHydro_UpdateMask.F90	2011-05-23 22:31:18 UTC (rev 248)
+++ trunk/src/GRHydro_UpdateMask.F90	2011-06-11 21:51:15 UTC (rev 249)
@@ -219,13 +219,19 @@
   CCTK_REAL :: det, psi4pt
CCTK_INT :: type_bits, atmosphere, not_atmosphere
+! begin EOS Omni vars
+  integer :: n,keytemp,anyerr,keyerr(1)
+  real*8  :: xpress,xeps,xtemp,xye
+  n=1;keytemp=0;anyerr=0;keyerr(1)=0
+  xpress=0.0d0;xye=0.0d0;xeps=0.0d0;xtemp=0.0d0
+! end EOS Omni vars
call SpaceMask_GetTypeBits(type_bits, "Hydro_Atmosphere")
   call SpaceMask_GetStateBits(atmosphere, "Hydro_Atmosphere",&
                               "in_atmosphere")
   call SpaceMask_GetStateBits(not_atmosphere, "Hydro_Atmosphere",&
                               "not_in_atmosphere")
-
+!$OMP PARALLEL DO PRIVATE(det,keytemp,i,j,k,anyerr,keyerr)
   do k = 1, cctk_lsh(3)
     do j = 1, cctk_lsh(2)
       do i = 1, cctk_lsh(1)
@@ -235,37 +241,52 @@
                     atmosphere)) &
              &) then
-!!$          write(*,*) 'Resetting at ',i,j,k, atmosphere_mask(i, j, k), &
-!!$             & (SpaceMask_CheckStateBitsF90(space_mask,i, j, k, type_bits,\
-!!$                    atmosphere)) 
-
           rho(i,j,k) = GRHydro_rho_min
           velx(i,j,k) = 0.0d0
           vely(i,j,k) = 0.0d0
           velz(i,j,k) = 0.0d0
           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 prim2conpolytype(GRHydro_polytrope_handle, &
-               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), scon(i,j,k,1), scon(i,j,k,2), scon(i,j,k,3), &
-               tau(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), w_lorentz(i,j,k))
-          if (wk_atmosphere .eq. 0) then
-            atmosphere_mask(i, j, k) = 0
-            SpaceMask_SetStateBitsF90(space_mask, i, j, k, type_bits,\
-                                      not_atmosphere)
-          end if
-        end if
+          if(evolve_temper.ne.0) then
+!             ! set the temperature to be relatively low
+             temperature(i,j,k) = grhydro_hot_atmo_temp
+             y_e(i,j,k) = grhydro_hot_atmo_Y_e
+             keytemp = 1
+             call EOS_Omni_press(GRHydro_eos_handle,keytemp,GRHydro_eos_rf_prec,n,&
+                  rho(i,j,k),eps(i,j,k),temperature(i,j,k),y_e(i,j,k),&
+                  press(i,j,k),keyerr,anyerr)
+             call prim2con_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),scon(i,j,k,1), scon(i,j,k,2), scon(i,j,k,3), &
+                  tau(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), w_lorentz(i,j,k),&
+                  temperature(i,j,k),y_e(i,j,k))
+             y_e_con(i,j,k) = dens(i,j,k) * y_e(i,j,k)
+          else
+             call prim2conpolytype(GRHydro_polytrope_handle, &
+                  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), scon(i,j,k,1), scon(i,j,k,2), scon(i,j,k,3), &
+                  tau(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), w_lorentz(i,j,k))
+             if (wk_atmosphere .eq. 0) then
+                atmosphere_mask(i, j, k) = 0
+                SpaceMask_SetStateBitsF90(space_mask, i, j, k, type_bits,\
+                not_atmosphere)
+             end if
+          endif
+
+       end if
+
       end do
     end do
   end do
+!$OMP END PARALLEL DO
-!!$  call GRHydro_Boundaries(CCTK_PASS_FTOF)
-
 end subroutine GRHydro_AtmosphereReset
subroutine GRHydro_InitialAtmosphereReset(CCTK_ARGUMENTS)
@@ -316,8 +337,9 @@
           velz(i,j,k) = 0.0d0
if(evolve_temper.ne.0) then
-             ! set the temperature to be relatively low
-             temperature(i,j,k) = 0.1d0
+!             ! set the temperature to be relatively low
+             temperature(i,j,k) = grhydro_hot_atmo_temp
+             y_e(i,j,k) = grhydro_hot_atmo_Y_e
              keytemp = 1
              call EOS_Omni_press(GRHydro_eos_handle,keytemp,GRHydro_eos_rf_prec,n,&
                   rho(i,j,k),eps(i,j,k),temperature(i,j,k),y_e(i,j,k),&
@@ -330,6 +352,7 @@
                   tau(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), w_lorentz(i,j,k),&
                   temperature(i,j,k),y_e(i,j,k))
+             y_e_con(i,j,k) = dens(i,j,k) * y_e(i,j,k)
           else
              call EOS_Omni_press(GRHydro_polytrope_handle,keytemp,GRHydro_eos_rf_prec,n,&
                   rho(i,j,k),eps(i,j,k),xtemp,xye,press(i,j,k),keyerr,anyerr)
@@ -355,7 +378,8 @@
if(evolve_temper.ne.0) then
              ! set the temperature to be relatively low
-             temperature_p(i,j,k) = 0.1d0
+             temperature_p(i,j,k) = grhydro_hot_atmo_temp
+             y_e_p(i,j,k) = grhydro_hot_atmo_Y_e
              keytemp = 1
              call EOS_Omni_press(GRHydro_eos_handle,keytemp,GRHydro_eos_rf_prec,n,&
                   rho_p(i,j,k),eps_p(i,j,k),temperature_p(i,j,k),y_e_p(i,j,k),&
@@ -368,6 +392,7 @@
                   tau_p(i,j,k), rho_p(i,j,k), velx_p(i,j,k), vely_p(i,j,k), &
                   velz_p(i,j,k), eps_p(i,j,k), press_p(i,j,k), w_lorentz_p(i,j,k),&
                   temperature_p(i,j,k),y_e_p(i,j,k))
+             y_e_con_p(i,j,k) = dens_p(i,j,k) * y_e_p(i,j,k)
             else
                call EOS_Omni_press(GRHydro_polytrope_handle,keytemp,GRHydro_eos_rf_prec,n,&
                     rho_p(i,j,k),eps_p(i,j,k),xtemp,xye,press_p(i,j,k),keyerr,anyerr)
@@ -394,7 +419,8 @@
if(evolve_temper.ne.0) then
              ! set the temperature to be relatively low
-             temperature_p_p(i,j,k) = 0.1d0
+             temperature_p_p(i,j,k) = grhydro_hot_atmo_temp
+             y_e_p_p(i,j,k) = grhydro_hot_atmo_Y_e
              keytemp = 1
              call EOS_Omni_press(GRHydro_eos_handle,keytemp,GRHydro_eos_rf_prec,n,&
                   rho_p_p(i,j,k),eps_p_p(i,j,k),temperature_p_p(i,j,k),y_e_p_p(i,j,k),&
@@ -406,6 +432,7 @@
                   tau_p_p(i,j,k), rho_p_p(i,j,k), velx_p_p(i,j,k), vely_p_p(i,j,k), &
                   velz_p_p(i,j,k), eps_p_p(i,j,k), press_p_p(i,j,k), w_lorentz_p_p(i,j,k),&
                   temperature_p_p(i,j,k),y_e_p_p(i,j,k))
+             y_e_con_p_p(i,j,k) = dens_p_p(i,j,k) * y_e_p_p(i,j,k)
             else
                call EOS_Omni_press(GRHydro_polytrope_handle,keytemp,GRHydro_eos_rf_prec,n,&
                     rho_p_p(i,j,k),eps_p_p(i,j,k),xtemp,xye,press_p_p(i,j,k),keyerr,anyerr)

    