IndexMap is a singleton class (#43)

Good to go

.gitcustom/git-done

@@ -5,4 +5,4 @@ currentBranch=$(git symbolic-ref --short -q HEAD) # Getting the current branch
 git checkout main
 git fetch
 git rebase origin/main
-git branch -D currentBranch
+git branch -D $currentBranch

.github/workflows/main.yml

@@ -144,9 +144,6 @@ jobs:
     runs-on: [self-hosted, linux, 4core]
     needs: build-non-hydrostatic
     steps:
-      - name: source
-        run: source ~/.bashrc
-
       - uses: actions/checkout@v3
 
       - name: create build directory
@@ -175,9 +172,6 @@ jobs:
     runs-on: [self-hosted, Linux, 8core]
     needs: build-non-hydrostatic
     steps:
-      - name: source
-        run: source ~/.bashrc
-
       - uses: actions/checkout@v3
 
       - name: create build directory
@@ -190,6 +184,8 @@ jobs:
         run: cmake --build ${{github.workspace}}/build --config ${{env.BUILD_TYPE}} -- -j8
 
       - name: run simulation
+        env:
+          PATH: $PATH:/usr/lib64/mpich/bin
         run: cd ${{github.workspace}}/build/bin && mpiexec -n 8 ./robert.release -i robert2d.inp
 
       - name: combile outputs
@@ -206,9 +202,6 @@ jobs:
     runs-on: [self-hosted, Linux, 8core]
     needs: build-non-hydrostatic
     steps:
-      - name: source
-        run: source ~/.bashrc
-
       - uses: actions/checkout@v3
 
       - name: create build directory
@@ -221,6 +214,8 @@ jobs:
         run: cmake --build ${{github.workspace}}/build --config ${{env.BUILD_TYPE}} -- -j8
 
       - name: run simulation
+        env:
+          PATH: $PATH:/usr/lib64/mpich/bin
         run: cd ${{github.workspace}}/build/bin && mpiexec -n 8 ./robert.release -i robert3d.inp
 
       - name: combile outputs

examples/2023-jupiter-mwr-eq/juno.cpp

@@ -147,6 +147,7 @@ void MeshBlock::ProblemGenerator(ParameterInput *pin)
   Real cp = gamma/(gamma - 1.)*Rd;
 
   // index
+  auto pindex = IndexMap::GetInstance();
   int iH2O = pindex->GetVaporId("H2O");
   int iNH3 = pindex->GetVaporId("NH3");
 
@@ -218,24 +219,6 @@ void MeshBlock::ProblemGenerator(ParameterInput *pin)
       }
     }
 
-  // set tracers, electron and Na
-  int ielec = pindex->GetTracerId("e-");
-  int iNa = pindex->GetTracerId("Na");
-  auto ptracer = pimpl->ptracer;
-
-  for (int k = ks; k <= ke; ++k)
-    for (int j = js; j <= je; ++j)
-      for (int i = is; i <= ie; ++i) {
-        Real temp = pthermo->GetTemp(this,k,j,i);
-        Real pH2S = xH2S*phydro->w(IPR,k,j,i);
-        Real pNa = xNa*phydro->w(IPR,k,j,i);
-        Real svp = sat_vapor_p_Na_H2S_Visscher(temp, pH2S);
-        pNa = std::min(svp, pNa);
-
-        ptracer->u(iNa,k,j,i) = pNa/(Constants::kBoltz*temp);
-        ptracer->u(ielec,k,j,i) = saha_ionization_electron_density(temp, ptracer->u(iNa,k,j,i), 5.14);
-      }
-
   // if requested
   // modify atmospheric stability
   Real adlnTdlnP = pin->GetOrAddReal("problem", "adlnTdlnP", 0.);
@@ -275,6 +258,24 @@ void MeshBlock::ProblemGenerator(ParameterInput *pin)
         }
   }
 
+  // set tracers, electron and Na
+  int ielec = pindex->GetTracerId("e-");
+  int iNa = pindex->GetTracerId("Na");
+  auto ptracer = pimpl->ptracer;
+
+  for (int k = ks; k <= ke; ++k)
+    for (int j = js; j <= je; ++j)
+      for (int i = is; i <= ie; ++i) {
+        Real temp = pthermo->GetTemp(this,k,j,i);
+        Real pH2S = xH2S*phydro->w(IPR,k,j,i);
+        Real pNa = xNa*phydro->w(IPR,k,j,i);
+        Real svp = sat_vapor_p_Na_H2S_Visscher(temp, pH2S);
+        pNa = std::min(svp, pNa);
+
+        ptracer->u(iNa,k,j,i) = pNa/(Constants::kBoltz*temp);
+        ptracer->u(ielec,k,j,i) = saha_ionization_electron_density(temp, ptracer->u(iNa,k,j,i), 5.14);
+      }
+
   // primitive to conserved conversion
   peos->PrimitiveToConserved(phydro->w, pfield->bcc,
     phydro->u, pcoord, is, ie, js, je, ks, ke);

examples/2023-jupiter-mwr-eq/juno.inp

@@ -116,7 +116,6 @@ xNa         = 3.97E-6
 xKCl        = 0.
 metallicity = -2
 
-#rdlnTdlnP   = 1.015
-rdlnTdlnP   = 0.985
-rdlnTdlnP.pmin = 1.E5
-rdlnTdlnP.pmax = 20.E5
+adlnTdlnP   = 0.015
+adlnTdlnP.pmin = 1.E5
+adlnTdlnP.pmax = 20.E5

src/harp/absorber.cpp

@@ -1,6 +1,5 @@
 // athena
 #include <athena/athena.hpp>
-#include <athena/mesh/mesh.hpp>
 
 // canoe
 #include <index_map.hpp>
@@ -17,15 +16,16 @@ Absorber::Absorber(std::string name) : name_(name) {
   app->Log("Create Absorber " + name_);
 }
 
-Absorber::Absorber(MeshBlock* pmb, std::string name,
-                   std::vector<std::string> species,
+Absorber::Absorber(std::string name, std::vector<std::string> species,
                    std::map<std::string, Real> params)
     : name_(name), params_(params) {
   Application::Logger app("harp");
   app->Log("Create Absorber " + name_);
 
+  auto pindex = IndexMap::GetInstance();
+
   for (auto s : species) {
-    imols_.push_back(pmb->pindex->GetSpeciesId(s));
+    imols_.push_back(pindex->GetSpeciesId(s));
   }
 
   app->Log("Dependent species ids", imols_);

src/harp/absorber.hpp

@@ -24,15 +24,14 @@
 #include <athena/athena.hpp>
 
 class Variable;
-class MeshBlock;
 
 using ParameterMap = std::map<std::string, Real>;
 
 class Absorber {
  public:
   explicit Absorber(std::string name);
 
-  Absorber(MeshBlock* pmb, std::string name, std::vector<std::string> species,
+  Absorber(std::string name, std::vector<std::string> species,
            ParameterMap params);
 
   std::string GetName() const { return name_; }

src/harp/add_absorber_giants.cpp

@@ -45,32 +45,32 @@ void RadiationBand::addAbsorberGiants(ParameterInput *pin, YAML::Node &node) {
   if (type_ == "radio") {
     if (aname == "NH3") {
       auto ab = std::make_unique<gp::MwrAbsorberNH3>(
-          pmy_block_, species, ToParameterMap(node["parameters"]));
+          species, ToParameterMap(node["parameters"]));
 
       absorbers_.push_back(std::move(ab));
     } else if (aname == "H2O") {
       auto ab = std::make_unique<gp::MwrAbsorberH2O>(
-          pmy_block_, species, ToParameterMap(node["parameters"]));
+          species, ToParameterMap(node["parameters"]));
 
       absorbers_.push_back(std::move(ab));
     } else if (aname == "H2S") {
       auto ab = std::make_unique<gp::MwrAbsorberH2S>(
-          pmy_block_, species, ToParameterMap(node["parameters"]));
+          species, ToParameterMap(node["parameters"]));
 
       absorbers_.push_back(std::move(ab));
     } else if (aname == "PH3") {
       auto ab = std::make_unique<gp::MwrAbsorberPH3>(
-          pmy_block_, species, ToParameterMap(node["parameters"]));
+          species, ToParameterMap(node["parameters"]));
 
       absorbers_.push_back(std::move(ab));
     } else if (aname == "CIA") {
       auto ab = std::make_unique<gp::MwrAbsorberCIA>(
-          pmy_block_, species, ToParameterMap(node["parameters"]));
+          species, ToParameterMap(node["parameters"]));
 
       absorbers_.push_back(std::move(ab));
     } else if (aname == "Electron") {
       auto ab = std::make_unique<gp::MwrAbsorberElectron>(
-          pmy_block_, species, ToParameterMap(node["parameters"]));
+          species, ToParameterMap(node["parameters"]));
 
       absorbers_.push_back(std::move(ab));
     } else {

src/impl.cpp

@@ -24,10 +24,14 @@
 
 // canoe
 #include "impl.hpp"
+#include "index_map.hpp"
 
 MeshBlock::Impl::Impl(MeshBlock *pmb, ParameterInput *pin) : pmy_block_(pmb) {
   du.NewAthenaArray(NHYDRO, pmb->ncells3, pmb->ncells2, pmb->ncells1);
 
+  // index map
+  IndexMap::InitFromAthenaInput(pin);
+
   // thermodynamics
   Thermodynamics::InitFromAthenaInput(pin);
 
@@ -64,96 +68,89 @@ MeshBlock::Impl::Impl(MeshBlock *pmb, ParameterInput *pin) : pmy_block_(pmb) {
 #endif  // HYDROSTATIC
 }
 
-MeshBlock::Impl::~Impl() { Thermodynamics::Destroy(); }
+MeshBlock::Impl::~Impl() {
+  Thermodynamics::Destroy();
+  IndexMap::Destroy();
+}
 
-void MeshBlock::Impl::GatherFromPrimitive(Variable *var, int k, int j, int i) const {
-  auto& pmb = pmy_block_;
+void MeshBlock::Impl::GatherFromPrimitive(Variable *var, int k, int j,
+                                          int i) const {
+  auto &pmb = pmy_block_;
   auto mytype = var->GetType();
   var->SetType(Variable::Type::MassFrac);
 
-  for (int n = 0; n < NHYDRO; ++n)
-    var->w[n] = pmb->phydro->w(n, k, j, i);
+  for (int n = 0; n < NHYDRO; ++n) var->w[n] = pmb->phydro->w(n, k, j, i);
 
-  for (int n = 0; n < NCLOUD; ++n)
-    var->c[n] = pcloud->w(n, k, j, i);
+  for (int n = 0; n < NCLOUD; ++n) var->c[n] = pcloud->w(n, k, j, i);
 
-  for (int n = 0; n < NTRACER; ++n)
-    var->x[n] = ptracer->w(n, k, j, i);
+  for (int n = 0; n < NTRACER; ++n) var->x[n] = ptracer->w(n, k, j, i);
 
-  for (int n = 0; n < NCHEMISTRY; ++n)
-    var->q[n] = pchem->w(n, k, j, i);
+  for (int n = 0; n < NCHEMISTRY; ++n) var->q[n] = pchem->w(n, k, j, i);
 
   var->ConvertTo(mytype);
 }
 
-void MeshBlock::Impl::GatherFromConserved(Variable *var, int k, int j, int i) const {
-  auto& pmb = pmy_block_;
+void MeshBlock::Impl::GatherFromConserved(Variable *var, int k, int j,
+                                          int i) const {
+  auto &pmb = pmy_block_;
   auto mytype = var->GetType();
   var->SetType(Variable::Type::MassConc);
 
-  for (int n = 0; n < NHYDRO; ++n)
-    var->w[n] = pmb->phydro->u(n, k, j, i);
+  for (int n = 0; n < NHYDRO; ++n) var->w[n] = pmb->phydro->u(n, k, j, i);
 
-  for (int n = 0; n < NCLOUD; ++n)
-    var->c[n] = pcloud->u(n, k, j, i);
+  for (int n = 0; n < NCLOUD; ++n) var->c[n] = pcloud->u(n, k, j, i);
 
-  for (int n = 0; n < NTRACER; ++n)
-    var->x[n] = ptracer->u(n, k, j, i);
+  for (int n = 0; n < NTRACER; ++n) var->x[n] = ptracer->u(n, k, j, i);
 
-  for (int n = 0; n < NCHEMISTRY; ++n)
-    var->q[n] = pchem->u(n, k, j, i);
+  for (int n = 0; n < NCHEMISTRY; ++n) var->q[n] = pchem->u(n, k, j, i);
 
   var->ConvertTo(mytype);
 }
 
-void MeshBlock::Impl::DistributeToPrimitive(Variable const& var_in, int k, int j, int i) {
+void MeshBlock::Impl::DistributeToPrimitive(Variable const &var_in, int k,
+                                            int j, int i) {
   Variable *var;
-  auto& pmb = pmy_block_;
+  auto &pmb = pmy_block_;
 
   if (var_in.GetType() != Variable::Type::MassFrac) {
     var = new Variable(var_in);
     var->ConvertToMassFraction();
   } else {
-    var = const_cast<Variable*>(&var_in);
+    var = const_cast<Variable *>(&var_in);
   }
 
   for (int n = 0; n < NHYDRO; ++n) pmb->phydro->w(n, k, j, i) = var->w[n];
 
-  for (int n = 0; n < NCLOUD; ++n)
-    pcloud->w(n, k, j, i) = var->c[n];
+  for (int n = 0; n < NCLOUD; ++n) pcloud->w(n, k, j, i) = var->c[n];
 
-  for (int n = 0; n < NCHEMISTRY; ++n)
-    pchem->w(n, k, j, i) = var->q[n];
+  for (int n = 0; n < NCHEMISTRY; ++n) pchem->w(n, k, j, i) = var->q[n];
 
-  for (int n = 0; n < NTRACER; ++n)
-    ptracer->w(n, k, j, i) = var->x[n];
+  for (int n = 0; n < NTRACER; ++n) ptracer->w(n, k, j, i) = var->x[n];
 
   if (var_in.GetType() != Variable::Type::MassFrac) {
     delete var;
   }
 }
 
-void MeshBlock::Impl::DistributeToConserved(Variable const& var_in, int k, int j, int i) {
+void MeshBlock::Impl::DistributeToConserved(Variable const &var_in, int k,
+                                            int j, int i) {
   Variable *var;
-  auto& pmb = pmy_block_;
+  auto &pmb = pmy_block_;
 
   if (var_in.GetType() != Variable::Type::MassConc) {
     var = new Variable(var_in);
     var->ConvertToMassConcentration();
   } else {
-    var = const_cast<Variable*>(&var_in);
+    var = const_cast<Variable *>(&var_in);
   }
 
   for (int n = 0; n < NHYDRO; ++n) pmb->phydro->u(n, k, j, i) = var->w[n];
 
-  for (int n = 0; n < NCLOUD; ++n)
-    pcloud->u(n, k, j, i) = var->c[n];
+  for (int n = 0; n < NCLOUD; ++n) pcloud->u(n, k, j, i) = var->c[n];
 
-  for (int n = 0; n < NCHEMISTRY; ++n)
-    pchem->u(n, k, j, i) = var->q[n];
+  for (int n = 0; n < NCHEMISTRY; ++n) pchem->u(n, k, j, i) = var->q[n];
 
-  for (int n = 0; n < NTRACER; ++n)
-    ptracer->u(n, k, j, i) = var->x[n];
+  for (int n = 0; n < NTRACER; ++n) ptracer->u(n, k, j, i) = var->x[n];
 
   if (var_in.GetType() != Variable::Type::MassConc) {
     delete var;