Move PeleRad source files to PelePhysics (#459)

* Move PeleRad source files to PelePhysics. --------- Co-authored-by: Bruce Perry <[email protected]> Co-authored-by: Wenjun Ge [email protected]
AMReX-Combustion · Dec 18, 2023 · 8e091fa · 8e091fa
1 parent 344c393
commit 8e091fa
Show file tree

Hide file tree

Showing 34 changed files with 4,667 additions and 0 deletions.
diff --git a/Source/Radiation/AMRParam.H b/Source/Radiation/AMRParam.H
@@ -0,0 +1,35 @@
+#ifndef AMR_PARAM_H
+#define AMR_PARAM_H
+
+#include <AMReX_ParmParse.H>
+#include <AMReX_Vector.H>
+#include <string>
+
+namespace PeleRad {
+
+struct AMRParam
+{
+public:
+  amrex::ParmParse pp_;
+
+  int max_level_;
+  int ref_ratio_;
+  int n_cell_;
+  int max_grid_size_;
+  int prob_type_;
+  std::string plot_file_name_;
+
+  AMREX_GPU_HOST
+  AMRParam(amrex::ParmParse const& pp) : pp_(pp)
+  {
+    pp_.query("max_level", max_level_);
+    pp_.query("ref_ratio", ref_ratio_);
+    pp_.query("n_cell", n_cell_);
+    pp_.query("max_grid_size", max_grid_size_);
+    pp_.query("prob_type", prob_type_);
+    pp.query("plot_file_name", plot_file_name_);
+  }
+};
+
+} // namespace PeleRad
+#endif
diff --git a/Source/Radiation/Constants.H b/Source/Radiation/Constants.H
@@ -0,0 +1,35 @@
+#ifndef CONSTANTS_H
+#define CONSTANTS_H
+
+namespace PeleRad {
+
+struct RadComps
+{
+  int co2Indx = -1;
+  int h2oIndx = -1;
+  int coIndx = -1;
+  int ch4Indx = -1;
+  int c2h4Indx = -1;
+
+  void checkIndices()
+  {
+    if (co2Indx > 0)
+      std::cout << " include radiative gas: CO2, the index is " << co2Indx
+                << std::endl;
+    if (h2oIndx > 0)
+      std::cout << " include radiative gas: H2O, the index is " << h2oIndx
+                << std::endl;
+    if (coIndx > 0)
+      std::cout << " include radiative gas: CO, the index is " << coIndx
+                << std::endl;
+    if (ch4Indx > 0)
+      std::cout << " include radiative gas: CH4, the index is " << coIndx
+                << std::endl;
+    if (c2h4Indx > 0)
+      std::cout << " include radiative gas: C2H4, the index is " << coIndx
+                << std::endl;
+  }
+};
+
+} // namespace PeleRad
+#endif
diff --git a/Source/Radiation/MLMGParam.H b/Source/Radiation/MLMGParam.H
@@ -0,0 +1,110 @@
+#ifndef MLMG_PARAM_H
+#define MLMG_PARAM_H
+
+#include <AMReX_ParmParse.H>
+
+namespace PeleRad {
+
+struct MLMGParam
+{
+public:
+  amrex::ParmParse pp_;
+
+  int verbose_;
+  int bottom_verbose_;
+  int max_iter_;
+  int max_fmg_iter_;
+  int max_bottom_iter_;
+  amrex::Real reltol_;
+  amrex::Real abstol_;
+  amrex::Real bottom_reltol_;
+  amrex::Real bottom_abstol_;
+  int linop_maxorder_;
+  int max_coarsening_level_;
+  int agg_grid_size_;
+  int con_grid_size_;
+  int agglomeration_;
+  int consolidation_;
+  int composite_solve_;
+  int fine_level_solve_only_;
+  bool use_hypre_;
+  amrex::Array<amrex::LinOpBCType, AMREX_SPACEDIM> lobc_;
+  amrex::Array<amrex::LinOpBCType, AMREX_SPACEDIM> hibc_;
+  int ebbc_type_;
+  std::string kppath_;
+
+  AMREX_GPU_HOST
+  MLMGParam()
+    : verbose_(0),
+      bottom_verbose_(0),
+      max_iter_(50),
+      max_bottom_iter_(20),
+      reltol_(1.e-4),
+      abstol_(1.e-4),
+      bottom_reltol_(1.e-6),
+      bottom_abstol_(1.e-6),
+      linop_maxorder_(3),
+      max_coarsening_level_(20),
+      agg_grid_size_(-1),
+      con_grid_size_(-1),
+      agglomeration_(0),
+      consolidation_(1),
+      composite_solve_(1),
+      fine_level_solve_only_(0),
+      use_hypre_(0),
+      ebbc_type_(2)
+  {
+  }
+
+  AMREX_GPU_HOST
+  MLMGParam(const amrex::ParmParse& pp) : pp_(pp)
+  {
+    pp_.query("kppath", kppath_);
+    pp_.query("verbose", verbose_);
+    pp_.query("bottom_verbose", bottom_verbose_);
+    pp_.query("max_iter", max_iter_);
+    pp_.query("max_fmg_iter", max_fmg_iter_);
+    pp_.query("max_bottom_iter", max_bottom_iter_);
+    pp_.query("reltol", reltol_);
+    pp_.query("abstol", abstol_);
+    pp_.query("bottom_reltol", bottom_reltol_);
+    pp_.query("bottom_abstol", bottom_abstol_);
+    pp_.query("linop_maxorder", linop_maxorder_);
+    pp_.query("max_coarsening_level", max_coarsening_level_);
+    pp_.query("agg_grid_size", agg_grid_size_);
+    pp_.query("con_grid_size", con_grid_size_);
+    pp_.query("agglomeration", agglomeration_);
+    pp_.query("consolidation", consolidation_);
+    pp_.query("composite_solve", composite_solve_);
+    pp_.query("fine_level_solve_only", fine_level_solve_only_);
+    pp_.query("use_hypre", use_hypre_);
+
+    amrex::Vector<std::string> lo_bc_char_(AMREX_SPACEDIM);
+    amrex::Vector<std::string> hi_bc_char_(AMREX_SPACEDIM);
+    pp_.getarr("lo_bc", lo_bc_char_, 0, AMREX_SPACEDIM);
+    pp_.getarr("hi_bc", hi_bc_char_, 0, AMREX_SPACEDIM);
+
+    std::unordered_map<std::string, amrex::LinOpBCType> bctype;
+    bctype["Dirichlet"] = amrex::LinOpBCType::Dirichlet;
+    bctype["Periodic"] = amrex::LinOpBCType::Periodic;
+    bctype["Neumann"] = amrex::LinOpBCType::Neumann;
+    bctype["Robin"] = amrex::LinOpBCType::Robin;
+
+    pp_.query("ebbc_type", ebbc_type_);
+
+    amrex::Print() << "Define radiation BC:"
+                   << "\n";
+
+    for (int idim = 0; idim < AMREX_SPACEDIM; idim++) {
+      lobc_[idim] = bctype[lo_bc_char_[idim]];
+      hibc_[idim] = bctype[hi_bc_char_[idim]];
+
+      amrex::Print() << "lobc[" << idim << "]=" << lobc_[idim] << ", ";
+      amrex::Print() << "hibc[" << idim << "]=" << hibc_[idim] << "\n";
+    }
+  }
+};
+
+} // namespace PeleRad
+
+#endif
diff --git a/Source/Radiation/Make.package b/Source/Radiation/Make.package
@@ -0,0 +1,4 @@
+CEXE_headers += AMRParam.H Constants.H MLMGParam.H POneMulti.H POneMultiEB.H POneMultiLevbyLev.H POneSingle.H POneSingleEB.H PeleCRad.H PeleLMRad.H PlanckMean.H SpectralModels.H
+
+VPATH_LOCATIONS += $(PELE_PHYSICS_HOME)/Source/Radiation
+INCLUDE_LOCATIONS += $(PELE_PHYSICS_HOME)/Source/Radiation
diff --git a/Source/Radiation/POneMulti.H b/Source/Radiation/POneMulti.H
@@ -0,0 +1,140 @@
+#ifndef PONEMULTI_H
+#define PONEMULTI_H
+
+#include <AMRParam.H>
+#include <AMReX_FArrayBox.H>
+#include <AMReX_MLABecLaplacian.H>
+#include <AMReX_MLMG.H>
+#include <AMReX_MultiFabUtil.H>
+#include <AMReX_ParmParse.H>
+#include <MLMGParam.H>
+
+namespace PeleRad {
+
+class POneMulti
+{
+private:
+  MLMGParam mlmgpp_;
+
+  amrex::LPInfo info_;
+
+public:
+  amrex::Vector<amrex::Geometry>& geom_;
+  amrex::Vector<amrex::BoxArray>& grids_;
+  amrex::Vector<amrex::DistributionMapping>& dmap_;
+
+  amrex::Vector<amrex::MultiFab>& solution_;
+  amrex::Vector<amrex::MultiFab> const& rhs_;
+  amrex::Vector<amrex::MultiFab> const& acoef_;
+  amrex::Vector<amrex::MultiFab> const& bcoef_;
+
+  amrex::Vector<amrex::MultiFab> const& robin_a_;
+  amrex::Vector<amrex::MultiFab> const& robin_b_;
+  amrex::Vector<amrex::MultiFab> const& robin_f_;
+
+  amrex::Real const ascalar_ = 1.0;
+  amrex::Real const bscalar_ = 1.0 / 3.0;
+
+  POneMulti() = delete;
+
+  // constructor
+  POneMulti(
+    MLMGParam const& mlmgpp,
+    amrex::Vector<amrex::Geometry>& geom,
+    amrex::Vector<amrex::BoxArray>& grids,
+    amrex::Vector<amrex::DistributionMapping>& dmap,
+    amrex::Vector<amrex::MultiFab>& solution,
+    amrex::Vector<amrex::MultiFab> const& rhs,
+    amrex::Vector<amrex::MultiFab> const& acoef,
+    amrex::Vector<amrex::MultiFab> const& bcoef,
+    amrex::Vector<amrex::MultiFab> const& robin_a,
+    amrex::Vector<amrex::MultiFab> const& robin_b,
+    amrex::Vector<amrex::MultiFab> const& robin_f)
+    : mlmgpp_(mlmgpp),
+      geom_(geom),
+      grids_(grids),
+      dmap_(dmap),
+      solution_(solution),
+      rhs_(rhs),
+      acoef_(acoef),
+      bcoef_(bcoef),
+      robin_a_(robin_a),
+      robin_b_(robin_b),
+      robin_f_(robin_f)
+  {
+    auto const max_coarsening_level = mlmgpp_.max_coarsening_level_;
+    auto const agglomeration = mlmgpp_.agglomeration_;
+    auto const consolidation = mlmgpp_.consolidation_;
+    info_.setAgglomeration(agglomeration);
+    info_.setConsolidation(consolidation);
+    info_.setMaxCoarseningLevel(max_coarsening_level);
+  }
+
+  void solve()
+  {
+    auto const nlevels = geom_.size();
+
+    auto const linop_maxorder = mlmgpp_.linop_maxorder_;
+
+    auto const& lobc = mlmgpp_.lobc_;
+    auto const& hibc = mlmgpp_.hibc_;
+
+    amrex::MLABecLaplacian mlabec(geom_, grids_, dmap_, info_);
+
+    mlabec.setDomainBC(lobc, hibc);
+    mlabec.setScalars(ascalar_, bscalar_);
+    mlabec.setMaxOrder(linop_maxorder);
+
+    auto const max_iter = mlmgpp_.max_iter_;
+    auto const max_fmg_iter = mlmgpp_.max_fmg_iter_;
+    auto const verbose = mlmgpp_.verbose_;
+    auto const bottom_verbose = mlmgpp_.bottom_verbose_;
+    auto const use_hypre = mlmgpp_.use_hypre_;
+    auto const bottom_reltol = mlmgpp_.bottom_reltol_;
+    auto const bottom_abstol = mlmgpp_.bottom_abstol_;
+
+    amrex::MLMG mlmg(mlabec);
+    mlmg.setMaxIter(max_iter);
+    mlmg.setMaxFmgIter(max_fmg_iter);
+    mlmg.setVerbose(verbose);
+    mlmg.setBottomVerbose(bottom_verbose);
+    mlmg.setBottomSolver(amrex::BottomSolver::bicgstab);
+    if (use_hypre)
+      mlmg.setBottomSolver(amrex::MLMG::BottomSolver::hypre);
+    mlmg.setBottomTolerance(bottom_reltol);
+    mlmg.setBottomToleranceAbs(bottom_abstol);
+
+    for (int ilev = 0; ilev < nlevels; ++ilev) {
+      auto const& geom = geom_[ilev];
+      auto& solution = solution_[ilev];
+      auto const& acoef = acoef_[ilev];
+      auto const& bcoef = bcoef_[ilev];
+      auto const& robin_a = robin_a_[ilev];
+      auto const& robin_b = robin_b_[ilev];
+      auto const& robin_f = robin_f_[ilev];
+
+      mlabec.setLevelBC(ilev, &solution, &robin_a, &robin_b, &robin_f);
+
+      mlabec.setACoeffs(ilev, acoef);
+
+      amrex::Array<amrex::MultiFab, AMREX_SPACEDIM> face_bcoef;
+      for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) {
+        amrex::BoxArray const& ba = amrex::convert(
+          bcoef.boxArray(), amrex::IntVect::TheDimensionVector(idim));
+        face_bcoef[idim].define(ba, bcoef.DistributionMap(), 1, 0);
+      }
+      amrex::average_cellcenter_to_face(GetArrOfPtrs(face_bcoef), bcoef, geom);
+      mlabec.setBCoeffs(ilev, amrex::GetArrOfConstPtrs(face_bcoef));
+    }
+
+    auto const tol_rel = mlmgpp_.reltol_;
+    auto const tol_abs = mlmgpp_.abstol_;
+
+    mlmg.solve(
+      GetVecOfPtrs(solution_), GetVecOfConstPtrs(rhs_), tol_rel, tol_abs);
+  }
+};
+
+} // namespace PeleRad
+
+#endif