thorns for general eos interface (moved from Whisky_Dev repositories)

git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinEOS/EOSG_Polytrope/trunk@2 a192deed-9f5d-41fb-bb60-7f59df5b1b8d

8 files changed, 492 insertions, 0 deletions

diff --git a/src/EOS_GP.F90 b/src/EOS_GP.F90
new file mode 100755
index 0000000..283b9ac
--- /dev/null
+++ b/src/EOS_GP.F90
@@ -0,0 +1,189 @@
+ /*@@
+   @file      EOS_GP.F90
+   @date      Mon Mar 14 16:42:15 2005
+   @author    Ian Hawke
+   @desc 
+   The functions that actually set the polytropic EOS
+   @enddesc 
+ @@*/
+
+#include "cctk.h"
+#include "cctk_Parameters.h"
+
+subroutine EOS_GP_Pressure(nelems, rho, press)
+
+  USE EOS_GP_Scalars
+
+  implicit none
+
+  CCTK_INT, intent(in) :: nelems
+  CCTK_REAL, dimension(nelems), intent(in) :: rho
+  CCTK_REAL, dimension(nelems), intent(out) :: press
+
+  press = p_geom_factor * eos_k_cgs * &
+       (rho * rho_geom_factor_inv) ** eos_gamma_local  
+
+end subroutine EOS_GP_Pressure
+
+subroutine EOS_GP_IntEn(nelems, rho, IntEn)
+
+  USE EOS_GP_Scalars
+
+  implicit none
+
+  CCTK_INT, intent(in) :: nelems
+  CCTK_REAL, dimension(nelems), intent(in) :: rho
+  CCTK_REAL, dimension(nelems), intent(out) :: IntEn
+
+  IntEn = p_geom_factor * eos_k_cgs * &
+       (rho * rho_geom_factor_inv) ** eos_gamma_local / &
+       (rho * (eos_gamma_local - 1.d0) )
+
+end subroutine EOS_GP_IntEn
+
+subroutine EOS_GP_DPDRho(nelems, rho, dpdrho)
+
+  USE EOS_GP_Scalars
+
+  implicit none
+
+  CCTK_INT, intent(in) :: nelems
+  CCTK_REAL, dimension(nelems), intent(in) :: rho
+  CCTK_REAL, dimension(nelems), intent(out) :: dpdrho
+
+  dpdrho = p_geom_factor * eos_k_cgs * &
+           eos_gamma_local * rho_geom_factor_inv * &
+       (rho * rho_geom_factor_inv) ** (eos_gamma_local - 1.d0)
+
+end subroutine EOS_GP_DPDRho
+
+subroutine EOS_GP_DPDIE(nelems, rho, dpdie)
+
+  USE EOS_GP_Scalars
+
+  implicit none
+
+  CCTK_INT, intent(in) :: nelems
+  CCTK_REAL, dimension(nelems), intent(in) :: rho
+  CCTK_REAL, dimension(nelems), intent(out) :: dpdie
+
+  dpdie = 0.d0
+
+end subroutine EOS_GP_DPDIE
+
+subroutine EOS_GP_cs2(nelems, rho, cs2)
+
+  USE EOS_GP_Scalars
+
+  implicit none
+
+  CCTK_INT, intent(in) :: nelems
+  CCTK_REAL, dimension(nelems), intent(in) :: rho
+  CCTK_REAL, dimension(nelems), intent(out) :: cs2
+
+  cs2 = (p_geom_factor * eos_k_cgs * eos_gamma_local * rho_geom_factor_inv * &
+         (rho * rho_geom_factor_inv) ** (eos_gamma_local - 1.d0)) / &
+       (1.d0 + p_geom_factor * eos_k_cgs * eos_gamma_local * &
+       rho_geom_factor_inv / (eos_gamma_local - 1.d0) * &
+       (rho * rho_geom_factor_inv) ** (eos_gamma_local - 1.d0))
+  
+end subroutine EOS_GP_cs2
+
+ /*@@
+   @routine    Inverse functions
+   @date       Mon Apr  4 11:03:04 2005
+   @author     Ian Hawke
+   @desc 
+   Give P find the results...
+   @enddesc 
+   @calls     
+   @calledby   
+   @history 
+ 
+   @endhistory 
+
+@@*/
+
+
+subroutine EOS_GP_Inv_Rho(nelems, press, rho)
+
+  USE EOS_GP_Scalars
+
+  implicit none
+
+  CCTK_INT, intent(in) :: nelems
+  CCTK_REAL, dimension(nelems), intent(in) :: press
+  CCTK_REAL, dimension(nelems), intent(out) :: rho
+
+  rho = rho_geom_factor * &
+       (press / p_geom_factor / eos_k_cgs) ** (1.d0 / eos_gamma_local)
+
+end subroutine EOS_GP_Inv_Rho
+
+subroutine EOS_GP_Inv_IntEn(nelems, press, IntEn)
+
+  USE EOS_GP_Scalars
+
+  implicit none
+
+  CCTK_INT, intent(in) :: nelems
+  CCTK_REAL, dimension(nelems), intent(in) :: press
+  CCTK_REAL, dimension(nelems), intent(out) :: IntEn
+
+  IntEn = press / &
+       ( (eos_gamma_local - 1.d0) * rho_geom_factor * &
+       (press / p_geom_factor / eos_k_cgs) ** (1.d0 / eos_gamma_local) )
+
+end subroutine EOS_GP_Inv_IntEn
+
+subroutine EOS_GP_Inv_DPDRho(nelems, press, dpdrho)
+
+  USE EOS_GP_Scalars
+
+  implicit none
+
+  CCTK_INT, intent(in) :: nelems
+  CCTK_REAL, dimension(nelems), intent(in) :: press
+  CCTK_REAL, dimension(nelems), intent(out) :: dpdrho
+
+  dpdrho = p_geom_factor * eos_k_cgs * &
+           eos_gamma_local * rho_geom_factor_inv * &
+       (press / (p_geom_factor * eos_k_cgs) ) ** &
+       ( (eos_gamma_local - 1.d0) / eos_gamma_local )
+
+end subroutine EOS_GP_Inv_DPDRho
+
+subroutine EOS_GP_Inv_DPDIE(nelems, press, dpdie)
+
+  USE EOS_GP_Scalars
+
+  implicit none
+
+  CCTK_INT, intent(in) :: nelems
+  CCTK_REAL, dimension(nelems), intent(in) :: press
+  CCTK_REAL, dimension(nelems), intent(out) :: dpdie
+
+  dpdie = 0.d0
+
+end subroutine EOS_GP_Inv_DPDIE
+
+subroutine EOS_GP_Inv_cs2(nelems, press, cs2)
+
+  USE EOS_GP_Scalars
+
+  implicit none
+
+  CCTK_INT, intent(in) :: nelems
+  CCTK_REAL, dimension(nelems), intent(in) :: press
+  CCTK_REAL, dimension(nelems), intent(out) :: cs2
+
+  cs2 = (p_geom_factor * eos_k_cgs * eos_gamma_local * rho_geom_factor_inv * &
+         (press / (p_geom_factor * eos_k_cgs) ) ** &
+         ( (eos_gamma_local - 1.d0) / eos_gamma_local ) ) / &
+       (1.d0 + p_geom_factor * eos_k_cgs * eos_gamma_local * &
+       rho_geom_factor_inv / (eos_gamma_local - 1.d0) * &
+       (press / (p_geom_factor * eos_k_cgs) ) ** &
+       ( (eos_gamma_local - 1.d0) / eos_gamma_local ) )
+  
+end subroutine EOS_GP_Inv_cs2
+
diff --git a/src/EOS_GP.c b/src/EOS_GP.c
new file mode 100755
index 0000000..a9ea0b6
--- /dev/null
+++ b/src/EOS_GP.c
@@ -0,0 +1,140 @@
+#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_GP.h"
+
+void CCTK_FCALL CCTK_FNAME(EOS_GP_Pressure) (const CCTK_INT* nelems,
+                                             const CCTK_REAL* rho,
+                                             const CCTK_REAL* press);
+void CCTK_FCALL CCTK_FNAME(EOS_GP_IntEn)    (const CCTK_INT* nelems,
+                                             const CCTK_REAL* const rho,
+                                             const CCTK_REAL* inten);
+void CCTK_FCALL CCTK_FNAME(EOS_GP_DPDRho)   (const CCTK_INT* nelems,
+                                             const CCTK_REAL* const rho,
+                                             const CCTK_REAL* dpdrho);
+void CCTK_FCALL CCTK_FNAME(EOS_GP_DPDIE)    (const CCTK_INT* nelems,
+                                             const CCTK_REAL* const rho,
+                                             const CCTK_REAL* dpdie);
+void CCTK_FCALL CCTK_FNAME(EOS_GP_cs2)      (const CCTK_INT* nelems,
+                                             const CCTK_REAL* const rho,
+                                             const CCTK_REAL* cs2);
+
+CCTK_INT EOS_GP_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_GP_Pressure) (&n_elems,
+                                       (const CCTK_REAL*)indep_vars[0], 
+                                       (const CCTK_REAL*)dep_vars[0]);
+          break;
+        case (1):
+          CCTK_FNAME(EOS_GP_IntEn)    (&n_elems,
+                                       (const CCTK_REAL*)indep_vars[0], 
+                                       (const CCTK_REAL*)dep_vars[1]);
+          break;
+        case (2):
+          CCTK_FNAME(EOS_GP_DPDRho)   (&n_elems,
+                                       (const CCTK_REAL*)indep_vars[0], 
+                                       (const CCTK_REAL*)dep_vars[2]);
+          break;
+        case (3):
+          CCTK_FNAME(EOS_GP_DPDIE)    (&n_elems,
+                                       (const CCTK_REAL*)indep_vars[0], 
+                                       (const CCTK_REAL*)dep_vars[3]);
+          break;
+        case (4):
+          CCTK_FNAME(EOS_GP_cs2)      (&n_elems,
+                                       (const CCTK_REAL*)indep_vars[0], 
+                                       (const CCTK_REAL*)dep_vars[4]);
+          break;
+      }
+    }
+  }
+  
+  return ierr;
+}
+
+void CCTK_FCALL CCTK_FNAME(EOS_GP_Inv_Rho)      (const CCTK_INT* nelems,
+                                                 const CCTK_REAL* rho,
+                                                 const CCTK_REAL* press);
+void CCTK_FCALL CCTK_FNAME(EOS_GP_Inv_IntEn)    (const CCTK_INT* nelems,
+                                                 const CCTK_REAL* const rho,
+                                                 const CCTK_REAL* inten);
+void CCTK_FCALL CCTK_FNAME(EOS_GP_Inv_DPDRho)   (const CCTK_INT* nelems,
+                                                 const CCTK_REAL* const rho,
+                                                 const CCTK_REAL* dpdrho);
+void CCTK_FCALL CCTK_FNAME(EOS_GP_Inv_DPDIE)    (const CCTK_INT* nelems,
+                                                 const CCTK_REAL* const rho,
+                                                 const CCTK_REAL* dpdie);
+void CCTK_FCALL CCTK_FNAME(EOS_GP_Inv_cs2)      (const CCTK_INT* nelems,
+                                                 const CCTK_REAL* const rho,
+                                                 const CCTK_REAL* cs2);
+
+CCTK_INT EOS_GP_Inverse_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, dim;
+
+  for (var = 0; var < N_DEPS; ++var)
+  {
+    if (which_deps_to_set[var])
+    {
+      switch (var)
+      {
+        case (0):
+          CCTK_FNAME(EOS_GP_Inv_Rho)      (&n_elems,
+                                           (const CCTK_REAL*)indep_vars[0], 
+                                           (const CCTK_REAL*)dep_vars[0]);
+          break;
+        case (1):
+          CCTK_FNAME(EOS_GP_Inv_IntEn)    (&n_elems,
+                                           (const CCTK_REAL*)indep_vars[0], 
+                                           (const CCTK_REAL*)dep_vars[1]);
+          break;
+        case (2):
+          CCTK_FNAME(EOS_GP_Inv_DPDRho)   (&n_elems,
+                                           (const CCTK_REAL*)indep_vars[0], 
+                                           (const CCTK_REAL*)dep_vars[2]);
+          break;
+        case (3):
+          CCTK_FNAME(EOS_GP_Inv_DPDIE)    (&n_elems,
+                                           (const CCTK_REAL*)indep_vars[0], 
+                                           (const CCTK_REAL*)dep_vars[3]);
+          break;
+        case (4):
+          CCTK_FNAME(EOS_GP_Inv_cs2)      (&n_elems,
+                                           (const CCTK_REAL*)indep_vars[0], 
+                                           (const CCTK_REAL*)dep_vars[4]);
+          break;
+      }
+    }
+  }
+  
+  return ierr;
+}
diff --git a/src/EOS_GP.h b/src/EOS_GP.h
new file mode 100755
index 0000000..6fb5eb4
--- /dev/null
+++ b/src/EOS_GP.h
@@ -0,0 +1,21 @@
+#ifndef EOS_GP_H
+#define EOS_GP_H
+
+#include "cctk.h"
+
+#define N_INDEPS 1
+#define N_DEPS 5
+
+CCTK_INT EOS_GP_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);
+
+CCTK_INT EOS_GP_Inverse_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_GP_Scalars.F90 b/src/EOS_GP_Scalars.F90
new file mode 100755
index 0000000..c82e1d6
--- /dev/null
+++ b/src/EOS_GP_Scalars.F90
@@ -0,0 +1,21 @@
+ /*@@
+   @file      EOS_GP_Scalars.F90
+   @date      Mon Mar 14 17:39:50 2005
+   @author    Ian Hawke
+   @desc 
+   Constants for the EOS routines.
+   @enddesc 
+ @@*/
+
+#include "cctk.h"
+
+module EOS_GP_Scalars
+
+  implicit none
+
+  CCTK_REAL :: p_geom_factor, rho_geom_factor, &
+       rho_geom_factor_inv, eos_k_cgs, eos_gamma_local
+
+  CCTK_REAL ::  m_solar_geom, m_solar_cgs, c_cgs, G_cgs
+  
+end module EOS_GP_Scalars
diff --git a/src/EOS_GP_Setup.F90 b/src/EOS_GP_Setup.F90
new file mode 100755
index 0000000..4601e92
--- /dev/null
+++ b/src/EOS_GP_Setup.F90
@@ -0,0 +1,46 @@
+ /*@@
+   @file      EOS_GP_Setup.F90
+   @date      Mon Mar 14 17:41:20 2005
+   @author    Ian Hawke
+   @desc 
+   Setup the scalar variables
+   @enddesc 
+ @@*/
+
+#include "cctk.h"
+#include "cctk_Arguments.h"
+#include "cctk_Parameters.h"
+
+subroutine EOS_GP_Setup()
+
+  USE EOS_GP_Scalars
+
+  implicit none
+
+  DECLARE_CCTK_PARAMETERS
+
+  m_solar_cgs = 1.987d33
+  c_cgs = 2.9979d10
+  G_cgs = 6.6732d-8
+  m_solar_geom = G_cgs / c_cgs**2 * m_solar_cgs
+
+  eos_gamma_local = eos_gamma
+  
+  if (use_cgs .ne. 0) then
+    
+    p_geom_factor = G_cgs / c_cgs**4 * m_solar_geom**2
+    rho_geom_factor = p_geom_factor * c_cgs**2
+    rho_geom_factor_inv = 1.d0 / rho_geom_factor 
+    
+    eos_k_cgs = eos_k * rho_geom_factor**gamma_ini / p_geom_factor
+    
+  else
+    
+    p_geom_factor = 1.d0
+    rho_geom_factor = 1.d0
+    rho_geom_factor_inv = 1.d0
+    eos_k_cgs = eos_k
+    
+  end if
+
+end subroutine EOS_GP_Setup
diff --git a/src/EOS_GP_Setup.c b/src/EOS_GP_Setup.c
new file mode 100755
index 0000000..7bc3aa2
--- /dev/null
+++ b/src/EOS_GP_Setup.c
@@ -0,0 +1,59 @@
+#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_GP.h"
+
+void EOS_GP_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",
+                      "Independent variable names");
+  Util_TableSetString(table_handle,
+                      "Pressure SpecificInternalEnergy DPressureDRho DPressureDSpecificInternalEnergy c_s^2",
+                      "Dependent variable names");
+  Util_TableSetString(table_handle,
+                      "Polytrope",
+                      "EOS Name");
+
+  ierr = EOS_RegisterCall(table_handle,
+                          EOS_GP_SetArray);
+  if (ierr)
+  {
+    CCTK_WARN(0, "Failed to set up EOS_Polytrope call");
+  }
+
+  Util_TableSetString(table_handle,
+                      "Pressure",
+                      "Independent variable names");
+  Util_TableSetString(table_handle,
+                      "Rho SpecificInternalEnergy DPressureDRho DPressureDSpecificInternalEnergy c_s^2",
+                      "Dependent variable names");
+  Util_TableSetString(table_handle,
+                      "Inverse Polytrope",
+                      "EOS Name");
+
+  ierr = EOS_RegisterCall(table_handle,
+                          EOS_GP_Inverse_SetArray);
+  if (ierr)
+  {
+    CCTK_WARN(0, "Failed to set up EOS_Polytrope (inverse) call");
+  }
+}
diff --git a/src/make.code.defn b/src/make.code.defn
new file mode 100755
index 0000000..ba6d0b4
--- /dev/null
+++ b/src/make.code.defn
@@ -0,0 +1,12 @@
+# Main make.code.defn file for thorn EOS_GeneralPolytrope
+# $Header$
+
+# Source files in this directory
+SRCS = 	EOS_GP_Setup.c \
+	EOS_GP.c \
+	EOS_GP.F90 \
+	EOS_GP_Scalars.F90 \
+	EOS_GP_Setup.F90
+
+# Subdirectories containing source files
+SUBDIRS = 
diff --git a/src/make.code.deps b/src/make.code.deps
new file mode 100755
index 0000000..111d0db
--- /dev/null
+++ b/src/make.code.deps
@@ -0,0 +1,4 @@
+# Module dependencies for thorn EOS_GeneralPolytrope
+
+EOS_GP.F90.o:		EOS_GP_Scalars.F90.o
+EOS_GP_Setup.F90.o:	EOS_GP_Scalars.F90.o