diff --git a/Source/Radiation/AMRParam.H b/Source/Radiation/AMRParam.H new file mode 100644 index 000000000..55546d309 --- /dev/null +++ b/Source/Radiation/AMRParam.H @@ -0,0 +1,35 @@ +#ifndef AMR_PARAM_H +#define AMR_PARAM_H + +#include +#include +#include + +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 new file mode 100644 index 000000000..507408b28 --- /dev/null +++ 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 new file mode 100644 index 000000000..422ca9361 --- /dev/null +++ b/Source/Radiation/MLMGParam.H @@ -0,0 +1,110 @@ +#ifndef MLMG_PARAM_H +#define MLMG_PARAM_H + +#include + +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 lobc_; + amrex::Array 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 lo_bc_char_(AMREX_SPACEDIM); + amrex::Vector 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 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 new file mode 100644 index 000000000..e71ea3e09 --- /dev/null +++ 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 new file mode 100644 index 000000000..3f7b2c959 --- /dev/null +++ b/Source/Radiation/POneMulti.H @@ -0,0 +1,140 @@ +#ifndef PONEMULTI_H +#define PONEMULTI_H + +#include +#include +#include +#include +#include +#include +#include + +namespace PeleRad { + +class POneMulti +{ +private: + MLMGParam mlmgpp_; + + amrex::LPInfo info_; + +public: + amrex::Vector& geom_; + amrex::Vector& grids_; + amrex::Vector& dmap_; + + amrex::Vector& solution_; + amrex::Vector const& rhs_; + amrex::Vector const& acoef_; + amrex::Vector const& bcoef_; + + amrex::Vector const& robin_a_; + amrex::Vector const& robin_b_; + amrex::Vector 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& geom, + amrex::Vector& grids, + amrex::Vector& dmap, + amrex::Vector& solution, + amrex::Vector const& rhs, + amrex::Vector const& acoef, + amrex::Vector const& bcoef, + amrex::Vector const& robin_a, + amrex::Vector const& robin_b, + amrex::Vector 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 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 diff --git a/Source/Radiation/POneMultiEB.H b/Source/Radiation/POneMultiEB.H new file mode 100644 index 000000000..803905f38 --- /dev/null +++ b/Source/Radiation/POneMultiEB.H @@ -0,0 +1,149 @@ +#ifndef PONEMULTIEB_H +#define PONEMULTIEB_H + +#include +#include +#include +#include +#include +#include +#include + +namespace PeleRad { + +class POneMultiEB +{ +private: + MLMGParam mlmgpp_; + + amrex::LPInfo info_; + +public: + amrex::Vector& geom_; + amrex::Vector& grids_; + amrex::Vector& dmap_; + + amrex::Vector& factory_; + + amrex::Vector& solution_; + amrex::Vector const& rhs_; + amrex::Vector const& acoef_; + amrex::Vector const& bcoef_; + + amrex::Vector const& robin_a_; + amrex::Vector const& robin_b_; + amrex::Vector const& robin_f_; + + amrex::Real const ascalar_ = 1.0; + amrex::Real const bscalar_ = 1.0 / 3.0; + + // constructor + POneMultiEB( + MLMGParam const& mlmgpp, + amrex::Vector& geom, + amrex::Vector& grids, + amrex::Vector& dmap, + amrex::Vector& factory, + amrex::Vector& solution, + amrex::Vector const& rhs, + amrex::Vector const& acoef, + amrex::Vector const& bcoef, + amrex::Vector const& robin_a, + amrex::Vector const& robin_b, + amrex::Vector const& robin_f) + : mlmgpp_(mlmgpp), + geom_(geom), + grids_(grids), + dmap_(dmap), + factory_(factory), + 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::MLEBABecLap mlabec( + // geom_, grids_, dmap_, info_, + // amrex::GetVecOfConstPtrs(factory_)); + + amrex::MLEBABecLap mlabec(geom_, grids_, dmap_, info_, factory_); + + 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 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)); + + if (mlmgpp_.ebbc_type_ == 1) { + mlabec.setEBDirichlet(ilev, solution, 1.0); + } + } + 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 diff --git a/Source/Radiation/POneMultiLevbyLev.H b/Source/Radiation/POneMultiLevbyLev.H new file mode 100644 index 000000000..acbd96aea --- /dev/null +++ b/Source/Radiation/POneMultiLevbyLev.H @@ -0,0 +1,154 @@ +#ifndef PONEMULTILEVBYLEV_H +#define PONEMULTILEVBYLEV_H + +#include +#include +#include +#include +#include +#include +#include + +namespace PeleRad { + +class POneMultiLevbyLev +{ +private: + MLMGParam mlmgpp_; + + amrex::LPInfo info_; + + int ref_ratio_; + +public: + amrex::Vector& geom_; + amrex::Vector& grids_; + amrex::Vector& dmap_; + + amrex::Vector& solution_; + amrex::Vector const& rhs_; + amrex::Vector const& acoef_; + amrex::Vector const& bcoef_; + + amrex::Vector const& robin_a_; + amrex::Vector const& robin_b_; + amrex::Vector const& robin_f_; + + amrex::Real const ascalar_ = 1.0; + amrex::Real const bscalar_ = 1.0 / 3.0; + + POneMultiLevbyLev() = delete; + + // constructor + POneMultiLevbyLev( + MLMGParam const& mlmgpp, + int const ref_ratio, + amrex::Vector& geom, + amrex::Vector& grids, + amrex::Vector& dmap, + amrex::Vector& solution, + amrex::Vector const& rhs, + amrex::Vector const& acoef, + amrex::Vector const& bcoef, + amrex::Vector const& robin_a, + amrex::Vector const& robin_b, + amrex::Vector const& robin_f) + : mlmgpp_(mlmgpp), + ref_ratio_(ref_ratio), + 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() + { + int const solver_level = 0; + auto const nlevels = geom_.size(); + auto const tol_rel = mlmgpp_.reltol_; + auto const tol_abs = mlmgpp_.abstol_; + + auto const linop_maxorder = mlmgpp_.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& lobc = mlmgpp_.lobc_; + auto const& hibc = mlmgpp_.hibc_; + + for (int ilev = 0; ilev < nlevels; ++ilev) { + auto const& geom = geom_[ilev]; + auto& solution = solution_[ilev]; + auto const& rhs = rhs_[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]; + + amrex::MLABecLaplacian mlabeclev( + {geom_[ilev]}, {grids_[ilev]}, {dmap_[ilev]}, info_); + + mlabeclev.setMaxOrder(linop_maxorder); + + mlabeclev.setDomainBC(lobc, hibc); + + if (ilev > 0) { + mlabeclev.setCoarseFineBC(&solution_[ilev - 1], ref_ratio_); + } + + mlabeclev.setLevelBC( + solver_level, &solution, &robin_a, &robin_b, &robin_f); + + mlabeclev.setScalars(ascalar_, bscalar_); + + mlabeclev.setACoeffs(solver_level, acoef); + + amrex::Array 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); + + mlabeclev.setBCoeffs(solver_level, amrex::GetArrOfConstPtrs(face_bcoef)); + + amrex::MLMG mlmglev(mlabeclev); + mlmglev.setMaxIter(max_iter); + mlmglev.setMaxFmgIter(max_fmg_iter); + mlmglev.setBottomMaxIter(100); + + mlmglev.setVerbose(verbose); + mlmglev.setBottomVerbose(bottom_verbose); + + mlmglev.setBottomSolver(amrex::BottomSolver::bicgstab); + + if (use_hypre) { + mlmglev.setBottomSolver(amrex::BottomSolver::hypre); + // mlmglev.setHypreInterface(hypre_interface); + } + + mlmglev.solve({&solution}, {&rhs}, tol_rel, tol_abs); + } + } +}; + +} // namespace PeleRad +#endif diff --git a/Source/Radiation/POneSingle.H b/Source/Radiation/POneSingle.H new file mode 100644 index 000000000..0197fa885 --- /dev/null +++ b/Source/Radiation/POneSingle.H @@ -0,0 +1,145 @@ +#ifndef PONESINGLE_H +#define PONESINGLE_H + +#include +#include +#include +#include +#include +#include +#include + +namespace PeleRad { + +class POneSingle +{ +private: + MLMGParam mlmgpp_; + +public: + amrex::Geometry const& geom_; + amrex::BoxArray const& grids_; + amrex::DistributionMapping const& dmap_; + + amrex::MultiFab& solution_; + amrex::MultiFab const& rhs_; + amrex::MultiFab const& acoef_; + amrex::MultiFab const& bcoef_; + + amrex::MultiFab const& robin_a_; + amrex::MultiFab const& robin_b_; + amrex::MultiFab const& robin_f_; + + amrex::Real const ascalar = 1.0; + amrex::Real const bscalar = 1.0 / 3.0; + + // constructor + POneSingle( + MLMGParam const& mlmgpp, + amrex::Geometry const& geom, + amrex::BoxArray const& grids, + amrex::DistributionMapping const& dmap, + amrex::MultiFab& solution, + amrex::MultiFab const& rhs, + amrex::MultiFab const& acoef, + amrex::MultiFab const& bcoef, + amrex::MultiFab const& robin_a, + amrex::MultiFab const& robin_b, + 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){}; + + void solve() + { + auto const max_coarsening_level = mlmgpp_.max_coarsening_level_; + 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 tol_rel = mlmgpp_.reltol_; + auto const tol_abs = mlmgpp_.abstol_; + auto const use_hypre = mlmgpp_.use_hypre_; + + auto const& lobc = mlmgpp_.lobc_; + auto const& hibc = mlmgpp_.hibc_; + + auto const& geom = geom_; + auto const& grids = grids_; + auto const& dmap = dmap_; + + auto& solution = solution_; + auto const& rhs = rhs_; + auto const& acoef = acoef_; + auto const& bcoef = bcoef_; + auto const& robin_a = robin_a_; + auto const& robin_b = robin_b_; + auto const& robin_f = robin_f_; + + amrex::MLABecLaplacian mlabec( + {geom}, {grids}, {dmap}, + amrex::LPInfo().setMaxCoarseningLevel(max_coarsening_level)); + + mlabec.setDomainBC(lobc, hibc); + + mlabec.setLevelBC(0, &solution, &robin_a, &robin_b, &robin_f); + + mlabec.setScalars(ascalar, bscalar); + + mlabec.setACoeffs(0, acoef); + + amrex::Array face_bcoef; + + for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) { + const amrex::BoxArray& 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(0, amrex::GetArrOfConstPtrs(face_bcoef)); + + amrex::MLMG mlmg(mlabec); + mlmg.setMaxIter(max_iter); + mlmg.setMaxFmgIter(max_fmg_iter); + mlmg.setVerbose(verbose); + mlmg.setBottomVerbose(bottom_verbose); + + if (use_hypre) + mlmg.setBottomSolver(amrex::MLMG::BottomSolver::hypre); + + mlmg.solve({&solution}, {&rhs}, tol_rel, tol_abs); + } + + void calcRadSource(amrex::MultiFab& rad_src) + { + for (amrex::MFIter mfi(rad_src); mfi.isValid(); ++mfi) { + amrex::Box const& bx = mfi.validbox(); + + auto const& rhsfab = rhs_.array(mfi); + auto const& solfab = solution_.array(mfi); + auto const& acfab = acoef_.array(mfi); + + auto radfab = rad_src.array(mfi); + + amrex::ParallelFor( + bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { + radfab(i, j, k, 4) = + acfab(i, j, k) * solfab(i, j, k) - rhsfab(i, j, k); + }); + } + } +}; + +} // namespace PeleRad + +#endif diff --git a/Source/Radiation/POneSingleEB.H b/Source/Radiation/POneSingleEB.H new file mode 100644 index 000000000..a74d9c455 --- /dev/null +++ b/Source/Radiation/POneSingleEB.H @@ -0,0 +1,158 @@ +#ifndef PONESINGLE_H +#define PONESINGLE_H + +#include +#include +#include +#include +#include +#include +#include +#include +#include + +namespace PeleRad { + +class POneSingleEB +{ +private: + MLMGParam mlmgpp_; + +public: + amrex::Geometry const& geom_; + amrex::BoxArray const& grids_; + amrex::DistributionMapping const& dmap_; + + std::unique_ptr const& factory_; + + amrex::MultiFab& solution_; + amrex::MultiFab const& rhs_; + amrex::MultiFab const& acoef_; + amrex::MultiFab const& bcoef_; + + amrex::MultiFab const& robin_a_; + amrex::MultiFab const& robin_b_; + amrex::MultiFab const& robin_f_; + + amrex::Real const ascalar = 1.0; + amrex::Real const bscalar = 1.0 / 3.0; + + // constructor + POneSingleEB( + MLMGParam const& mlmgpp, + amrex::Geometry const& geom, + amrex::BoxArray const& grids, + amrex::DistributionMapping const& dmap, + std::unique_ptr const& factory, + amrex::MultiFab& solution, + amrex::MultiFab const& rhs, + amrex::MultiFab const& acoef, + amrex::MultiFab const& bcoef, + amrex::MultiFab const& robin_a, + amrex::MultiFab const& robin_b, + amrex::MultiFab const& robin_f) + : mlmgpp_(mlmgpp), + geom_(geom), + grids_(grids), + dmap_(dmap), + factory_(factory), + solution_(solution), + rhs_(rhs), + acoef_(acoef), + bcoef_(bcoef), + robin_a_(robin_a), + robin_b_(robin_b), + robin_f_(robin_f){}; + + void solve() + { + amrex::MultiFab const& vfrc = factory_->getVolFrac(); + + auto const max_coarsening_level = mlmgpp_.max_coarsening_level_; + 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 tol_rel = mlmgpp_.reltol_; + auto const tol_abs = mlmgpp_.abstol_; + auto const use_hypre = mlmgpp_.use_hypre_; + + auto const& lobc = mlmgpp_.lobc_; + auto const& hibc = mlmgpp_.hibc_; + + auto const& geom = geom_; + auto const& grids = grids_; + auto const& dmap = dmap_; + auto const& factory = factory_; + + auto& solution = solution_; + auto const& rhs = rhs_; + auto const& acoef = acoef_; + auto const& bcoef = bcoef_; + auto const& robin_a = robin_a_; + auto const& robin_b = robin_b_; + auto const& robin_f = robin_f_; + + amrex::MLEBABecLap mlabec( + {geom}, {grids}, {dmap}, + amrex::LPInfo().setMaxCoarseningLevel(max_coarsening_level), + {factory.get()}); + + mlabec.setDomainBC(lobc, hibc); + + mlabec.setLevelBC(0, &solution, &robin_a, &robin_b, &robin_f); + + mlabec.setScalars(ascalar, bscalar); + + mlabec.setACoeffs(0, acoef); + + amrex::Array 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(0, amrex::GetArrOfConstPtrs(face_bcoef)); + + if (mlmgpp_.ebbc_type_ == 1) { + mlabec.setEBDirichlet(0, solution, 1.0); + } + + amrex::MLMG mlmg(mlabec); + mlmg.setMaxIter(max_iter); + mlmg.setMaxFmgIter(max_fmg_iter); + mlmg.setVerbose(verbose); + mlmg.setBottomVerbose(bottom_verbose); + + if (use_hypre) + mlmg.setBottomSolver(amrex::MLMG::BottomSolver::hypre); + + mlmg.solve({&solution}, {&rhs}, tol_rel, tol_abs); + } + + void calcRadSource(amrex::MultiFab& rad_src) + { + for (amrex::MFIter mfi(rad_src); mfi.isValid(); ++mfi) { + amrex::Box const& bx = mfi.validbox(); + + auto const& rhsfab = rhs_.array(mfi); + auto const& solfab = solution_.array(mfi); + auto const& acfab = acoef_.array(mfi); + + auto radfab = rad_src.array(mfi); + + amrex::ParallelFor( + bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { + radfab(i, j, k, 4) = + acfab(i, j, k) * solfab(i, j, k) - rhsfab(i, j, k); + }); + } + } +}; + +} // namespace PeleRad + +#endif diff --git a/Source/Radiation/PeleCRad.H b/Source/Radiation/PeleCRad.H new file mode 100644 index 000000000..740fcb577 --- /dev/null +++ b/Source/Radiation/PeleCRad.H @@ -0,0 +1,201 @@ +#ifndef RADIATION_H +#define RADIATION_H + +#include +#include +#include +#include +#include +#include + +namespace PeleRad { + +class Radiation +{ +private: + // AMRParam amrpp_; + // MLMGParam mlmgpp_; + + // amrex::ParmParse const& pp_; + + PlanckMean radprop; + + amrex::Geometry const& geom_; + amrex::BoxArray const& grids_; + amrex::DistributionMapping const& dmap_; + + amrex::MultiFab solution_; + amrex::MultiFab rhs_; + amrex::MultiFab acoef_; + amrex::MultiFab bcoef_; + amrex::MultiFab robin_a_; + amrex::MultiFab robin_b_; + amrex::MultiFab robin_f_; + + amrex::MultiFab absc_; + + RadComps rc_; + +public: + AMREX_GPU_HOST + Radiation( + amrex::Geometry const& geom, + amrex::BoxArray const& grids, + amrex::DistributionMapping const& dmap, + RadComps rc) + : geom_(geom), grids_(grids), dmap_(dmap), rc_(rc) + { + if (amrex::ParallelDescriptor::IOProcessor()) + rc_.checkIndices(); + solution_.define( + grids_, dmap_, 1, 1); // one ghost cell to store boundary conditions + rhs_.define(grids_, dmap_, 1, 0); + acoef_.define(grids_, dmap_, 1, 0); + bcoef_.define(grids_, dmap_, 1, 1); // one ghost cell for averaging to faces + robin_a_.define(grids_, dmap_, 1, 1); + robin_b_.define(grids_, dmap_, 1, 1); + robin_f_.define(grids_, dmap_, 1, 1); + absc_.define(grids_, dmap_, 1, 0); + + solution_.setVal(0.0, 0, 1, amrex::IntVect(0)); + + loadSpecModel(); + } + + AMREX_GPU_HOST + void loadSpecModel() + { + std::string data_path; + radprop.load(data_path); + amrex::Print() << "The radiative property database is loaded" + << "\n"; + } + + void readRadParams(amrex::ParmParse const& pp) {} + + void updateSpecProp( + amrex::MFIter const& mfi, + amrex::Array4 const& Yco2, + amrex::Array4 const& Yh2o, + amrex::Array4 const& Yco, + amrex::Array4 const& T, + amrex::Array4 const& P +#ifdef PELEC_USE_SOOT + , + amrex::Array4 const& fv +#endif + ) + { + // std::cout << "update radiative properties" << std::endl; + auto const& kpco2 = radprop.kpco2(); + auto const& kph2o = radprop.kph2o(); + auto const& kpco = radprop.kpco(); + + amrex::Box const& bx = mfi.validbox(); + amrex::Box const& gbx = mfi.growntilebox(1); + + auto const dlo = amrex::lbound(geom_.Domain()); + auto const dhi = amrex::ubound(geom_.Domain()); + + auto kappa = absc_.array(mfi); + auto rhsfab = rhs_.array(mfi); + auto alphafab = acoef_.array(mfi); + auto betafab = bcoef_.array(mfi); + auto robin_a_fab = robin_a_.array(mfi); + auto robin_b_fab = robin_b_.array(mfi); + auto robin_f_fab = robin_f_.array(mfi); + + amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { + RadProp::getRadPropGas( + i, j, k, Yco2, Yh2o, Yco, T, P, kappa, kpco2, kph2o, kpco); + kappa(i, j, k) *= 0.1; // correction for P in cgs + }); + +#ifdef PELEC_USE_SOOT + auto const& kpsoot = radprop.kpsoot(); + amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { + RadProp::getRadPropSoot(i, j, k, fv, T, kappa, kpsoot); + }); +#endif + + amrex::ParallelFor(gbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { + betafab(i, j, k) = 1.0 / 0.001; + + if (bx.contains(i, j, k)) { + double ka = std::max(0.001, kappa(i, j, k)); + betafab(i, j, k) = 1.0 / ka; + // rhsfab(i, j, k) + // = 4.0 * ka * 5.67e-8 * std::pow(T(i, j, k), 4.0); //si + rhsfab(i, j, k) = 4.0 * ka * 5.67e-5 * std::pow(T(i, j, k), 4.0); // cgs + alphafab(i, j, k) = ka; + } + + // Robin BC + bool robin_cell = false; + if (j >= dlo.y && j <= dhi.y && k >= dlo.z && k <= dhi.z) { + if (i > dhi.x || i < dlo.x) { + robin_cell = true; + } + } else if (i >= dlo.x && i <= dhi.x && k >= dlo.z && k <= dhi.z) { + if (j > dhi.y || j < dlo.y) { + robin_cell = true; + } + } else if (i >= dlo.x && i <= dhi.x && j >= dlo.y && j <= dhi.y) { + if (k > dhi.z || k < dlo.z) { + robin_cell = true; + } + } + + if (robin_cell) { + robin_a_fab(i, j, k) = -1.0 / betafab(i, j, k); + robin_b_fab(i, j, k) = -2.0 / 3.0; + robin_f_fab(i, j, k) = 0.0; + } + }); + } + + void evaluateRad(amrex::MultiFab& rad_src) + { + // std::cout << "evaluateRad() is called" << std::endl; + amrex::Array lobc{AMREX_D_DECL( + amrex::LinOpBCType::Robin, amrex::LinOpBCType::Periodic, + amrex::LinOpBCType::Periodic)}; + amrex::Array hibc{AMREX_D_DECL( + amrex::LinOpBCType::Neumann, amrex::LinOpBCType::Periodic, + amrex::LinOpBCType::Periodic)}; + + amrex::ParmParse pp("pelerad"); + MLMGParam mlmgpp(pp); + + bcoef_.FillBoundary(); + + // std::cout << "before the rte constructor" << std::endl; + POneSingle rte( + mlmgpp, geom_, grids_, dmap_, solution_, rhs_, acoef_, bcoef_, lobc, hibc, + robin_a_, robin_b_, robin_f_); + + rte.solve(); + + rte.calcRadSource(rad_src); + // write(); + // amrex::Abort(); + } + + RadComps const readRadIndices() const { return rc_; } + + void write() + { + amrex::MultiFab plotmf(grids_, dmap_, 4, 0); + amrex::MultiFab::Copy(plotmf, solution_, 0, 0, 1, 0); + amrex::MultiFab::Copy(plotmf, rhs_, 0, 1, 1, 0); + amrex::MultiFab::Copy(plotmf, acoef_, 0, 2, 1, 0); + amrex::MultiFab::Copy(plotmf, bcoef_, 0, 3, 1, 0); + + amrex::WriteSingleLevelPlotfile( + "plotrad", plotmf, {"solution", "rhs", "acoef", "bcoef"}, geom_, 0.0, 0); + } +}; + +} // namespace PeleRad + +#endif diff --git a/Source/Radiation/PeleLMRad.H b/Source/Radiation/PeleLMRad.H new file mode 100644 index 000000000..2d13107da --- /dev/null +++ b/Source/Radiation/PeleLMRad.H @@ -0,0 +1,381 @@ +#ifndef RADIATION_H +#define RADIATION_H + +#include +#include +#include + +#ifdef AMREX_USE_EB +#include +#else +#include +#endif + +#include +#include + +namespace PeleRad { + +class Radiation +{ +private: + PlanckMean radprop; + + amrex::Vector& geom_; + amrex::Vector& grids_; + amrex::Vector& dmap_; + + RadComps rc_; + + MLMGParam mlmgpp_; + + amrex::Vector solution_; + amrex::Vector rhs_; + amrex::Vector acoef_; + amrex::Vector bcoef_; + amrex::Vector robin_a_; + amrex::Vector robin_b_; + amrex::Vector robin_f_; + + amrex::Vector absc_; + + // bool composite_solve_; + +#ifdef AMREX_USE_EB + std::unique_ptr rte_; +#else + std::unique_ptr rte_; +#endif + + // std::unique_ptr rtelevbylev_; + +#ifdef AMREX_USE_EB + amrex::Vector ebfactVec_; +#endif + +public: + AMREX_GPU_HOST + Radiation( + amrex::Vector& geom, + amrex::Vector& grids, + amrex::Vector& dmap, + RadComps rc, + amrex::ParmParse const& mlmgpp +#ifdef AMREX_USE_EB + , + amrex::Vector>> const& + factory +#endif + ) + : geom_(geom), grids_(grids), dmap_(dmap), rc_(rc), mlmgpp_(mlmgpp) + { + if (amrex::ParallelDescriptor::IOProcessor()) + rc_.checkIndices(); + + auto const nlevels = geom_.size(); + + solution_.resize(nlevels); + rhs_.resize(nlevels); + acoef_.resize(nlevels); + bcoef_.resize(nlevels); + + robin_a_.resize(nlevels); + robin_b_.resize(nlevels); + robin_f_.resize(nlevels); + absc_.resize(nlevels); + +#ifdef AMREX_USE_EB + for (int ilev = 0; ilev < nlevels; ++ilev) { + ebfactVec_.push_back( + &(static_cast(*factory[ilev]))); + } + + initVars(grids, dmap, ebfactVec_); +#else + initVars(grids, dmap); +#endif + + loadSpecModel(); + + // composite_solve_ = mlmgpp_.composite_solve_; + + // if (composite_solve_) + // { + +#ifdef AMREX_USE_EB + rte_ = std::make_unique( + mlmgpp_, geom_, grids_, dmap_, ebfactVec_, solution_, rhs_, acoef_, + bcoef_, robin_a_, robin_b_, robin_f_); +#else + rte_ = std::make_unique( + mlmgpp_, geom_, grids_, dmap_, solution_, rhs_, acoef_, bcoef_, robin_a_, + robin_b_, robin_f_); +#endif + // } + // else + // { + // level by level option is not ready + // rtelevbylev_ = + // std::make_unique(mlmgpp_, 2, + // geom_, grids_, dmap_, solution_, rhs_, acoef_, bcoef_, + // robin_a_, robin_b_, robin_f_); + // } + } + + AMREX_GPU_HOST + void loadSpecModel() + { + // path to spectral database + auto data_path = mlmgpp_.kppath_; + + radprop.load(data_path); + + amrex::Print() << "The radiative property database is loaded." + << "\n" + << "kp path: " << data_path << "\n"; + } + + void updateSpecProp( + amrex::MFIter const& mfi, + amrex::Array4 const& Yco2, + amrex::Array4 const& Yh2o, + amrex::Array4 const& Yco, + amrex::Array4 const& T, + amrex::Array4 const& P +#ifdef PELELM_USE_SOOT + , + amrex::Array4 const& fv +#endif + , + int ilev) + { + // amrex::Print() << "update radiative properties \n"; + + auto const& kpco2 = radprop.kpco2(); + auto const& kph2o = radprop.kph2o(); + auto const& kpco = radprop.kpco(); + + amrex::Box const& bx = mfi.validbox(); + amrex::Box const& gbx = amrex::grow(bx, 1); + + auto const dlo = amrex::lbound(geom_[ilev].Domain()); + auto const dhi = amrex::ubound(geom_[ilev].Domain()); + + auto const& kappa = absc_[ilev].array(mfi); + auto const& rhsfab = rhs_[ilev].array(mfi); + auto const& alphafab = acoef_[ilev].array(mfi); + auto const& betafab = bcoef_[ilev].array(mfi); + auto const& robin_a_fab = robin_a_[ilev].array(mfi); + auto const& robin_b_fab = robin_b_[ilev].array(mfi); + auto const& robin_f_fab = robin_f_[ilev].array(mfi); + + amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { + RadProp::getRadPropGas( + i, j, k, Yco2, Yh2o, Yco, T, P, kappa, kpco2, kph2o, kpco); + }); + +#ifdef PELELM_USE_SOOT + auto const& kpsoot = radprop.kpsoot(); + amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { + RadProp::getRadPropSoot(i, j, k, fv, T, kappa, kpsoot); + }); +#endif + + amrex::ParallelFor(gbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { + betafab(i, j, k) = 1.0; + + if (bx.contains(i, j, k)) { + double ka = std::max(0.01, kappa(i, j, k) * 100); + betafab(i, j, k) = 1.0 / ka; + + rhsfab(i, j, k) = 4.0 * ka * 5.67e-8 * std::pow(T(i, j, k), + 4.0); // SI + + /*rhsfab(i, j, k) = 4.0 * ka * 5.67e-5 + * std::pow(T(i, j, k), + 4.0);*/ // cgs + alphafab(i, j, k) = ka; + } + + // Robin BC + bool robin_cell = false; + if (j >= dlo.y && j <= dhi.y && k >= dlo.z && k <= dhi.z) { + int jj, kk; + + if (j < bx.loVect3d()[1]) + jj = bx.loVect3d()[1]; + else if (j > bx.hiVect3d()[1]) + jj = bx.hiVect3d()[1]; + else + jj = j; + if (k < bx.loVect3d()[2]) + kk = bx.loVect3d()[2]; + else if (k > bx.hiVect3d()[2]) + kk = bx.hiVect3d()[2]; + else + kk = k; + + if (i > dhi.x) { + robin_cell = true; + betafab(i, j, k) = 1.0 / std::max(1.0, kappa(dhi.x, jj, kk) * 100); + } + if (i < dlo.x) { + robin_cell = true; + betafab(i, j, k) = 1.0 / std::max(1.0, kappa(dlo.x, jj, kk) * 100); + } + } else if (i >= dlo.x && i <= dhi.x && k >= dlo.z && k <= dhi.z) { + int ii, kk; + + if (i < bx.loVect3d()[0]) + ii = bx.loVect3d()[0]; + else if (i > bx.hiVect3d()[0]) + ii = bx.hiVect3d()[0]; + else + ii = i; + if (k < bx.loVect3d()[2]) + kk = bx.loVect3d()[2]; + else if (k > bx.hiVect3d()[2]) + kk = bx.hiVect3d()[2]; + else + kk = k; + + if (j > dhi.y) { + robin_cell = true; + betafab(i, j, k) = 1.0 / std::max(1.0, kappa(ii, dhi.y, kk) * 100); + } + if (j < dlo.y) { + robin_cell = true; + betafab(i, j, k) = 1.0 / std::max(1.0, kappa(ii, dlo.y, kk) * 100); + } + } else if (i >= dlo.x && i <= dhi.x && j >= dlo.y && j <= dhi.y) { + int ii, jj; + + if (i < bx.loVect3d()[0]) + ii = bx.loVect3d()[0]; + else if (i > bx.hiVect3d()[0]) + ii = bx.hiVect3d()[0]; + else + ii = i; + if (j < bx.loVect3d()[1]) + jj = bx.loVect3d()[1]; + else if (j > bx.hiVect3d()[1]) + jj = bx.hiVect3d()[1]; + else + jj = j; + + if (k > dhi.z) { + robin_cell = true; + betafab(i, j, k) = 1.0 / std::max(1.0, kappa(ii, jj, dhi.z) * 100); + } + if (k < dlo.z) { + robin_cell = true; + betafab(i, j, k) = 1.0 / std::max(1.0, kappa(ii, jj, dlo.z) * 100); + } + } + + if (robin_cell) { + robin_a_fab(i, j, k) = -1.0 / betafab(i, j, k); + robin_b_fab(i, j, k) = -2.0 / 3.0; + robin_f_fab(i, j, k) = 0.0; + } + }); + } + + void initVars( + amrex::Vector const& grids, + amrex::Vector const& dmap) + { + amrex::IntVect ng = amrex::IntVect{1}; + grids_ = grids; + dmap_ = dmap; + + for (int ilev = 0; ilev < grids.size(); ++ilev) { + solution_[ilev].define(grids[ilev], dmap[ilev], 1, ng); + rhs_[ilev].define(grids[ilev], dmap[ilev], 1, 0); + acoef_[ilev].define(grids[ilev], dmap[ilev], 1, 0); + bcoef_[ilev].define(grids[ilev], dmap[ilev], 1, ng); + robin_a_[ilev].define(grids[ilev], dmap[ilev], 1, ng); + robin_b_[ilev].define(grids[ilev], dmap[ilev], 1, ng); + robin_f_[ilev].define(grids[ilev], dmap[ilev], 1, ng); + absc_[ilev].define(grids[ilev], dmap[ilev], 1, 0); + + solution_[ilev].setVal(0.0, 0, 1, ng); + bcoef_[ilev].setVal(1.0, 0, 1, ng); + } + } + +#ifdef AMREX_USE_EB + void initVars( + amrex::Vector const& grids, + amrex::Vector const& dmap, + amrex::Vector& factory) + { + amrex::IntVect ng = amrex::IntVect{1}; + grids_ = grids; + dmap_ = dmap; + + for (int ilev = 0; ilev < grids.size(); ++ilev) { + solution_[ilev].define( + grids[ilev], dmap[ilev], 1, ng, amrex::MFInfo(), *factory[ilev]); + rhs_[ilev].define( + grids[ilev], dmap[ilev], 1, 0, amrex::MFInfo(), *factory[ilev]); + acoef_[ilev].define( + grids[ilev], dmap[ilev], 1, 0, amrex::MFInfo(), *factory[ilev]); + bcoef_[ilev].define( + grids[ilev], dmap[ilev], 1, ng, amrex::MFInfo(), *factory[ilev]); + robin_a_[ilev].define( + grids[ilev], dmap[ilev], 1, ng, amrex::MFInfo(), *factory[ilev]); + robin_b_[ilev].define( + grids[ilev], dmap[ilev], 1, ng, amrex::MFInfo(), *factory[ilev]); + robin_f_[ilev].define( + grids[ilev], dmap[ilev], 1, ng, amrex::MFInfo(), *factory[ilev]); + absc_[ilev].define( + grids[ilev], dmap[ilev], 1, 0, amrex::MFInfo(), *factory[ilev]); + + solution_[ilev].setVal(0.0, 0, 1, ng); + bcoef_[ilev].setVal(1.0, 0, 1, ng); + } + } +#endif + + void evaluateRad() + { + // std::cout << "begin of evaluateRad() \n"; + + for (int ilev = 0; ilev < grids_.size(); ++ilev) + bcoef_[ilev].FillBoundary(); + // if (composite_solve_) + rte_->solve(); + // else + // rtelevbylev_->solve(); + } + + void calcRadSource( + amrex::MFIter const& mfi, + amrex::Array4 const& radfab, + int ilev) + { + amrex::Box const& bx = mfi.validbox(); + auto const& rhsfab = rhs_[ilev].array(mfi); + auto const& solfab = solution_[ilev].array(mfi); + auto const& acfab = acoef_[ilev].array(mfi); + + amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { + radfab(i, j, k) += acfab(i, j, k) * solfab(i, j, k) - rhsfab(i, j, k); + }); + } + + RadComps const readRadIndices() const { return rc_; } + + amrex::Vector const& G() { return solution_; } + + amrex::Vector const& kappa() { return acoef_; } + + amrex::Vector const& emis() { return rhs_; } + + amrex::Vector const& grids() { return grids_; } +}; + +} // namespace PeleRad + +#endif diff --git a/Source/Radiation/PlanckMean.H b/Source/Radiation/PlanckMean.H new file mode 100644 index 000000000..049537a64 --- /dev/null +++ b/Source/Radiation/PlanckMean.H @@ -0,0 +1,96 @@ +#ifndef PLANCK_MEAN_H +#define PLANCK_MEAN_H + +#include + +#include +#include + +namespace PeleRad { + +class PlanckMean +{ +private: + amrex::GpuArray kpco2_; + amrex::GpuArray kph2o_; + amrex::GpuArray kpco_; + amrex::GpuArray kpch4_; + amrex::GpuArray kpc2h4_; + amrex::GpuArray kpsoot_; + + AMREX_GPU_HOST + void read_kp_file( + std::string file_path, amrex::GpuArray& file_data) + { + std::ifstream datastream(file_path); + if (!datastream.is_open()) { + amrex::Abort("PeleRad: Failed to open data file: " + file_path); + } + size_t idx = 0; + for (amrex::Real T_temp, kp_temp; datastream >> T_temp >> kp_temp;) { + if (!amrex::almostEqual(T_temp, 300.0 + 20.0 * idx)) { + amrex::Abort( + "PeleRad: Invalid temperatures in data file: " + file_path); + } + if (idx > 125) { + amrex::Abort( + "PeleRad: Invalid number of entries in data file, must have 126: " + + file_path); + } + file_data[idx] = kp_temp; + idx++; + } + if (idx < 126) { + amrex::Abort( + "PeleRad: Invalid number of entries in data file, must have 126: " + + file_path); + } + datastream.close(); + } + +public: + AMREX_GPU_HOST + PlanckMean() = default; + + AMREX_GPU_HOST + PlanckMean(std::string data_path) { load(data_path); } + + AMREX_GPU_HOST + PlanckMean(PlanckMean const&) = delete; + + AMREX_GPU_HOST + PlanckMean& operator=(PlanckMean const&) = delete; + + AMREX_GPU_HOST + void load(std::string data_path) + { + read_kp_file(data_path + "/kpl_co2.dat", kpco2_); + read_kp_file(data_path + "/kpl_h2o.dat", kph2o_); + read_kp_file(data_path + "/kpl_co.dat", kpco_); + read_kp_file(data_path + "/kpl_ch4.dat", kpch4_); + read_kp_file(data_path + "/kpl_c2h4.dat", kpc2h4_); + read_kp_file(data_path + "/kpl_soot.dat", kpsoot_); + } + + AMREX_GPU_HOST_DEVICE + const amrex::GpuArray& kpco2() const { return kpco2_; } + + AMREX_GPU_HOST_DEVICE + const amrex::GpuArray& kph2o() const { return kph2o_; } + + AMREX_GPU_HOST_DEVICE + const amrex::GpuArray& kpco() const { return kpco_; } + + AMREX_GPU_HOST_DEVICE + const amrex::GpuArray& kpch4() const { return kpch4_; } + + AMREX_GPU_HOST_DEVICE + const amrex::GpuArray& kpc2h4() const { return kpc2h4_; } + + AMREX_GPU_HOST_DEVICE + const amrex::GpuArray& kpsoot() const { return kpsoot_; } +}; + +} // namespace PeleRad + +#endif diff --git a/Source/Radiation/SpectralModels.H b/Source/Radiation/SpectralModels.H new file mode 100644 index 000000000..26646edb6 --- /dev/null +++ b/Source/Radiation/SpectralModels.H @@ -0,0 +1,107 @@ +#ifndef SPECTRAL_MODELS_H +#define SPECTRAL_MODELS_H + +#include +#include +#include + +namespace PeleRad { + +namespace RadProp { +AMREX_GPU_HOST_DEVICE +AMREX_FORCE_INLINE +void +interpT(amrex::Real const& T, int& TindexL, amrex::Real& weight) +{ + if (T < 300) { + TindexL = 0; + weight = 0.0; + return; + } + if (T > 2800) { + TindexL = 125; + weight = 1.0; + return; + } + + amrex::Real TindexReal = (T - 300.0) / 20.0; + amrex::Real TindexInte = floor((T - 300.0) / 20.0); + TindexL = static_cast(TindexInte); + weight = TindexReal - TindexInte; + + assert(weight <= 1 && weight >= 0); +} + +AMREX_GPU_HOST_DEVICE +AMREX_FORCE_INLINE +amrex::Real +interpk( + int const& TindexL, + amrex::Real const& weight, + amrex::GpuArray const& k) +{ + return (1.0 - weight) * k[TindexL] + weight * k[TindexL + 1]; +} + +AMREX_GPU_HOST_DEVICE +AMREX_FORCE_INLINE +void +getRadPropGas( + int i, + int j, + int k, + amrex::Array4 const& yco2, + amrex::Array4 const& yh2o, + amrex::Array4 const& yco, + amrex::Array4 const& temp, + amrex::Array4 const& pressure, + amrex::Array4 const& absc, + amrex::GpuArray const& kdataco2, + amrex::GpuArray const& kdatah2o, + amrex::GpuArray const& kdataco) +{ + int TindexL = 0; + amrex::Real weight = 1.0; + interpT(temp(i, j, k), TindexL, weight); + + amrex::Real kp_co2 = interpk(TindexL, weight, kdataco2); + amrex::Real kp_h2o = interpk(TindexL, weight, kdatah2o); + amrex::Real kp_co = interpk(TindexL, weight, kdataco); + + absc(i, j, k) = + yco2(i, j, k) * kp_co2 + yh2o(i, j, k) * kp_h2o + yco(i, j, k) * kp_co; + + // absc(i, j, k) *= pressure(i, j, k) / 1.0e5 * 100.0; //si, in m-1 + absc(i, j, k) *= + pressure(i, j, k) / + 1.0e5; // cgs, in cm-1; if P is in bar(PeleLM), no correction + // is needed, otherwise the absc needs to be corrected. +} + +AMREX_GPU_HOST_DEVICE +AMREX_FORCE_INLINE +void +getRadPropSoot( + int i, + int j, + int k, + amrex::Array4 const& fv, + amrex::Array4 const& temp, + amrex::Array4 const& absc, + amrex::GpuArray const& kdatasoot) +{ + int TindexL = 0; + amrex::Real weight = 1.0; + + interpT(temp(i, j, k), TindexL, weight); + + amrex::Real kp_soot = interpk(TindexL, weight, kdatasoot); + + // absc(i, j, k) += fv(i, j, k) * kp_soot * 100.0; //si, in m-1 + absc(i, j, k) += fv(i, j, k) * kp_soot; // cgs, in cm-1 +} + +} // namespace RadProp + +} // namespace PeleRad +#endif diff --git a/Testing/Exec/Radiation/CMakeLists.txt b/Testing/Exec/Radiation/CMakeLists.txt new file mode 100644 index 000000000..e9aad6c36 --- /dev/null +++ b/Testing/Exec/Radiation/CMakeLists.txt @@ -0,0 +1,62 @@ +add_library(PeleRad OBJECT) + +set(RADIATION_INCLUDE ${CMAKE_SOURCE_DIR}/Submodules/PelePhysics/Source/Radiation) +target_sources(PeleRad + PRIVATE + ${RADIATION_INCLUDE}/AMRParam.H + ${RADIATION_INCLUDE}/Constants.H + ${RADIATION_INCLUDE}/MLMGParam.H + ${RADIATION_INCLUDE}/POneMulti.H + ${RADIATION_INCLUDE}/POneMultiEB.H + ${RADIATION_INCLUDE}/POneMultiLevbyLev.H + ${RADIATION_INCLUDE}/POneSingle.H + ${RADIATION_INCLUDE}/POneSingleEB.H + ${RADIATION_INCLUDE}/PeleCRad.H + ${RADIATION_INCLUDE}/PeleLMRad.H + ${RADIATION_INCLUDE}/PlanckMean.H + ${RADIATION_INCLUDE}/SpectralModels.H +) +target_include_directories(PeleRad PUBLIC ${RADIATION_INCLUDE}) +target_link_libraries(PeleRad PUBLIC AMReX::amrex) + +file(COPY kpDB DESTINATION "${CMAKE_CURRENT_BINARY_DIR}/") + +if(PELELMEX_ENABLE_EB) + add_executable(PeleRad_POneSingleEB.exe tstPOneSingleEB.cpp) + target_link_libraries(PeleRad_POneSingleEB.exe PRIVATE PeleRad) + add_test(NAME PeleRad_POneSingleEB_Test COMMAND ${MPIEXEC_EXECUTABLE} ${MPIEXEC_NUMPROC_FLAG} ${MPIEXEC_MAX_NUMPROCS} ${MPIEXEC_PREFLAGS} ${CMAKE_CURRENT_BINARY_DIR}/PeleRad_POneSingleEB.exe ${CMAKE_CURRENT_SOURCE_DIR}/inputs/inputs.tstPOneSingleEB) + + add_executable(PeleRad_POneMultiEB.exe tstPOneMultiEB.cpp) + target_link_libraries(PeleRad_POneMultiEB.exe PRIVATE PeleRad) + add_test(NAME PeleRad_POneMultiEB_Test COMMAND ${MPIEXEC_EXECUTABLE} ${MPIEXEC_NUMPROC_FLAG} ${MPIEXEC_MAX_NUMPROCS} ${MPIEXEC_PREFLAGS} ${CMAKE_CURRENT_BINARY_DIR}/PeleRad_POneMultiEB.exe ${CMAKE_CURRENT_SOURCE_DIR}/inputs/inputs.tstPOneMultiEB) +else() + add_executable(PeleRad_POneSingle.exe tstPOneSingle.cpp) + target_link_libraries(PeleRad_POneSingle.exe PRIVATE PeleRad AMReX::amrex) + add_test(NAME PeleRad_POneSingle_Test COMMAND ${MPIEXEC_EXECUTABLE} ${MPIEXEC_NUMPROC_FLAG} ${MPIEXEC_MAX_NUMPROCS} ${MPIEXEC_PREFLAGS} ${CMAKE_CURRENT_BINARY_DIR}/PeleRad_POneSingle.exe ${CMAKE_CURRENT_SOURCE_DIR}/inputs/inputs.tstPOneSingle) + + add_executable(PeleRad_POneMulti.exe tstPOneMulti.cpp) + target_link_libraries(PeleRad_POneMulti.exe PRIVATE PeleRad) + add_test(NAME PeleRad_POneMulti_Test_Composite COMMAND ${MPIEXEC_EXECUTABLE} ${MPIEXEC_NUMPROC_FLAG} ${MPIEXEC_MAX_NUMPROCS} ${MPIEXEC_PREFLAGS} ${CMAKE_CURRENT_BINARY_DIR}/PeleRad_POneMulti.exe ${CMAKE_CURRENT_SOURCE_DIR}/inputs/inputs.tstPOneMulti) + add_test(NAME PeleRad_POneMulti_Test_LevbyLev COMMAND ${MPIEXEC_EXECUTABLE} ${MPIEXEC_NUMPROC_FLAG} ${MPIEXEC_MAX_NUMPROCS} ${MPIEXEC_PREFLAGS} ${CMAKE_CURRENT_BINARY_DIR}/PeleRad_POneMulti.exe ${CMAKE_CURRENT_SOURCE_DIR}/inputs/inputs.tstPOneMultiLevbyLev) + + add_executable(PeleRad_POneSingleAF.exe tstPOneSingleAF.cpp) + target_link_libraries(PeleRad_POneSingleAF.exe PRIVATE PeleRad) + add_test(NAME PeleRad_POneSingleAF_Test COMMAND ${MPIEXEC_EXECUTABLE} ${MPIEXEC_NUMPROC_FLAG} ${MPIEXEC_MAX_NUMPROCS} ${MPIEXEC_PREFLAGS} ${CMAKE_CURRENT_BINARY_DIR}/PeleRad_POneSingleAF.exe ${CMAKE_CURRENT_SOURCE_DIR}/inputs/inputs.tstPOneSingleAF) + + add_executable(PeleRad_POneMultiAF.exe tstPOneMultiAF.cpp) + target_link_libraries(PeleRad_POneMultiAF.exe PRIVATE PeleRad) + add_test(NAME PeleRad_POneMultiAF_Test_Composite COMMAND ${MPIEXEC_EXECUTABLE} ${MPIEXEC_NUMPROC_FLAG} ${MPIEXEC_MAX_NUMPROCS} ${MPIEXEC_PREFLAGS} ${CMAKE_CURRENT_BINARY_DIR}/PeleRad_POneMultiAF.exe ${CMAKE_CURRENT_SOURCE_DIR}/inputs/inputs.tstPOneMultiAF) + add_test(NAME PeleRad_POneMultiAF_Test_LevbyLev COMMAND ${MPIEXEC_EXECUTABLE} ${MPIEXEC_NUMPROC_FLAG} ${MPIEXEC_MAX_NUMPROCS} ${MPIEXEC_PREFLAGS} ${CMAKE_CURRENT_BINARY_DIR}/PeleRad_POneMultiAF.exe ${CMAKE_CURRENT_SOURCE_DIR}/inputs/inputs.tstPOneMultiAFLevbyLev) +endif() + +if(PELELMEX_ENABLE_CUDA) + if(PELELMEX_ENABLE_EB) + setup_target_for_cuda_compilation(PeleRad_POneSingleEB.exe) + setup_target_for_cuda_compilation(PeleRad_POneMultiEB.exe) + else() + setup_target_for_cuda_compilation(PeleRad_POneSingle.exe) + setup_target_for_cuda_compilation(PeleRad_POneMulti.exe) + setup_target_for_cuda_compilation(PeleRad_POneSingleAF.exe) + setup_target_for_cuda_compilation(PeleRad_POneMultiAF.exe) + endif() +endif() diff --git a/Testing/Exec/Radiation/inputs/inputs.tstPOneMulti b/Testing/Exec/Radiation/inputs/inputs.tstPOneMulti new file mode 100644 index 000000000..5415e809c --- /dev/null +++ b/Testing/Exec/Radiation/inputs/inputs.tstPOneMulti @@ -0,0 +1,26 @@ +#amr +max_level = 2 +ref_ratio = 2 +n_cell = 64 +max_grid_size = 32 +plot_file_name = "pltRobinMultiComp" + +#mlmg +composite_solve = 1 +use_hypre = 0 +verbose = 0 +bottom_verbose = 0 +bottom_reltol = 1.e-6 +bottom_abstol = 1.e-6 +max_iter = 20 +max_fmg_iter = 0 +max_coarsening_level = 10 +reltol = 1.e-6 +abstol = 1.e-6 +linop_maxorder = 2 +agglomeration = 1 +consolidation = 0 + +#bc +lo_bc = Robin Dirichlet Neumann +hi_bc = Robin Dirichlet Neumann diff --git a/Testing/Exec/Radiation/inputs/inputs.tstPOneMultiAF b/Testing/Exec/Radiation/inputs/inputs.tstPOneMultiAF new file mode 100644 index 000000000..fe6f696c1 --- /dev/null +++ b/Testing/Exec/Radiation/inputs/inputs.tstPOneMultiAF @@ -0,0 +1,25 @@ +#amr +max_level = 2 +ref_ratio = 2 +n_cell = 32 +max_grid_size = 32 +plot_file_name = "pltRobinAFMultiComp" + +#mlmg +composite_solve = 1 +use_hypre = 0 +verbose = 0 +bottom_verbose = 0 +bottom_reltol = 1.0e-5 +bottom_abstol = 1.0e-5 +max_iter = 200 +max_fmg_iter = 0 +max_coarsening_level = 10 +reltol = 1.0e-5 +abstol = 1.0e-5 +linop_maxorder = 2 +agglomeration = 1 +consolidation = 0 + +lo_bc = Robin Robin Neumann +hi_bc = Robin Robin Neumann diff --git a/Testing/Exec/Radiation/inputs/inputs.tstPOneMultiAFLevbyLev b/Testing/Exec/Radiation/inputs/inputs.tstPOneMultiAFLevbyLev new file mode 100644 index 000000000..94f60e9b8 --- /dev/null +++ b/Testing/Exec/Radiation/inputs/inputs.tstPOneMultiAFLevbyLev @@ -0,0 +1,23 @@ +#amr +max_level = 2 +ref_ratio = 2 +n_cell = 32 +max_grid_size = 32 +plot_file_name = "pltRobinAFMultiLevbyLev" + +#mlmg +composite_solve = 0 +use_hypre = 0 +verbose = 0 +bottom_verbose = 0 +max_iter = 200 +max_fmg_iter = 0 +max_coarsening_level = 10 +reltol = 1.0e-5 +abstol = 1.0e-5 +linop_maxorder = 2 +agglomeration = 1 +consolidation = 1 + +lo_bc = Robin Robin Neumann +hi_bc = Robin Robin Neumann diff --git a/Testing/Exec/Radiation/inputs/inputs.tstPOneMultiEB b/Testing/Exec/Radiation/inputs/inputs.tstPOneMultiEB new file mode 100644 index 000000000..b2966a19a --- /dev/null +++ b/Testing/Exec/Radiation/inputs/inputs.tstPOneMultiEB @@ -0,0 +1,28 @@ +#amr +max_level = 1 +ref_ratio = 2 +n_cell = 64 +max_grid_size = 32 + +plot_file_name = "pltMultiEB" + +#mlmg +composite_solve = 1 +use_hypre = 0 +verbose = 0 +bottom_verbose = 0 +max_iter = 20 +max_fmg_iter = 0 +max_coarsening_level = 30 +reltol = 1.e-6 +abstol = 1.e-6 +bottom_reltol = 1.e-6 +bottom_abstol = 1.e-6 +linop_maxorder = 2 +agglomeration = 1 +consolidation = 0 + +lo_bc = Periodic Periodic Periodic +hi_bc = Periodic Periodic Periodic + +ebbc_type = 1 diff --git a/Testing/Exec/Radiation/inputs/inputs.tstPOneMultiLevbyLev b/Testing/Exec/Radiation/inputs/inputs.tstPOneMultiLevbyLev new file mode 100644 index 000000000..d9569eb8b --- /dev/null +++ b/Testing/Exec/Radiation/inputs/inputs.tstPOneMultiLevbyLev @@ -0,0 +1,23 @@ +#amr +max_level = 2 +ref_ratio = 2 +n_cell = 64 +max_grid_size = 32 +plot_file_name = "pltRobinMultiLevbyLev" + +#mlmg +composite_solve = 0 +use_hypre = 0 +verbose = 0 +bottom_verbose = 0 +max_iter = 20 +max_fmg_iter = 0 +max_coarsening_level = 10 +reltol = 1.e-6 +abstol = 1.e-6 +linop_maxorder = 2 +agglomeration = 1 +consolidation = 1 + +lo_bc = Robin Dirichlet Neumann +hi_bc = Robin Dirichlet Neumann diff --git a/Testing/Exec/Radiation/inputs/inputs.tstPOneSingle b/Testing/Exec/Radiation/inputs/inputs.tstPOneSingle new file mode 100644 index 000000000..5b8e177e5 --- /dev/null +++ b/Testing/Exec/Radiation/inputs/inputs.tstPOneSingle @@ -0,0 +1,23 @@ +#amr +max_level = 0 +n_cell = 64 +max_grid_size = 32 +plot_file_name = "pltRobinSingle" + +#mlmg +use_hypre = 0 +verbose = 0 +bottom_verbose = 0 +max_iter = 10 +max_fmg_iter = 0 +max_coarsening_level = 10 +reltol = 1.e-5 +abstol = 1.e-5 +bottom_reltol = 1.e-5 +bottom_abstol = 1.e-5 +linop_maxorder = 2 +agglomeration = 1 +consolidation = 0 + +lo_bc = Robin Dirichlet Neumann +hi_bc = Robin Dirichlet Neumann diff --git a/Testing/Exec/Radiation/inputs/inputs.tstPOneSingleAF b/Testing/Exec/Radiation/inputs/inputs.tstPOneSingleAF new file mode 100644 index 000000000..a935f3abb --- /dev/null +++ b/Testing/Exec/Radiation/inputs/inputs.tstPOneSingleAF @@ -0,0 +1,24 @@ +#amr +max_level = 0 +n_cell = 32 +max_grid_size = 32 +plot_file_name = "pltRobinAF" + +#mlmg +use_hypre = 0 +verbose = 0 +bottom_verbose = 0 +max_iter = 200 +max_fmg_iter = 0 +max_coarsening_level = 10 +reltol = 1.e-6 +abstol = 1.e-6 +bottom_reltol = 1.e-6 +bottom_abstol = 1.e-6 +linop_maxorder = 2 +agglomeration = 1 +consolidation = 0 + +#bc +lo_bc = Robin Robin Neumann +hi_bc = Robin Robin Neumann diff --git a/Testing/Exec/Radiation/inputs/inputs.tstPOneSingleEB b/Testing/Exec/Radiation/inputs/inputs.tstPOneSingleEB new file mode 100644 index 000000000..ffa6bdba7 --- /dev/null +++ b/Testing/Exec/Radiation/inputs/inputs.tstPOneSingleEB @@ -0,0 +1,25 @@ +#amr +max_level = 0 +n_cell = 64 +max_grid_size = 32 +plot_file_name = "pltSingleEB" + +#mlmg +use_hypre = 0 +verbose = 0 +bottom_verbose = 0 +max_iter = 10 +max_fmg_iter = 0 +max_coarsening_level = 10 +reltol = 1.e-5 +abstol = 1.e-5 +bottom_reltol = 1.e-5 +bottom_abstol = 1.e-5 +linop_maxorder = 2 +agglomeration = 1 +consolidation = 0 + +lo_bc = Periodic Periodic Periodic +hi_bc = Periodic Periodic Periodic + +ebbc_type = 1 diff --git a/Testing/Exec/Radiation/kpDB/kpl_c2h4.dat b/Testing/Exec/Radiation/kpDB/kpl_c2h4.dat new file mode 100755 index 000000000..e82ea7c40 --- /dev/null +++ b/Testing/Exec/Radiation/kpDB/kpl_c2h4.dat @@ -0,0 +1,126 @@ + 300.0 0.22175E+00 + 320.0 0.20916E+00 + 340.0 0.19394E+00 + 360.0 0.17738E+00 + 380.0 0.16047E+00 + 400.0 0.14388E+00 + 420.0 0.12809E+00 + 440.0 0.11338E+00 + 460.0 0.99894E-01 + 480.0 0.87694E-01 + 500.0 0.76764E-01 + 520.0 0.67049E-01 + 540.0 0.58469E-01 + 560.0 0.50927E-01 + 580.0 0.44322E-01 + 600.0 0.38556E-01 + 620.0 0.33533E-01 + 640.0 0.29163E-01 + 660.0 0.25367E-01 + 680.0 0.22072E-01 + 700.0 0.19212E-01 + 720.0 0.16731E-01 + 740.0 0.14577E-01 + 760.0 0.12709E-01 + 780.0 0.11086E-01 + 800.0 0.96764E-02 + 820.0 0.84515E-02 + 840.0 0.73863E-02 + 860.0 0.64594E-02 + 880.0 0.56524E-02 + 900.0 0.49493E-02 + 920.0 0.43363E-02 + 940.0 0.38015E-02 + 960.0 0.33346E-02 + 980.0 0.29268E-02 + 1000.0 0.25704E-02 + 1020.0 0.22586E-02 + 1040.0 0.19857E-02 + 1060.0 0.17468E-02 + 1080.0 0.15374E-02 + 1100.0 0.13539E-02 + 1120.0 0.11929E-02 + 1140.0 0.10516E-02 + 1160.0 0.92752E-03 + 1180.0 0.81850E-03 + 1200.0 0.72266E-03 + 1220.0 0.63837E-03 + 1240.0 0.56420E-03 + 1260.0 0.49889E-03 + 1280.0 0.44137E-03 + 1300.0 0.39068E-03 + 1320.0 0.34598E-03 + 1340.0 0.30654E-03 + 1360.0 0.27174E-03 + 1380.0 0.24101E-03 + 1400.0 0.21386E-03 + 1420.0 0.18987E-03 + 1440.0 0.16864E-03 + 1460.0 0.14987E-03 + 1480.0 0.13325E-03 + 1500.0 0.11854E-03 + 1520.0 0.10550E-03 + 1540.0 0.93939E-04 + 1560.0 0.83688E-04 + 1580.0 0.74593E-04 + 1600.0 0.66519E-04 + 1620.0 0.59349E-04 + 1640.0 0.52977E-04 + 1660.0 0.47313E-04 + 1680.0 0.42275E-04 + 1700.0 0.37793E-04 + 1720.0 0.33802E-04 + 1740.0 0.30247E-04 + 1760.0 0.27079E-04 + 1780.0 0.24255E-04 + 1800.0 0.21736E-04 + 1820.0 0.19488E-04 + 1840.0 0.17481E-04 + 1860.0 0.15688E-04 + 1880.0 0.14086E-04 + 1900.0 0.12654E-04 + 1920.0 0.11372E-04 + 1940.0 0.10225E-04 + 1960.0 0.91986E-05 + 1980.0 0.82789E-05 + 2000.0 0.74546E-05 + 2020.0 0.67155E-05 + 2040.0 0.60525E-05 + 2060.0 0.54575E-05 + 2080.0 0.49232E-05 + 2100.0 0.44433E-05 + 2120.0 0.40120E-05 + 2140.0 0.36242E-05 + 2160.0 0.32754E-05 + 2180.0 0.29615E-05 + 2200.0 0.26789E-05 + 2220.0 0.24243E-05 + 2240.0 0.21949E-05 + 2260.0 0.19881E-05 + 2280.0 0.18016E-05 + 2300.0 0.16332E-05 + 2320.0 0.14813E-05 + 2340.0 0.13440E-05 + 2360.0 0.12200E-05 + 2380.0 0.11079E-05 + 2400.0 0.10066E-05 + 2420.0 0.91488E-06 + 2440.0 0.83187E-06 + 2460.0 0.75672E-06 + 2480.0 0.68864E-06 + 2500.0 0.62694E-06 + 2520.0 0.57100E-06 + 2540.0 0.52027E-06 + 2560.0 0.47423E-06 + 2580.0 0.43245E-06 + 2600.0 0.39450E-06 + 2620.0 0.36002E-06 + 2640.0 0.32869E-06 + 2660.0 0.30020E-06 + 2680.0 0.27428E-06 + 2700.0 0.25071E-06 + 2720.0 0.22924E-06 + 2740.0 0.20969E-06 + 2760.0 0.19189E-06 + 2780.0 0.17566E-06 + 2800.0 0.16086E-06 diff --git a/Testing/Exec/Radiation/kpDB/kpl_ch4.dat b/Testing/Exec/Radiation/kpDB/kpl_ch4.dat new file mode 100755 index 000000000..35a1ba67c --- /dev/null +++ b/Testing/Exec/Radiation/kpDB/kpl_ch4.dat @@ -0,0 +1,126 @@ + 300. 0.04497709 + 320. 0.04844904 + 340. 0.05092656 + 360. 0.05253431 + 380. 0.05341675 + 400. 0.05371585 + 420. 0.05355871 + 440. 0.05305221 + 460. 0.05228213 + 480. 0.05131487 + 500. 0.05020037 + 520. 0.04897542 + 540. 0.04766674 + 560. 0.04629374 + 580. 0.04487073 + 600. 0.04340859 + 620. 0.04191607 + 640. 0.04040063 + 660. 0.03886905 + 680. 0.03732774 + 700. 0.03578290 + 720. 0.03424064 + 740. 0.03270691 + 760. 0.03118750 + 780. 0.02968791 + 800. 0.02821332 + 820. 0.02676851 + 840. 0.02535779 + 860. 0.02398498 + 880. 0.02265336 + 900. 0.02136569 + 920. 0.02012418 + 940. 0.01893054 + 960. 0.01778599 + 980. 0.01669127 + 1000. 0.01564672 + 1020. 0.01465228 + 1040. 0.01370758 + 1060. 0.01281190 + 1080. 0.01196432 + 1100. 0.01116367 + 1120. 0.01040860 + 1140. 0.00969764 + 1160. 0.00902918 + 1180. 0.00840154 + 1200. 0.00781298 + 1220. 0.00726174 + 1240. 0.00674601 + 1260. 0.00626402 + 1280. 0.00581399 + 1300. 0.00539419 + 1320. 0.00500292 + 1340. 0.00463852 + 1360. 0.00429939 + 1380. 0.00398399 + 1400. 0.00369085 + 1420. 0.00341854 + 1440. 0.00316573 + 1460. 0.00293114 + 1480. 0.00271353 + 1500. 0.00251178 + 1520. 0.00232479 + 1540. 0.00215154 + 1560. 0.00199108 + 1580. 0.00184249 + 1600. 0.00170494 + 1620. 0.00157764 + 1640. 0.00145984 + 1660. 0.00135085 + 1680. 0.00125003 + 1700. 0.00115678 + 1720. 0.00107054 + 1740. 0.00099079 + 1760. 0.00091704 + 1780. 0.00084886 + 1800. 0.00078582 + 1820. 0.00072753 + 1840. 0.00067364 + 1860. 0.00062381 + 1880. 0.00057774 + 1900. 0.00053514 + 1920. 0.00049576 + 1940. 0.00045934 + 1960. 0.00042565 + 1980. 0.00039450 + 2000. 0.00036569 + 2020. 0.00033904 + 2040. 0.00031439 + 2060. 0.00029157 + 2080. 0.00027047 + 2100. 0.00025093 + 2120. 0.00023286 + 2140. 0.00021612 + 2160. 0.00020063 + 2180. 0.00018628 + 2200. 0.00017300 + 2220. 0.00016069 + 2240. 0.00014929 + 2260. 0.00013873 + 2280. 0.00012894 + 2300. 0.00011987 + 2320. 0.00011146 + 2340. 0.00010366 + 2360. 0.00009643 + 2380. 0.00008972 + 2400. 0.00008350 + 2420. 0.00007773 + 2440. 0.00007237 + 2460. 0.00006739 + 2480. 0.00006278 + 2500. 0.00005849 + 2520. 0.00005450 + 2540. 0.00005080 + 2560. 0.00004736 + 2580. 0.00004417 + 2600. 0.00004120 + 2620. 0.00003843 + 2640. 0.00003587 + 2660. 0.00003348 + 2680. 0.00003125 + 2700. 0.00002918 + 2720. 0.00002726 + 2740. 0.00002546 + 2760. 0.00002379 + 2780. 0.00002224 + 2800. 0.00002079 diff --git a/Testing/Exec/Radiation/kpDB/kpl_co.dat b/Testing/Exec/Radiation/kpDB/kpl_co.dat new file mode 100755 index 000000000..a068b0ecf --- /dev/null +++ b/Testing/Exec/Radiation/kpDB/kpl_co.dat @@ -0,0 +1,126 @@ + 300. 0.00683693 + 320. 0.00939837 + 340. 0.01222332 + 360. 0.01519314 + 380. 0.01819187 + 400. 0.02111634 + 420. 0.02388186 + 440. 0.02642425 + 460. 0.02869933 + 480. 0.03068088 + 500. 0.03235780 + 520. 0.03373101 + 540. 0.03481054 + 560. 0.03561292 + 580. 0.03615891 + 600. 0.03647178 + 620. 0.03657587 + 640. 0.03649554 + 660. 0.03625444 + 680. 0.03587501 + 700. 0.03537812 + 720. 0.03478290 + 740. 0.03410669 + 760. 0.03336501 + 780. 0.03257164 + 800. 0.03173871 + 820. 0.03087683 + 840. 0.02999518 + 860. 0.02910166 + 880. 0.02820304 + 900. 0.02730505 + 920. 0.02641251 + 940. 0.02552943 + 960. 0.02465915 + 980. 0.02380436 + 1000. 0.02296724 + 1020. 0.02214950 + 1040. 0.02135246 + 1060. 0.02057710 + 1080. 0.01982410 + 1100. 0.01909392 + 1120. 0.01838679 + 1140. 0.01770276 + 1160. 0.01704177 + 1180. 0.01640360 + 1200. 0.01578794 + 1220. 0.01519443 + 1240. 0.01462260 + 1260. 0.01407197 + 1280. 0.01354199 + 1300. 0.01303209 + 1320. 0.01254169 + 1340. 0.01207019 + 1360. 0.01161697 + 1380. 0.01118143 + 1400. 0.01076296 + 1420. 0.01036094 + 1440. 0.00997478 + 1460. 0.00960390 + 1480. 0.00924771 + 1500. 0.00890566 + 1520. 0.00857719 + 1540. 0.00826178 + 1560. 0.00795890 + 1580. 0.00766805 + 1600. 0.00738875 + 1620. 0.00712054 + 1640. 0.00686295 + 1660. 0.00661556 + 1680. 0.00637794 + 1700. 0.00614970 + 1720. 0.00593044 + 1740. 0.00571979 + 1760. 0.00551741 + 1780. 0.00532293 + 1800. 0.00513604 + 1820. 0.00495642 + 1840. 0.00478377 + 1860. 0.00461779 + 1880. 0.00445821 + 1900. 0.00430477 + 1920. 0.00415721 + 1940. 0.00401529 + 1960. 0.00387877 + 1980. 0.00374743 + 2000. 0.00362106 + 2020. 0.00349946 + 2040. 0.00338242 + 2060. 0.00326976 + 2080. 0.00316130 + 2100. 0.00305687 + 2120. 0.00295631 + 2140. 0.00285946 + 2160. 0.00276617 + 2180. 0.00267630 + 2200. 0.00258971 + 2220. 0.00250626 + 2240. 0.00242583 + 2260. 0.00234831 + 2280. 0.00227357 + 2300. 0.00220151 + 2320. 0.00213202 + 2340. 0.00206500 + 2360. 0.00200035 + 2380. 0.00193798 + 2400. 0.00187780 + 2420. 0.00181973 + 2440. 0.00176368 + 2460. 0.00170958 + 2480. 0.00165734 + 2500. 0.00160691 + 2520. 0.00155821 + 2540. 0.00151118 + 2560. 0.00146574 + 2580. 0.00142185 + 2600. 0.00137945 + 2620. 0.00133847 + 2640. 0.00129886 + 2660. 0.00126058 + 2680. 0.00122357 + 2700. 0.00118779 + 2720. 0.00115320 + 2740. 0.00111974 + 2760. 0.00108738 + 2780. 0.00105607 + 2800. 0.00102579 diff --git a/Testing/Exec/Radiation/kpDB/kpl_co2.dat b/Testing/Exec/Radiation/kpDB/kpl_co2.dat new file mode 100755 index 000000000..79d64ff0c --- /dev/null +++ b/Testing/Exec/Radiation/kpDB/kpl_co2.dat @@ -0,0 +1,126 @@ + 300. 0.26371273 + 320. 0.25467201 + 340. 0.25050793 + 360. 0.25078330 + 380. 0.25474820 + 400. 0.26151696 + 420. 0.27019220 + 440. 0.27994417 + 460. 0.29005530 + 480. 0.29993911 + 500. 0.30914198 + 520. 0.31733408 + 540. 0.32429445 + 560. 0.32989369 + 580. 0.33407636 + 600. 0.33684464 + 620. 0.33824375 + 640. 0.33834966 + 660. 0.33725884 + 680. 0.33508018 + 700. 0.33192858 + 720. 0.32792021 + 740. 0.32316893 + 760. 0.31778387 + 780. 0.31186776 + 800. 0.30551590 + 820. 0.29881570 + 840. 0.29184644 + 860. 0.28467945 + 880. 0.27737837 + 900. 0.26999954 + 920. 0.26259251 + 940. 0.25520055 + 960. 0.24786115 + 980. 0.24060659 + 1000. 0.23346441 + 1020. 0.22645790 + 1040. 0.21960657 + 1060. 0.21292650 + 1080. 0.20643078 + 1100. 0.20012981 + 1120. 0.19403165 + 1140. 0.18814224 + 1160. 0.18246571 + 1180. 0.17700455 + 1200. 0.17175985 + 1220. 0.16673146 + 1240. 0.16191814 + 1260. 0.15731770 + 1280. 0.15292716 + 1300. 0.14874281 + 1320. 0.14476035 + 1340. 0.14097494 + 1360. 0.13738134 + 1380. 0.13397392 + 1400. 0.13074675 + 1420. 0.12769365 + 1440. 0.12480828 + 1460. 0.12208412 + 1480. 0.11951458 + 1500. 0.11709298 + 1520. 0.11481262 + 1540. 0.11266682 + 1560. 0.11064891 + 1580. 0.10875228 + 1600. 0.10697042 + 1620. 0.10529690 + 1640. 0.10372540 + 1660. 0.10224977 + 1680. 0.10086398 + 1700. 0.09956216 + 1720. 0.09833864 + 1740. 0.09718789 + 1760. 0.09610460 + 1780. 0.09508366 + 1800. 0.09412012 + 1820. 0.09320926 + 1840. 0.09234658 + 1860. 0.09152775 + 1880. 0.09074867 + 1900. 0.09000545 + 1920. 0.08929439 + 1940. 0.08861200 + 1960. 0.08795499 + 1980. 0.08732026 + 2000. 0.08670491 + 2020. 0.08610624 + 2040. 0.08552170 + 2060. 0.08494896 + 2080. 0.08438582 + 2100. 0.08383027 + 2120. 0.08328047 + 2140. 0.08273471 + 2160. 0.08219145 + 2180. 0.08164927 + 2200. 0.08110692 + 2220. 0.08056323 + 2240. 0.08001719 + 2260. 0.07946790 + 2280. 0.07891457 + 2300. 0.07835649 + 2320. 0.07779309 + 2340. 0.07722386 + 2360. 0.07664837 + 2380. 0.07606631 + 2400. 0.07547741 + 2420. 0.07488147 + 2440. 0.07427838 + 2460. 0.07366806 + 2480. 0.07305050 + 2500. 0.07242575 + 2520. 0.07179388 + 2540. 0.07115502 + 2560. 0.07050933 + 2580. 0.06985702 + 2600. 0.06919830 + 2620. 0.06853345 + 2640. 0.06786272 + 2660. 0.06718642 + 2680. 0.06650488 + 2700. 0.06581842 + 2720. 0.06512739 + 2740. 0.06443214 + 2760. 0.06373305 + 2780. 0.06303048 + 2800. 0.06232481 diff --git a/Testing/Exec/Radiation/kpDB/kpl_h2o.dat b/Testing/Exec/Radiation/kpDB/kpl_h2o.dat new file mode 100755 index 000000000..14689e21f --- /dev/null +++ b/Testing/Exec/Radiation/kpDB/kpl_h2o.dat @@ -0,0 +1,126 @@ + 300. 0.52025064 + 320. 0.44995562 + 340. 0.39460901 + 360. 0.35021992 + 380. 0.31397267 + 400. 0.28386284 + 420. 0.25845138 + 440. 0.23669601 + 460. 0.21783442 + 480. 0.20130280 + 500. 0.18667882 + 520. 0.17364155 + 540. 0.16194326 + 560. 0.15138949 + 580. 0.14182497 + 600. 0.13312362 + 620. 0.12518140 + 640. 0.11791116 + 660. 0.11123899 + 680. 0.10510144 + 700. 0.09944356 + 720. 0.09421739 + 740. 0.08938076 + 760. 0.08489645 + 780. 0.08073146 + 800. 0.07685641 + 820. 0.07324513 + 840. 0.06987420 + 860. 0.06672266 + 880. 0.06377169 + 900. 0.06100440 + 920. 0.05840562 + 940. 0.05596165 + 960. 0.05366016 + 980. 0.05149002 + 1000. 0.04944115 + 1020. 0.04750445 + 1040. 0.04567165 + 1060. 0.04393525 + 1080. 0.04228844 + 1100. 0.04072501 + 1120. 0.03923931 + 1140. 0.03782616 + 1160. 0.03648083 + 1180. 0.03519900 + 1200. 0.03397669 + 1220. 0.03281024 + 1240. 0.03169630 + 1260. 0.03063174 + 1280. 0.02961372 + 1300. 0.02863956 + 1320. 0.02770683 + 1340. 0.02681323 + 1360. 0.02595665 + 1380. 0.02513510 + 1400. 0.02434676 + 1420. 0.02358990 + 1440. 0.02286291 + 1460. 0.02216431 + 1480. 0.02149267 + 1500. 0.02084668 + 1520. 0.02022510 + 1540. 0.01962676 + 1560. 0.01905057 + 1580. 0.01849549 + 1600. 0.01796056 + 1620. 0.01744484 + 1640. 0.01694748 + 1660. 0.01646766 + 1680. 0.01600459 + 1700. 0.01555755 + 1720. 0.01512585 + 1740. 0.01470881 + 1760. 0.01430583 + 1780. 0.01391630 + 1800. 0.01353966 + 1820. 0.01317539 + 1840. 0.01282297 + 1860. 0.01248193 + 1880. 0.01215180 + 1900. 0.01183214 + 1920. 0.01152255 + 1940. 0.01122263 + 1960. 0.01093200 + 1980. 0.01065030 + 2000. 0.01037720 + 2020. 0.01011236 + 2040. 0.00985548 + 2060. 0.00960625 + 2080. 0.00936440 + 2100. 0.00912965 + 2120. 0.00890176 + 2140. 0.00868046 + 2160. 0.00846553 + 2180. 0.00825673 + 2200. 0.00805386 + 2220. 0.00785671 + 2240. 0.00766508 + 2260. 0.00747878 + 2280. 0.00729763 + 2300. 0.00712146 + 2320. 0.00695010 + 2340. 0.00678339 + 2360. 0.00662119 + 2380. 0.00646334 + 2400. 0.00630970 + 2420. 0.00616014 + 2440. 0.00601453 + 2460. 0.00587274 + 2480. 0.00573466 + 2500. 0.00560018 + 2520. 0.00546917 + 2540. 0.00534154 + 2560. 0.00521718 + 2580. 0.00509600 + 2600. 0.00497790 + 2620. 0.00486279 + 2640. 0.00475058 + 2660. 0.00464118 + 2680. 0.00453453 + 2700. 0.00443053 + 2720. 0.00432911 + 2740. 0.00423020 + 2760. 0.00413373 + 2780. 0.00403963 + 2800. 0.00394784 diff --git a/Testing/Exec/Radiation/kpDB/kpl_soot.dat b/Testing/Exec/Radiation/kpDB/kpl_soot.dat new file mode 100755 index 000000000..31b3cfa93 --- /dev/null +++ b/Testing/Exec/Radiation/kpDB/kpl_soot.dat @@ -0,0 +1,126 @@ + 300.000000000000 3077.80995450365 + 320.000000000000 3389.14347543128 + 340.000000000000 3703.94292476713 + 360.000000000000 4021.28926738427 + 380.000000000000 4340.37377312421 + 400.000000000000 4660.48805179194 + 420.000000000000 4981.01430327583 + 440.000000000000 5301.41605415981 + 460.000000000000 5621.22952568398 + 480.000000000000 5940.05569737670 + 500.000000000000 6257.55308013873 + 520.000000000000 6573.43118162418 + 540.000000000000 6887.44462864955 + 560.000000000000 7199.38790160765 + 580.000000000000 7509.09063151284 + 600.000000000000 7816.41340939527 + 620.000000000000 8121.24405898481 + 640.000000000000 8423.49432611040 + 660.000000000000 8723.09694140918 + 680.000000000000 9020.00301641339 + 700.000000000000 9314.17973661617 + 720.000000000000 9605.60831855520 + 740.000000000000 9894.28220120525 + 760.000000000000 10180.2054449843 + 780.000000000000 10463.3913144345 + 800.000000000000 10743.8610231357 + 820.000000000000 11021.6426216528 + 840.000000000000 11296.7700113372 + 860.000000000000 11569.2820686035 + 880.000000000000 11839.2218659302 + 900.000000000000 12106.6359772945 + 920.000000000000 12371.5738570765 + 940.000000000000 12634.0872826849 + 960.000000000000 12894.2298522639 + 980.000000000000 13152.0565298739 + 1000.00000000000 13407.6232314890 + 1020.00000000000 13660.9864460464 + 1040.00000000000 13912.2028866051 + 1060.00000000000 14161.3291674378 + 1080.00000000000 14408.4215035906 + 1100.00000000000 14653.5354300888 + 1120.00000000000 14896.7255385594 + 1140.00000000000 15138.0452295668 + 1160.00000000000 15377.5464794256 + 1180.00000000000 15615.2796206587 + 1200.00000000000 15851.2931356124 + 1220.00000000000 16085.6334630245 + 1240.00000000000 16318.3448175691 + 1260.00000000000 16549.4690225711 + 1280.00000000000 16779.0453562091 + 1300.00000000000 17007.1104115972 + 1320.00000000000 17233.6979711719 + 1340.00000000000 17458.8388958079 + 1360.00000000000 17682.5610290531 + 1380.00000000000 17904.8891168135 + 1400.00000000000 18125.8447427395 + 1420.00000000000 18345.4462794709 + 1440.00000000000 18563.7088557896 + 1460.00000000000 18780.6443396192 + 1480.00000000000 18996.2613366927 + 1500.00000000000 19210.5652045966 + 1520.00000000000 19423.5580817869 + 1540.00000000000 19635.2389310691 + 1560.00000000000 19845.6035969333 + 1580.00000000000 20054.6448760524 + 1600.00000000000 20262.3526001688 + 1620.00000000000 20468.7137305313 + 1640.00000000000 20673.7124629890 + 1660.00000000000 20877.3303428061 + 1680.00000000000 21079.5463882273 + 1700.00000000000 21280.3372218081 + 1720.00000000000 21479.6772085100 + 1740.00000000000 21677.5385995640 + 1760.00000000000 21873.8916811141 + 1780.00000000000 22068.7049266707 + 1800.00000000000 22261.9451524291 + 1820.00000000000 22453.5776745416 + 1840.00000000000 22643.5664674663 + 1860.00000000000 22831.8743225626 + 1880.00000000000 23018.4630061449 + 1900.00000000000 23203.2934162589 + 1920.00000000000 23386.3257374954 + 1940.00000000000 23567.5195932084 + 1960.00000000000 23746.8341945606 + 1980.00000000000 23924.2284858717 + 2000.00000000000 24099.6612857996 + 2020.00000000000 24273.0914239397 + 2040.00000000000 24444.4778724753 + 2060.00000000000 24613.7798725677 + 2080.00000000000 24780.9570552203 + 2100.00000000000 24945.9695563975 + 2120.00000000000 25108.7781262268 + 2140.00000000000 25269.3442321509 + 2160.00000000000 25427.6301559373 + 2180.00000000000 25583.5990844919 + 2200.00000000000 25737.2151944548 + 2220.00000000000 25888.4437305875 + 2240.00000000000 26037.2510779941 + 2260.00000000000 26183.6048282401 + 2280.00000000000 26327.4738394596 + 2300.00000000000 26468.8282905621 + 2320.00000000000 26607.6397296692 + 2340.00000000000 26743.8811169272 + 2360.00000000000 26877.5268618584 + 2380.00000000000 27008.5528554240 + 2400.00000000000 27136.9364969825 + 2420.00000000000 27262.6567163377 + 2440.00000000000 27385.6939910754 + 2460.00000000000 27506.0303593938 + 2480.00000000000 27623.6494286378 + 2500.00000000000 27738.5363797468 + 2520.00000000000 27850.6779678317 + 2540.00000000000 27960.0625190912 + 2560.00000000000 28066.6799242814 + 2580.00000000000 28170.5216289484 + 2600.00000000000 28271.5806206317 + 2620.00000000000 28369.8514132427 + 2640.00000000000 28465.3300288191 + 2660.00000000000 28558.0139768496 + 2680.00000000000 28647.9022313617 + 2700.00000000000 28734.9952059569 + 2720.00000000000 28819.2947269720 + 2740.00000000000 28900.8040049417 + 2760.00000000000 28979.5276045287 + 2780.00000000000 29055.4714130815 + 2800.00000000000 29128.6426079759 diff --git a/Testing/Exec/Radiation/tstPOneMulti.cpp b/Testing/Exec/Radiation/tstPOneMulti.cpp new file mode 100644 index 000000000..5295d19ba --- /dev/null +++ b/Testing/Exec/Radiation/tstPOneMulti.cpp @@ -0,0 +1,321 @@ +#include +#include +#include +#include + +AMREX_GPU_DEVICE AMREX_FORCE_INLINE void +actual_init_coefs( + int i, + int j, + int k, + amrex::Array4 const& rhs, + amrex::Array4 const& sol, + amrex::Array4 const& alpha, + amrex::Array4 const& beta, + amrex::Array4 const& robin_a, + amrex::Array4 const& robin_b, + amrex::Array4 const& robin_f, + amrex::GpuArray const& prob_lo, + amrex::GpuArray const& prob_hi, + amrex::GpuArray const& dx, + amrex::Dim3 const& dlo, + amrex::Dim3 const& dhi, + amrex::Box const& vbx) +{ + amrex::Real const L = 2.0; + amrex::Real const n = 3.0; + amrex::Real const npioverL = n * M_PI / L; + + amrex::Real x = prob_lo[0] + dx[0] * (i + 0.5); + amrex::Real y = prob_lo[1] + dx[1] * (j + 0.5); + amrex::Real z = prob_lo[2] + dx[2] * (k + 0.5); + + beta(i, j, k) = 1.0; + + x = amrex::min(amrex::max(x, prob_lo[0]), prob_hi[0]); + y = amrex::min(amrex::max(y, prob_lo[1]), prob_hi[1]); + z = amrex::min(amrex::max(z, prob_lo[2]), prob_hi[2]); + + amrex::Real sincossin = + std::sin(npioverL * x) * std::cos(npioverL * y) * std::sin(npioverL * z); + + amrex::Real coscossin = + std::cos(npioverL * x) * std::cos(npioverL * y) * std::sin(npioverL * z); + + sol(i, j, k) = sincossin; + + if (vbx.contains(i, j, k)) { + rhs(i, j, k) = (1.0 + npioverL * npioverL) * beta(i, j, k) * sincossin; + alpha(i, j, k) = 1.0; + } + + // Robin BC + if (j >= dlo.y && j <= dhi.y && k >= dlo.z && k <= dhi.z) { + if (i > dhi.x || i < dlo.x) { + robin_a(i, j, k) = -1.0; + robin_b(i, j, k) = -2.0 / 3.0; + + robin_f(i, j, k) = robin_a(i, j, k) * sol(i, j, k) + + robin_b(i, j, k) * npioverL * coscossin; + } + } else if (i >= dlo.x && i <= dhi.x && k >= dlo.z && k <= dhi.z) { + if (j > dhi.y || j < dlo.y) { + robin_a(i, j, k) = -1.0; + robin_b(i, j, k) = -2.0 / 3.0; + + amrex::Real msinsinsin = -std::sin(npioverL * x) * + std::sin(npioverL * y) * std::sin(npioverL * z); + + robin_f(i, j, k) = robin_a(i, j, k) * sol(i, j, k) + + robin_b(i, j, k) * npioverL * msinsinsin; + } + } else if (i >= dlo.x && i <= dhi.x && j >= dlo.y && j <= dhi.y) { + if (k > dhi.z || k < dlo.z) { + robin_a(i, j, k) = -1.0; + robin_b(i, j, k) = -2.0 / 3.0; + + amrex::Real sincoscos = -std::sin(npioverL * x) * std::cos(npioverL * y) * + std::cos(npioverL * z); + + robin_f(i, j, k) = robin_a(i, j, k) * sol(i, j, k) + + robin_b(i, j, k) * npioverL * sincoscos; + } + } +} + +void +initProbABecLaplacian( + amrex::Vector& geom, + amrex::Vector& solution, + amrex::Vector& rhs, + amrex::Vector& exact_solution, + amrex::Vector& acoef, + amrex::Vector& bcoef, + amrex::Vector& robin_a, + amrex::Vector& robin_b, + amrex::Vector& robin_f) +{ + auto nlevels = geom.size(); + for (int ilev = 0; ilev < nlevels; ++ilev) { + auto const prob_lo = geom[ilev].ProbLoArray(); + auto const prob_hi = geom[ilev].ProbHiArray(); + auto const dx = geom[ilev].CellSizeArray(); + auto const dlo = amrex::lbound(geom[ilev].Domain()); + auto const dhi = amrex::ubound(geom[ilev].Domain()); + for (amrex::MFIter mfi(rhs[ilev]); mfi.isValid(); ++mfi) { + amrex::Box const& bx = mfi.validbox(); + amrex::Box const& gbx = amrex::grow(bx, 1); + auto const& rhsfab = rhs[ilev].array(mfi); + auto const& solfab = solution[ilev].array(mfi); + auto const& acfab = acoef[ilev].array(mfi); + auto const& bcfab = bcoef[ilev].array(mfi); + auto const& rafab = robin_a[ilev].array(mfi); + auto const& rbfab = robin_b[ilev].array(mfi); + auto const& rffab = robin_f[ilev].array(mfi); + amrex::ParallelFor( + gbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { + actual_init_coefs( + i, j, k, rhsfab, solfab, acfab, bcfab, rafab, rbfab, rffab, prob_lo, + prob_hi, dx, dlo, dhi, bx); + }); + } + + amrex::MultiFab::Copy(exact_solution[ilev], solution[ilev], 0, 0, 1, 0); + solution[ilev].setVal(0.0, 0, 1, amrex::IntVect(0)); + } +} + +void +initMeshandData( + PeleRad::AMRParam const& amrpp, + amrex::Vector& geom, + amrex::Vector& grids, + amrex::Vector& dmap, + amrex::Vector& solution, + amrex::Vector& rhs, + amrex::Vector& exact_solution, + amrex::Vector& acoef, + amrex::Vector& bcoef, + amrex::Vector& robin_a, + amrex::Vector& robin_b, + amrex::Vector& robin_f) +{ + int const nlevels = amrpp.max_level_ + 1; + int const ref_ratio = amrpp.ref_ratio_; + int const n_cell = amrpp.n_cell_; + int const max_grid_size = amrpp.max_grid_size_; + + // initialize mesh + // std::cout << "initialize the mesh" << std::endl; + geom.resize(nlevels); + grids.resize(nlevels); + dmap.resize(nlevels); + + amrex::RealBox rb( + {AMREX_D_DECL(-1.0, -1.0, -1.0)}, {AMREX_D_DECL(1.0, 1.0, 1.0)}); + amrex::Array is_periodic{AMREX_D_DECL(0, 0, 0)}; + amrex::Geometry::Setup(&rb, 0, is_periodic.data()); + amrex::Box domain0( + amrex::IntVect{AMREX_D_DECL(0, 0, 0)}, + amrex::IntVect{AMREX_D_DECL(n_cell - 1, n_cell - 1, n_cell - 1)}); + + amrex::Box domain = domain0; + + for (int ilev = 0; ilev < nlevels; ++ilev) { + geom[ilev].define(domain); + domain.refine(ref_ratio); + } + + domain = domain0; + for (int ilev = 0; ilev < nlevels; ++ilev) { + grids[ilev].define(domain); + grids[ilev].maxSize(max_grid_size); + domain.refine(ref_ratio); + } + + // initialize variables + // std::cout << "initialize the data" << std::endl; + + solution.resize(nlevels); + rhs.resize(nlevels); + exact_solution.resize(nlevels); + acoef.resize(nlevels); + bcoef.resize(nlevels); + robin_a.resize(nlevels); + robin_b.resize(nlevels); + robin_f.resize(nlevels); + + amrex::IntVect ng = amrex::IntVect{1}; + + for (int ilev = 0; ilev < nlevels; ++ilev) { + dmap[ilev].define(grids[ilev]); + solution[ilev].define(grids[ilev], dmap[ilev], 1, ng); + rhs[ilev].define(grids[ilev], dmap[ilev], 1, 0); + exact_solution[ilev].define(grids[ilev], dmap[ilev], 1, ng); + acoef[ilev].define(grids[ilev], dmap[ilev], 1, 0); + bcoef[ilev].define(grids[ilev], dmap[ilev], 1, ng); + robin_a[ilev].define(grids[ilev], dmap[ilev], 1, ng); + robin_b[ilev].define(grids[ilev], dmap[ilev], 1, ng); + robin_f[ilev].define(grids[ilev], dmap[ilev], 1, ng); + } + + initProbABecLaplacian( + geom, solution, rhs, exact_solution, acoef, bcoef, robin_a, robin_b, + robin_f); +} + +amrex::Real +check_norm(amrex::MultiFab const& phi, amrex::MultiFab const& exact) +{ + amrex::MultiFab mf(phi.boxArray(), phi.DistributionMap(), 1, 0); + amrex::MultiFab::Copy(mf, phi, 0, 0, 1, 0); + + amrex::MultiFab::Subtract(mf, exact, 0, 0, 1, 0); + + amrex::Real L0norm = mf.norm0(); + amrex::Real L1norm = mf.norm1(); + std::cout << " L0 norm:" << L0norm << ", L1 norm:" << L1norm << std::endl; + + return L1norm; +} + +int +main(int argc, char* argv[]) +{ + amrex::Initialize(argc, argv); + amrex::ParmParse pp; + PeleRad::AMRParam amrpp(pp); + PeleRad::MLMGParam mlmgpp(pp); + + bool const write = false; + int const nlevels = amrpp.max_level_ + 1; + int const ref_ratio = amrpp.ref_ratio_; + int const composite_solve = mlmgpp.composite_solve_; + + amrex::Vector geom; + amrex::Vector grids; + amrex::Vector dmap; + + amrex::Vector solution; + amrex::Vector rhs; + amrex::Vector exact_solution; + + amrex::Vector acoef; + amrex::Vector bcoef; + amrex::Vector robin_a; + amrex::Vector robin_b; + amrex::Vector robin_f; + + // std::cout << "initialize data ... \n"; + initMeshandData( + amrpp, geom, grids, dmap, solution, rhs, exact_solution, acoef, bcoef, + robin_a, robin_b, robin_f); + // std::cout << "construct the PDE ... \n"; + if (composite_solve) { + + std::cout << "composite solve ... \n"; + PeleRad::POneMulti rte( + mlmgpp, geom, grids, dmap, solution, rhs, acoef, bcoef, robin_a, robin_b, + robin_f); + rte.solve(); + } else { + std::cout << "level by level solve ... \n"; + PeleRad::POneMultiLevbyLev rte( + mlmgpp, ref_ratio, geom, grids, dmap, solution, rhs, acoef, bcoef, + robin_a, robin_b, robin_f); + rte.solve(); + } + + amrex::Real eps = 0.0; + amrex::Real eps_max = 0.0; + for (int ilev = 0; ilev < nlevels; ++ilev) { + auto dx = geom[ilev].CellSize(); + amrex::Real dvol = AMREX_D_TERM(dx[0], *dx[1], *dx[2]); + eps = check_norm(solution[ilev], exact_solution[ilev]); + eps *= dvol; + std::cout << "Level=" << ilev << ", normalized L1 norm:" << eps + << std::endl; + if (eps > eps_max) + eps_max = eps; + } + + // plot results + if (write) { + amrex::Vector plotmf(nlevels); + + for (int ilev = 0; ilev < nlevels; ++ilev) { + // std::cout << "write the results ... \n"; + plotmf[ilev].define(grids[ilev], dmap[ilev], 4, 0); + amrex::MultiFab::Copy(plotmf[ilev], solution[ilev], 0, 0, 1, 0); + amrex::MultiFab::Copy(plotmf[ilev], rhs[ilev], 0, 1, 1, 0); + amrex::MultiFab::Copy(plotmf[ilev], exact_solution[ilev], 0, 2, 1, 0); + amrex::MultiFab::Copy(plotmf[ilev], solution[ilev], 0, 3, 1, 0); + amrex::MultiFab::Subtract(plotmf[ilev], plotmf[ilev], 2, 3, 1, 0); + + // For amrvis + /* + amrex::writeFabs(solution[ilev], + "solution_lev"+std::to_string(ilev)); amrex::writeFabs(bcoef[ilev], + "bcoef_lev"+std::to_string(ilev)); amrex::writeFabs(robin_a[ilev], + "robin_a_lev"+std::to_string(ilev)); amrex::writeFabs(robin_b[ilev], + "robin_b_lev"+std::to_string(ilev)); amrex::writeFabs(robin_f[ilev], + "robin_f_lev"+std::to_string(ilev)); + */ + } + + auto const plot_file_name = amrpp.plot_file_name_; + amrex::WriteMultiLevelPlotfile( + plot_file_name, nlevels, amrex::GetVecOfConstPtrs(plotmf), + {"phi", "rhs", "exact", "error"}, geom, 0.0, + amrex::Vector(nlevels, 0), + amrex::Vector(nlevels, amrex::IntVect{ref_ratio})); + } + + amrex::Finalize(); + + if (eps_max < 1e-2) { + return (0); + } else { + return (-1); + } +} diff --git a/Testing/Exec/Radiation/tstPOneMultiAF.cpp b/Testing/Exec/Radiation/tstPOneMultiAF.cpp new file mode 100644 index 000000000..ac851032e --- /dev/null +++ b/Testing/Exec/Radiation/tstPOneMultiAF.cpp @@ -0,0 +1,396 @@ +#include +#include + +#include +#include +#include +#include + +AMREX_GPU_DEVICE AMREX_FORCE_INLINE void +initGasField( + int i, + int j, + int k, + amrex::Array4 const& y_co2, + amrex::Array4 const& y_h2o, + amrex::Array4 const& y_co, + amrex::Array4 const& fv_soot, + amrex::Array4 const& temp, + amrex::Array4 const& pressure, + amrex::GpuArray const& dx, + amrex::GpuArray const& plo, + amrex::GpuArray const& phi) +{ + amrex::Real coef = 100; + amrex::Real xc = (phi[0] + plo[0]) * 0.5; + amrex::Real yc = (phi[1] + plo[1]) * 0.5; + + amrex::Real x = plo[0] + (i + 0.5) * dx[0]; + amrex::Real y = plo[1] + (j + 0.5) * dx[1]; + amrex::Real z = plo[2] + (k + 0.5) * dx[2]; + + amrex::Real r = std::sqrt((x - xc) * (x - xc) + (y - yc) * (y - yc)); + + z /= coef; + r /= coef; + + amrex::Real expr = std::exp( + -(4.0 * r / (0.05 + 0.1 * 4.0 * z)) * (4.0 * r / (0.05 + 0.1 * 4.0 * z))); + amrex::Real expTz = std::exp( + -((4.0 * z - 1.3) / (0.7 + 0.5 * 4.0 * z)) * + ((4.0 * z - 1.3) / (0.7 + 0.5 * 4.0 * z))); + + temp(i, j, k) = 300.0 + 1700.0 * expr * expTz; + + pressure(i, j, k) = 1.0e5; // Pa + + amrex::Real expSoot = + std::exp(-((4.0 * z - 1.0) / 0.7) * ((4.0 * z - 1.0) / 0.7)); + + fv_soot(i, j, k) = 1e-6 * expr * expSoot; + + amrex::Real expCO2z = std::exp( + -((4.0 * z - 1.1) / (0.6 + 0.5 * 4.0 * z)) * + ((4.0 * z - 1.1) / (0.6 + 0.5 * 4.0 * z))); + + y_co2(i, j, k) = 0.1 * expr * expCO2z; + + amrex::Real expH2Oz = std::exp( + -((4.0 * z - 1.0) / (0.7 + 0.5 * 4.0 * z)) * + ((4.0 * z - 1.0) / (0.7 + 0.5 * 4.0 * z))); + + y_h2o(i, j, k) = 0.2 * expr * expH2Oz; + + amrex::Real expCOz = + std::exp(-((4.0 * z - 1.0) / 0.7) * ((4.0 * z - 1.0) / 0.7)); + + y_co(i, j, k) = 0.09 * expr * expCOz; +} + +AMREX_GPU_DEVICE AMREX_FORCE_INLINE void +actual_init_coefs( + int i, + int j, + int k, + amrex::Array4 const& rhs, + amrex::Array4 const& alpha, + amrex::Array4 const& beta, + amrex::Dim3 const& dlo, + amrex::Dim3 const& dhi, + amrex::Box const& vbx, + amrex::Array4 const& absc, + amrex::Array4 const& T) +{ + beta(i, j, k) = 1.0; + if (vbx.contains(i, j, k)) { + amrex::Real ka = std::max(0.01, absc(i, j, k) * 100); + beta(i, j, k) = 1.0 / ka; + rhs(i, j, k) = 4.0 * ka * 5.67e-8 * std::pow(T(i, j, k), 4.0); // SI + alpha(i, j, k) = ka; + } +} + +AMREX_GPU_DEVICE AMREX_FORCE_INLINE void +actual_init_bc_coefs( + int i, + int j, + int k, + amrex::Array4 const& alpha, + amrex::Array4 const& beta, + amrex::Array4 const& robin_a, + amrex::Array4 const& robin_b, + amrex::Array4 const& robin_f, + amrex::Dim3 const& dlo, + amrex::Dim3 const& dhi, + amrex::Array4 const& T) +{ + // Robin BC + bool robin_cell = false; + if (j >= dlo.y && j <= dhi.y && k >= dlo.z && k <= dhi.z) { + if (i > dhi.x || i < dlo.x) { + robin_cell = true; + } + } else if (i >= dlo.x && i <= dhi.x && k >= dlo.z && k <= dhi.z) { + if (j > dhi.y || j < dlo.y) { + robin_cell = true; + } + } + + if (robin_cell) { + robin_a(i, j, k) = -1.0 / beta(i, j, k); + robin_b(i, j, k) = -2.0 / 3.0; + robin_f(i, j, k) = 0.0; + } +} + +void +initProbABecLaplacian( + amrex::Vector& geom, + amrex::Vector& solution, + amrex::Vector& rhs, + amrex::Vector& acoef, + amrex::Vector& bcoef, + amrex::Vector& robin_a, + amrex::Vector& robin_b, + amrex::Vector& robin_f, + amrex::Vector& y_co2, + amrex::Vector& y_h2o, + amrex::Vector& y_co, + amrex::Vector& soot_fv_rad, + amrex::Vector& temperature, + amrex::Vector& pressure, + amrex::Vector& absc) +{ + using PeleRad::PlanckMean; + using PeleRad::RadProp::getRadPropGas; + using PeleRad::RadProp::getRadPropSoot; + + std::string data_path("kpDB/"); + PlanckMean radprop(data_path); + auto const& kpco2 = radprop.kpco2(); + auto const& kph2o = radprop.kph2o(); + auto const& kpco = radprop.kpco(); + auto const& kpsoot = radprop.kpsoot(); + + auto nlevels = geom.size(); + + for (int ilev = 0; ilev < nlevels; ++ilev) { + auto const prob_lo = geom[ilev].ProbLoArray(); + auto const prob_hi = geom[ilev].ProbHiArray(); + auto const dx = geom[ilev].CellSizeArray(); + auto const dlo = amrex::lbound(geom[ilev].Domain()); + auto const dhi = amrex::ubound(geom[ilev].Domain()); + + for (amrex::MFIter mfi(rhs[ilev]); mfi.isValid(); ++mfi) { + amrex::Box const& bx = mfi.validbox(); + amrex::Box const& gbx = amrex::grow(bx, 1); + auto const& Yco2 = y_co2[ilev].array(mfi); + auto const& Yh2o = y_h2o[ilev].array(mfi); + auto const& Yco = y_co[ilev].array(mfi); + auto const& fv = soot_fv_rad[ilev].array(mfi); + auto const& T = temperature[ilev].array(mfi); + auto const& P = pressure[ilev].array(mfi); + auto const& kappa = absc[ilev].array(mfi); + + amrex::ParallelFor( + bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { + initGasField( + i, j, k, Yco2, Yh2o, Yco, fv, T, P, dx, prob_lo, prob_hi); + getRadPropGas( + i, j, k, Yco2, Yh2o, Yco, T, P, kappa, kpco2, kph2o, kpco); + }); + + // if soot exists + amrex::ParallelFor( + bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { + getRadPropSoot(i, j, k, fv, T, kappa, kpsoot); + }); + + auto const& rhsfab = rhs[ilev].array(mfi); + auto const& acfab = acoef[ilev].array(mfi); + auto const& bcfab = bcoef[ilev].array(mfi); + auto const& rafab = robin_a[ilev].array(mfi); + auto const& rbfab = robin_b[ilev].array(mfi); + auto const& rffab = robin_f[ilev].array(mfi); + amrex::ParallelFor( + gbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { + actual_init_coefs( + i, j, k, rhsfab, acfab, bcfab, dlo, dhi, bx, kappa, T); + }); + + amrex::ParallelFor( + gbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { + actual_init_bc_coefs( + i, j, k, acfab, bcfab, rafab, rbfab, rffab, dlo, dhi, T); + }); + } + solution[ilev].setVal(0.0, 0, 1, amrex::IntVect(0)); + } +} + +void +initMeshandData( + PeleRad::AMRParam const& amrpp, + amrex::Vector& geom, + amrex::Vector& grids, + amrex::Vector& dmap, + amrex::Vector& solution, + amrex::Vector& rhs, + amrex::Vector& acoef, + amrex::Vector& bcoef, + amrex::Vector& robin_a, + amrex::Vector& robin_b, + amrex::Vector& robin_f, + amrex::Vector& y_co2, + amrex::Vector& y_h2o, + amrex::Vector& y_co, + amrex::Vector& soot_fv_rad, + amrex::Vector& temperature, + amrex::Vector& pressure, + amrex::Vector& absc) +{ + int const nlevels = amrpp.max_level_ + 1; + int const ref_ratio = amrpp.ref_ratio_; + int const n_cell = amrpp.n_cell_; + int const max_grid_size = amrpp.max_grid_size_; + + // initialize mesh + geom.resize(nlevels); + grids.resize(nlevels); + dmap.resize(nlevels); + + amrex::RealBox rb( + {AMREX_D_DECL(0.0, 0.0, 0.0)}, {AMREX_D_DECL(12.5, 12.5, 75)}); + amrex::Array is_periodic{AMREX_D_DECL(0, 0, 0)}; + amrex::Geometry::Setup(&rb, 0, is_periodic.data()); + + std::vector npts{n_cell, n_cell, 6 * n_cell}; + amrex::Box domain0( + amrex::IntVect{AMREX_D_DECL(0, 0, 0)}, + amrex::IntVect{AMREX_D_DECL(npts[0] - 1, npts[1] - 1, npts[2] - 1)}); + + amrex::Box domain = domain0; + + for (int ilev = 0; ilev < nlevels; ++ilev) { + geom[ilev].define(domain); + domain.refine(ref_ratio); + } + + domain = domain0; + + for (int ilev = 0; ilev < nlevels; ++ilev) { + grids[ilev].define(domain); + grids[ilev].maxSize(max_grid_size); + domain.refine(ref_ratio); + } + + // initialize variables + solution.resize(nlevels); + rhs.resize(nlevels); + acoef.resize(nlevels); + bcoef.resize(nlevels); + robin_a.resize(nlevels); + robin_b.resize(nlevels); + robin_f.resize(nlevels); + + y_co2.resize(nlevels); + y_h2o.resize(nlevels); + y_co.resize(nlevels); + soot_fv_rad.resize(nlevels); + temperature.resize(nlevels); + pressure.resize(nlevels); + absc.resize(nlevels); + + amrex::IntVect ng = amrex::IntVect{1}; + + for (int ilev = 0; ilev < nlevels; ++ilev) { + dmap[ilev].define(grids[ilev]); + solution[ilev].define(grids[ilev], dmap[ilev], 1, ng); + rhs[ilev].define(grids[ilev], dmap[ilev], 1, 0); + acoef[ilev].define(grids[ilev], dmap[ilev], 1, 0); + bcoef[ilev].define(grids[ilev], dmap[ilev], 1, ng); + robin_a[ilev].define(grids[ilev], dmap[ilev], 1, ng); + robin_b[ilev].define(grids[ilev], dmap[ilev], 1, ng); + robin_f[ilev].define(grids[ilev], dmap[ilev], 1, ng); + + y_co2[ilev].define(grids[ilev], dmap[ilev], 1, 0); + y_h2o[ilev].define(grids[ilev], dmap[ilev], 1, 0); + y_co[ilev].define(grids[ilev], dmap[ilev], 1, 0); + soot_fv_rad[ilev].define(grids[ilev], dmap[ilev], 1, 0); + temperature[ilev].define(grids[ilev], dmap[ilev], 1, 0); + pressure[ilev].define(grids[ilev], dmap[ilev], 1, 0); + absc[ilev].define(grids[ilev], dmap[ilev], 1, 0); + } + + initProbABecLaplacian( + geom, solution, rhs, acoef, bcoef, robin_a, robin_b, robin_f, y_co2, y_h2o, + y_co, soot_fv_rad, temperature, pressure, absc); + + for (int ilev = 0; ilev < nlevels; ++ilev) { + bcoef[ilev].FillBoundary(); + } +} + +int +main(int argc, char* argv[]) +{ + amrex::Initialize(argc, argv); + amrex::ParmParse pp; + PeleRad::AMRParam amrpp(pp); + PeleRad::MLMGParam mlmgpp(pp); + + bool const write = false; + int const nlevels = amrpp.max_level_ + 1; + int const ref_ratio = amrpp.ref_ratio_; + int const composite_solve = mlmgpp.composite_solve_; + + amrex::Vector geom; + amrex::Vector grids; + amrex::Vector dmap; + + amrex::Vector solution; + amrex::Vector rhs; + + amrex::Vector acoef; + amrex::Vector bcoef; + amrex::Vector robin_a; + amrex::Vector robin_b; + amrex::Vector robin_f; + + amrex::Vector y_co2; + amrex::Vector y_h2o; + amrex::Vector y_co; + amrex::Vector soot_fv_rad; + amrex::Vector temperature; + amrex::Vector pressure; + amrex::Vector absc; + + std::cout << "initialize data ... \n"; + initMeshandData( + amrpp, geom, grids, dmap, solution, rhs, acoef, bcoef, robin_a, robin_b, + robin_f, y_co2, y_h2o, y_co, soot_fv_rad, temperature, pressure, absc); + + std::cout << "construct the PDE ... \n"; + + if (composite_solve) { + PeleRad::POneMulti rte( + mlmgpp, geom, grids, dmap, solution, rhs, acoef, bcoef, robin_a, robin_b, + robin_f); + rte.solve(); + } else { + PeleRad::POneMultiLevbyLev rte( + mlmgpp, ref_ratio, geom, grids, dmap, solution, rhs, acoef, bcoef, + robin_a, robin_b, robin_f); + rte.solve(); + } + + // plot results + if (write) { + amrex::Vector plotmf(nlevels); + + for (int ilev = 0; ilev < nlevels; ++ilev) { + std::cout << "write the results ... \n"; + plotmf[ilev].define(grids[ilev], dmap[ilev], 9, 0); + amrex::MultiFab::Copy(plotmf[ilev], solution[ilev], 0, 0, 1, 0); + amrex::MultiFab::Copy(plotmf[ilev], rhs[ilev], 0, 1, 1, 0); + amrex::MultiFab::Copy(plotmf[ilev], acoef[ilev], 0, 2, 1, 0); + amrex::MultiFab::Copy(plotmf[ilev], bcoef[ilev], 0, 3, 1, 0); + + amrex::MultiFab::Copy(plotmf[ilev], y_co2[ilev], 0, 4, 1, 0); + amrex::MultiFab::Copy(plotmf[ilev], y_h2o[ilev], 0, 5, 1, 0); + amrex::MultiFab::Copy(plotmf[ilev], y_co[ilev], 0, 6, 1, 0); + amrex::MultiFab::Copy(plotmf[ilev], soot_fv_rad[ilev], 0, 7, 1, 0); + amrex::MultiFab::Copy(plotmf[ilev], temperature[ilev], 0, 8, 1, 0); + } + + auto const plot_file_name = amrpp.plot_file_name_; + amrex::WriteMultiLevelPlotfile( + plot_file_name, nlevels, amrex::GetVecOfConstPtrs(plotmf), + {"solution", "rhs", "acoef", "bcoef", "Y_co2", "Y_h2o", "Y_co", + "Soot_fv_rad", "Temperature"}, + geom, 0.0, amrex::Vector(nlevels, 0), + amrex::Vector(nlevels, amrex::IntVect{ref_ratio})); + } + amrex::Finalize(); + return (0); +} diff --git a/Testing/Exec/Radiation/tstPOneMultiEB.cpp b/Testing/Exec/Radiation/tstPOneMultiEB.cpp new file mode 100644 index 000000000..6390a8391 --- /dev/null +++ b/Testing/Exec/Radiation/tstPOneMultiEB.cpp @@ -0,0 +1,325 @@ +#include +#include +#include +#include +#include + +#include +#include +#include +#include + +AMREX_GPU_DEVICE AMREX_FORCE_INLINE void +actual_init_coefs_eb( + int i, + int j, + int k, + amrex::Array4 const& phi, + amrex::Array4 const& phi_exact, + amrex::Array4 const& rhs, + amrex::Array4 const& acoef, + amrex::Array4 const& bx, + amrex::Array4 const& flag, + amrex::GpuArray const& dx, + amrex::Box const& vbx) +{ + amrex::Real const L = std::sqrt(2.0); + amrex::Real const n = 3.0; + amrex::Real npioverL = n * M_PI / L; + + amrex::Real pioverfour = M_PI / 4.0; + amrex::Real cospioverfour = std::cos(pioverfour); + amrex::Real sinpioverfour = std::sin(pioverfour); + + if (vbx.contains(i, j, k)) { + bx(i, j, k) = 1.0; + phi(i, j, k) = 0.0; + + if (flag(i, j, k).isCovered()) { + rhs(i, j, k) = 0.0; + phi_exact(i, j, k) = 0.0; + } else { + amrex::Real x = dx[0] * (i + 0.5) - 0.5; + amrex::Real y = dx[1] * (j + 0.5) - 0.5; + amrex::Real z = dx[2] * k; + + amrex::Real xp = x * cospioverfour + y * sinpioverfour; + amrex::Real yp = -x * sinpioverfour + y * cospioverfour; + + amrex::Real sincossin = std::sin(npioverL * xp) * + std::cos(npioverL * yp) * std::sin(npioverL * z); + + rhs(i, j, k) = (1.0 + npioverL * npioverL) * bx(i, j, k) * sincossin; + acoef(i, j, k) = 1.0; + + phi_exact(i, j, k) = sincossin; + } + } +} + +void +initProbABecLaplacian( + amrex::Vector& geom, + amrex::Vector>& factory, + amrex::Vector& solution, + amrex::Vector& rhs, + amrex::Vector& exact_solution, + amrex::Vector& acoef, + amrex::Vector& bcoef) +{ + int const nlevels = geom.size(); + for (int ilev = 0; ilev < nlevels; ++ilev) { + amrex::FabArray const& flags = + factory[ilev]->getMultiEBCellFlagFab(); + auto const dx = geom[ilev].CellSizeArray(); + amrex::Box const& domainbox = geom[ilev].Domain(); + + for (amrex::MFIter mfi(rhs[ilev]); mfi.isValid(); ++mfi) { + amrex::Box const& bx = mfi.validbox(); + amrex::Box const& nbx = amrex::surroundingNodes(bx); + // amrex::Box const& gbx = amrex::grow(bx, 1); + amrex::Array4 const& phi_arr = solution[ilev].array(mfi); + amrex::Array4 const& phi_ex_arr = + exact_solution[ilev].array(mfi); + amrex::Array4 const& rhs_arr = rhs[ilev].array(mfi); + + amrex::Array4 const& b_arr = bcoef[ilev].array(mfi); + + auto fabtyp = flags[mfi].getType(bx); + if (fabtyp == amrex::FabType::covered) { + std::cout << " amrex::FabType::covered == fabtyp \n"; + } else if (fabtyp == amrex::FabType::regular) { + std::cout << " amrex::FabType::regular == fabtyp \n"; + } else { + // std::cout << " amrex::FabType else \n"; + amrex::Array4 const& acoef_arr = acoef[ilev].array(mfi); + amrex::Array4 const& flag_arr = + flags.const_array(mfi); + // amrex::Array4 const& cent_arr + // = cent.const_array(mfi); + // amrex::Array4 const& bcent_arr + // = bcent.const_array(mfi); + + amrex::ParallelFor( + nbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { + actual_init_coefs_eb( + i, j, k, phi_arr, phi_ex_arr, rhs_arr, acoef_arr, b_arr, flag_arr, + dx, bx); + }); + } + } + } +} + +amrex::Real +check_norm( + amrex::MultiFab const& phi, + amrex::MultiFab const& exact, + amrex::MultiFab const& vfrc) +{ + amrex::MultiFab mf(phi.boxArray(), phi.DistributionMap(), 1, 0); + amrex::MultiFab::Copy(mf, phi, 0, 0, 1, 0); + + amrex::MultiFab::Subtract(mf, exact, 0, 0, 1, 0); + amrex::MultiFab::Multiply(mf, vfrc, 0, 0, 1, 0); + + amrex::Real L0norm = mf.norm0(); + amrex::Real L1norm = mf.norm1(); + std::cout << " L0 norm:" << L0norm << ", L1 norm:" << L1norm << std::endl; + + return L1norm; +} + +int +main(int argc, char* argv[]) +{ + amrex::Initialize(argc, argv); + amrex::ParmParse pp; + PeleRad::AMRParam amrpp(pp); + PeleRad::MLMGParam mlmgpp(pp); + + bool const write = false; + int const n_cell = amrpp.n_cell_; + int const nlevels = amrpp.max_level_ + 1; + int const max_level = amrpp.max_level_; + int const ref_ratio = amrpp.ref_ratio_; + int const composite_solve = mlmgpp.composite_solve_; + + int const max_grid_size = amrpp.max_grid_size_; + + amrex::Vector geom; + amrex::Vector grids; + amrex::Vector dmap; + + amrex::Vector solution; + amrex::Vector rhs; + amrex::Vector exact_solution; + + amrex::Vector acoef; + amrex::Vector bcoef; + amrex::Vector robin_a; + amrex::Vector robin_b; + amrex::Vector robin_f; + + amrex::Vector> factory; + + geom.resize(nlevels); + grids.resize(nlevels); + + amrex::RealBox rb( + {AMREX_D_DECL(-1.0, -1.0, -1.0)}, {AMREX_D_DECL(1.0, 1.0, 1.0)}); + amrex::Array is_periodic{AMREX_D_DECL(1, 1, 1)}; + amrex::Geometry::Setup(&rb, 0, is_periodic.data()); + amrex::Box domain0( + amrex::IntVect{AMREX_D_DECL(0, 0, 0)}, + amrex::IntVect{AMREX_D_DECL(n_cell - 1, n_cell - 1, n_cell - 1)}); + amrex::Box domain = domain0; + for (int ilev = 0; ilev < nlevels; ++ilev) { + geom[ilev].define(domain); + domain.refine(ref_ratio); + } + + grids[0].define(domain0); + grids[0].maxSize(max_grid_size); + + // rotated box + int const max_coarsening_level = mlmgpp.max_coarsening_level_; + amrex::Real const la = std::sqrt(2.0) / 2.0; + amrex::Real const shift = std::sqrt(2.0) / 3.0; + amrex::EB2::BoxIF box( + {AMREX_D_DECL(-la, -la, -1.0)}, {AMREX_D_DECL(la, la, -1.0 + 4.0 * shift)}, + true); + auto gshop = amrex::EB2::makeShop(amrex::EB2::translate( + amrex::EB2::rotate( + amrex::EB2::translate(box, {AMREX_D_DECL(-0.0, -0.0, -0.0)}), + std::atan(1.0) * 1.0, 2), + {AMREX_D_DECL(0.0, 0.0, 0.0)})); + amrex::EB2::Build( + gshop, geom.back(), max_level, max_level + max_coarsening_level); + + // refine grid + for (int ilev = 1; ilev < nlevels; ++ilev) { + amrex::EB2::IndexSpace const& eb_is = amrex::EB2::IndexSpace::top(); + amrex::EB2::Level const& eb_level = eb_is.getLevel(geom[ilev]); + amrex::BoxList bl = eb_level.boxArray().boxList(); + amrex::Box const& domain = geom[ilev].Domain(); + for (auto& b : bl) { + b &= domain; + } + grids[ilev].define(bl); + } + + dmap.resize(nlevels); + factory.resize(nlevels); + solution.resize(nlevels); + exact_solution.resize(nlevels); + rhs.resize(nlevels); + acoef.resize(nlevels); + bcoef.resize(nlevels); + robin_a.resize(nlevels); + robin_b.resize(nlevels); + robin_f.resize(nlevels); + + for (int ilev = 0; ilev < nlevels; ++ilev) { + dmap[ilev].define(grids[ilev]); + amrex::EB2::IndexSpace const& eb_is = amrex::EB2::IndexSpace::top(); + amrex::EB2::Level const& eb_level = eb_is.getLevel(geom[ilev]); + factory[ilev] = std::make_unique( + eb_level, geom[ilev], grids[ilev], dmap[ilev], + amrex::Vector{2, 2, 2}, amrex::EBSupport::full); + + solution[ilev].define( + grids[ilev], dmap[ilev], 1, 1, amrex::MFInfo(), *factory[ilev]); + exact_solution[ilev].define( + grids[ilev], dmap[ilev], 1, 0, amrex::MFInfo(), *factory[ilev]); + rhs[ilev].define( + grids[ilev], dmap[ilev], 1, 0, amrex::MFInfo(), *factory[ilev]); + acoef[ilev].define( + grids[ilev], dmap[ilev], 1, 0, amrex::MFInfo(), *factory[ilev]); + bcoef[ilev].define( + grids[ilev], dmap[ilev], 1, 1, amrex::MFInfo(), *factory[ilev]); + robin_a[ilev].define( + grids[ilev], dmap[ilev], 1, 0, amrex::MFInfo(), *factory[ilev]); + robin_b[ilev].define( + grids[ilev], dmap[ilev], 1, 0, amrex::MFInfo(), *factory[ilev]); + robin_f[ilev].define( + grids[ilev], dmap[ilev], 1, 0, amrex::MFInfo(), *factory[ilev]); + + solution[ilev].setVal(0.0); + rhs[ilev].setVal(0.0); + acoef[ilev].setVal(1.0); + bcoef[ilev].setVal(1.0); + robin_a[ilev].setVal(0.0); + robin_b[ilev].setVal(0.0); + robin_f[ilev].setVal(0.0); + } + + initProbABecLaplacian( + geom, factory, solution, rhs, exact_solution, acoef, bcoef); + + // Get vector of EB factory + std::cout << "construct the PDE ... \n"; + amrex::Vector ebfactVec; + for (int ilev = 0; ilev < nlevels; ++ilev) { + ebfactVec.push_back(factory[ilev].get()); + } + PeleRad::POneMultiEB rte( + mlmgpp, geom, grids, dmap, ebfactVec, solution, rhs, acoef, bcoef, robin_a, + robin_b, robin_f); + + std::cout << "solve the PDE ... \n"; + rte.solve(); + + amrex::Real eps = 0.0; + amrex::Real eps_max = 0.0; + for (int ilev = 0; ilev < nlevels; ++ilev) { + auto const& vfrc = factory[ilev]->getVolFrac(); + auto dx = geom[ilev].CellSize(); + amrex::Real dvol = AMREX_D_TERM(dx[0], *dx[1], *dx[2]); + eps = check_norm(solution[ilev], exact_solution[ilev], vfrc); + eps *= dvol / 8.0; // total volume/cell volume + std::cout << "Level=" << ilev << ", normalized L1 norm:" << eps + << std::endl; + if (eps > eps_max) + eps_max = eps; + } + + // plot results + if (write) { + amrex::Vector plotmf(nlevels); + std::cout << "write the results ... \n"; + for (int ilev = 0; ilev < nlevels; ++ilev) { + auto const& vfrc = factory[ilev]->getVolFrac(); + plotmf[ilev].define(grids[ilev], dmap[ilev], 6, 0); + amrex::MultiFab::Copy(plotmf[ilev], solution[ilev], 0, 0, 1, 0); + amrex::MultiFab::Copy(plotmf[ilev], exact_solution[ilev], 0, 1, 1, 0); + amrex::MultiFab::Copy(plotmf[ilev], rhs[ilev], 0, 2, 1, 0); + amrex::MultiFab::Copy(plotmf[ilev], acoef[ilev], 0, 3, 1, 0); + amrex::MultiFab::Copy(plotmf[ilev], bcoef[ilev], 0, 4, 1, 0); + amrex::MultiFab::Copy(plotmf[ilev], vfrc, 0, 5, 1, 0); + } + auto const plot_file_name = amrpp.plot_file_name_; + amrex::WriteMultiLevelPlotfile( + plot_file_name, nlevels, amrex::GetVecOfConstPtrs(plotmf), + {"phi", "exact", "rhs", "acoef", "bcoef", "vfrac"}, geom, 0.0, + amrex::Vector(nlevels, 0), + amrex::Vector(nlevels, amrex::IntVect{ref_ratio})); + + // for amrvis + /* + amrex::writeFabs(solution, "solution"); + amrex::writeFabs(bcoef, "bcoef"); + amrex::writeFabs(robin_a, "robin_a"); + amrex::writeFabs(robin_b, "robin_b"); + amrex::writeFabs(robin_f, "robin_f"); + */ + } + + amrex::Finalize(); + + if (eps < 1e-2) { + return (0); + } else { + return (-1); + } +} diff --git a/Testing/Exec/Radiation/tstPOneSingle.cpp b/Testing/Exec/Radiation/tstPOneSingle.cpp new file mode 100644 index 000000000..bc5b9d1b1 --- /dev/null +++ b/Testing/Exec/Radiation/tstPOneSingle.cpp @@ -0,0 +1,257 @@ +#include +#include +#include +#include + +AMREX_GPU_DEVICE AMREX_FORCE_INLINE void +actual_init_coefs( + int i, + int j, + int k, + amrex::Array4 const& rhs, + amrex::Array4 const& sol, + amrex::Array4 const& alpha, + amrex::Array4 const& beta, + amrex::Array4 const& robin_a, + amrex::Array4 const& robin_b, + amrex::Array4 const& robin_f, + amrex::GpuArray const& prob_lo, + amrex::GpuArray const& prob_hi, + amrex::GpuArray const& dx, + amrex::Dim3 const& dlo, + amrex::Dim3 const& dhi, + amrex::Box const& vbx) +{ + amrex::Real const L = 2.0; + amrex::Real const n = 3.0; + amrex::Real const npioverL = n * M_PI / L; + + amrex::Real x = prob_lo[0] + dx[0] * (i + 0.5); + amrex::Real y = prob_lo[1] + dx[1] * (j + 0.5); + amrex::Real z = prob_lo[2] + dx[2] * (k + 0.5); + + beta(i, j, k) = 1.0; + + x = amrex::min(amrex::max(x, prob_lo[0]), prob_hi[0]); + y = amrex::min(amrex::max(y, prob_lo[1]), prob_hi[1]); + z = amrex::min(amrex::max(z, prob_lo[2]), prob_hi[2]); + + amrex::Real sincossin = + std::sin(npioverL * x) * std::cos(npioverL * y) * std::sin(npioverL * z); + + amrex::Real coscossin = + std::cos(npioverL * x) * std::cos(npioverL * y) * std::sin(npioverL * z); + + sol(i, j, k) = sincossin; + + if (vbx.contains(i, j, k)) { + rhs(i, j, k) = (1.0 + npioverL * npioverL) * beta(i, j, k) * sincossin; + alpha(i, j, k) = 1.0; + } + + // Robin BC + if (j >= dlo.y && j <= dhi.y && k >= dlo.z && k <= dhi.z) { + if (i > dhi.x || i < dlo.x) { + robin_a(i, j, k) = -1.0; + robin_b(i, j, k) = -2.0 / 3.0; + robin_f(i, j, k) = robin_a(i, j, k) * sol(i, j, k) + + robin_b(i, j, k) * npioverL * coscossin; + } + } else if (i >= dlo.x && i <= dhi.x && k >= dlo.z && k <= dhi.z) { + if (j > dhi.y || j < dlo.y) { + robin_a(i, j, k) = -1.0; + robin_b(i, j, k) = -2.0 / 3.0; + + amrex::Real msinsinsin = -std::sin(npioverL * x) * + std::sin(npioverL * y) * std::sin(npioverL * z); + robin_f(i, j, k) = robin_a(i, j, k) * sol(i, j, k) + + robin_b(i, j, k) * npioverL * msinsinsin; + } + } else if (i >= dlo.x && i <= dhi.x && j >= dlo.y && j <= dhi.y) { + if (k > dhi.z || k < dlo.z) { + robin_a(i, j, k) = -1.0; + robin_b(i, j, k) = -2.0 / 3.0; + + amrex::Real sincoscos = -std::sin(npioverL * x) * std::cos(npioverL * y) * + std::cos(npioverL * z); + + robin_f(i, j, k) = robin_a(i, j, k) * sol(i, j, k) + + robin_b(i, j, k) * npioverL * sincoscos; + } + } +} + +void +initProbABecLaplacian( + amrex::Geometry& geom, + amrex::MultiFab& solution, + amrex::MultiFab& rhs, + amrex::MultiFab& exact_solution, + amrex::MultiFab& acoef, + amrex::MultiFab& bcoef, + amrex::MultiFab& robin_a, + amrex::MultiFab& robin_b, + amrex::MultiFab& robin_f) +{ + auto const prob_lo = geom.ProbLoArray(); + auto const prob_hi = geom.ProbHiArray(); + auto const dx = geom.CellSizeArray(); + auto const dlo = amrex::lbound(geom.Domain()); + auto const dhi = amrex::ubound(geom.Domain()); + for (amrex::MFIter mfi(rhs); mfi.isValid(); ++mfi) { + amrex::Box const& bx = mfi.validbox(); + amrex::Box const& gbx = amrex::grow(bx, 1); + auto const& rhsfab = rhs.array(mfi); + auto const& solfab = solution.array(mfi); + auto const& acfab = acoef.array(mfi); + auto const& bcfab = bcoef.array(mfi); + auto const& rafab = robin_a.array(mfi); + auto const& rbfab = robin_b.array(mfi); + auto const& rffab = robin_f.array(mfi); + amrex::ParallelFor(gbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { + actual_init_coefs( + i, j, k, rhsfab, solfab, acfab, bcfab, rafab, rbfab, rffab, prob_lo, + prob_hi, dx, dlo, dhi, bx); + }); + } + + amrex::MultiFab::Copy(exact_solution, solution, 0, 0, 1, 0); + solution.setVal(0.0, 0, 1, amrex::IntVect(0)); +} + +void +initMeshandData( + PeleRad::AMRParam const& amrpp, + amrex::Geometry& geom, + amrex::BoxArray& grids, + amrex::DistributionMapping& dmap, + amrex::MultiFab& solution, + amrex::MultiFab& rhs, + amrex::MultiFab& exact_solution, + amrex::MultiFab& acoef, + amrex::MultiFab& bcoef, + amrex::MultiFab& robin_a, + amrex::MultiFab& robin_b, + amrex::MultiFab& robin_f) +{ + int const n_cell = amrpp.n_cell_; + int const max_grid_size = amrpp.max_grid_size_; + + amrex::RealBox rb( + {AMREX_D_DECL(-1.0, -1.0, -1.0)}, {AMREX_D_DECL(1.0, 1.0, 1.0)}); + amrex::Array is_periodic{AMREX_D_DECL(0, 0, 0)}; + amrex::Geometry::Setup(&rb, 0, is_periodic.data()); + amrex::Box domain0( + amrex::IntVect{AMREX_D_DECL(0, 0, 0)}, + amrex::IntVect{AMREX_D_DECL(n_cell - 1, n_cell - 1, n_cell - 1)}); + geom.define(domain0); + + grids.define(domain0); + grids.maxSize(max_grid_size); + + amrex::IntVect ng = amrex::IntVect{1}; + + dmap.define(grids); + solution.define(grids, dmap, 1, ng); + rhs.define(grids, dmap, 1, 0); + exact_solution.define(grids, dmap, 1, ng); + acoef.define(grids, dmap, 1, 0); + bcoef.define(grids, dmap, 1, ng); + robin_a.define(grids, dmap, 1, ng); + robin_b.define(grids, dmap, 1, ng); + robin_f.define(grids, dmap, 1, ng); + + initProbABecLaplacian( + geom, solution, rhs, exact_solution, acoef, bcoef, robin_a, robin_b, + robin_f); +} + +amrex::Real +check_norm(amrex::MultiFab const& phi, amrex::MultiFab const& exact) +{ + amrex::MultiFab mf(phi.boxArray(), phi.DistributionMap(), 1, 0); + amrex::MultiFab::Copy(mf, phi, 0, 0, 1, 0); + + amrex::MultiFab::Subtract(mf, exact, 0, 0, 1, 0); + + amrex::Real L0norm = mf.norm0(); + amrex::Real L1norm = mf.norm1(); + std::cout << " L0 norm:" << L0norm << ", L1 norm:" << L1norm << std::endl; + + return L1norm; +} + +int +main(int argc, char* argv[]) +{ + amrex::Initialize(argc, argv); + amrex::ParmParse pp; + PeleRad::AMRParam amrpp(pp); + PeleRad::MLMGParam mlmgpp(pp); + + bool const write = false; + int const n_cell = amrpp.n_cell_; + + amrex::Geometry geom; + amrex::BoxArray grids; + amrex::DistributionMapping dmap; + + amrex::MultiFab solution; + amrex::MultiFab rhs; + amrex::MultiFab exact_solution; + + amrex::MultiFab acoef; + amrex::MultiFab bcoef; + amrex::MultiFab robin_a; + amrex::MultiFab robin_b; + amrex::MultiFab robin_f; + + // std::cout << "initialize data ... \n"; + initMeshandData( + amrpp, geom, grids, dmap, solution, rhs, exact_solution, acoef, bcoef, + robin_a, robin_b, robin_f); + + // std::cout << "construct the PDE ... \n"; + PeleRad::POneSingle rte( + mlmgpp, geom, grids, dmap, solution, rhs, acoef, bcoef, robin_a, robin_b, + robin_f); + + // std::cout << "solve the PDE ... \n"; + rte.solve(); + + auto eps = check_norm(solution, exact_solution); + eps /= static_cast(n_cell * n_cell * n_cell); + std::cout << "n_cell=" << n_cell << ", normalized L1 norm:" << eps + << std::endl; + + // plot results + if (write) { + std::cout << "write the results ... \n"; + amrex::MultiFab plotmf(grids, dmap, 4, 0); + amrex::MultiFab::Copy(plotmf, solution, 0, 0, 1, 0); + amrex::MultiFab::Copy(plotmf, rhs, 0, 1, 1, 0); + amrex::MultiFab::Copy(plotmf, exact_solution, 0, 2, 1, 0); + amrex::MultiFab::Copy(plotmf, solution, 0, 3, 1, 0); + amrex::MultiFab::Subtract(plotmf, plotmf, 2, 3, 1, 0); + + auto const plot_file_name = amrpp.plot_file_name_; + amrex::WriteSingleLevelPlotfile( + plot_file_name, plotmf, {"phi", "rhs", "exact", "error"}, geom, 0.0, 0); + + // for amrvis + /* + amrex::writeFabs(solution, "solution"); + amrex::writeFabs(bcoef, "bcoef"); + amrex::writeFabs(robin_a, "robin_a"); + amrex::writeFabs(robin_b, "robin_b"); + amrex::writeFabs(robin_f, "robin_f"); + */ + } + + amrex::Finalize(); + if (eps < 1e-3) { + return (0); + } else { + return (-1); + } +} diff --git a/Testing/Exec/Radiation/tstPOneSingleAF.cpp b/Testing/Exec/Radiation/tstPOneSingleAF.cpp new file mode 100644 index 000000000..a91ff381e --- /dev/null +++ b/Testing/Exec/Radiation/tstPOneSingleAF.cpp @@ -0,0 +1,365 @@ +#include +#include +#include + +#include +#include +AMREX_GPU_DEVICE AMREX_FORCE_INLINE void +initGasField( + int i, + int j, + int k, + amrex::Array4 const& y_co2, + amrex::Array4 const& y_h2o, + amrex::Array4 const& y_co, + amrex::Array4 const& fv_soot, + amrex::Array4 const& temp, + amrex::Array4 const& pressure, + amrex::GpuArray const& dx, + amrex::GpuArray const& plo, + amrex::GpuArray const& phi) +{ + amrex::Real coef = 100; + amrex::Real xc = (phi[0] + plo[0]) * 0.5; + amrex::Real yc = (phi[1] + plo[1]) * 0.5; + + amrex::Real x = plo[0] + (i + 0.5) * dx[0]; + amrex::Real y = plo[1] + (j + 0.5) * dx[1]; + amrex::Real z = plo[2] + (k + 0.5) * dx[2]; + + amrex::Real r = std::sqrt((x - xc) * (x - xc) + (y - yc) * (y - yc)); + + z /= coef; + r /= coef; + + amrex::Real expr = std::exp( + -(4.0 * r / (0.05 + 0.1 * 4.0 * z)) * (4.0 * r / (0.05 + 0.1 * 4.0 * z))); + amrex::Real expTz = std::exp( + -((4.0 * z - 1.3) / (0.7 + 0.5 * 4.0 * z)) * + ((4.0 * z - 1.3) / (0.7 + 0.5 * 4.0 * z))); + + temp(i, j, k) = 300.0 + 1700.0 * expr * expTz; + + pressure(i, j, k) = 1.0e5; // Pa + + amrex::Real expSoot = + std::exp(-((4.0 * z - 1.0) / 0.7) * ((4.0 * z - 1.0) / 0.7)); + + fv_soot(i, j, k) = 1e-6 * expr * expSoot; + + amrex::Real expCO2z = std::exp( + -((4.0 * z - 1.1) / (0.6 + 0.5 * 4.0 * z)) * + ((4.0 * z - 1.1) / (0.6 + 0.5 * 4.0 * z))); + + y_co2(i, j, k) = 0.1 * expr * expCO2z; + + amrex::Real expH2Oz = std::exp( + -((4.0 * z - 1.0) / (0.7 + 0.5 * 4.0 * z)) * + ((4.0 * z - 1.0) / (0.7 + 0.5 * 4.0 * z))); + + y_h2o(i, j, k) = 0.2 * expr * expH2Oz; + + amrex::Real expCOz = + std::exp(-((4.0 * z - 1.0) / 0.7) * ((4.0 * z - 1.0) / 0.7)); + + y_co(i, j, k) = 0.09 * expr * expCOz; +} + +AMREX_GPU_DEVICE AMREX_FORCE_INLINE void +actual_init_coefs( + int i, + int j, + int k, + amrex::Array4 const& rhs, + amrex::Array4 const& alpha, + amrex::Array4 const& beta, + amrex::Dim3 const& dlo, + amrex::Dim3 const& dhi, + amrex::Box const& vbx, + amrex::Array4 const& absc, + amrex::Array4 const& T) +{ + beta(i, j, k) = 1.0; + if (vbx.contains(i, j, k)) { + amrex::Real ka = std::max(0.001, absc(i, j, k) * 100); + beta(i, j, k) = 1.0 / ka; + rhs(i, j, k) = 4.0 * ka * 5.67e-8 * std::pow(T(i, j, k), 4.0); + alpha(i, j, k) = ka; + } +} + +AMREX_GPU_DEVICE AMREX_FORCE_INLINE void +actual_init_bc_coefs( + int i, + int j, + int k, + amrex::Array4 const& alpha, + amrex::Array4 const& beta, + amrex::Array4 const& robin_a, + amrex::Array4 const& robin_b, + amrex::Array4 const& robin_f, + amrex::Dim3 const& dlo, + amrex::Dim3 const& dhi, + amrex::Array4 const& T) +{ + // Robin BC + bool robin_cell = false; + + if (j >= dlo.y && j <= dhi.y && k >= dlo.z && k <= dhi.z) { + if (i > dhi.x || i < dlo.x) { + robin_cell = true; + } + } else if (i >= dlo.x && i <= dhi.x && k >= dlo.z && k <= dhi.z) { + if (j > dhi.y || j < dlo.y) { + robin_cell = true; + } + } + + if (robin_cell) { + robin_a(i, j, k) = -1.0 / beta(i, j, k); + robin_b(i, j, k) = -2.0 / 3.0; + robin_f(i, j, k) = 0.0; + } +} + +void +initProbABecLaplacian( + amrex::Geometry& geom, + amrex::MultiFab& solution, + amrex::MultiFab& rhs, + amrex::MultiFab& acoef, + amrex::MultiFab& bcoef, + amrex::MultiFab& robin_a, + amrex::MultiFab& robin_b, + amrex::MultiFab& robin_f, + amrex::MultiFab& y_co2, + amrex::MultiFab& y_h2o, + amrex::MultiFab& y_co, + amrex::MultiFab& soot_fv_rad, + amrex::MultiFab& temperature, + amrex::MultiFab& pressure, + amrex::MultiFab& absc) +{ + using PeleRad::PlanckMean; + using PeleRad::RadProp::getRadPropGas; + using PeleRad::RadProp::getRadPropSoot; + + std::string data_path("kpDB/"); + PlanckMean radprop(data_path); + auto const& kpco2 = radprop.kpco2(); + auto const& kph2o = radprop.kph2o(); + auto const& kpco = radprop.kpco(); + auto const& kpsoot = radprop.kpsoot(); + + auto const prob_lo = geom.ProbLoArray(); + auto const prob_hi = geom.ProbHiArray(); + auto const dx = geom.CellSizeArray(); + auto const dlo = amrex::lbound(geom.Domain()); + auto const dhi = amrex::ubound(geom.Domain()); + + for (amrex::MFIter mfi(rhs); mfi.isValid(); ++mfi) { + amrex::Box const& bx = mfi.validbox(); + amrex::Box const& gbx = amrex::grow(bx, 1); + auto const& Yco2 = y_co2.array(mfi); + auto const& Yh2o = y_h2o.array(mfi); + auto const& Yco = y_co.array(mfi); + auto const& fv = soot_fv_rad.array(mfi); + auto const& T = temperature.array(mfi); + auto const& P = pressure.array(mfi); + auto const& kappa = absc.array(mfi); + + amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { + initGasField(i, j, k, Yco2, Yh2o, Yco, fv, T, P, dx, prob_lo, prob_hi); + getRadPropGas(i, j, k, Yco2, Yh2o, Yco, T, P, kappa, kpco2, kph2o, kpco); + }); + + // if soot exists + // cudaDeviceSynchronize(); + amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { + getRadPropSoot(i, j, k, fv, T, kappa, kpsoot); + }); + + auto const& rhsfab = rhs.array(mfi); + auto const& acfab = acoef.array(mfi); + auto const& bcfab = bcoef.array(mfi); + auto const& rafab = robin_a.array(mfi); + auto const& rbfab = robin_b.array(mfi); + auto const& rffab = robin_f.array(mfi); + + amrex::ParallelFor(gbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { + actual_init_coefs(i, j, k, rhsfab, acfab, bcfab, dlo, dhi, bx, kappa, T); + }); + + // cudaDeviceSynchronize(); + + amrex::ParallelFor(gbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { + actual_init_bc_coefs( + i, j, k, acfab, bcfab, rafab, rbfab, rffab, dlo, dhi, T); + }); + } + + solution.setVal(0.0, 0, 1, amrex::IntVect(0)); +} + +void +initMeshandData( + PeleRad::AMRParam const& amrpp, + amrex::Geometry& geom, + amrex::BoxArray& grids, + amrex::DistributionMapping& dmap, + amrex::MultiFab& solution, + amrex::MultiFab& rhs, + amrex::MultiFab& rad_src, + amrex::MultiFab& acoef, + amrex::MultiFab& bcoef, + amrex::MultiFab& robin_a, + amrex::MultiFab& robin_b, + amrex::MultiFab& robin_f, + amrex::MultiFab& y_co2, + amrex::MultiFab& y_h2o, + amrex::MultiFab& y_co, + amrex::MultiFab& soot_fv_rad, + amrex::MultiFab& temperature, + amrex::MultiFab& pressure, + amrex::MultiFab& absc) +{ + int const n_cell = amrpp.n_cell_; + int const max_grid_size = amrpp.max_grid_size_; + + amrex::RealBox rb( + {AMREX_D_DECL(0.0, 0.0, 0.0)}, {AMREX_D_DECL(12.5, 12.5, 75)}); + amrex::Array is_periodic{AMREX_D_DECL(0, 0, 0)}; + amrex::Geometry::Setup(&rb, 0, is_periodic.data()); + + std::vector npts{n_cell, n_cell, 6 * n_cell}; + amrex::Box domain0( + amrex::IntVect{AMREX_D_DECL(0, 0, 0)}, + amrex::IntVect{AMREX_D_DECL(npts[0] - 1, npts[1] - 1, npts[2] - 1)}); + geom.define(domain0); + + grids.define(domain0); + grids.maxSize(max_grid_size); + + amrex::IntVect ng = amrex::IntVect{1}; + + dmap.define(grids); + solution.define(grids, dmap, 1, ng); + rhs.define(grids, dmap, 1, 0); + rad_src.define(grids, dmap, 1, 0); + acoef.define(grids, dmap, 1, 0); + bcoef.define(grids, dmap, 1, ng); + robin_a.define(grids, dmap, 1, ng); + robin_b.define(grids, dmap, 1, ng); + robin_f.define(grids, dmap, 1, ng); + + y_co2.define(grids, dmap, 1, 0); + y_h2o.define(grids, dmap, 1, 0); + y_co.define(grids, dmap, 1, 0); + soot_fv_rad.define(grids, dmap, 1, 0); + temperature.define(grids, dmap, 1, 0); + pressure.define(grids, dmap, 1, 0); + absc.define(grids, dmap, 1, 0); + + initProbABecLaplacian( + geom, solution, rhs, acoef, bcoef, robin_a, robin_b, robin_f, y_co2, y_h2o, + y_co, soot_fv_rad, temperature, pressure, absc); + + bcoef.FillBoundary(); +} + +int +main(int argc, char* argv[]) +{ + amrex::Initialize(argc, argv); + amrex::ParmParse pp; + PeleRad::AMRParam amrpp(pp); + PeleRad::MLMGParam mlmgpp(pp); + + bool const write = false; + + amrex::Geometry geom; + amrex::BoxArray grids; + amrex::DistributionMapping dmap; + + amrex::MultiFab solution; + amrex::MultiFab rhs; + amrex::MultiFab rad_src; + + amrex::MultiFab acoef; + amrex::MultiFab bcoef; + amrex::MultiFab robin_a; + amrex::MultiFab robin_b; + amrex::MultiFab robin_f; + + amrex::MultiFab y_co2; + amrex::MultiFab y_h2o; + amrex::MultiFab y_co; + amrex::MultiFab soot_fv_rad; + amrex::MultiFab temperature; + amrex::MultiFab pressure; + amrex::MultiFab absc; + + // std::cout << "initialize data ... \n"; + initMeshandData( + amrpp, geom, grids, dmap, solution, rhs, rad_src, acoef, bcoef, robin_a, + robin_b, robin_f, y_co2, y_h2o, y_co, soot_fv_rad, temperature, pressure, + absc); + + // std::cout << "construct the PDE ... \n"; + + PeleRad::POneSingle rte( + mlmgpp, geom, grids, dmap, solution, rhs, acoef, bcoef, robin_a, robin_b, + robin_f); + // std::cout << "solve the PDE ... \n"; + rte.solve(); + + // calculate radiative heat source + for (amrex::MFIter mfi(rad_src); mfi.isValid(); ++mfi) { + amrex::Box const& bx = mfi.validbox(); + + auto const& rhsfab = rhs.array(mfi); + auto const& solfab = solution.array(mfi); + auto const& acfab = acoef.array(mfi); + + auto radfab = rad_src.array(mfi); + + amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { + radfab(i, j, k) = acfab(i, j, k) * solfab(i, j, k) - rhsfab(i, j, k); + }); + } + + // plot results + if (write) { + // std::cout << "write the results ... \n"; + amrex::MultiFab plotmf(grids, dmap, 10, 0); + amrex::MultiFab::Copy(plotmf, solution, 0, 0, 1, 0); + amrex::MultiFab::Copy(plotmf, rhs, 0, 1, 1, 0); + amrex::MultiFab::Copy(plotmf, acoef, 0, 2, 1, 0); + amrex::MultiFab::Copy(plotmf, bcoef, 0, 3, 1, 0); + amrex::MultiFab::Copy(plotmf, rad_src, 0, 4, 1, 0); + + amrex::MultiFab::Copy(plotmf, y_co2, 0, 5, 1, 0); + amrex::MultiFab::Copy(plotmf, y_h2o, 0, 6, 1, 0); + amrex::MultiFab::Copy(plotmf, y_co, 0, 7, 1, 0); + amrex::MultiFab::Copy(plotmf, soot_fv_rad, 0, 8, 1, 0); + amrex::MultiFab::Copy(plotmf, temperature, 0, 9, 1, 0); + + auto const plot_file_name = amrpp.plot_file_name_; + amrex::WriteSingleLevelPlotfile( + plot_file_name, plotmf, + {"G", "Emis", "kappa", "bcoef", "RadSrc", "Y_co2", "Y_h2o", "Y_co", + "Soot_fv", "Temperature"}, + geom, 0.0, 0); + + // for amrvis + /* + amrex::writeFabs(solution, "solution"); + amrex::writeFabs(bcoef, "bcoef"); + amrex::writeFabs(robin_a, "robin_a"); + amrex::writeFabs(robin_b, "robin_b"); + amrex::writeFabs(robin_f, "robin_f"); + */ + } + amrex::Finalize(); + return (0); +} diff --git a/Testing/Exec/Radiation/tstPOneSingleEB.cpp b/Testing/Exec/Radiation/tstPOneSingleEB.cpp new file mode 100644 index 000000000..5ec82c0a4 --- /dev/null +++ b/Testing/Exec/Radiation/tstPOneSingleEB.cpp @@ -0,0 +1,273 @@ +#include +#include +#include +#include +#include + +#include +#include +#include +#include + +AMREX_GPU_DEVICE AMREX_FORCE_INLINE void +actual_init_coefs_eb( + int i, + int j, + int k, + amrex::Array4 const& phi, + amrex::Array4 const& phi_exact, + amrex::Array4 const& rhs, + amrex::Array4 const& acoef, + amrex::Array4 const& bx, + amrex::Array4 const& flag, + amrex::GpuArray const& dx, + amrex::Box const& vbx) +{ + amrex::Real const L = std::sqrt(2.0); + amrex::Real const n = 3.0; + amrex::Real npioverL = n * M_PI / L; + + amrex::Real pioverfour = M_PI / 4.0; + amrex::Real cospioverfour = std::cos(pioverfour); + amrex::Real sinpioverfour = std::sin(pioverfour); + + if (vbx.contains(i, j, k)) { + bx(i, j, k) = 1.0; + phi(i, j, k) = 0.0; + + if (flag(i, j, k).isCovered()) { + rhs(i, j, k) = 0.0; + phi_exact(i, j, k) = 0.0; + } else { + amrex::Real x = dx[0] * (i + 0.5) - 0.5; + amrex::Real y = dx[1] * (j + 0.5) - 0.5; + amrex::Real z = dx[2] * k; + + // rotation + amrex::Real xp = x * cospioverfour + y * sinpioverfour; + amrex::Real yp = -x * sinpioverfour + y * cospioverfour; + + amrex::Real sincossin = std::sin(npioverL * xp) * + std::cos(npioverL * yp) * std::sin(npioverL * z); + + rhs(i, j, k) = (1.0 + npioverL * npioverL) * bx(i, j, k) * sincossin; + acoef(i, j, k) = 1.0; + + phi_exact(i, j, k) = sincossin; + } + } +} + +void +initProbABecLaplacian( + amrex::Geometry& geom, + std::unique_ptr& factory, + amrex::MultiFab& solution, + amrex::MultiFab& rhs, + amrex::MultiFab& exact_solution, + amrex::MultiFab& acoef, + amrex::MultiFab& bcoef) +{ + + amrex::FabArray const& flags = + factory->getMultiEBCellFlagFab(); + // amrex::MultiCutFab const& bcent = factory->getBndryCent(); + // amrex::MultiCutFab const& cent = factory->getCentroid(); + + auto const dx = geom.CellSizeArray(); + amrex::Box const& domainbox = geom.Domain(); + + for (amrex::MFIter mfi(rhs); mfi.isValid(); ++mfi) { + amrex::Box const& bx = mfi.validbox(); + amrex::Box const& nbx = amrex::surroundingNodes(bx); + // amrex::Box const& gbx = amrex::grow(bx, 1); + amrex::Array4 const& phi_arr = solution.array(mfi); + amrex::Array4 const& phi_ex_arr = exact_solution.array(mfi); + amrex::Array4 const& rhs_arr = rhs.array(mfi); + + amrex::Array4 const& bx_arr = bcoef.array(mfi); + + auto fabtyp = flags[mfi].getType(bx); + if (fabtyp == amrex::FabType::covered) { + std::cout << " amrex::FabType::covered == fabtyp \n"; + } else if (fabtyp == amrex::FabType::regular) { + std::cout << " amrex::FabType::regular == fabtyp \n"; + } else { + // std::cout << " amrex::FabType else \n"; + amrex::Array4 const& acoef_arr = acoef.array(mfi); + amrex::Array4 const& flag_arr = + flags.const_array(mfi); + // amrex::Array4 const& cent_arr + // = cent.const_array(mfi); + // amrex::Array4 const& bcent_arr + // = bcent.const_array(mfi); + + amrex::ParallelFor( + nbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { + actual_init_coefs_eb( + i, j, k, phi_arr, phi_ex_arr, rhs_arr, acoef_arr, bx_arr, flag_arr, + dx, bx); + }); + } + } +} + +amrex::Real +check_norm( + amrex::MultiFab const& phi, + amrex::MultiFab const& exact, + amrex::MultiFab const& vfrc) +{ + amrex::MultiFab mf(phi.boxArray(), phi.DistributionMap(), 1, 0); + amrex::MultiFab::Copy(mf, phi, 0, 0, 1, 0); + + amrex::MultiFab::Subtract(mf, exact, 0, 0, 1, 0); + amrex::MultiFab::Multiply(mf, vfrc, 0, 0, 1, 0); + + amrex::Real L0norm = mf.norm0(); + amrex::Real L1norm = mf.norm1(); + std::cout << " L0 norm:" << L0norm << ", L1 norm:" << L1norm << std::endl; + + return L1norm; +} + +int +main(int argc, char* argv[]) +{ + amrex::Initialize(argc, argv); + amrex::ParmParse pp; + PeleRad::AMRParam amrpp(pp); + PeleRad::MLMGParam mlmgpp(pp); + + bool const write = false; + int const n_cell = amrpp.n_cell_; + + amrex::Geometry geom; + amrex::BoxArray grids; + amrex::DistributionMapping dmap; + + amrex::MultiFab solution; + amrex::MultiFab exact_solution; + amrex::MultiFab rhs; + amrex::MultiFab acoef; + amrex::MultiFab bcoef; + + amrex::MultiFab robin_a; + amrex::MultiFab robin_b; + amrex::MultiFab robin_f; + + std::unique_ptr factory; + + int const max_grid_size = amrpp.max_grid_size_; + + amrex::RealBox rb( + {AMREX_D_DECL(-1.0, -1.0, -1.0)}, {AMREX_D_DECL(1.0, 1.0, 1.0)}); + amrex::Array is_periodic{AMREX_D_DECL(1, 1, 1)}; + amrex::Geometry::Setup(&rb, 0, is_periodic.data()); + amrex::Box domain0( + amrex::IntVect{AMREX_D_DECL(0, 0, 0)}, + amrex::IntVect{AMREX_D_DECL(n_cell - 1, n_cell - 1, n_cell - 1)}); + // geom.define(domain0); + geom.define(domain0, rb, amrex::CoordSys::cartesian, is_periodic); + + grids.define(domain0); + grids.maxSize(max_grid_size); + + // rotated box + int const max_coarsening_level = mlmgpp.max_coarsening_level_; + amrex::Real const la = std::sqrt(2.0) / 2.0; + amrex::Real const shift = std::sqrt(2.0) / 3.0; + amrex::EB2::BoxIF box( + {AMREX_D_DECL(-la, -la, -1.0)}, {AMREX_D_DECL(la, la, -1 + 4.0 * shift)}, + true); + auto gshop = amrex::EB2::makeShop(amrex::EB2::translate( + amrex::EB2::rotate( + amrex::EB2::translate(box, {AMREX_D_DECL(-0.0, -0.0, -0.0)}), + std::atan(1.0) * 1.0, 2), + {AMREX_D_DECL(0.0, 0.0, 0.0)})); + amrex::EB2::Build(gshop, geom, 0, max_coarsening_level); + + amrex::IntVect ng = amrex::IntVect{1}; + + dmap.define(grids); + + amrex::EB2::IndexSpace const& eb_is = amrex::EB2::IndexSpace::top(); + amrex::EB2::Level const& eb_level = eb_is.getLevel(geom); + factory = std::make_unique( + eb_level, geom, grids, dmap, amrex::Vector{2, 2, 2}, + amrex::EBSupport::full); + + amrex::MultiFab const& vfrc = factory->getVolFrac(); + + solution.define(grids, dmap, 1, 1, amrex::MFInfo(), *factory); + exact_solution.define(grids, dmap, 1, 0, amrex::MFInfo(), *factory); + rhs.define(grids, dmap, 1, 0, amrex::MFInfo(), *factory); + acoef.define(grids, dmap, 1, 0, amrex::MFInfo(), *factory); + bcoef.define(grids, dmap, 1, 1, amrex::MFInfo(), *factory); + robin_a.define(grids, dmap, 1, 0, amrex::MFInfo(), *factory); + robin_b.define(grids, dmap, 1, 0, amrex::MFInfo(), *factory); + robin_f.define(grids, dmap, 1, 0, amrex::MFInfo(), *factory); + + solution.setVal(0.0); + rhs.setVal(0.0); + acoef.setVal(1.0); + bcoef.setVal(1.0); + robin_a.setVal(0.0); + robin_b.setVal(0.0); + robin_f.setVal(0.0); + + initProbABecLaplacian( + geom, factory, solution, rhs, exact_solution, acoef, bcoef); + + // std::cout << "initialize data ... \n"; + // initMeshandData(amrpp, mlmgpp, geom, grids, dmap, factory, solution, + // rhs, + // exact_solution, acoef, bcoef); + + std::cout << "construct the PDE ... \n"; + PeleRad::POneSingleEB rte( + mlmgpp, geom, grids, dmap, factory, solution, rhs, acoef, bcoef, robin_a, + robin_b, robin_f); + + // std::cout << "solve the PDE ... \n"; + rte.solve(); + + auto eps = check_norm(solution, exact_solution, vfrc); + eps /= static_cast(n_cell * n_cell * n_cell); + std::cout << "n_cell=" << n_cell << ", normalized L1 norm:" << eps + << std::endl; + + // plot results + if (write) { + std::cout << "write the results ... \n"; + amrex::MultiFab plotmf(grids, dmap, 6, 0); + amrex::MultiFab::Copy(plotmf, solution, 0, 0, 1, 0); + amrex::MultiFab::Copy(plotmf, exact_solution, 0, 1, 1, 0); + amrex::MultiFab::Copy(plotmf, rhs, 0, 2, 1, 0); + amrex::MultiFab::Copy(plotmf, acoef, 0, 3, 1, 0); + amrex::MultiFab::Copy(plotmf, bcoef, 0, 4, 1, 0); + amrex::MultiFab::Copy(plotmf, vfrc, 0, 5, 1, 0); + + auto const plot_file_name = amrpp.plot_file_name_; + amrex::WriteSingleLevelPlotfile( + plot_file_name, plotmf, + {"phi", "exact", "rhs", "acoef", "bcoefx", "vfrac"}, geom, 0.0, 0); + + // for amrvis + /* + amrex::writeFabs(solution, "solution"); + amrex::writeFabs(bcoef, "bcoef"); + amrex::writeFabs(robin_a, "robin_a"); + amrex::writeFabs(robin_b, "robin_b"); + amrex::writeFabs(robin_f, "robin_f"); + */ + } + + amrex::Finalize(); + + if (eps < 5e-3) { + return (0); + } else { + return (-1); + } +}