From ccf11e833de2c2b4a6f3619be3862d1ebdf96ff6 Mon Sep 17 00:00:00 2001
From: Rene Gassmoeller <rene.gassmoeller@mailbox.org>
Date: Thu, 18 Jul 2019 11:42:15 -0700
Subject: [PATCH] Allow different ANEOS EoS

---
 src/EoS.h | 7 ++++++-
 src/GI.h  | 8 +++++---
 2 files changed, 11 insertions(+), 4 deletions(-)

diff --git a/src/EoS.h b/src/EoS.h
index 1658b9b..f429d7b 100644
--- a/src/EoS.h
+++ b/src/EoS.h
@@ -269,9 +269,14 @@ namespace EoS{
 	};
 }
 
+// These are hard-coded equations of state of certain materials. They
+// are useful for benchmarking purposes, or if you need exactly these EoS.
 static const EoS::IdealGas<PS::F64>  Monoatomic(5./3.);
 static const EoS::IdealGas<PS::F64>  Diatomic  (1.4);
 static const EoS::Tillotson<PS::F64> Granite   (2680.0, 16.0e+6, 3.5e+6, 18.00e+6,  18.0e+9,  18.0e+9, 0.5, 1.3, 5.0, 5.0);
 static const EoS::Tillotson<PS::F64> Iron      (7800.0,  9.5e+6, 2.4e+6 , 8.67e+6, 128.0e+9, 105.0e+9, 0.5, 1.5, 5.0, 5.0);
-static const EoS::ANEOS<PS::F64> AGranite      ("eos/granite.rho_u.txt");
+
+// This vector can hold any number of ANEOS equations of state.
+// The actual equations of state are later read from files by the specific problem setup.
+static std::vector<EoS::ANEOS<PS::F64>> ANEOS_EoS;
 
diff --git a/src/GI.h b/src/GI.h
index fbc05c4..143b2e5 100644
--- a/src/GI.h
+++ b/src/GI.h
@@ -31,6 +31,8 @@ template <class Ptcl> class GI : public Problem<Ptcl>{
         end_time = parameter_file.getValueOf("end_time",1.0e+4);
         damping = parameter_file.getValueOf("damping",1.);
         
+        ANEOS_EoS.push_back(EoS::ANEOS<PS::F64>("eos/granite.rho_u.txt"));
+
 		const PS::F64 Expand = 1.1;
 		const PS::F64 tarMass = UnitMass;
 		const PS::F64 tarRadi = UnitRadi;
@@ -190,7 +192,7 @@ template <class Ptcl> class GI : public Problem<Ptcl>{
                             ith.eng  = 0.1 * Grav * tarMass / tarRadi;
                             ith.id   = id++;
                             // TODO: Modify this line for all particles that need new EoS
-                            ith.setPressure(&AGranite);
+                            ith.setPressure(&ANEOS_EoS[0]);
                             ith.tag = 0;
                             if(ith.id / NptclIn1Node == PS::Comm::getRank()) tar.push_back(ith);
                         }
@@ -240,7 +242,7 @@ template <class Ptcl> class GI : public Problem<Ptcl>{
                             ith.eng  = 0.1 * Grav * tarMass / tarRadi;
                             ith.id   = id++;
                             // TODO: Modify this line for all particles that need new EoS
-                            ith.setPressure(&AGranite);
+                            ith.setPressure(&ANEOS_EoS[0]);
                             ith.tag = 2;
                             if(ith.id / NptclIn1Node == PS::Comm::getRank()) imp.push_back(ith);
                         }
@@ -289,7 +291,7 @@ template <class Ptcl> class GI : public Problem<Ptcl>{
 		for(PS::U64 i = 0 ; i < sph_system.getNumberOfParticleLocal() ; ++ i){
 			// TODO: Modify the lines below for all particles that need new EoS
 			if(sph_system[i].tag % 2 == 0){
-				sph_system[i].setPressure(&AGranite);
+				sph_system[i].setPressure(&ANEOS_EoS[0]);
 			}else{
 				sph_system[i].setPressure(&Iron);
 			}