move development of EOS_Polytrope and EOS_Hybrid from Whisky to Cactus

git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinEOS/EOS_Hybrid/trunk@2 9c56c927-b35f-42c5-a3f6-df2e9f55b230

7 files changed, 344 insertions, 0 deletions

diff --git a/src/EOS_Hybrid.F b/src/EOS_Hybrid.F
new file mode 100644
index 0000000..c74af75
--- /dev/null
+++ b/src/EOS_Hybrid.F
@@ -0,0 +1,144 @@
+ /*@@
+   @file      EOS_Hybrid.F
+   @date      Wed Mar 20 14:56:35 2002
+   @author    Ian Hawke
+   @desc 
+      Routines to calculate the EOS used by 
+      Dimmelmeier Novak Font Ibanez Mueller PRD71 064023 (2005)
+      in supernova core collapse simulations.
+   @enddesc 
+ @@*/
+
+#include "cctk.h"
+#include "cctk_Parameters.h"
+
+      CCTK_REAL function EOS_Hybrid_Pressure(rho, eps)
+
+      USE EOS_Polytrope_Scalars   
+      USE EOS_Hybrid_Scalars
+
+      implicit none
+
+      DECLARE_CCTK_PARAMETERS
+
+      CCTK_REAL rho, eps, local_eos_gamma, p_poly, p_th, local_eos_k_cgs, zero
+
+      zero = 0.d0
+
+      if (rho > rho_nuc) then
+        local_eos_gamma = eos_gamma_supernuclear
+        local_eos_k_cgs = eos_k_supernuclear_cgs
+      else
+        local_eos_gamma = eos_gamma
+        local_eos_k_cgs = eos_k_cgs
+      end if
+
+      p_poly = p_geom_factor * local_eos_k_cgs * 
+     .     (rho * rho_geom_factor_inv)**local_eos_gamma
+
+      p_th = - p_geom_factor * local_eos_k_cgs * (eos_gamma_th - 1.d0) / 
+     .     (local_eos_gamma - 1.d0) * (rho * rho_geom_factor_inv)**local_eos_gamma + 
+     .     (eos_gamma_th - 1.d0) * rho * eps - 
+     .     (eos_gamma_th - 1.d0) * (local_eos_gamma - eos_gamma) / 
+     .     (eos_gamma - 1.d0) / (eos_gamma_supernuclear - 1.d0) * 
+     .     p_geom_factor * eos_k_cgs * rho_geom_factor_inv**eos_gamma * 
+     .     rho_nuc**(eos_gamma - 1.d0) * rho
+
+      p_th = max(zero, p_th)
+
+      EOS_Hybrid_Pressure = p_poly + p_th
+
+      end function EOS_Hybrid_Pressure
+
+c     The specific internal energy isn''t correct yet
+
+      CCTK_REAL function EOS_Hybrid_SpecificIE(rho, press)
+
+      USE EOS_Polytrope_Scalars   
+      USE EOS_Hybrid_Scalars
+
+      implicit none
+      DECLARE_CCTK_PARAMETERS
+
+      CCTK_REAL rho, press
+
+c$$$      if (rho > rho_nuc) then
+c$$$        EOS_Hybrid_SpecificIE = eos_k_cgs / (eos_gamma_2 - 1.d0) * 
+c$$$     .       (rho / rho_nuc) ** eos_gamma_2 * rho_nuc ** eos_gamma_1 / 
+c$$$     .       rho + 
+c$$$     .       eos_k_cgs * (eos_gamma_2 - eos_gamma_1) / 
+c$$$     .       (eos_gamma_2 - 1.d0) / (eos_gamma_1 - 1.d0) * 
+c$$$     .       rho_nuc ** (eos_gamma_1 - 1.d0)
+c$$$      else
+c$$$        EOS_Hybrid_SpecificIE = eos_k_cgs * (rho ** (eos_gamma_1 - 1.d0)) /
+c$$$     .       (eos_gamma_1 - 1.d0)
+c$$$      end if
+
+      EOS_Hybrid_SpecificIE = 2.34567890d0
+
+      end function EOS_Hybrid_SpecificIE
+
+c     The rest mass density isn''t correct yet
+
+      CCTK_REAL function EOS_Hybrid_RestMassDens(eps, press)
+
+      implicit none
+      DECLARE_CCTK_PARAMETERS
+
+      CCTK_REAL eps, press
+
+      EOS_Hybrid_RestMassDens = 1.23456789d0
+      
+      end function EOS_Hybrid_RestMassDens
+      
+      CCTK_REAL function EOS_Hybrid_DPressByDRho(rho, eps)
+
+      USE EOS_Polytrope_Scalars   
+      USE EOS_Hybrid_Scalars   
+
+      implicit none
+      DECLARE_CCTK_PARAMETERS
+
+      CCTK_REAL rho, eps, local_eos_gamma, local_eos_k_cgs, d_p_poly, d_p_th_1, 
+     .     d_p_th_2, zero
+
+      zero = 0.d0
+
+      if (rho > rho_nuc) then
+        local_eos_gamma = eos_gamma_supernuclear
+        local_eos_k_cgs = eos_k_supernuclear_cgs
+      else
+        local_eos_gamma = eos_gamma
+        local_eos_k_cgs = eos_k_cgs
+      end if
+
+      d_p_poly = local_eos_gamma * p_geom_factor * local_eos_k_cgs * 
+     .     rho**(local_eos_gamma - 1.d0) * rho_geom_factor_inv**local_eos_gamma
+
+      d_p_th_1 = - local_eos_gamma * p_geom_factor * local_eos_k_cgs * 
+     .     (eos_gamma_th - 1.d0) / (local_eos_gamma - 1.d0) * 
+     .     rho**(local_eos_gamma - 1.d0) * rho_geom_factor_inv**local_eos_gamma
+
+      d_p_th_2 = (eos_gamma_th - 1.d0) * eps
+     .     - (eos_gamma_th - 1.d0) * (local_eos_gamma - eos_gamma) / 
+     .     (eos_gamma - 1.d0) / (eos_gamma_supernuclear - 1.d0) * 
+     .     p_geom_factor * eos_k_cgs * rho_geom_factor_inv**eos_gamma * 
+     .     rho_nuc**(eos_gamma - 1.d0)
+
+!      d_p_th_1 = max(d_p_th_1, zero)
+!      d_p_th_2 = max(d_p_th_2, zero)
+      
+      EOS_Hybrid_DPressByDRho = d_p_poly + d_p_th_1 + d_p_th_2
+
+      end function EOS_Hybrid_DPressByDRho
+
+      CCTK_REAL function EOS_Hybrid_DPressByDEps(rho, eps)
+
+      implicit none
+      DECLARE_CCTK_PARAMETERS
+
+      CCTK_REAL rho, eps
+
+      EOS_Hybrid_DPressByDEps = (eos_gamma_th - 1.d0) * rho
+
+      end function EOS_Hybrid_DPressByDEps
diff --git a/src/EOS_Hybrid_Analysis.F b/src/EOS_Hybrid_Analysis.F
new file mode 100644
index 0000000..9f82a41
--- /dev/null
+++ b/src/EOS_Hybrid_Analysis.F
@@ -0,0 +1,76 @@
+ /*@@
+   @file      EOS_Hybrid_Analysis.F
+   @date      Fri Apr 26 16:14:47 2002
+   @author    Harry Dimmelmeier, Ian Hawke, Christian Ott
+   @desc 
+      Calculates the polytropic and thermal contributions to the 
+      total pressure.
+   @enddesc 
+ @@*/
+
+#include "cctk.h"
+#include "cctk_Arguments.h"
+#include "cctk_Parameters.h"
+
+ /*@@
+   @routine    Check_Poly_Contrib
+   @date       Fri Apr 26 16:25:35 2002
+   @author     Harry Dimmelmeier, Ian Hawke, Christian Ott
+   @desc 
+      The routine that calculates the contributions.
+   @enddesc 
+   @calls     
+   @calledby   
+   @history 
+ 
+   @endhistory 
+@@*/
+
+      subroutine Check_Poly_Contrib(CCTK_ARGUMENTS)
+
+      USE EOS_Polytrope_Scalars
+      USE EOS_Hybrid_Scalars
+
+      implicit none
+
+      DECLARE_CCTK_ARGUMENTS
+      DECLARE_CCTK_PARAMETERS
+
+      CCTK_INT :: i,j,k
+
+
+      CCTK_REAL local_eos_gamma, local_eos_k_cgs, d_p_poly, d_p_th_1,
+     .     d_p_th_2, zero
+
+      do k = 1, cctk_lsh(3)
+        do j = 1, cctk_lsh(2)
+          do i = 1, cctk_lsh(1)
+
+            if (rho(i,j,k) > rho_nuc) then
+              local_eos_gamma = eos_gamma_supernuclear
+              local_eos_k_cgs = eos_k_supernuclear_cgs
+            else
+              local_eos_gamma = eos_gamma
+              local_eos_k_cgs = eos_k_cgs
+            end if
+
+            pressure_poly(i,j,k) = p_geom_factor * local_eos_k_cgs * 
+     .           (rho(i,j,k) * rho_geom_factor_inv)**local_eos_gamma
+            
+            pressure_th(i,j,k) = - p_geom_factor * local_eos_k_cgs * 
+     .           (eos_gamma_th - 1.d0) / (local_eos_gamma - 1.d0) * 
+     .           (rho(i,j,k) * rho_geom_factor_inv)**local_eos_gamma + 
+     .           (eos_gamma_th - 1.d0) * rho(i,j,k) * eps(i,j,k) - 
+     .           (eos_gamma_th - 1.d0) * (local_eos_gamma - eos_gamma) / 
+     .           (eos_gamma - 1.d0) / (eos_gamma_supernuclear - 1.d0) * 
+     .           p_geom_factor * eos_k_cgs * 
+     .           rho_geom_factor_inv**eos_gamma * 
+     .           rho_nuc**(eos_gamma - 1.d0) * rho(i,j,k)
+
+c$$$            pressure_th(i,j,k) = dmax1(0.d0, pressure_th(i,j,k))
+
+          end do
+        end do
+      end do
+
+      end subroutine Check_Poly_Contrib
diff --git a/src/EOS_Hybrid_Scalars.F b/src/EOS_Hybrid_Scalars.F
new file mode 100644
index 0000000..3e85c39
--- /dev/null
+++ b/src/EOS_Hybrid_Scalars.F
@@ -0,0 +1,18 @@
+ /*@@
+   @file      EOS_Hybrid_Scalars.F
+   @date      Fri Apr 26 13:14:18 2002
+   @author    Harry Dimmelmeier, Ian Hawke, Christian Ott
+   @desc 
+      Constant for the EOS routines.
+   @enddesc 
+ @@*/
+
+#include "cctk.h"
+
+      module EOS_Hybrid_Scalars
+
+      implicit none
+
+      CCTK_REAL :: eos_k_supernuclear_cgs
+   
+      end module EOS_Hybrid_Scalars
diff --git a/src/EOS_Hybrid_Startup.F b/src/EOS_Hybrid_Startup.F
new file mode 100644
index 0000000..24111d4
--- /dev/null
+++ b/src/EOS_Hybrid_Startup.F
@@ -0,0 +1,84 @@
+ /*@@
+   @file      EOS_Hybrid_Startup.F
+   @date      Wed Mar 20 14:50:22 2002
+   @author    Ian Hawke
+   @desc 
+   Startup for EOS_Hybrid. 
+   @enddesc 
+ @@*/
+
+#include "cctk.h"
+#include "cctk_Parameters.h"
+
+ /*@@
+   @routine    EOS_Hybrid_Startup
+   @date       Wed Mar 20 14:52:41 2002
+   @author     Ian Hawke
+   @desc 
+   Startup for EOS_Hybrid
+   @enddesc 
+   @calls     
+   @calledby   
+   @history 
+ 
+   @endhistory 
+@@*/
+
+      integer function EOS_Hybrid_Startup()
+
+      USE EOS_Polytrope_Scalars
+      USE EOS_Hybrid_Scalars
+
+      implicit none
+
+      DECLARE_CCTK_PARAMETERS
+
+      integer handle, ierr
+      
+      external EOS_Hybrid_Pressure
+      external EOS_Hybrid_SpecificIE
+      external EOS_Hybrid_RestMassDens
+      external EOS_Hybrid_DPressByDRho
+      external EOS_Hybrid_DPressByDEps
+
+      call EOS_RegisterMethod(handle,"EOS_Hybrid")
+
+      if (handle .ge. 0) then
+
+        call EOS_RegisterPressure(ierr, handle, EOS_Hybrid_Pressure)
+        call EOS_RegisterSpecificIntEnergy(ierr, handle, 
+     .                                     EOS_Hybrid_SpecificIE)
+        call EOS_RegisterRestMassDens(ierr, handle, 
+     .                                EOS_Hybrid_RestMassDens)
+        call EOS_RegisterDPressByDRho(ierr, handle, 
+     .                                EOS_Hybrid_DPressByDRho)
+        call EOS_RegisterDPressByDEps(ierr, handle, 
+     .                                EOS_Hybrid_DPressByDEps)
+
+      else
+
+        call CCTK_WARN(0, "Unable to register the EOS method!")
+
+      end if
+
+c$$$      all these are now set by the polytropic EOS
+
+c$$$      These constants match the values in RNSID: consts.h
+
+c$$$      m_solar_cgs = 1.987d33
+c$$$      c_cgs = 2.9979d10
+c$$$      G_cgs = 6.6732d-8
+c$$$
+c$$$      m_solar_geom = G_cgs / c_cgs**2 * m_solar_cgs
+c$$$      p_geom_factor = G_cgs / c_cgs**4 * m_solar_geom**2
+c$$$      rho_geom_factor = p_geom_factor * c_cgs**2
+c$$$      rho_geom_factor_inv = 1.d0 / rho_geom_factor 
+c$$$
+c$$$      eos_k_cgs = eos_k * rho_geom_factor**gamma_ini / p_geom_factor 
+
+      eos_k_supernuclear_cgs = eos_k_cgs * (rho_nuc * rho_geom_factor_inv)**
+     .     (eos_gamma - eos_gamma_supernuclear)
+
+      EOS_Hybrid_Startup = 0
+        
+      end function EOS_Hybrid_Startup
diff --git a/src/make.code.defn b/src/make.code.defn
new file mode 100644
index 0000000..103c9f0
--- /dev/null
+++ b/src/make.code.defn
@@ -0,0 +1,12 @@
+# Main make.code.defn file for thorn EOS_Hybrid
+# $Header$
+
+# Source files in this directory
+SRCS = 	EOS_Hybrid.F \
+	EOS_Hybrid_Analysis.F \
+	EOS_Hybrid_Scalars.F \
+	EOS_Hybrid_Startup.F
+
+# Subdirectories containing source files
+SUBDIRS = 
+
diff --git a/src/make.code.deps b/src/make.code.deps
new file mode 100644
index 0000000..76300b4
--- /dev/null
+++ b/src/make.code.deps
@@ -0,0 +1,5 @@
+# Module dependencies for thorn EOS_Hybrid
+
+EOS_Hybrid.F.o:			EOS_Hybrid_Scalars.F.o
+EOS_Hybrid_Analysis.F.o:	EOS_Hybrid_Scalars.F.o
+EOS_Hybrid_Startup.F.o:		EOS_Hybrid_Scalars.F.o
\ No newline at end of file
diff --git a/src/make.configuration.deps b/src/make.configuration.deps
new file mode 100644
index 0000000..a1d3ce6
--- /dev/null
+++ b/src/make.configuration.deps
@@ -0,0 +1,5 @@
+# Make sure that EOS_Polytrope is compiled first
+
+ifneq (,$(findstring EUHydro/EOS_Polytrope,$(THORNS)))
+$(CCTK_LIBDIR)$(DIRSEP)libEOS_Hybrid.a : $(CCTK_LIBDIR)$(DIRSEP)libEOS_Polytrope.a
+endif