Refactoring Impl (#97)

Clean up Impl

examples/2018-Li-harp/giants_re_pgen.cpp

@@ -19,6 +19,7 @@
 
 // canoe
 #include <configure.hpp>
+#include <dirty.hpp>
 #include <impl.hpp>
 #include <index_map.hpp>
 
@@ -91,7 +92,7 @@ void MeshBlock::ProblemGenerator(ParameterInput *pin) {
   // mesh limits
   Real x1min = pmy_mesh->mesh_size.x1min;
   Real x1max = pmy_mesh->mesh_size.x1max;
-  Real H0 = pimpl->GetPressureScaleHeight();
+  Real H0 = pcoord->GetPressureScaleHeight();
 
   // request temperature and pressure
   app->Log("request T", T0);

examples/2023-Li-saturn-vla/saturn_radio.cpp

@@ -24,6 +24,7 @@
 #include <air_parcel.hpp>
 #include <configure.hpp>
 #include <constants.hpp>
+#include <dirty.hpp>
 #include <impl.hpp>
 #include <index_map.hpp>
 
@@ -45,8 +46,12 @@
 #include <astro/planets.hpp>
 
 // harp
+#include <harp/radiation.hpp>
 #include <harp/radiation_band.hpp>
 
+// tracer
+#include <tracer/tracer.hpp>
+
 // opacity
 #include <opacity/Giants/microwave/mwr_absorbers.hpp>
 
@@ -134,7 +139,7 @@ void MeshBlock::ProblemGenerator(ParameterInput *pin) {
   // mesh limits
   Real x1min = pmy_mesh->mesh_size.x1min;
   Real x1max = pmy_mesh->mesh_size.x1max;
-  Real H0 = pimpl->GetPressureScaleHeight();
+  Real H0 = pcoord->GetPressureScaleHeight();
 
   // request temperature and pressure
   app->Log("request T", T0);

examples/2023-Li-uranus/uranus_mwr.cpp

@@ -24,6 +24,7 @@
 #include <air_parcel.hpp>
 #include <configure.hpp>
 #include <constants.hpp>
+#include <dirty.hpp>
 #include <impl.hpp>
 #include <index_map.hpp>
 
@@ -45,8 +46,12 @@
 #include <astro/planets.hpp>
 
 // harp
+#include <harp/radiation.hpp>
 #include <harp/radiation_band.hpp>
 
+// tracer
+#include <tracer/tracer.hpp>
+
 // opacity
 #include <opacity/Giants/microwave/mwr_absorbers.hpp>
 
@@ -130,7 +135,7 @@ void MeshBlock::ProblemGenerator(ParameterInput *pin) {
   // mesh limits
   Real x1min = pmy_mesh->mesh_size.x1min;
   Real x1max = pmy_mesh->mesh_size.x1max;
-  Real H0 = pimpl->GetPressureScaleHeight();
+  Real H0 = pcoord->GetPressureScaleHeight();
 
   // request temperature and pressure
   app->Log("request T", T0);

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

@@ -24,6 +24,7 @@
 #include <air_parcel.hpp>
 #include <configure.hpp>
 #include <constants.hpp>
+#include <dirty.hpp>
 #include <impl.hpp>
 #include <index_map.hpp>
 
@@ -44,8 +45,12 @@
 #include <astro/planets.hpp>
 
 // harp
+#include <harp/radiation.hpp>
 #include <harp/radiation_band.hpp>
 
+// tracer
+#include <tracer/tracer.hpp>
+
 // opacity
 #include <opacity/Giants/microwave/mwr_absorbers.hpp>
 
@@ -119,7 +124,7 @@ void MeshBlock::ProblemGenerator(ParameterInput *pin) {
   // mesh limits
   Real x1min = pmy_mesh->mesh_size.x1min;
   Real x1max = pmy_mesh->mesh_size.x1max;
-  Real H0 = pimpl->GetPressureScaleHeight();
+  Real H0 = pcoord->GetPressureScaleHeight();
 
   // request temperature and pressure
   app->Log("request T", T0);

patches/19.calculate_fluxes.patch

@@ -1,31 +1,32 @@
 diff --git a/src/hydro/calculate_fluxes.cpp b/src/hydro/calculate_fluxes.cpp
-index eeaf381b..679e6401 100644
+index eeaf381b..3b5d89b6 100644
 --- a/src/hydro/calculate_fluxes.cpp
 +++ b/src/hydro/calculate_fluxes.cpp
-@@ -24,6 +24,9 @@
+@@ -24,6 +24,10 @@
  #include "hydro.hpp"
  #include "hydro_diffusion/hydro_diffusion.hpp"
-
+ 
 +// snap injection
 +#include <impl.hpp>
++#include <snap/decomposition/decomposition.hpp>
 +
  // OpenMP header
  #ifdef OPENMP_PARALLEL
  #include <omp.h>
-@@ -74,6 +77,9 @@ void Hydro::CalculateFluxes(AthenaArray<Real> &w, FaceField &b,
+@@ -74,6 +78,9 @@ void Hydro::CalculateFluxes(AthenaArray<Real> &w, FaceField &b,
      }
    }
-
+ 
 +  // decompose pressure to pertubation pressure and hydrostatic pressure
 +  pmb->pimpl->pdec->ChangeToEntropy(w, kl, ku, jl, ju);
 +
    for (int k=kl; k<=ku; ++k) {
      for (int j=jl; j<=ju; ++j) {
        // reconstruct L/R states
-@@ -89,6 +95,9 @@ void Hydro::CalculateFluxes(AthenaArray<Real> &w, FaceField &b,
+@@ -89,6 +96,9 @@ void Hydro::CalculateFluxes(AthenaArray<Real> &w, FaceField &b,
          pmb->precon->PiecewiseParabolicX1(k, j, is-1, ie+1, w, bcc, wl_, wr_);
        }
-
+ 
 +      // assemble pressure pertubation
 +      pmb->pimpl->pdec->RestoreFromEntropy(w, wl_, wr_, k, j, is, ie + 1);
 +

src/air_parcel.cpp

@@ -11,6 +11,12 @@
 // microphysics
 #include "microphysics/microphysics.hpp"
 
+// chemistry
+#include "c3m/chemistry.hpp"
+
+// tracer
+#include "tracer/tracer.hpp"
+
 // snap
 #include "snap/thermodynamics/thermodynamics.hpp"
 

src/dirty.hpp

@@ -0,0 +1,16 @@
+//! \file dirty.hpp
+//! \brief dirty.hpp contains dirty functions that should go somewhere else
+
+#ifndef SRC_DIRTY_HPP_
+#define SRC_DIRTY_HPP_
+
+// helper functions, will be moved in the future
+int find_pressure_level_lesser(Real pres, AthenaArray<Real> const &w, int k,
+                               int j, int is, int ie) {
+  for (int i = is; i <= ie; ++i)
+    if (w(IPR, k, j, i) < pres) return i;
+
+  return ie + 1;
+}
+
+#endif  // SRC_DIRTY_HPP_

src/harp/set_spectral_properties.cpp

@@ -84,7 +84,7 @@ void RadiationBand::SetSpectralProperties(int k, int j, int il, int iu) {
 #ifdef HYDROSTATIC
       auto phydro = pmb->phydro;
       Real grav = -phydro->hsrc.GetG1();
-      Real H0 = pmb->pimpl->GetPressureScaleHeight();
+      Real H0 = pmb->pcoord->GetPressureScaleHeight();
       // TODO(cli) check this
       // \delta z = \delt Z * P/(\rho g H)
       tau_(m, i) *= pmb->pcoord->dx1f(i) * phydro->w(IPR, k, j, i) /

src/impl.cpp

@@ -15,7 +15,6 @@
 
 // inversion
 #include "inversion/inversion.hpp"
-#include "inversion/inversion_helper.hpp"
 
 // snap
 #include "snap/decomposition/decomposition.hpp"
@@ -26,6 +25,12 @@
 // microphysics
 #include "microphysics/microphysics.hpp"
 
+// c3m
+#include "c3m/chemistry.hpp"
+
+// tracer
+#include "tracer/tracer.hpp"
+
 // n-body
 #include "nbody/particles.hpp"
 
@@ -45,54 +50,27 @@ MeshBlock::Impl::Impl(MeshBlock *pmb, ParameterInput *pin) : pmy_block_(pmb) {
   // microphysics
   pmicro = std::make_shared<Microphysics>(pmb, pin);
 
+  // radiation
+  prad = std::make_shared<Radiation>(pmb, pin);
+
   // chemistry
   pchem = std::make_shared<Chemistry>(pmb, pin);
 
   // tracer
   ptracer = std::make_shared<Tracer>(pmb, pin);
 
   // turbulence
-  if (pin->DoesParameterExist("hydro", "turbulence")) {
-    std::string turbulence_model = pin->GetString("hydro", "turbulence");
-    if (turbulence_model == "none") {
-      pturb = std::make_shared<TurbulenceModel>(pmb, pin);
-    } else if (turbulence_model == "kepsilon") {
-      pturb = std::make_shared<KEpsilonTurbulence>(pmb, pin);
-      if (NTURBULENCE < 2)
-        throw NotImplementedError(
-            "NTURBULENCE must be at least 2 for k-epsilon model");
-    } else {
-      throw NotImplementedError(turbulence_model);
-    }
-  }
-
-  // radiation
-  prad = std::make_shared<Radiation>(pmb, pin);
+  pturb = TurbulenceFactory::CreateTurbulenceModel(pmb, pin);
 
   // inversion queue
-  // all_fits = Inversion::NewInversionQueue(pmb, pin);
-
-  // particles
-  all_particles = ParticlesFactory::create_all_particles(pmb, pin);
+  all_fits = InversionsFactory::CreateAllInversions(pmb, pin);
 
-#ifdef HYDROSTATIC
-  // reference pressure
-  reference_pressure_ = pin->GetReal("mesh", "ReferencePressure");
-
-  // pressure scale height
-  pressure_scale_height_ = pin->GetReal("mesh", "PressureScaleHeight");
-#else
-  reference_pressure_ = 1.0;
-  pressure_scale_height_ = 1.0;
-#endif  // HYDROSTATIC
+  // particle queue
+  all_particles = ParticlesFactory::CreateAllParticles(pmb, pin);
 }
 
 MeshBlock::Impl::~Impl() {}
 
-int find_pressure_level_lesser(Real pres, AthenaArray<Real> const &w, int k,
-                               int j, int is, int ie) {
-  for (int i = is; i <= ie; ++i)
-    if (w(IPR, k, j, i) < pres) return i;
-
-  return ie + 1;
+void MeshBlock::Impl::MapScalarsConserved(AthenaArray<Real> &s) {
+  if (NCLOUD > 0) pmicro->u.InitWithShallowSlice(s, 4, 0, NCLOUD);
 }

src/impl.hpp

@@ -11,31 +11,16 @@
 #include <athena/athena.hpp>
 #include <athena/mesh/mesh.hpp>
 
-// canoe
-#include "air_parcel.hpp"
-
-// snap
-#include "snap/decomposition/decomposition.hpp"
-#include "snap/implicit/implicit_solver.hpp"
-#include "snap/thermodynamics/thermodynamics.hpp"
-
-// harp
-#include "harp/radiation.hpp"
-
-// microphysics
-#include "microphysics/microphysics.hpp"
-
-// tracer
-#include "tracer/tracer.hpp"
-
-// chemistry
-#include "c3m/chemistry.hpp"
-
-// turbulence
-#include "snap/turbulence/turbulence_model.hpp"
-
 class ParameterInput;
 
+class Decomposition;
+class ImplicitSolver;
+class Microphysics;
+class Radiation;
+class Chemistry;
+class Tracer;
+class TurbulenceModel;
+
 class Inversion;
 class ParticleBase;
 
@@ -50,19 +35,17 @@ class MeshBlock::Impl {
   /// public data
   AthenaArray<Real> du;  // stores tendency
 
-  DecompositionPtr pdec;
-  ImplicitSolverPtr phevi;
+  std::shared_ptr<Decomposition> pdec;
+  std::shared_ptr<ImplicitSolver> phevi;
+  std::shared_ptr<Microphysics> pmicro;
+  std::shared_ptr<Radiation> prad;
+  std::shared_ptr<Chemistry> pchem;
+  std::shared_ptr<Tracer> ptracer;
+  std::shared_ptr<TurbulenceModel> pturb;
 
-  MicrophysicsPtr pmicro;
-  ChemistryPtr pchem;
-  TracerPtr ptracer;
   // StaticVariablePtr pstatic;
 
-  TurbulenceModelPtr pturb;
-
-  RadiationPtr prad;
-
-  // std::vector<std::shared_ptr<Inversion>> all_fits;
+  std::vector<std::shared_ptr<Inversion>> all_fits;
   std::vector<std::shared_ptr<ParticleBase>> all_particles;
 
   /// constructor and destructor
@@ -74,34 +57,21 @@ class MeshBlock::Impl {
     return boundary_exchanger_.at(name);
   }
 
-  auto GetMeshOutputGroups() const { return mesh_outputs_; }
-  auto GetFITSOutputGroups() const { return fits_outputs_; }
-
-  Real GetReferencePressure() const { return reference_pressure_; }
-  Real GetPressureScaleHeight() const { return pressure_scale_height_; }
+  auto &GetMeshOutputGroups() const { return mesh_outputs_; }
+  auto &GetFITSOutputGroups() const { return fits_outputs_; }
 
   // TODO(cli) : more needs to be changed
   // called in task_list/time_integration.cpp
-  void MapScalarsConserved(AthenaArray<Real> &s) {
-    if (NCLOUD > 0) pmicro->u.InitWithShallowSlice(s, 4, 0, NCLOUD);
-  }
+  void MapScalarsConserved(AthenaArray<Real> &s);
 
  protected:
   std::map<std::string, void *> boundary_exchanger_;
-
   std::vector<std::weak_ptr<MeshOutputGroup>> mesh_outputs_;
   std::vector<std::weak_ptr<FITSOutputGroup>> fits_outputs_;
 
-  Real reference_pressure_;
-  Real pressure_scale_height_;
-
  private:
-  MeshBlock *pmy_block_;
+  MeshBlock const *pmy_block_;
   int my_stage_;
 };
 
-// helper functions
-int find_pressure_level_lesser(Real pmax, AthenaArray<Real> const &w, int k,
-                               int j, int is, int ie);
-
 #endif  // SRC_IMPL_HPP_

src/index_map.hpp

@@ -2,7 +2,6 @@
 #define SRC_INDEX_MAP_HPP_
 
 // C/C++
-#include <iostream>
 #include <map>
 #include <string>