thorns for general eos interface (moved from Whisky_Dev repositories)

git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinEOS/EOSG_Hybrid/trunk@2 0a999221-0efe-4d2d-b270-2ec6803ffb49

10 files changed, 514 insertions, 0 deletions

diff --git a/src/EOS_GeneralHybrid.F90 b/src/EOS_GeneralHybrid.F90
new file mode 100755
index 0000000..e842600
--- /dev/null
+++ b/src/EOS_GeneralHybrid.F90
@@ -0,0 +1,228 @@
+ /*@@
+   @file      EOS_GeneralHybrid.F90
+   @date      
+   @author    Ian Hawke, Christian D. Ott
+   @desc 
+      Routines to calculate the EOS used by Dimmelmeier et al. 
+      in supernova core collapse simulations.
+   @enddesc 
+ @@*/
+
+#include "cctk.h"
+#include "cctk_Parameters.h"
+
+subroutine EOS_GeneralHybrid_Pressure(nelems, rho, eps, press)
+
+  USE EOS_GP_Scalars   
+  USE EOS_GeneralHybrid_Scalars
+
+  implicit none
+
+  DECLARE_CCTK_PARAMETERS
+
+  CCTK_INT, intent(in) :: nelems
+  CCTK_REAL, dimension(nelems), intent(in) :: rho
+  CCTK_REAL, dimension(nelems), intent(in) :: eps
+  CCTK_REAL, dimension(nelems), intent(out) :: press
+
+  CCTK_REAL, allocatable, dimension(:) :: p_poly
+  CCTK_REAL, allocatable, dimension(:) :: p_th
+  CCTK_REAL, allocatable, dimension(:) :: local_eos_gamma
+  CCTK_REAL, allocatable, dimension(:) :: local_eos_k_cgs
+
+  CCTK_REAL zero
+
+  allocate(p_poly(1:nelems))
+  allocate(p_th(1:nelems))
+  allocate(local_eos_gamma(1:nelems))
+  allocate(local_eos_k_cgs(1:nelems))
+
+  zero = 0.d0
+
+  where (rho > rho_nuc) 
+     local_eos_gamma = eos_gamma_supernuclear
+     local_eos_k_cgs = eos_k_supernuclear_cgs
+  elsewhere
+     local_eos_gamma = eos_gamma
+     local_eos_k_cgs = eos_k_cgs
+  end where
+
+  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)
+
+  press = p_poly + p_th
+
+  deallocate(p_poly)
+  deallocate(p_th)
+  deallocate(local_eos_gamma)
+  deallocate(local_eos_k_cgs)
+
+
+end subroutine EOS_GeneralHybrid_Pressure
+
+subroutine EOS_GeneralHybrid_DPDRho(nelems, rho, eps, dpdrho)
+
+  USE EOS_GP_Scalars
+  USE EOS_GeneralHybrid_Scalars
+
+  implicit none
+  DECLARE_CCTK_PARAMETERS
+
+  CCTK_INT, intent(in) :: nelems
+  CCTK_REAL, dimension(nelems), intent(in) :: rho
+  CCTK_REAL, dimension(nelems), intent(in) :: eps
+  CCTK_REAL, dimension(nelems), intent(out) :: dpdrho
+
+  CCTK_REAL, allocatable,dimension(:) :: d_p_poly
+  CCTK_REAL, allocatable,dimension(:) :: d_p_th_1
+  CCTK_REAL, allocatable,dimension(:) :: d_p_th_2
+  CCTK_REAL, allocatable,dimension(:) :: local_eos_gamma
+  CCTK_REAL, allocatable,dimension(:) :: local_eos_k_cgs
+
+  CCTK_REAL zero
+
+  allocate(d_p_poly(1:nelems))
+  allocate(d_p_th_1(1:nelems))
+  allocate(d_p_th_2(1:nelems))
+  allocate(local_eos_gamma(1:nelems))
+  allocate(local_eos_k_cgs(1:nelems))
+
+  
+  zero = 0.d0
+
+  where (rho > rho_nuc) 
+     local_eos_gamma = eos_gamma_supernuclear
+     local_eos_k_cgs = eos_k_supernuclear_cgs
+  elsewhere
+     local_eos_gamma = eos_gamma
+     local_eos_k_cgs = eos_k_cgs
+  end where
+
+  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)
+      
+  dpdrho = d_p_poly + d_p_th_1 + d_p_th_2
+
+  deallocate(d_p_poly)
+  deallocate(d_p_th_1)
+  deallocate(d_p_th_2)
+  deallocate(local_eos_gamma)
+  deallocate(local_eos_k_cgs)
+
+
+end subroutine EOS_GeneralHybrid_DPDRho
+
+
+subroutine EOS_GeneralHybrid_cs2(nelems, rho, eps, cs2)
+
+! This is an ugly and quick hack.
+! Uses way to many operations and memory.
+
+  USE EOS_GP_Scalars
+  USE EOS_GeneralHybrid_Scalars
+
+  implicit none
+  DECLARE_CCTK_PARAMETERS
+
+  CCTK_INT, intent(in) :: nelems
+  CCTK_REAL, dimension(nelems), intent(in) :: rho
+  CCTK_REAL, dimension(nelems), intent(in) :: eps
+  CCTK_REAL, dimension(nelems), intent(out) :: cs2
+  CCTK_REAL, allocatable,dimension(:) :: p_poly
+  CCTK_REAL, allocatable,dimension(:) :: p_th
+  CCTK_REAL, allocatable,dimension(:) :: pressure
+
+  CCTK_REAL, allocatable,dimension(:) :: local_eos_gamma
+  CCTK_REAL, allocatable,dimension(:) :: local_eos_k_cgs
+
+  CCTK_REAL zero
+
+  allocate(p_poly(1:nelems))
+  allocate(p_th(1:nelems))
+  allocate(pressure(1:nelems))
+  allocate(local_eos_gamma(1:nelems))
+  allocate(local_eos_k_cgs(1:nelems))
+
+
+  zero = 0.d0
+
+  where (rho > rho_nuc) 
+     local_eos_gamma = eos_gamma_supernuclear
+     local_eos_k_cgs = eos_k_supernuclear_cgs
+  elsewhere
+     local_eos_gamma = eos_gamma
+     local_eos_k_cgs = eos_k_cgs
+  end where
+
+! First calculate the pressure 
+  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
+
+  pressure = p_poly + p_th
+  
+  p_th = max(zero, p_th)
+
+! This may look incorrect; It's in Harry's code; in his thesis and I worked it out as well.
+! Should be okay...
+  cs2 = (local_eos_gamma * p_poly + eos_gamma_th * p_th) / &
+       rho / (1.0d0 + eps + pressure/rho)
+
+  deallocate(pressure)
+  deallocate(p_poly)
+  deallocate(p_th)
+  deallocate(local_eos_gamma)
+  deallocate(local_eos_k_cgs)
+
+end subroutine EOS_GeneralHybrid_cs2
+
+
+
+subroutine EOS_GeneralHybrid_DPDIE(nelems, rho, eps, dpdie)
+
+
+  USE EOS_GP_Scalars
+  USE EOS_GeneralHybrid_Scalars
+
+  implicit none
+  DECLARE_CCTK_PARAMETERS
+
+  CCTK_INT, intent(in) :: nelems
+  CCTK_REAL, dimension(nelems), intent(in) :: rho
+  CCTK_REAL, dimension(nelems), intent(in) :: eps
+  CCTK_REAL, dimension(nelems), intent(out) :: dpdie
+
+  dpdie = (eos_gamma_th - 1.d0) * rho
+
+end subroutine EOS_GeneralHybrid_DPDIE
+
diff --git a/src/EOS_GeneralHybrid.c b/src/EOS_GeneralHybrid.c
new file mode 100755
index 0000000..9e4d443
--- /dev/null
+++ b/src/EOS_GeneralHybrid.c
@@ -0,0 +1,76 @@
+#include <ctype.h>
+#include <string.h>
+#include <stdlib.h>
+#include <stdio.h>
+
+#include "cctk.h"
+#include "cctk_Parameters.h"
+#include "cctk_Arguments.h"
+
+#include "EOS_GeneralHybrid.h"
+
+void CCTK_FCALL CCTK_FNAME(EOS_GeneralHybrid_Pressure) (const CCTK_INT* nelems,
+                                             const CCTK_REAL* const rho,
+                                             const CCTK_REAL* const eps,
+                                             const CCTK_REAL* press);
+void CCTK_FCALL CCTK_FNAME(EOS_GeneralHybrid_DPDRho)   (const CCTK_INT* nelems,
+                                             const CCTK_REAL* const rho,
+                                             const CCTK_REAL* const eps,
+                                             const CCTK_REAL* dpdrho);
+void CCTK_FCALL CCTK_FNAME(EOS_GeneralHybrid_DPDIE)    (const CCTK_INT* nelems,
+                                             const CCTK_REAL* const rho,
+                                             const CCTK_REAL* const eps,
+                                             const CCTK_REAL* dpdie);
+void CCTK_FCALL CCTK_FNAME(EOS_GeneralHybrid_cs2)      (const CCTK_INT* nelems,
+                                             const CCTK_REAL* const rho,
+                                             const CCTK_REAL* const eps,
+                                             const CCTK_REAL* cs2);
+
+CCTK_INT EOS_GeneralHybrid_SetArray(const CCTK_INT param_table,
+                                    const CCTK_INT n_elems,
+                                    const CCTK_POINTER* indep_vars,
+                                    const CCTK_INT* which_deps_to_set,
+                                    CCTK_POINTER* dep_vars)
+{
+
+  DECLARE_CCTK_PARAMETERS;
+  
+  CCTK_INT ierr = 0, var;
+
+  for (var = 0; var < N_DEPS; ++var)
+  {
+    if (which_deps_to_set[var])
+    {
+      switch (var)
+      {
+        case (0):
+          CCTK_FNAME(EOS_GeneralHybrid_Pressure) (&n_elems,
+                                                  (const CCTK_REAL*)indep_vars[0],
+						  (const CCTK_REAL*)indep_vars[1], 
+                                                  (const CCTK_REAL*)dep_vars[0]);
+          break;
+        case (1):
+          CCTK_FNAME(EOS_GeneralHybrid_DPDRho)   (&n_elems,
+                                       
+                                                  (const CCTK_REAL*)indep_vars[0], 
+                                                  (const CCTK_REAL*)indep_vars[1], 
+                                                  (const CCTK_REAL*)dep_vars[1]);
+          break;
+        case (2):
+          CCTK_FNAME(EOS_GeneralHybrid_DPDIE)    (&n_elems,
+                                                  (const CCTK_REAL*)indep_vars[0], 
+                                                  (const CCTK_REAL*)indep_vars[1], 
+                                                  (const CCTK_REAL*)dep_vars[2]);
+          break;
+        case (3):
+          CCTK_FNAME(EOS_GeneralHybrid_cs2)      (&n_elems,
+                                                  (const CCTK_REAL*)indep_vars[0], 
+                                                  (const CCTK_REAL*)indep_vars[1], 
+                                                  (const CCTK_REAL*)dep_vars[3]);
+          break;
+      }
+    }
+  }
+  
+  return ierr;
+}
diff --git a/src/EOS_GeneralHybrid.h b/src/EOS_GeneralHybrid.h
new file mode 100755
index 0000000..8c98260
--- /dev/null
+++ b/src/EOS_GeneralHybrid.h
@@ -0,0 +1,15 @@
+#ifndef EOS_GENERALHYBRID_H
+#define EOS_GENERALHYBRID_H
+
+#include "cctk.h"
+
+#define N_INDEPS 2
+#define N_DEPS 4
+
+CCTK_INT EOS_GeneralHybrid_SetArray(const CCTK_INT param_table,
+			 const CCTK_INT n_elems,
+                         const CCTK_POINTER* indep_vars,
+                         const CCTK_INT* which_deps_to_set,
+                         CCTK_POINTER* dep_vars);
+
+#endif
diff --git a/src/EOS_GeneralHybrid_Analysis.F90 b/src/EOS_GeneralHybrid_Analysis.F90
new file mode 100755
index 0000000..3fea75e
--- /dev/null
+++ b/src/EOS_GeneralHybrid_Analysis.F90
@@ -0,0 +1,74 @@
+ /*@@
+   @file      EOS_GeneralHybrid_Analysis.F
+   @date      
+   @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_General
+   @date       
+   @author     Harry Dimmelmeier, Ian Hawke, Christian Ott
+   @desc 
+      The routine that calculates the contributions.
+   @enddesc 
+   @calls     
+   @calledby   
+   @history 
+ 
+   @endhistory 
+@@*/
+
+subroutine Check_Poly_Contrib_General(CCTK_ARGUMENTS)
+
+  USE EOS_GP_Scalars
+  USE EOS_GeneralHybrid_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)
+
+         end do
+      end do
+   end do
+
+ end subroutine Check_Poly_Contrib_General
diff --git a/src/EOS_GeneralHybrid_Scalars.F90 b/src/EOS_GeneralHybrid_Scalars.F90
new file mode 100755
index 0000000..382e70e
--- /dev/null
+++ b/src/EOS_GeneralHybrid_Scalars.F90
@@ -0,0 +1,18 @@
+ /*@@
+   @file      EOS_GeneralHybrid_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_GeneralHybrid_Scalars
+
+      implicit none
+
+      CCTK_REAL :: eos_k_supernuclear_cgs
+   
+      end module EOS_GeneralHybrid_Scalars
diff --git a/src/EOS_GeneralHybrid_Setup.F90 b/src/EOS_GeneralHybrid_Setup.F90
new file mode 100755
index 0000000..1bb8439
--- /dev/null
+++ b/src/EOS_GeneralHybrid_Setup.F90
@@ -0,0 +1,40 @@
+ /*@@
+   @file      EOS_GeneralHybrid_Setup.F90
+   @date      
+   @author    Ian Hawke, Christian D. Ott
+   @desc 
+   Startup for EOS_Hybrid. 
+   @enddesc 
+ @@*/
+
+#include "cctk.h"
+#include "cctk_Parameters.h"
+
+ /*@@
+   @routine    EOS_GeneralHybrid_Setup
+   @date       
+   @author     Ian Hawkem, Christian d. Ott
+   @desc 
+   Setup EOS_HybridScalars
+   @enddesc 
+   @calls     
+   @calledby   
+   @history 
+ 
+   @endhistory 
+@@*/
+
+subroutine EOS_GeneralHybrid_Setup()
+
+  USE EOS_GP_Scalars
+  USE EOS_GeneralHybrid_Scalars
+
+  implicit none
+
+  DECLARE_CCTK_PARAMETERS
+  
+  
+  eos_k_supernuclear_cgs = eos_k_cgs * (rho_nuc * rho_geom_factor_inv)** &
+       (eos_gamma - eos_gamma_supernuclear)
+
+end subroutine EOS_GeneralHybrid_Setup
diff --git a/src/EOS_GeneralHybrid_Startup.c b/src/EOS_GeneralHybrid_Startup.c
new file mode 100755
index 0000000..034aa11
--- /dev/null
+++ b/src/EOS_GeneralHybrid_Startup.c
@@ -0,0 +1,42 @@
+#include "cctk.h"
+#include "cctk_Arguments.h"
+#include "cctk_Parameters.h"
+
+#include "util_String.h"
+#include "util_ErrorCodes.h"
+#include "util_Table.h"
+
+#include "EOS_GeneralHybrid.h"
+
+void EOS_GeneralHybrid_Startup (CCTK_ARGUMENTS)
+{
+  DECLARE_CCTK_ARGUMENTS;
+  DECLARE_CCTK_PARAMETERS;
+
+  CCTK_INT ierr, table_handle;
+
+  table_handle = Util_TableCreate(UTIL_TABLE_FLAGS_CASE_INSENSITIVE);
+
+  Util_TableSetInt(table_handle,
+                   N_INDEPS,
+                   "N independent variables");
+  Util_TableSetInt(table_handle,
+                   N_DEPS,
+                   "N dependent variables");
+  Util_TableSetString(table_handle,
+                      "Rho SpecificInternalEnergy",
+                      "Independent variable names");
+  Util_TableSetString(table_handle,
+                      "Pressure DPressureDRho DPressureDSpecificInternalEnergy c_s^2",
+                      "Dependent variable names");
+  Util_TableSetString(table_handle,
+                      "Hybrid",
+                      "EOS Name");
+
+  ierr = EOS_RegisterCall(table_handle,
+                          EOS_GeneralHybrid_SetArray);
+  if (ierr)
+    {
+      CCTK_WARN(0, "Failed to set up EOS_GeneralHybrid call");
+    }
+}
diff --git a/src/make.code.defn b/src/make.code.defn
new file mode 100755
index 0000000..f5cab0a
--- /dev/null
+++ b/src/make.code.defn
@@ -0,0 +1,11 @@
+# Main make.code.defn file for thorn EOS_GeneralHybrid
+# $Header$
+
+# Source files in this directory
+SRCS = EOS_GeneralHybrid.F90 EOS_GeneralHybrid_Scalars.F90 \
+	EOS_GeneralHybrid_Setup.F90 EOS_GeneralHybrid_Analysis.F90 \
+	EOS_GeneralHybrid.c EOS_GeneralHybrid_Startup.c \
+
+
+# Subdirectories containing source files
+SUBDIRS = 
diff --git a/src/make.code.deps b/src/make.code.deps
new file mode 100755
index 0000000..b6bdb4b
--- /dev/null
+++ b/src/make.code.deps
@@ -0,0 +1,5 @@
+# Module dependencies for thorn EOS_GeneralHybrid
+
+EOS_GeneralHybrid.F90.o:		EOS_GeneralHybrid_Scalars.F90.o
+EOS_GeneralHybrid_Analysis.F90.o:	EOS_GeneralHybrid_Scalars.F90.o
+EOS_GeneralHybrid_Setup.F90.o:		EOS_GeneralHybrid_Scalars.F90.o
diff --git a/src/make.configuration.deps b/src/make.configuration.deps
new file mode 100755
index 0000000..61740e1
--- /dev/null
+++ b/src/make.configuration.deps
@@ -0,0 +1,5 @@
+# Make sure that EOS_Polytrope is compiled first
+
+ifneq (,$(findstring Whisky_Dev/EOS_GeneralPolytrope,$(THORNS)))
+$(CCTK_LIBDIR)$(DIRSEP)libEOS_GeneralHybrid.a : $(CCTK_LIBDIR)$(DIRSEP)libEOS_GeneralPolytrope.a
+endif