microphysical_schemes.hpp

Go to the documentation of this file.

#ifndef SRC_MICROPHYSICS_MICROPHYSICAL_SCHEMES_HPP_
#define SRC_MICROPHYSICS_MICROPHYSICAL_SCHEMES_HPP_
 
// C/C++
#include <map>
#include <memory>
#include <string>
#include <vector>
 
// external
#include <yaml-cpp/yaml.h>
 
#include <Eigen/Core>
#include <Eigen/Dense>
 
// athena
#include <athena/athena.hpp>
#include <athena/mesh/mesh.hpp>
 
// canoe
#include <air_parcel.hpp>
#include <virtual_groups.hpp>
 
// microphysics
#include "chemistry_solver.hpp"
 
class MeshBlock;
class ParameterInput;
 
class MicrophysicalSchemeBase : public NamedGroup {
 public:  
  MicrophysicalSchemeBase(std::string name) : NamedGroup(name) {}
  virtual ~MicrophysicalSchemeBase() {}
 
 public:  
  virtual void AssembleReactionMatrix(AirParcel const &air, Real time) = 0;
 
  virtual void EvolveOneStep(AirParcel *air, Real time, Real dt) = 0;
 
 public:  
  virtual void SetVsedFromConserved(AthenaArray<Real> vsed[3],
                                    Hydro const *phydro, int kl, int ku, int jl,
                                    int ju, int il, int iu) = 0;
};
 
using MicrophysicalSchemePtr = std::shared_ptr<MicrophysicalSchemeBase>;
using AllMicrophysicalSchemes = std::vector<MicrophysicalSchemePtr>;
 
class MicrophysicalSchemesFactory {
 public:
  static AllMicrophysicalSchemes Create(MeshBlock *pmb, ParameterInput *pin);
};
 
template <int D>
class MicrophysicalScheme : public MicrophysicalSchemeBase,
                            public ParameterGroup,
                            public SpeciesIndexGroup {
 public:
  // needed for Eigen small matrix
  EIGEN_MAKE_ALIGNED_OPERATOR_NEW
 
  enum { Size = D };
  using RealVector = Eigen::Matrix<Real, D, 1>;
  using RealMatrix = Eigen::Matrix<Real, D, D>;
 
 public:  
  MicrophysicalScheme(std::string name, YAML::Node const &node)
      : MicrophysicalSchemeBase(name) {
    if (node["parameters"]) SetRealsFrom(node["parameters"]);
 
    if (node["dependent-species"]) {
      auto species = node["dependent-species"].as<std::vector<std::string>>();
      SetSpeciesIndex(species);
    }
  }
 
  virtual ~MicrophysicalScheme() {}
 
 public:  
  Real const *GetRatePtr() const { return rate_.data(); }
  Real const *GetJacobianPtr() const { return jacb_.data(); }
 
 public:  
  void SetVsedFromConserved(AthenaArray<Real> vsed[3], Hydro const *phydro,
                            int kl, int ku, int jl, int ju, int il,
                            int iu) override {
    for (auto n : GetCloudIndexArray())
      for (int k = kl; k <= ku; ++k)
        for (int j = jl; j <= ju; ++j)
          for (int i = il; i <= iu; ++i) {
            vsed[0](n, k, j, i) = 0.0;
            vsed[1](n, k, j, i) = 0.0;
            vsed[2](n, k, j, i) = 0.0;
          }
  }
 
 protected:
  Eigen::Matrix<Real, D, 1> rate_;
  Eigen::Matrix<Real, D, D> jacb_;
};
 
class Kessler94 : public MicrophysicalScheme<3> {
 public:
  using Base = MicrophysicalScheme<3>;
  enum { Size = Base::Size };
 
 public:  
  Kessler94(std::string name, YAML::Node const &node);
  ~Kessler94();
 
 public:  
  void AssembleReactionMatrix(AirParcel const &air, Real time) override;
  void EvolveOneStep(AirParcel *air, Real time, Real dt) override;
 
  void SetVsedFromConserved(AthenaArray<Real> vsed[3], Hydro const *phydro,
                            int kl, int ku, int jl, int ju, int il,
                            int iu) override;
 
 protected:
  ChemistrySolver<Size> solver_;
};
 
#endif  // SRC_MICROPHYSICS_MICROPHYSICAL_SCHEMES_HPP_