From 5d8b01c7487bb6c4b5a2ec36d0fad16769903d6b Mon Sep 17 00:00:00 2001 From: Xi Zhang Date: Sun, 8 Dec 2024 12:32:28 -0800 Subject: [PATCH 01/10] modified CK and add Toon solver --- cmake/examples/rtdry.cmake | 7 +- examples/2024-XZhang-cloud-rt/rtdry.cpp | 6 +- examples/2024-XZhang-cloud-rt/rtdry.yaml | 95 +++ src/harp/radiation_band.cpp | 298 +++++++++ src/harp/radiation_band.hpp | 246 +++++++ src/harp/rt_solver_disort.cpp | 425 ++++++++++++ src/harp/rt_solver_toon.cpp | 797 +++++++++++++++++++++++ src/harp/rt_solvers.hpp | 155 +++++ src/harp/set_spectral_properties.cpp | 143 ++++ src/opacity/absorber.hpp | 114 ++++ src/opacity/absorber_ck.cpp | 170 +++++ src/opacity/absorber_ck.hpp | 42 ++ src/opacity/absorber_factory.cpp | 147 +++++ 13 files changed, 2637 insertions(+), 8 deletions(-) create mode 100644 examples/2024-XZhang-cloud-rt/rtdry.yaml create mode 100644 src/harp/radiation_band.cpp create mode 100644 src/harp/radiation_band.hpp create mode 100644 src/harp/rt_solver_disort.cpp create mode 100644 src/harp/rt_solver_toon.cpp create mode 100644 src/harp/rt_solvers.hpp create mode 100644 src/harp/set_spectral_properties.cpp create mode 100644 src/opacity/absorber.hpp create mode 100644 src/opacity/absorber_ck.cpp create mode 100644 src/opacity/absorber_ck.hpp create mode 100644 src/opacity/absorber_factory.cpp diff --git a/cmake/examples/rtdry.cmake b/cmake/examples/rtdry.cmake index 3c039389..f29191fd 100644 --- a/cmake/examples/rtdry.cmake +++ b/cmake/examples/rtdry.cmake @@ -3,11 +3,8 @@ # athena variables set_if_empty(NUMBER_GHOST_CELLS 3) -set(CUBED_SPHERE ON) -set(NBLOCKS 864) -# set(COORDINATE_SYSTEM "spherical_polar") -set(NETCDF OFF) -set(PNETCDF ON) +set(NETCDF ON) +set(PNETCDF OFF) set(MPI ON) set(DISORT ON) set(PYTHON_BINDINGS ON) diff --git a/examples/2024-XZhang-cloud-rt/rtdry.cpp b/examples/2024-XZhang-cloud-rt/rtdry.cpp index 2b29b303..495d8af1 100644 --- a/examples/2024-XZhang-cloud-rt/rtdry.cpp +++ b/examples/2024-XZhang-cloud-rt/rtdry.cpp @@ -203,9 +203,9 @@ void MeshBlock::ProblemGenerator(ParameterInput *pin) { for (; i <= ie; ++i) { air.w[IVX] = 0.001 * (1. * rand() / RAND_MAX - 0.5); AirParcelHelper::distribute_to_conserved(this, k, j, i, air); - // pthermo->Extrapolate(&air, pcoord->dx1f(i), "isothermal", - // grav); - pthermo->Extrapolate(&air, pcoord->dx1f(i), "dry", grav, 1.e-3); + pthermo->Extrapolate(&air, pcoord->dx1f(i), "isothermal", + grav); + //pthermo->Extrapolate(&air, pcoord->dx1f(i), "dry", grav, 1.e-4); } } } diff --git a/examples/2024-XZhang-cloud-rt/rtdry.yaml b/examples/2024-XZhang-cloud-rt/rtdry.yaml new file mode 100644 index 00000000..7f5f5bbc --- /dev/null +++ b/examples/2024-XZhang-cloud-rt/rtdry.yaml @@ -0,0 +1,95 @@ +opacity-sources: + - name: H2-vis + class: FreedmanMean + parameters: {scale: 0.5, met: 1.7} + + - name: H2-ir + class: FreedmanMean + parameters: {scale: 1., met: 1.7} + + - name: premix + class: AbsorberCK + data: ../data/ExoFMS/output.nc + + - name: heliosmix + class: HeliosCK + data: ../data/ck_data_01242024/ck/PM_ck_HELIOSK_cond_11_nOPT_wcia.txt + parameters: {band_id: 1} + +bands: [vistck] +#bands: [ird, irt] + #bands: [visd, vist] + +ird: + units: cm-1 + grid-type: regular + wavenumber-range: [100., 10000.] + num-bins: 1 + opacity: [H2-ir] + rt-solver: Disort + parameters: {btemp: 0., ttemp: 0} #assuming Tint = 350. K + flags: [broad_band, thermal_emission] + +irt: + units: cm-1 + grid-type: regular + wavenumber-range: [100., 10000.] + num-bins: 1 + opacity: [H2-ir] + rt-solver: Toon + flags: [broad_band, thermal_emission] + +irtck: + units: cm-1 + wavenumber-range: [100., 10000.] + grid-type: cktable + opacity: [premix] + rt-solver: Toon + flags: [thermal_emission] + +visd: + units: cm-1 + grid-type: regular + wavenumber-range: [10000., 50000.] + num-bins: 1 + opacity: [H2-vis] + rt-solver: Disort + parameters: {temp0: 1000., dist_au: 1., umu0: 1., phi0: 0.} + flags: [broad_band, stellar_beam] + +vist: + units: cm-1 + grid-type: regular + wavenumber-range: [10000., 50000.] + num-bins: 1 + opacity: [H2-vis] + rt-solver: Toon + parameters: {temp0: 1000., dist_au: 1., umu0: 1., phi0: 0.} + flags: [broad_band, stellar_beam] + +vistck: + units: cm-1 + wavenumber-range: [100., 10000.] + grid-type: cktable + opacity: [premix] + rt-solver: Toon + parameters: {temp0: 1000., dist_au: 1., umu0: 1., phi0: 0.} + flags: [broad_band, stellar_beam] + +Disort-flags: + ibcnd: false + usrtau: false + usrang: false + onlyfl: true + spher: false + intensity_correction: true + old_intensity_correction: false + general_source: false + output_uum: false + quiet: true + print-input: false + print-fluxes: false + print-intensity: false + print-transmissivity: false + print-phase-function: true + diff --git a/src/harp/radiation_band.cpp b/src/harp/radiation_band.cpp new file mode 100644 index 00000000..f1eddd6d --- /dev/null +++ b/src/harp/radiation_band.cpp @@ -0,0 +1,298 @@ +// C/C++ +#include +#include +#include +#include + +// external +#include + +// athena +#include +#include +#include +#include + +// canoe +#include + +// climath +#include + +// application +#include +#include + +// utils +#include +#include +#include + +// opacity +#include + +// harp +#include "radiation.hpp" +#include "radiation_band.hpp" +#include "rt_solvers.hpp" + +RadiationBand::RadiationBand(std::string myname, YAML::Node const &rad, + bool load_opacity) + : NamedGroup(myname) { + Application::Logger app("harp"); + app->Log("Initialize RadiationBand " + myname); + + if (!rad[myname]) { + throw NotFoundError("RadiationBand", myname); + } + + auto my = rad[myname]; + + if (my["parameters"]) { + SetRealsFrom(my["parameters"]); + } + + pgrid_ = SpectralGridFactory::CreateFrom(my); + + wrange_ = RadiationHelper::parse_wave_range(my); + + if (my["outdir"]) { + if (!my["outdir"].IsSequence()) { + throw RuntimeError("RadiationBand", "outdir must be a sequence"); + } + + for (const auto &item : my["outdir"]) { + rayOutput_.push_back( + RadiationHelper::parse_radiation_direction(item.as())); + } + } + + // set absorbers + if (my["opacity"]) { + if (!my["opacity"].IsSequence()) { + throw RuntimeError("RadiationBand", "opacity must be a sequence"); + } + + auto names = my["opacity"].as>(); + absorbers_ = AbsorberFactory::CreateFrom(names, GetName(), rad); + + if (load_opacity) { + for (auto &ab : absorbers_) { + ab->LoadOpacity(RadiationBandsFactory::GetBandId(myname)); + // Correlated-k absorbers need to modify the spectral grid + ab->ModifySpectralGrid(pgrid_->spec); + } + } + } + + // set flags + if (my["flags"]) { + if (!my["flags"].IsSequence()) { + throw RuntimeError("RadiationBand", "flags must be a sequence"); + } + + auto flag_strs = my["flags"].as>(); + for (auto flag : flag_strs) { + SetFlag(RadiationHelper::parse_radiation_flags(flag)); + } + } + + // set rt solver + if (my["rt-solver"]) { + psolver_ = CreateRTSolverFrom(my["rt-solver"].as(), rad); + } else { + psolver_ = nullptr; + } + + char buf[80]; + snprintf(buf, sizeof(buf), "%.2f - %.2f", wrange_.first, wrange_.second); + app->Log("Spectral range", buf); + app->Log("Number of spectral bins", pgrid_->spec.size()); +} + +RadiationBand::~RadiationBand() { + Application::Logger app("harp"); + app->Log("Destroy RadiationBand " + GetName()); +} + +void RadiationBand::Resize(int nc1, int nc2, int nc3, int nstr, + MeshBlock const *pmb) { + // allocate memory for spectral properties + tem_.resize(nc1); + temf_.resize(nc1 + 1); + + tau_.NewAthenaArray(pgrid_->spec.size(), nc1); + tau_.ZeroClear(); + + ssa_.NewAthenaArray(pgrid_->spec.size(), nc1); + ssa_.ZeroClear(); + + pmom_.NewAthenaArray(pgrid_->spec.size(), nc1, nstr + 1); + pmom_.ZeroClear(); + + // spectral grids properties + toa_.NewAthenaArray(pgrid_->spec.size(), rayOutput_.size(), nc3, nc2); + toa_.ZeroClear(); + + flxup_.NewAthenaArray(pgrid_->spec.size(), nc3, nc2, nc1); + flxup_.ZeroClear(); + + flxdn_.NewAthenaArray(pgrid_->spec.size(), nc3, nc2, nc1); + flxdn_.ZeroClear(); + + // band properties + btau.NewAthenaArray(nc3, nc2, nc1); + bssa.NewAthenaArray(nc3, nc2, nc1); + bpmom.NewAthenaArray(nstr + 1, nc3, nc2, nc1); + + btoa.NewAthenaArray(rayOutput_.size(), nc3, nc2); + bflxup.NewAthenaArray(nc3, nc2, nc1 + 1); + bflxdn.NewAthenaArray(nc3, nc2, nc1 + 1); + + // exchange buffer + pexv = std::make_shared>(GetName()); + + int nlayers = GetNumLayers(); + int npmom = GetNumPhaseMoments(); + pexv->send_buffer[0].resize(temf_.size()); + pexv->send_buffer[1].resize(pgrid_->spec.size()*nlayers * (npmom + 3)); + + pexv->Regroup(pmb, X1DIR); + int nblocks = pexv->GetGroupSize(); + pexv->recv_buffer[0].resize(nblocks * pexv->send_buffer[0].size()); + pexv->recv_buffer[1].resize(nblocks * pexv->send_buffer[1].size()); + + if (psolver_ != nullptr) { + psolver_->Resize(nblocks * (nc1 - 2 * NGHOST), nstr); + } +} + +AbsorberPtr RadiationBand::GetAbsorberByName(std::string const &name) { + for (auto &absorber : absorbers_) { + if (absorber->GetName() == name) { + return absorber; + } + } + + throw NotFoundError("RadiationBand", "Absorber " + name); + + return nullptr; +} + +RadiationBand const *RadiationBand::CalBandFlux(MeshBlock const *pmb, int k, + int j) { + // reset flux of this column + for (int i = pmb->is; i <= pmb->ie + 1; ++i) { + bflxup(k, j, i) = 0.; + bflxdn(k, j, i) = 0.; + } + + psolver_->Prepare(pmb, k, j); + + psolver_->CalBandFlux(pmb, k, j); + + return this; +} + +RadiationBand const *RadiationBand::CalBandRadiance(MeshBlock const *pmb, int k, + int j) { + // reset radiance of this column + for (int n = 0; n < GetNumOutgoingRays(); ++n) { + btoa(n, k, j) = 0.; + } + + psolver_->Prepare(pmb, k, j); + psolver_->CalBandRadiance(pmb, k, j); + + return this; +} + +void RadiationBand::WriteAsciiHeader(OutputParameters const *pout) const { + if (!TestFlag(RadiationFlags::WriteBinRadiance)) return; + + char fname[80], number[6]; + snprintf(number, sizeof(number), "%05d", pout->file_number); + snprintf(fname, sizeof(fname), "%s.radiance.%s.txt", GetName().c_str(), + number); + FILE *pfile = fopen(fname, "w"); + + fprintf(pfile, "# Bin Radiances of Band %s: %.3g - %.3g\n", GetName().c_str(), + wrange_.first, wrange_.second); + fprintf(pfile, "# Ray output size: %lu\n", rayOutput_.size()); + + fprintf(pfile, "# Polar angles: "); + for (auto &ray : rayOutput_) { + fprintf(pfile, "%.3f", rad2deg(acos(ray.mu))); + } + fprintf(pfile, "\n"); + + fprintf(pfile, "# Azimuthal angles: "); + for (auto &ray : rayOutput_) { + fprintf(pfile, "%.3f", rad2deg(ray.phi)); + } + fprintf(pfile, "\n"); + + fprintf(pfile, "#%12s%12s", "wave1", "wave2"); + for (size_t j = 0; j < rayOutput_.size(); ++j) { + fprintf(pfile, "%12s%lu", "Radiance", j + 1); + } + fprintf(pfile, "%12s\n", "weight"); + + fclose(pfile); +} + +void RadiationBand::WriteAsciiData(OutputParameters const *pout) const { + if (!TestFlag(RadiationFlags::WriteBinRadiance)) return; + + char fname[80], number[6]; + snprintf(number, sizeof(number), "%05d", pout->file_number); + snprintf(fname, sizeof(fname), "%s.radiance.%s.txt", GetName().c_str(), + number); + FILE *pfile = fopen(fname, "w"); + + for (size_t i = 0; i < pgrid_->spec.size(); ++i) { + fprintf(pfile, "%13.3g%12.3g", pgrid_->spec[i].wav1, pgrid_->spec[i].wav2); + for (size_t j = 0; j < rayOutput_.size(); ++j) { + fprintf(pfile, "%12.3f", toa_(i, j)); + } + if (TestFlag(RadiationFlags::Normalize) && + (wrange_.first != wrange_.second)) { + fprintf(pfile, "%12.3g\n", + pgrid_->spec[i].wght / (wrange_.second - wrange_.first)); + } else { + fprintf(pfile, "%12.3g\n", pgrid_->spec[i].wght); + } + } + + fclose(pfile); +} + +std::string RadiationBand::ToString() const { + std::stringstream ss; + ss << "RadiationBand: " << GetName() << std::endl; + ss << "Absorbers: "; + for (auto &ab : absorbers_) { + ss << ab->GetName() << ", "; + } + ss << std::endl << "RT-Solver: " << psolver_->GetName(); + return ss.str(); +} + +std::shared_ptr RadiationBand::CreateRTSolverFrom( + std::string const &rt_name, YAML::Node const &rad) { + std::shared_ptr psolver; + + if (rt_name == "Lambert") { + psolver = std::make_shared(this, rad); +#ifdef RT_DISORT + } else if (rt_name == "Disort") { + psolver = std::make_shared(this, rad); + } else if (rt_name == "Toon") { + psolver = std::make_shared(this, rad); +#endif // RT_DISORT + } else { + throw NotFoundError("RadiationBand", rt_name); + } + + return psolver; +} diff --git a/src/harp/radiation_band.hpp b/src/harp/radiation_band.hpp new file mode 100644 index 00000000..76a34337 --- /dev/null +++ b/src/harp/radiation_band.hpp @@ -0,0 +1,246 @@ +#ifndef SRC_HARP_RADIATION_BAND_HPP_ +#define SRC_HARP_RADIATION_BAND_HPP_ + +// C/C++ headers +#include +#include +#include +#include + +// athena +#include // Real + +// canoe +#include +#include + +// exchanger +#include + +// harp +#include "spectral_grid.hpp" + +namespace YAML { +class Node; +} + +class OutputParameters; +class Absorber; +class SpectralGridBase; +struct Direction; + +class RadiationBand : public NamedGroup, + public FlagGroup, + public ParameterGroup, + public ASCIIOutputGroup, + public StringReprGroup { + public: // public access data + // implementation of RT Solver + class RTSolver; + class RTSolverLambert; + class RTSolverDisort; + class RTSolverToon; + + std::shared_ptr> pexv; + + //! band optical depth + AthenaArray btau; + + //! band single scattering albedo + AthenaArray bssa; + + //! band phase function moments + AthenaArray bpmom; + + //! band upward flux (shallow reference to fluxup) + AthenaArray bflxup; + + //! band downward flux (shallow reference to flxdn) + AthenaArray bflxdn; + + //! \brief band top-of-the-atmosphere radiance (shallow reference to radiance) + AthenaArray btoa; + + public: // constructor and destructor + RadiationBand(std::string name, YAML::Node const &rad, bool load_opacity); + virtual ~RadiationBand(); + + public: // member functions + //! \brief Allocate memory for radiation band + void Resize(int nc1, int nc2, int nc3, int nstr, MeshBlock const *pmb); + + //! \brief Create radiative transfer solver from YAML node + //! + //! \param[in] rad YAML node containing whole radiation configuration + std::shared_ptr CreateRTSolverFrom(std::string const &name, + YAML::Node const &rad); + + //! Get number of spectral grids + size_t GetNumSpecGrids() const { return pgrid_->spec.size(); } + + //! Get number of absorbers + size_t GetNumAbsorbers() const { return absorbers_.size(); } + + std::vector> const &Absorbers() const { + return absorbers_; + } + + //! Get number of phase function moments + size_t GetNumPhaseMoments() const { return pmom_.GetDim1() - 1; } + + //! Get number of phase function moments + size_t GetNumLayers() const { return tem_.size() - 2 * NGHOST; } + + //! Get range of wavenumbers + std::pair GetRange() const { return wrange_; } + + //! Get an individual absorber + std::shared_ptr GetAbsorber(int i) { return absorbers_[i]; } + + //! Get an individual absorber by name + std::shared_ptr GetAbsorberByName(std::string const &name); + + //! Get number of outgoing rays + size_t GetNumOutgoingRays() { return rayOutput_.size(); } + + //! Get cosine polar angle of the outgoing ray + Real GetCosinePolarAngle(int n) const { return rayOutput_[n].mu; } + + //! Get azimuthal angle of the outgoing ray + Real GetAzimuthalAngle(int n) const { return rayOutput_[n].phi; } + + public: // python bindings + Real const *_GetToaPtr() const { return toa_.data(); } + Real const *_GetTauPtr() const { return tau_.data(); } + + public: // inbound functions + //! \brief Set spectral properties for an air column + //! + //! \param[in] air mole fraction representation of air column + //! \param[in] pcoord coordinates + //! \param[in] gH0 grav * H0 + //! \param[in] k horizontal index + //! \param[in] j horizontal index + void SetSpectralProperties(AirColumn &air, Real const *x1f, Real gH0 = 0, + int k = 0, int j = 0); + + //! \brief Calculate band radiative fluxes + RadiationBand const *CalBandFlux(MeshBlock const *pmb, int k, int j); + + //! \brief Calculate band radiances + RadiationBand const *CalBandRadiance(MeshBlock const *pmb, int k, int j); + + //! \brief Set outgoing ray directions + //! + //! \param[in] rayOutput outgoing ray directions + void SetOutgoingRays(std::vector const &ray) { rayOutput_ = ray; } + + public: // ASCIIOutputGroup functions + void WriteAsciiHeader(OutputParameters const *) const override; + void WriteAsciiData(OutputParameters const *) const override; + + public: // StringReprGroup functions + std::string ToString() const override; + + public: // Exchanger functions + //! \brief Pack temperature at cell face into send buffer 0 + void packTemperature(); + + //! \brief Unpack temperature at cell face from receive buffer 0 + bool unpackTemperature(void *arg); + + //! \brief Pack data in all spectral grids into send buffer 1 + //! + //! \param[in] b spectral bin index + void packSpectralProperties(); + + //! \brief Unpack data from receive buffer 1 into spectral grid b + void unpackSpectralProperties(int b, void *arg); + + protected: + //! radiative transfer solver + std::shared_ptr psolver_; + + //! spectral grid + std::shared_ptr pgrid_; + + //! \brief identify wave- number (wavelength, frequency) range + std::pair wrange_; + + //! all absorbers + std::vector> absorbers_; + + //! outgoing rays + std::vector rayOutput_; + + //! \brief spectral bin optical depth + //! + //! This is a two-dimensional array. The first dimension is the size of the + //! spectral bin and the second dimension is the size of the vertical + //! dimension. + AthenaArray tau_; + + //! \brief spectral bin single scattering albedo + //! + //! This is a two-dimensional array. The first dimension is the size of the + //! spectral bin and the second dimension is the size of the vertical + //! dimension. + AthenaArray ssa_; + + //! \brief spectral bin top-of-the-atmosphere radiance + //! + //! This is a two-dimensional array. The first dimension is the size of the + //! spectral bin and the second dimension is the size of the vertical + //! dimension. + AthenaArray toa_; + + //! \brief spectral bin upward flux + //! + //! This is a two-dimensional array. The first dimension is the size of the + //! spectral bin and the second dimension is the size of the vertical + //! dimension. + AthenaArray flxup_; + + //! \brief spectral bin downward flux + //! + //! This is a two-dimensional array. The first dimension is the size of the + //! spectral bin and the second dimension is the size of the vertical + //! dimension. + AthenaArray flxdn_; + + //! \brief spectral bin phase function moments + //! + //! This is a three-dimensional array. The first dimension is the size of the + //! spectral bin. The second dimension is the size of the vertical dimension. + //! The third dimension is the size of the phase function moments (1 + npmom). + AthenaArray pmom_; + + //! temperature at cell center + std::vector tem_; + + //! temperature at cell boundary (face) + std::vector temf_; +}; + +using RadiationBandPtr = std::shared_ptr; +using RadiationBandContainer = std::vector; + +class RadiationBandsFactory { + public: + static RadiationBandContainer CreateFrom(ParameterInput *pin, + std::string key); + static RadiationBandContainer CreateFrom(std::string filename); + + static int GetBandId(std::string const &bname) { + if (band_id_.find(bname) == band_id_.end()) { + return -1; + } + return band_id_.at(bname); + } + + protected: + static std::map band_id_; + static int next_band_id_; +}; + +#endif // SRC_HARP_RADIATION_BAND_HPP_ diff --git a/src/harp/rt_solver_disort.cpp b/src/harp/rt_solver_disort.cpp new file mode 100644 index 00000000..6762d84e --- /dev/null +++ b/src/harp/rt_solver_disort.cpp @@ -0,0 +1,425 @@ +//! \file rt_solver_disort.cpp +//! \brief Call DISORT to perform radiative transfer calculation + +// C/C++ +#include +#include + +// external +#include + +// athena +#include +#include + +// application +#include +#include + +// climath +#include + +// canoe +#include +#include + +// astro +#include + +// exo3 +#include +#include + +// harp +#include "radiation.hpp" +#include "rt_solvers.hpp" + +#ifdef RT_DISORT + +std::map to_map_bool(YAML::Node const &node) { + std::map flags; + + for (auto it = node.begin(); it != node.end(); ++it) { + flags[it->first.as()] = it->second.as(); + } + + return flags; +} + +RadiationBand::RTSolverDisort::RTSolverDisort(RadiationBand *pmy_band, + YAML::Node const &rad) + : RTSolver(pmy_band, "Disort") { + if (rad["Disort-flags"]) { + // SetFlags is a member function of DisortWrapper + SetFlags(to_map_bool(rad["Disort-flags"])); + } + + if (pmy_band->HasPar("btemp")) { + ds_.bc.btemp = pmy_band->GetPar("btemp"); + } + + if (pmy_band->HasPar("albedo")) { + ds_.bc.albedo = pmy_band->GetPar("albedo"); + } else { + // hard code hack to fix the surface albedo to 0, emissivity to 1, since it + // won't pickup from the yaml file + ds_.bc.albedo = 0.; + } + + if (pmy_band->HasPar("ttemp")) { + ds_.bc.ttemp = pmy_band->GetPar("ttemp"); + } + + if (pmy_band->HasPar("temis")) { + ds_.bc.temis = pmy_band->GetPar("temis"); + } else { + ds_.bc.temis = 0.; + } + + if (pmy_band->HasPar("fluor")) { + ds_.bc.fluor = pmy_band->GetPar("fluor"); + } else { + ds_.bc.fluor = 0.; + } + + if (pmy_band->HasPar("fisot")) { + ds_.bc.fisot = pmy_band->GetPar("fisot"); + } else { + ds_.bc.fisot = 0.; + } + + Application::Logger app("harp"); + + // override disort planck flag + if (pmy_band->TestFlag(RadiationFlags::ThermalEmission)) { + ds_.flag.planck = true; + app->Log("Planck function is enabled for band " + pmy_band->GetName()); + } else { + ds_.flag.planck = false; + app->Log("Planck function is disabled for band " + pmy_band->GetName()); + } + + SetHeader("Disort solving band " + pmy_band_->GetName()); +} + +//! \todo update based on band outdir +void RadiationBand::RTSolverDisort::Resize(int nlyr, int nstr) { + RadiationBand::RTSolver::Resize(nlyr, nstr); + Unseal(); + + auto &rayout = pmy_band_->rayOutput_; + auto &&uphi_umu = RadiationHelper::get_direction_grids(rayout); + + SetAtmosphereDimension(nlyr, nstr, nstr); + + dir_dim_[0] = uphi_umu.second.size(); // umu + dir_dim_[1] = uphi_umu.first.size(); // uphi + dir_axis_.resize(dir_dim_[0] + dir_dim_[1]); + + SetIntensityDimension(std::max(dir_dim_[1], 1lu), 1, + std::max(dir_dim_[0], 1lu)); + Seal(); +} + +//! \note Counting Disort Index +//! Example, il = 0, iu = 2, ds_.nlyr = 6, partition in to 3 blocks +//! face id -> 0 - 1 - 2 - 3 - 4 - 5 - 6 +//! cell id -> | 0 | 1 | 2 | 3 | 4 | 5 | +//! disort id -> 6 - 5 - 4 - 3 - 2 - 1 - 0 +//! blocks -> --------- * * +//! -> r = 0 * * * +//! -> --------- * +//! -> r = 1 * * +//! -> --------- +//! -> r = 2 +//! block r = 0 gets, 6 - 5 - 4 +//! block r = 1 gets, 4 - 3 - 2 +//! block r = 2 gets, 2 - 1 - 0 + +void RadiationBand::RTSolverDisort::Prepare(MeshBlock const *pmb, int k, + int j) { + auto &wmin = pmy_band_->wrange_.first; + auto &wmax = pmy_band_->wrange_.second; + + Real dist_au = 1.; + Direction ray = pmb->pimpl->prad->GetRayInput(0); + auto planet = pmb->pimpl->planet; + + if (planet && pmy_band_->TestFlag(RadiationFlags::TimeDependent)) { + Real time = pmb->pmy_mesh->time; + Real lat, lon; + + CubedSphereUtility::get_latlon_on_sphere(&lat, &lon, pmb, k, j, pmb->is); + + ray = planet->ParentZenithAngle(time, lat, lon); + dist_au = planet->ParentDistanceInAu(time); + } else { // constant zenith angle + if (pmy_band_->HasPar("umu0")) { + ray.mu = pmy_band_->GetPar("umu0"); + } + + if (pmy_band_->HasPar("phi0")) { + ray.phi = pmy_band_->GetPar("phi0"); + } + + if (pmy_band_->HasPar("dist_au")) { + dist_au = pmy_band_->GetPar("dist_au"); + } + } + + if (ds_.flag.ibcnd != 0) { + throw ValueError("RTSolverDisort::CalRadtranFlux", "ibcnd", ds_.flag.ibcnd, + 0); + } + + // pack temperature + if (ds_.flag.planck) { + pmy_band_->packTemperature(); + } + + // pack spectral properties + pmy_band_->packSpectralProperties(); + + ds_.bc.umu0 = ray.mu > 1.E-3 ? ray.mu : 1.E-3; + ds_.bc.phi0 = ray.phi; + + ds_.wvnmlo = wmin; + ds_.wvnmhi = wmax; + + if (pmy_band_->TestFlag(RadiationFlags::BroadBand)) { + // stellar source function overrides fbeam + if (pmy_band_->HasPar("S0")) { + ds_.bc.fbeam = pmy_band_->GetPar("S0"); + } else if (pmy_band_->HasPar("temp0")) { + Real temp0 = pmy_band_->GetPar("temp0"); + ds_.bc.fbeam = Constants::stefanBoltzmann * pow(temp0, 4); + } else if (planet && planet->HasParentFlux()) { + ds_.bc.fbeam = planet->ParentInsolationFlux(wmin, wmax, 1.); + } else { + ds_.bc.fbeam = 0.; + } + ds_.bc.fbeam /= dist_au * dist_au; + } + + pmb->pcoord->Face1Area(k, j, pmb->is, pmb->ie + 1, farea_); + pmb->pcoord->CellVolume(k, j, pmb->is, pmb->ie, vol_); + + auto &&uphi_umu = RadiationHelper::get_direction_grids(pmy_band_->rayOutput_); + auto &uphi = uphi_umu.first; + auto &umu = uphi_umu.second; + + if (umu.size() > 0) { + SetUserCosinePolarAngle(umu); + + for (int i = 0; i < umu.size(); ++i) { + dir_axis_[i] = umu[i]; + } + } + + if (uphi.size() > 0) { + SetUserAzimuthalAngle(uphi); + + for (int i = 0; i < uphi.size(); ++i) { + dir_axis_[ds_.numu + i] = uphi[i]; + } + } + + // set the surface temperature from the user input variables + // FIXME(cli) + // AthenaArray &ts = pmb->ruser_meshblock_data[1]; + // ds_.bc.btemp = ts(j); +} + +void RadiationBand::RTSolverDisort::CalBandFlux(MeshBlock const *pmb, int k, + int j) { + Real dist_au = 1.; + auto planet = pmb->pimpl->planet; + + if (planet && pmy_band_->TestFlag(RadiationFlags::TimeDependent)) { + dist_au = planet->ParentDistanceInAu(pmb->pmy_mesh->time); + } else if (pmy_band_->HasPar("dist_au")) { + dist_au = pmy_band_->GetPar("dist_au"); + } + + // bflxup and bflxdn has been reset in RadiationBand::CalBandFlux + + // loop over spectral grids in the band + bool override_with_stellar_spectra = false; + if (!pmy_band_->TestFlag(RadiationFlags::BroadBand) && + !pmy_band_->HasPar("S0") && !pmy_band_->HasPar("temp0") && planet && + planet->HasParentFlux()) { + override_with_stellar_spectra = true; + } + + pmy_band_->pexv->GatherAll(pmb); + if (ds_.flag.planck) { + pmy_band_->unpackTemperature(&ds_); + } + + int b = 0; + for (auto &spec : pmy_band_->pgrid_->spec) { + if (override_with_stellar_spectra) { + // stellar source function + ds_.bc.fbeam = + planet->ParentInsolationFlux(spec.wav1, spec.wav2, dist_au); + // planck source function + ds_.wvnmlo = spec.wav1; + ds_.wvnmhi = spec.wav2; + } + + // transfer spectral grid data + pmy_band_->unpackSpectralProperties(b, &ds_); + + if (ds_.nmom_nstr == 2) { + int nlay = ds_.nlyr; + int ierror[TWOSTR_NERR]; + std::vector gg(nlay, 0.); + for (int i = 0; i < nlay; ++i) gg[i] = ds_.pmom[i * (ds_.nmom_nstr + 1) + 1]; + c_twostr(&ds_, &ds_out_, 0, gg.data(), ierror, 1000.); + } else { + c_disort(&ds_, &ds_out_); + } + + // add spectral bin flux + addDisortFlux(pmb->pcoord, b++, k, j, pmb->is, pmb->ie + 1); + } +} + +void RadiationBand::RTSolverDisort::addDisortFlux(Coordinates const *pcoord, + int b, int k, int j, int il, + int iu) { + auto &bflxup = pmy_band_->bflxup; + auto &bflxdn = pmy_band_->bflxdn; + + auto &flxup = pmy_band_->flxup_; + auto &flxdn = pmy_band_->flxdn_; + auto const &spec = pmy_band_->pgrid_->spec; + + int rank_in_column = pmy_band_->pexv->GetRankInGroup(); + + /// accumulate flux from spectral bins + for (int i = il; i <= iu; ++i) { + int m = ds_.nlyr - (rank_in_column * (iu - il) + i - il); + //! \bug does not work for spherical geometry, need to scale area using + //! farea(il)/farea(i) + // flux up + flxup(b, k, j, i) = ds_out_.rad[m].flup; + + //! \bug does not work for spherical geomtry, need to scale area using + //! farea(il)/farea(i) + // flux down + flxdn(b, k, j, i) = ds_out_.rad[m].rfldir + ds_out_.rad[m].rfldn; + + bflxup(k, j, i) += spec[b].wght * flxup(b, k, j, i); + bflxdn(k, j, i) += spec[b].wght * flxdn(b, k, j, i); + } + + //! \note Spherical correction by XIZ + //! xiz 2022 flux scaling so that the heating rate is the same as the + //! plane-parallel scheme volheating scaling: first calculate flux divergence + //! from DISORT using Plane-parallel in a cell then mulitpled by the cell + //! volume divided by dx1f then solve for F using F1*S1-F2*S2 = volheating + //! the top fluxes are the still the same as the plane-paralell values + Real volh; + Real bflxup_iu = bflxup(k, j, iu); + Real bflxdn_iu = bflxdn(k, j, iu); + + for (int i = iu - 1; i >= il; --i) { + // upward + volh = (bflxup_iu - bflxup(k, j, i)) / pcoord->dx1f(i) * vol_(i); + bflxup_iu = bflxup(k, j, i); + bflxup(k, j, i) = (bflxup(k, j, i + 1) * farea_(i + 1) - volh) / farea_(i); + + // downward + volh = (bflxdn_iu - bflxdn(k, j, i)) / pcoord->dx1f(i) * vol_(i); + bflxdn_iu = bflxdn(k, j, i); + bflxdn(k, j, i) = (bflxdn(k, j, i + 1) * farea_(i + 1) - volh) / farea_(i); + } + + /* + for (int i = iu; i >= il; --i) { + std::cout << "i: " << iu-i+1 <<" flxup: " << bflxup(k, j, i) << " flxdn: " << bflxdn(k, j, i) << " fluxdiff: " << bflxup(k, j, i) - bflxdn(k, j, i) << std::endl; + } +*/ + +} + +void RadiationBand::RTSolverDisort::CalBandRadiance(MeshBlock const *pmb, int k, + int j) { + if (ds_.flag.onlyfl) { + throw RuntimeError("RTSolverDisort::CalBandRadiance", + "Radiance calculation disabled"); + } + + if (ds_.ntau != 1) { + throw RuntimeError("RTSolverDisort::CalBandRadiance", + "Only toa radiance (ds.ntau = 1) is supported"); + } + + int nrays = ds_.nphi * ds_.numu; + + if (nrays < pmy_band_->GetNumOutgoingRays()) { + throw RuntimeError("RTSolverDisort::CalBandRadiance", + "Number of outgoing rays more than DISORT can host"); + } + + // toa has been reset in RadiationBand::CalBandRadiance + + Real dist_au = 1.; + auto planet = pmb->pimpl->planet; + + if (planet && pmy_band_->TestFlag(RadiationFlags::TimeDependent)) { + dist_au = pmb->pimpl->planet->ParentDistanceInAu(pmb->pmy_mesh->time); + } else if (pmy_band_->HasPar("dist_au")) { + dist_au = pmy_band_->GetPar("dist_au"); + } + + pmy_band_->pexv->GatherAll(pmb); + if (ds_.flag.planck) { + pmy_band_->unpackTemperature(&ds_); + } + + // loop over spectral grids in the band + int b = 0; + for (auto &spec : pmy_band_->pgrid_->spec) { + // override source function for non-broadband radiation + if (!(pmy_band_->TestFlag(RadiationFlags::BroadBand))) { + // stellar source function + if (planet) { + ds_.bc.fbeam = + planet->ParentInsolationFlux(spec.wav1, spec.wav2, dist_au); + } else { + ds_.bc.fbeam = 0.0; + } + // planck source function + ds_.wvnmlo = spec.wav1; + ds_.wvnmhi = spec.wav2; + } + + // transfer spectral grid data + pmy_band_->unpackSpectralProperties(b, &ds_); + + // run disort + c_disort(&ds_, &ds_out_); + + // add spectral bin radiance + addDisortRadiance(b++, k, j); + } +} + +void RadiationBand::RTSolverDisort::addDisortRadiance(int b, int k, int j) { + auto &toa = pmy_band_->toa_; + auto &btoa = pmy_band_->btoa; + auto &spec = pmy_band_->pgrid_->spec; + auto &rayout = pmy_band_->rayOutput_; + + for (int n = 0; n < pmy_band_->GetNumOutgoingRays(); ++n) { + Real val; + Real coor[2] = {rayout[n].mu, rayout[n].phi}; + interpn(&val, coor, ds_out_.uu, dir_axis_.data(), dir_dim_, 2, 1); + toa(b, n, k, j) = val; + btoa(n, k, j) += spec[b].wght * toa(b, n, k, j); + } +} + +#endif // RT_DISORT diff --git a/src/harp/rt_solver_toon.cpp b/src/harp/rt_solver_toon.cpp new file mode 100644 index 00000000..dd3c3235 --- /dev/null +++ b/src/harp/rt_solver_toon.cpp @@ -0,0 +1,797 @@ +// RT solvers based on Toon 1989 method by Xi Zhang +// Reference: Toon, O.B., 1989, JGR, 94, 16287-16301. + +#include +#include +#include +#include +#include + +// external +#include + +// athena +#include +#include + +// application +#include +#include + +// climath +#include + +// canoe +#include +#include + +// astro +#include + +// exo3 +#include +#include + +// harp +#include "radiation.hpp" +#include "rt_solvers.hpp" + +#ifdef RT_DISORT + +RadiationBand::RTSolverToon::RTSolverToon(RadiationBand *pmy_band, YAML::Node const &rad) + : RTSolver(pmy_band, "Toon") { + Application::Logger app("harp"); + app->Log("Toon solver initialized for band " + pmy_band_->GetName()); +} + + +//! \todo update based on band outdir +void RadiationBand::RTSolverToon::Resize(int nlyr, int nstr) { + RadiationBand::RTSolver::Resize(nlyr, nstr); + Unseal(); + SetAtmosphereDimension(nlyr, nstr, nstr); + Seal(); +} + +void RadiationBand::RTSolverToon::Prepare(MeshBlock const *pmb, int k, int j) { + auto &wmin = pmy_band_->wrange_.first; + auto &wmax = pmy_band_->wrange_.second; + + Real dist_au = 1.0; + Direction ray = pmb->pimpl->prad->GetRayInput(0); + auto planet = pmb->pimpl->planet; + + if (planet && pmy_band_->TestFlag(RadiationFlags::TimeDependent)) { + Real time = pmb->pmy_mesh->time; + Real lat, lon; + + CubedSphereUtility::get_latlon_on_sphere(&lat, &lon, pmb, k, j, pmb->is); + + ray = planet->ParentZenithAngle(time, lat, lon); + dist_au = planet->ParentDistanceInAu(time); + } else { + if (pmy_band_->HasPar("umu0")) { + ray.mu = pmy_band_->GetPar("umu0"); + } + + if (pmy_band_->HasPar("phi0")) { + ray.phi = pmy_band_->GetPar("phi0"); + } + + if (pmy_band_->HasPar("dist_au")) { + dist_au = pmy_band_->GetPar("dist_au"); + } + } + + // pack temperature + if (pmy_band_->TestFlag(RadiationFlags::ThermalEmission)) { + pmy_band_->packTemperature(); + } + + // pack spectral properties + pmy_band_->packSpectralProperties(); + ds_.bc.umu0 = ray.mu > 1.E-3 ? ray.mu : 1.E-3; + + if (pmy_band_->TestFlag(RadiationFlags::BroadBand)) { + // stellar source function overrides fbeam + if (pmy_band_->HasPar("S0")) { + ds_.bc.fbeam = pmy_band_->GetPar("S0"); + } else if (pmy_band_->HasPar("temp0")) { + Real temp0 = pmy_band_->GetPar("temp0"); + ds_.bc.fbeam = Constants::stefanBoltzmann * pow(temp0, 4); + } else if (planet && planet->HasParentFlux()) { + ds_.bc.fbeam = planet->ParentInsolationFlux(wmin, wmax, 1.); + } else { + ds_.bc.fbeam = 0.; + } + ds_.bc.fbeam /= dist_au * dist_au; + } + + pmb->pcoord->Face1Area(k, j, pmb->is, pmb->ie + 1, farea_); + pmb->pcoord->CellVolume(k, j, pmb->is, pmb->ie, vol_); +} + +void RadiationBand::RTSolverToon::CalBandFlux(MeshBlock const *pmb, int k, int j) { + Real dist_au = 1.0; + auto planet = pmb->pimpl->planet; + + if (planet && pmy_band_->TestFlag(RadiationFlags::TimeDependent)) { + dist_au = planet->ParentDistanceInAu(pmb->pmy_mesh->time); + } else if (pmy_band_->HasPar("dist_au")) { + dist_au = pmy_band_->GetPar("dist_au"); + } + + // loop over spectral grids in the band + bool override_with_stellar_spectra = false; + if (!pmy_band_->TestFlag(RadiationFlags::BroadBand) && + !pmy_band_->HasPar("S0") && !pmy_band_->HasPar("temp0") && planet && + planet->HasParentFlux()) { + override_with_stellar_spectra = true; + } + + pmy_band_->pexv->GatherAll(pmb); + if (pmy_band_->TestFlag(RadiationFlags::ThermalEmission)) { + pmy_band_->unpackTemperature(&ds_); + } + + int b = 0; + for (auto &spec : pmy_band_->pgrid_->spec) { + + if (override_with_stellar_spectra) { + // stellar source function + ds_.bc.fbeam = + planet->ParentInsolationFlux(spec.wav1, spec.wav2, dist_au); + } + + // Transfer spectral grid data + pmy_band_->unpackSpectralProperties(b, &ds_); + + // Determine flux for each spectral bin using Toon method + int nlay = ds_.nlyr; + Eigen::VectorXd w0(nlay), g(nlay); + Eigen::VectorXd tau_cum(nlay+1); + tau_cum.setZero(); + + for (int i = 0; i < nlay; ++i) { + tau_cum(i+1) = tau_cum(i) + ds_.dtauc[i]; + w0(i) = ds_.ssalb[i]; + g(i) = ds_.pmom[i * (ds_.nmom_nstr + 1) + 1]; + } + + Eigen::VectorXd flux_up(nlay), flux_down(nlay); + double mu0 = ds_.bc.umu0; + double surface_albedo = ds_.bc.albedo; + double Finc = ds_.bc.fbeam; + +// Call shortwave or longwave solver based on the type of radiation + if (pmy_band_->TestFlag(RadiationFlags::StellarBeam)) { + Eigen::VectorXd mu0_in(nlay+1); + mu0_in.setConstant(mu0); + toonShortwaveSolver(nlay, Finc, mu0_in, tau_cum, w0, g, surface_albedo, flux_up, flux_down); + + } else if (pmy_band_->TestFlag(RadiationFlags::ThermalEmission)) { + Eigen::VectorXd temp(nlay + 1); + Eigen::VectorXd be(nlay + 1); + for (int i = 0; i < nlay + 1; ++i) { + temp(i) = ds_.temper[i]; + be(i) = BB_integrate(ds_.temper[i], spec.wav1, spec.wav2); + } + double surface_emissivity = ds_.bc.temis; + toonLongwaveSolver(nlay, be, tau_cum, w0, g, surface_emissivity, flux_up, flux_down); + } + + // add spectral bin flux + addToonFlux(pmb->pcoord, b++, k, j, pmb->is, pmb->ie + 1, flux_up, flux_down); + } +} + +void RadiationBand::RTSolverToon::addToonFlux(Coordinates const *pcoord, int b, int k, int j, int il, int iu, + const Eigen::VectorXd &flux_up, const Eigen::VectorXd &flux_down) { + auto &bflxup = pmy_band_->bflxup; + auto &bflxdn = pmy_band_->bflxdn; + + auto &flxup = pmy_band_->flxup_; + auto &flxdn = pmy_band_->flxdn_; + auto const &spec = pmy_band_->pgrid_->spec; + + int rank_in_column = pmy_band_->pexv->GetRankInGroup(); + + // Accumulate flux from spectral bins + for (int i = il; i <= iu; ++i) { + int m = ds_.nlyr - (rank_in_column * (iu - il) + i - il); + // Flux up + flxup(b, k, j, i) = flux_up(m); + // Flux down + flxdn(b, k, j, i) = flux_down(m); + + bflxup(k, j, i) += spec[b].wght * flxup(b, k, j, i); + bflxdn(k, j, i) += spec[b].wght * flxdn(b, k, j, i); + } + + // Spherical correction + Real volh; + Real bflxup_iu = bflxup(k, j, iu); + Real bflxdn_iu = bflxdn(k, j, iu); + + for (int i = iu - 1; i >= il; --i) { + // Upward + volh = (bflxup_iu - bflxup(k, j, i)) / pcoord->dx1f(i) * vol_(i); + bflxup_iu = bflxup(k, j, i); + bflxup(k, j, i) = (bflxup(k, j, i + 1) * farea_(i + 1) - volh) / farea_(i); + + // Downward + volh = (bflxdn_iu - bflxdn(k, j, i)) / pcoord->dx1f(i) * vol_(i); + bflxdn_iu = bflxdn(k, j, i); + bflxdn(k, j, i) = (bflxdn(k, j, i + 1) * farea_(i + 1) - volh) / farea_(i); + } + + /* + for (int i = iu; i >= il; --i) { + std::cout << "i: " << iu-i+1 <<" flxup: " << bflxup(k, j, i) << " flxdn: " << bflxdn(k, j, i) << " fluxdiff: " << bflxup(k, j, i) - bflxdn(k, j, i) << std::endl; + } +*/ +} + + +// Toon 1989 shortwave solver, based on Elsie Lee's implementation, which was based on CHIMERA code by Mike Line. Reference: Toon, O.B., 1989, JGR, 94, 16287-16301. +void RadiationBand::RTSolverToon::toonShortwaveSolver(int nlay, double F0_in, + const Eigen::VectorXd& mu_in, + const Eigen::VectorXd& tau_in, + const Eigen::VectorXd& w_in, + const Eigen::VectorXd& g_in, + double w_surf_in, + Eigen::VectorXd& flx_up, + Eigen::VectorXd& flx_down) + { + int nlev = nlay + 1; + // Input validation + if(mu_in.size() != static_cast(nlev) || + tau_in.size() != static_cast(nlev) || + w_in.size() != static_cast(nlay) || + g_in.size() != static_cast(nlay)) + { + throw std::invalid_argument("Input vectors have incorrect sizes."); + } + + // Initialize output flux arrays + flx_down = Eigen::VectorXd::Zero(nlev); + flx_up = Eigen::VectorXd::Zero(nlev); + + // Constants + const double sqrt3 = std::sqrt(3.0); + const double sqrt3d2 = sqrt3 / 2.0; + const double bsurf = 0.0; + const double btop = 0.0; + + // Check if all albedos are effectively zero + bool all_w0_zero = (w_in.array() <= 1.0e-12).all(); + + if(all_w0_zero){ + // Direct beam only + double mu_top = mu_in(nlev - 1); + double mu_first = mu_in(0); + + if(mu_top > 0.0){ + if(std::abs(mu_top - mu_first) < 1e-12){ + // No zenith correction, use regular method + flx_down = F0_in * mu_top * (-tau_in.array() / mu_top).exp(); + } + else{ + // Zenith angle correction using cumulative transmission + Eigen::VectorXd cum_trans(nlev); + cum_trans(0) = tau_in(0) / mu_in(0); + for(int k = 1; k < nlev; ++k){ + cum_trans(k) = cum_trans(k - 1) + (tau_in(k) - tau_in(k - 1)) / mu_in(k); + } + flx_down = F0_in * mu_top * (-cum_trans.array()).exp(); + } + // Adjust the downward flux at the surface layer for surface albedo + flx_down(nlev - 1) *= (1.0 - w_surf_in); + } + // Upward flux remains zero + return; + } + + // Delta Eddington scaling + Eigen::VectorXd w0 = ((1.0 - g_in.array().square()) * w_in.array()) / (1.0 - w_in.array() * g_in.array().square()); + Eigen::VectorXd dtau = (1.0 - w_in.array() * g_in.array().square()) + .cwiseProduct( (tau_in.segment(1, nlay) - tau_in.head(nlay)).array() ) + .matrix(); + Eigen::VectorXd hg = g_in.array() / (1.0 + g_in.array()); + + // Initialize tau_total + Eigen::VectorXd tau_total(nlev); + tau_total(0) = 0.0; + for(int k = 0; k < nlay; ++k){ + tau_total(k + 1) = tau_total(k) + dtau(k); + } + + // Compute g1, g2, g3, g4 + Eigen::VectorXd g1 = sqrt3d2 * (2.0 - w0.array() * (1.0 + hg.array())); + Eigen::VectorXd g2 = (sqrt3d2 * w0.array()) * (1.0 - hg.array()); + // Prevent division by zero + for(int i = 0; i < nlay; ++i){ + if(std::abs(g2(i)) < 1.0e-10){ + g2(i) = 1.0e-10; + } + } + // Compute mu_zm at midpoints + Eigen::VectorXd mu_zm(nlay); + mu_zm = (mu_in.head(nlay) + mu_in.tail(nlay)) / 2.0; + Eigen::VectorXd g3 = (1.0 - sqrt3 * hg.array() * mu_zm.array()) / 2.0; + Eigen::VectorXd g4 = 1.0 - g3.array(); + + // Compute lam and gam + Eigen::VectorXd lam = (g1.array().square() - g2.array().square()).sqrt(); + Eigen::VectorXd gam = (g1.array() - lam.array()) / g2.array(); + + // Compute denom and handle denom == 0 + Eigen::VectorXd denom = lam.array().square() - (1.0 / (mu_in(nlev - 1) * mu_in(nlev - 1))); + for(int i = 0; i < nlay; ++i){ + if(std::abs(denom(i)) < 1e-10){ + denom(i) = 1.0e-10; + } + } + + // Compute Am and Ap + Eigen::VectorXd Am = F0_in * w0.array() * (g4.array() * (g1.array() + 1.0 / mu_in(nlev - 1)) + g2.array() * g3.array()) / denom.array(); + Eigen::VectorXd Ap = F0_in * w0.array() * (g3.array() * (g1.array() - 1.0 / mu_in(nlev - 1)) + g2.array() * g4.array()) / denom.array(); + + // Compute Cpm1 and Cmm1 at the top of the layer + Eigen::VectorXd Cpm1 = Ap.array() * (-tau_total.head(nlay).array() / mu_in(nlev - 1)).exp(); + Eigen::VectorXd Cmm1 = Am.array() * (-tau_total.head(nlay).array() / mu_in(nlev - 1)).exp(); + + // Compute Cp and Cm at the bottom of the layer + Eigen::VectorXd Cp = Ap.array() * (-tau_total.segment(1, nlay).array() / mu_in(nlev - 1)).exp(); + Eigen::VectorXd Cm = Am.array() * (-tau_total.segment(1, nlay).array() / mu_in(nlev - 1)).exp(); + + // Compute exponential terms, clamped to prevent overflow + Eigen::VectorXd exptrm = (lam.array() * dtau.array()).min(35.0); + Eigen::VectorXd Ep = exptrm.array().exp(); + Eigen::VectorXd Em = 1.0 / Ep.array(); + Eigen::VectorXd E1 = (Ep.array() + gam.array() * Em.array()).matrix(); + Eigen::VectorXd E2 = (Ep.array() - gam.array() * Em.array()).matrix(); + Eigen::VectorXd E3 = (gam.array() * Ep.array() + Em.array()).matrix(); + Eigen::VectorXd E4 = (gam.array() * Ep.array() - Em.array()).matrix(); + + // Initialize Af, Bf, Cf, Df + int l = 2 * nlay; + Eigen::VectorXd Af_vec = Eigen::VectorXd::Zero(l); + Eigen::VectorXd Bf_vec = Eigen::VectorXd::Zero(l); + Eigen::VectorXd Cf_vec = Eigen::VectorXd::Zero(l); + Eigen::VectorXd Df_vec = Eigen::VectorXd::Zero(l); + + // Boundary conditions at the top + Af_vec(0) = 0.0; + Bf_vec(0) = gam(0) + 1.0; + Cf_vec(0) = gam(0) - 1.0; + Df_vec(0) = btop - Cmm1(0); + for(int i =1, n=1; i < l -1; i +=2, ++n){ + if(n >= nlay){ + throw std::out_of_range("Index out of range in sw_Toon89 Af, Bf, Cf, Df population."); + } + Af_vec(i) = (E1(n-1) + E3(n-1)) * (gam(n) - 1.0); + Bf_vec(i) = (E2(n-1) + E4(n-1)) * (gam(n) - 1.0); + Cf_vec(i) = 2.0 * (1.0 - gam(n) * gam(n)); + Df_vec(i) = (gam(n) - 1.0) * (Cpm1(n) - Cp(n - 1)) + (1.0 - gam(n)) * (Cm(n - 1) - Cmm1(n)); + } + + // Populate Af, Bf, Cf, Df for even indices + // Start from n=1 to avoid negative indexing (Cp(n-1) when n=0) + for(int i = 2, n=1; i < l -1; i +=2, ++n){ + if(n >= nlay){ + throw std::out_of_range("Index out of range in sw_Toon89 Af, Bf, Cf, Df population."); + } + Af_vec(i) = 2.0 * (1.0 - gam(n) * gam(n)); + Bf_vec(i) = (E1(n-1) - E3(n-1)) * (1.0 + gam(n)); + Cf_vec(i) = (E1(n-1) + E3(n-1)) * (gam(n) - 1.0); + Df_vec(i) = E3(n-1) * (Cpm1(n) - Cp(n - 1)) + E1(n-1) * (Cm(n - 1) - Cmm1(n)); + } + + // Boundary conditions at l (last index) + Af_vec(l - 1) = E1(nlay - 1) - w_surf_in * E3(nlay - 1); + Bf_vec(l - 1) = E2(nlay - 1) - w_surf_in * E4(nlay - 1); + Cf_vec(l - 1) = 0.0; + Df_vec(l - 1) = bsurf - Cp(nlay - 1) + w_surf_in * Cm(nlay - 1); + + // Prepare a, b, c, d for the solver + Eigen::VectorXd a_tridiag = Af_vec.segment(1, l - 1); + Eigen::VectorXd b_tridiag = Bf_vec; + Eigen::VectorXd c_tridiag = Cf_vec.segment(0, l - 1); + Eigen::VectorXd d_tridiag = Df_vec; + + // Solve the tridiagonal system + Eigen::VectorXd xk = tridiagonal_solver(a_tridiag, b_tridiag, c_tridiag, d_tridiag); + + // Compute xk1 and xk2 from xk + Eigen::VectorXd xk1(nlay); + Eigen::VectorXd xk2(nlay); + for(int idx = 0; idx < nlay; ++idx){ + int two_n = 2 * idx; + if(two_n + 1 >= xk.size()){ + throw std::out_of_range("Index out of range when accessing xk."); + } + xk1(idx) = xk(two_n) + xk(two_n + 1); + xk2(idx) = xk(two_n) - xk(two_n + 1); + if(std::abs(xk2(idx) / xk(two_n)) < 1e-30){ + xk2(idx) = 0.0; + } + } + + // Populate flx_up and flx_down for layers 1 to nlay + flx_up.head(nlay) = (xk1.array() + gam.array() * xk2.array() + Cpm1.array()).matrix(); + flx_down.head(nlay) = (xk1.array() * gam.array() + xk2.array() + Cmm1.array()).matrix(); + + // Compute flx_up and flx_down at level nlev + flx_up(nlev - 1) = xk1(nlay - 1) * std::exp(1.0) + gam(nlay - 1) * xk2(nlay - 1) * std::exp(-1.0) + Cp(nlay - 1); + flx_down(nlev - 1) = xk1(nlay - 1) * std::exp(1.0) * gam(nlay - 1) + xk2(nlay - 1) * std::exp(-1.0) + Cm(nlay - 1); + + // Compute dir flux + Eigen::VectorXd dir = Eigen::VectorXd::Zero(nlev); + double mu_top_nonzero = mu_in(nlev - 1); + double mu_first_nonzero = mu_in(0); + if(std::abs(mu_top_nonzero - mu_first_nonzero) < 1e-12){ + // No zenith correction + dir = F0_in * mu_top_nonzero * (-tau_in.array() / mu_top_nonzero).exp(); + } + else{ + // Zenith angle correction + Eigen::VectorXd cum_trans(nlev); + cum_trans(0) = tau_total(0) / mu_in(0); + for(int k = 1; k < nlev; ++k){ + cum_trans(k) = cum_trans(k - 1) + (tau_total(k) - tau_total(k - 1)) / mu_in(k); + } + dir = F0_in * mu_top_nonzero * (-cum_trans.array()).exp(); + } + + // Adjust the downward flux at the surface layer for surface albedo + dir(nlev - 1) *= (1.0 - w_surf_in); + +//for(int i=0; i nlay){ + throw std::out_of_range("Index out of range in lw_Toon89 Af, Bf, Cf, Df population."); + } + Af_vec(i) = (E1(n-1) + E3(n-1)) * (gam(n) - 1.0); + Bf_vec(i) = (E2(n-1) + E4(n-1)) * (gam(n) - 1.0); + Cf_vec(i) = 2.0 * (1.0 - gam(n) * gam(n)); + Df_vec(i) = (gam(n) - 1.0) * (Cpm1(n) - Cp(n-1)) + (1.0 - gam(n)) * (Cm(n-1) - Cmm1(n)); + } + + // Populate Af, Bf, Cf, Df for even indices (i=2,4,...) + n = 1; // Reset n to 1 for even indices + for(int i=2; i < l -1; i +=2, ++n){ + if(n > nlay){ + throw std::out_of_range("Index out of range in lw_Toon89 Af, Bf, Cf, Df population."); + } + Af_vec(i) = 2.0 * (1.0 - gam(n-1) * gam(n-1)); + Bf_vec(i) = (E1(n-1) - E3(n-1)) * (1.0 + gam(n)); + Cf_vec(i) = (E1(n-1) + E3(n-1)) * (gam(n) - 1.0); + Df_vec(i) = E3(n-1) * (Cpm1(n) - Cp(n-1)) + E1(n-1) * (Cm(n-1) - Cmm1(n)); + } + + // Boundary conditions at the last index + Af_vec(l -1) = E1(nlay -1) - a_surf_in * E3(nlay -1); + Bf_vec(l -1) = E2(nlay -1) - a_surf_in * E4(nlay -1); + Cf_vec(l -1) = 0.0; + Df_vec(l -1) = Bsurf - Cp(nlay -1) + a_surf_in * Cm(nlay -1); + + // Prepare a, b, c, d for the tridiagonal solver + Eigen::VectorXd a_tridiag = Af_vec.segment(1, l -1); + Eigen::VectorXd b_tridiag = Bf_vec; + Eigen::VectorXd c_tridiag = Cf_vec.segment(0, l -1); + Eigen::VectorXd d_tridiag = Df_vec; + + // Solve the tridiagonal system + Eigen::VectorXd xkk = tridiagonal_solver(a_tridiag, b_tridiag, c_tridiag, d_tridiag); + + // Compute xk1 and xk2 from xkk + Eigen::VectorXd xk1 = Eigen::VectorXd::Zero(nlay); + Eigen::VectorXd xk2 = Eigen::VectorXd::Zero(nlay); + for(int idx=0; idx < nlay; ++idx){ + int two_n = 2 * idx; + if(two_n +1 >= xkk.size()){ + throw std::out_of_range("Index out of range when accessing xk."); + } + xk1(idx) = xkk(two_n) + xkk(two_n +1); + xk2(idx) = xkk(two_n) - xkk(two_n +1); + if(std::abs(xk2(idx) / xkk(two_n)) < 1e-30){ + xk2(idx) = 0.0; + } + } + + // Apply conditional "where (w0 <= 1e-4)" + // Using Eigen's select function for element-wise conditional assignments + Eigen::ArrayXd mask = (w0.array() <= 1e-4).cast(); + + Eigen::VectorXd g_var_final = mask * 0.0 + (1.0 - mask) * (twopi * w0.array() * xk1.array() * (1.0 + hg.array() * alp.array()) / (1.0 + alp.array())); + Eigen::VectorXd h_var_final = mask * 0.0 + (1.0 - mask) * (twopi * w0.array() * xk2.array() * (1.0 - hg.array() * alp.array()) / (1.0 + alp.array())); + Eigen::VectorXd xj_final = mask * 0.0 + (1.0 - mask) * (twopi * w0.array() * xk1.array() * (1.0 - hg.array() * alp.array()) / (1.0 + alp.array())); + Eigen::VectorXd xk_final = mask * 0.0 + (1.0 - mask) * (twopi * w0.array() * xk2.array() * (1.0 + hg.array() * alp.array()) / (1.0 + alp.array())); + Eigen::VectorXd alpha1_vec_final = mask * (twopi * B0.array()) + (1.0 - mask) * (twopi * (B0.array() + B1.array() * (ubari * w0.array() * hg.array() / (1.0 - w0.array() * hg.array())))); + Eigen::VectorXd alpha2_vec_final = twopi * B1.array(); + Eigen::VectorXd sigma1_final = mask * (twopi * B0.array()) + (1.0 - mask) * (twopi * (B0.array() - B1.array() * (ubari * w0.array() * hg.array() / (1.0 - w0.array() * hg.array())))); + Eigen::VectorXd sigma2_final = alpha2_vec_final; + + // Compute obj, epp, em1, obj2, epp2 + Eigen::VectorXd obj_vec = (lam.array() * dtau.array()).min(35.0); + Eigen::VectorXd epp_vec_final = obj_vec.array().exp(); + Eigen::VectorXd em1_final = 1.0 / epp_vec_final.array(); + Eigen::VectorXd obj2_vec = (0.5 * lam.array() * dtau.array()).min(35.0); + Eigen::VectorXd epp2_vec_final = obj2_vec.array().exp(); + + // Initialize and resize flux arrays to ensure correct dimensions + flx_up = Eigen::VectorXd::Zero(nlev); + flx_down = Eigen::VectorXd::Zero(nlev); + + // Initialize lw_up_g and lw_down_g + Eigen::VectorXd lw_up_g = Eigen::VectorXd::Zero(nlev); + Eigen::VectorXd lw_down_g = Eigen::VectorXd::Zero(nlev); + + // Loop over m=1 to nmu (0 to nmu-1 in C++) + for(int m=0; m < nmu; ++m){ + // Compute em2 and em3 + Eigen::VectorXd em2_vec = (-dtau.array() / uarr(m)).exp(); + Eigen::VectorXd em3_vec = em1_final.array() * em2_vec.array(); + + // Downward loop + lw_down_g(0) = twopi * (1.0 - std::exp(-tautop / uarr(m))) * be(0); + for(int k=0; k < nlay; ++k){ + lw_down_g(k+1) = lw_down_g(k) * em2_vec(k) + + (xj_final(k) / (lam(k) * uarr(m) + 1.0)) * (epp_vec_final(k) - em2_vec(k)) + + (xk_final(k) / (lam(k) * uarr(m) - 1.0)) * (em2_vec(k) - em1_final(k)) + + sigma1_final(k) * (1.0 - em2_vec(k)) + + sigma2_final(k) * (uarr(m) * em2_vec(k) + dtau(k) - uarr(m)); + } + + // Upward loop + lw_up_g(nlev -1) = twopi * (Bsurf + B1(nlay -1) * uarr(m)); + for(int k = nlay -1; k >=0; --k){ + lw_up_g(k) = lw_up_g(k+1) * em2_vec(k) + + (g_var_final(k) / (lam(k) * uarr(m) - 1.0)) * (epp_vec_final(k) * em2_vec(k) - 1.0) + + (h_var_final(k) / (lam(k) * uarr(m) + 1.0)) * (1.0 - em3_vec(k)) + + alpha1_vec_final(k) * (1.0 - em2_vec(k)) + + alpha2_vec_final(k) * (uarr(m) - (dtau(k) + uarr(m)) * em2_vec(k)); + } + + // Sum up flux arrays with Gaussian quadrature weights + flx_down += (lw_down_g.array() * wuarr(m)).matrix(); + flx_up += (lw_up_g.array() * wuarr(m)).matrix(); + } +//for(int i=0; i (m) * v[i]; + exm *= ex; + d[i] += exm * (6.0 + mv * (6.0 + mv * (3.0 + mv))) / (m * m); + } + d[i] *= conc; + } + } + + double ans; + if (smallv == 2) { + // Both wavenumbers are small + ans = p[1] - p[0]; + } else if (smallv == 1) { + // One wavenumber is small, the other is large + ans = 1.0 - p[0] - d[1]; + } else { + // Both wavenumbers are large + ans = d[0] - d[1]; + } + + ans *= sigdpi * T * T * T * T; + + if (ans == 0.0) { + std::cerr << "BB_integrate: Warning - result is zero; possible underflow" << std::endl; + } + + return ans; +} + + +//Tridiagonal Solver using the Thomas Algorithm +inline Eigen::VectorXd RadiationBand::RTSolverToon::tridiagonal_solver( + const Eigen::VectorXd& a, const Eigen::VectorXd& b, + const Eigen::VectorXd& c, const Eigen::VectorXd& d) +{ + int l = b.size(); + if(a.size() != static_cast(l - 1) || c.size() != static_cast(l - 1) || d.size() != static_cast(l)){ + throw std::invalid_argument("Incorrect vector sizes for tridiagonal_solver."); + } + + Eigen::VectorXd c_prime(l - 1); + Eigen::VectorXd d_prime(l); + + // Forward sweep + c_prime(0) = c(0) / b(0); + d_prime(0) = d(0) / b(0); + for(int i = 1; i < l - 1; ++i){ + double denom = b(i) - a(i - 1) * c_prime(i - 1); + if(std::abs(denom) < 1e-12){ + throw std::runtime_error("Tridiagonal solver failed: near-zero denominator."); + } + c_prime(i) = c(i) / denom; + d_prime(i) = (d(i) - a(i - 1) * d_prime(i - 1)) / denom; + } + + // Last equation + double denom_last = b(l - 1) - a(l - 2) * c_prime(l - 2); + if(std::abs(denom_last) < 1e-12){ + throw std::runtime_error("Tridiagonal solver failed: near-zero denominator at last equation."); + } + d_prime(l - 1) = (d(l - 1) - a(l - 2) * d_prime(l - 2)) / denom_last; + + // Back substitution + Eigen::VectorXd x(l); + x(l - 1) = d_prime(l - 1); + for(int i = l - 2; i >= 0; --i){ + x(i) = d_prime(i) - c_prime(i) * x(i + 1); + } + + return x; +} + +#endif + diff --git a/src/harp/rt_solvers.hpp b/src/harp/rt_solvers.hpp new file mode 100644 index 00000000..5f069d39 --- /dev/null +++ b/src/harp/rt_solvers.hpp @@ -0,0 +1,155 @@ +#ifndef SRC_HARP_RT_SOLVERS_HPP_ +#define SRC_HARP_RT_SOLVERS_HPP_ + +// C/C++ +#include +#include // Needed for Eigen::VectorXd +#include +#include +#include +#include +#include +#include + +// athena +#include + +// application +#include +#include + +// canoe +#include +#include + +// cppdisort +#ifdef RT_DISORT +#include +#endif + +// harp +#include "radiation_band.hpp" + +class RadiationBand::RTSolver : public NamedGroup { + public: // constructor and destructor + RTSolver(RadiationBand *pmy_band, std::string name) + : NamedGroup(name), pmy_band_(pmy_band) { + Application::Logger app("harp"); + app->Log("Initialize RTSolver " + GetName()); + } + + virtual ~RTSolver() { + Application::Logger app("harp"); + app->Log("Destroy RTSolver " + GetName()); + } + + public: // member functions + //! \brief Prepare and seal the solver for the current column + virtual void Prepare(MeshBlock const *pmb, int k, int j) {} + + //! \brief Allocate memory for radiation solver + virtual void Resize(int nlyr, int nstr) { + farea_.DeleteAthenaArray(); + vol_.DeleteAthenaArray(); + + farea_.NewAthenaArray(nlyr + 2 * NGHOST); + vol_.NewAthenaArray(nlyr + 2 * NGHOST); + } + + public: // inbound functions + virtual void CalBandFlux(MeshBlock const *pmb, int k, int j) { + throw NotImplementedError("CalBandFlux not implemented."); + } + virtual void CalBandRadiance(MeshBlock const *pmb, int k, int j) { + throw NotImplementedError("CalBandRadiance not implemented."); + } + + protected: + RadiationBand *pmy_band_; + AthenaArray farea_, vol_; +}; + +class RadiationBand::RTSolverLambert : public RadiationBand::RTSolver { + public: // constructor and destructor + RTSolverLambert(RadiationBand *pmy_band, YAML::Node const &rad) + : RTSolver(pmy_band, "Lambert") {} + ~RTSolverLambert() {} + + public: // inbound functions + void CalBandRadiance(MeshBlock const *pmb, int k, int j) override; +}; + +#ifdef RT_DISORT +class RadiationBand::RTSolverDisort : public RadiationBand::RTSolver, + protected DisortWrapper { + public: // constructor and destructor + RTSolverDisort(RadiationBand *pmy_band, YAML::Node const &rad); + ~RTSolverDisort() {} + + public: // member functions + void Prepare(MeshBlock const *pmb, int k, int j) override; + void Resize(int nlyr, int nstr) override; + + public: // inbound functions + void CalBandFlux(MeshBlock const *pmb, int k, int j) override; + void CalBandRadiance(MeshBlock const *pmb, int k, int j) override; + + protected: + size_t dir_dim_[2]; + std::vector dir_axis_; + + void addDisortFlux(Coordinates const *pcoord, int n, int k, int j, int il, + int iu); + + void addDisortRadiance(int n, int k, int j); +}; + + +class RadiationBand::RTSolverToon : public RadiationBand::RTSolver, + protected DisortWrapper { + public: // Constructor and Destructor + RTSolverToon(RadiationBand *pmy_band, YAML::Node const &rad); + ~RTSolverToon() {} + + public: // Member Functions + void Prepare(MeshBlock const *pmb, int k, int j) override; + void Resize(int nlyr, int nstr) override; + + public: // Inbound Functions + void CalBandFlux(MeshBlock const *pmb, int k, int j) override; + + protected: + void addToonFlux(Coordinates const *pcoord, int b, int k, int j, int il, int iu, + const Eigen::VectorXd &flux_up, const Eigen::VectorXd &flux_down); + + void toonShortwaveSolver(int nlay, double F0_in, + const Eigen::VectorXd& mu_in, + const Eigen::VectorXd& tau_in, + const Eigen::VectorXd& w_in, + const Eigen::VectorXd& g_in, + double w_surf_in, + Eigen::VectorXd& flx_up, + Eigen::VectorXd& flx_down); + + void toonLongwaveSolver(int nlay, + const Eigen::VectorXd& be, + const Eigen::VectorXd& tau_in, + const Eigen::VectorXd& w_in, + const Eigen::VectorXd& g_in, + double a_surf_in, + Eigen::VectorXd& flx_up, + Eigen::VectorXd& flx_down); + + double BB_integrate(double T, double wn1, double wn2); + + private: + // Private Tridiagonal Solver using the Thomas Algorithm + inline Eigen::VectorXd tridiagonal_solver(const Eigen::VectorXd& a, // Sub-diagonal (size n-1) + const Eigen::VectorXd& b, // Main diagonal (size n) + const Eigen::VectorXd& c, // Super-diagonal (size n-1) + const Eigen::VectorXd& d); // Right-hand side (size n) +}; + +#endif + +#endif // SRC_HARP_RT_SOLVERS_HPP_ diff --git a/src/harp/set_spectral_properties.cpp b/src/harp/set_spectral_properties.cpp new file mode 100644 index 00000000..b5b2b51e --- /dev/null +++ b/src/harp/set_spectral_properties.cpp @@ -0,0 +1,143 @@ +// cnaoe +#include + +// athena +#include +#include +#include +#include +#include + +// canoe +#include +#include +#include + +// snap +#include +#include + +// opacity +#include +#include + +// harp +#include "radiation.hpp" +#include "radiation_band.hpp" + +// setting optical properties +void RadiationBand::SetSpectralProperties(AirColumn& ac, Real const* x1f, + Real gH0, int k, int j) { + int nspec = GetNumSpecGrids(); + int npmom = GetNumPhaseMoments(); + + // set tau, ssalb, pmom, etc... + tau_.ZeroClear(); + ssa_.ZeroClear(); + pmom_.ZeroClear(); + + std::vector mypmom(1 + npmom); + + for (int i = 0; i < ac.size(); ++i) { + auto& air = ac[i]; + air.ToMoleFraction(); + tem_[i] = air.w[IDN]; + + for (auto& a : absorbers_) { + for (int m = 0; m < nspec; ++m) { + auto& spec = pgrid_->spec[m]; + Real kcoeff; // Variable to store attenuation coefficient + if (auto* ck = dynamic_cast(a.get())) { + kcoeff = ck->GetAttenuation(m, air); + } else { + kcoeff = a->GetAttenuation(spec.wav1, spec.wav2, air); + } + Real dssalb = + a->GetSingleScatteringAlbedo(spec.wav1, spec.wav2, air) * kcoeff; + // tau + tau_(m, i) += kcoeff; + + // ssalb + ssa_(m, i) += dssalb; + // pmom + a->GetPhaseMomentum(mypmom.data(), spec.wav1, spec.wav2, air, npmom); + for (int p = 0; p <= npmom; ++p) pmom_(m, i, p) += mypmom[p] * dssalb; + + } + } + } + + // set temperature at cell interface + int il = NGHOST, iu = ac.size() - 1 - NGHOST; + temf_[il] = (3. * tem_[il] - tem_[il + 1]) / 2.; + temf_[il + 1] = (tem_[il] + tem_[il + 1]) / 2.; + for (int i = il + 2; i <= iu - 1; ++i) + temf_[i] = interp_cp4(tem_[i - 2], tem_[i - 1], tem_[i], tem_[i + 1]); + temf_[iu] = (tem_[iu] + tem_[iu - 1]) / 2.; + // temf_[iu + 1] = (3. * tem_[iu] - tem_[iu - 1]) / 2.; + temf_[iu + 1] = tem_[iu]; // isothermal top boundary + + for (int i = 0; i < il; ++i) temf_[i] = tem_[il]; + for (int i = iu + 2; i < ac.size(); ++i) temf_[i] = tem_[iu + 1]; + + bool error = false; + for (int i = 0; i < ac.size(); ++i) { + if (temf_[i] < 0.) { + temf_[i] = tem_[i]; + // error = true; + } + } + for (int i = il; i <= iu; ++i) { + if (tem_[i] < 0.) error = true; + } + if (error) { + for (int i = il; i <= iu; ++i) { + std::cout << "--- temf[" << i << "] = " << temf_[i] << std::endl; + std::cout << "tem[" << i << "] = " << tem_[i] << std::endl; + } + std::cout << "--- temf[" << iu + 1 << "] = " << temf_[iu + 1] << std::endl; + throw std::runtime_error("Negative temperature at cell interface"); + } + + // absorption coefficients -> optical thickness + for (int m = 0; m < nspec; ++m) { + for (int i = 0; i < ac.size(); ++i) { + if (tau_(m, i) > 1e-10 && ssa_(m, i) > 1e-10) { // has scattering + for (int p = 0; p <= npmom; ++p) pmom_(m, i, p) /= ssa_(m, i); + ssa_(m, i) /= tau_(m, i); + } else { + ssa_(m, i) = 0.; + pmom_(m, i, 0) = 1.; + for (int p = 1; p <= npmom; ++p) pmom_(m, i, p) = 0.; + } + +#ifdef HYDROSTATIC + auto pthermo = Thermodynamics::GetInstance(); + Real Rgas = get_rovrd(ac[i], pthermo->GetMuRatio()) * pthermo->GetRd(); + // TODO(cli) check this + // \delta z = \delt Z * (R T)/(g H0) + tau_(m, i) *= (x1f[i + 1] - x1f[i]) * (Rgas * tem_[i]) / (gH0); +#else + tau_(m, i) *= x1f[i + 1] - x1f[i]; +#endif + } + } + + // aggregated band properties + for (int i = 0; i < ac.size(); ++i) { + btau(k, j, i) = 0; + bssa(k, j, i) = 0; + for (int p = 0; p <= npmom; ++p) bpmom(p, k, j, i) = 0.; + + for (int m = 0; m < nspec; ++m) { + btau(k, j, i) += tau_(m, i); + bssa(k, j, i) += ssa_(m, i) * tau_(m, i); + for (int p = 0; p <= npmom; ++p) + bpmom(p, k, j, i) += pmom_(m, i, p) * ssa_(m, i); + } + + for (int p = 0; p <= npmom; ++p) bpmom(p, k, j, i) /= bssa(k, j, i); + bssa(k, j, i) /= btau(k, j, i); + btau(k, j, i) /= nspec; + } +} diff --git a/src/opacity/absorber.hpp b/src/opacity/absorber.hpp new file mode 100644 index 00000000..91067da0 --- /dev/null +++ b/src/opacity/absorber.hpp @@ -0,0 +1,114 @@ +#ifndef SRC_OPACITY_ABSORBER_HPP_ +#define SRC_OPACITY_ABSORBER_HPP_ + +// C/C++ +#include +#include +#include +#include + +// athena +#include + +// canoe +#include + +class AirParcel; +class SpectralBin; +namespace YAML { +class Node; +} + +//! \brief base class of all absorbers +class Absorber : public NamedGroup, + public ParameterGroup, + public SpeciesIndexGroup, + public StringReprGroup { + public: // constructor and destructor + Absorber(std::string name); + virtual ~Absorber(); + + public: // member functions + //! Set absorption model + void SetModel(std::string name) { model_name_ = name; } + + //! Combines SetOpacityFile() and LoadOpacity() + void LoadOpacityFromFile(std::string filename); + + //! Set opacity filename to internal variable, does not load opacity + void SetOpacityFile(std::string filename); + + //! Get opacity filename + std::string GetOpacityFile() const { return opacity_filename_; } + + //! Load opacity from internal variable + void LoadOpacity(int bid); + + //! Load absorption coefficient from file + virtual void LoadCoefficient(std::string fname, int bid) {} + + //! Get attenuation coefficient [1/m] + virtual Real GetAttenuation(Real wave1, Real wave2, + AirParcel const& var) const { + return 0.; + } + + //! Get attenuation coefficient for CK + virtual Real GetAttenuation(int m, AirParcel const& var) const { + return 0.; + } + + //! Get single scattering albedo [1] + virtual Real GetSingleScatteringAlbedo(Real wave1, Real wave2, + AirParcel const& var) const { + return 0.; + } + + //! Get phase function [1] + virtual void GetPhaseMomentum(Real* pp, Real wave1, Real wave2, + AirParcel const& var, int np) const {} + + virtual void CheckFail() const {} + + virtual void ModifySpectralGrid(std::vector& spec) const {} + + public: // StringRepr + std::string ToString() const override; + + protected: + //! absorption model model + std::string model_name_; + + //! opacity filename + std::string opacity_filename_; +}; + +using AbsorberPtr = std::shared_ptr; +using AbsorberContainer = std::vector; + +class AbsorberFactory { + public: + //! \brief Create an absorber from YAML node + //! + //! \param[in] my YAML node containing the current absorber + //! \param[in] band_name name of the radiation band + static AbsorberPtr CreateFrom(YAML::Node const& my, std::string band_name); + + //! \brief Create absorbers from YAML node + //! + //! \param[in] names names of absorbers + //! \param[in] band_name name of the radiation band + //! \param[in] rad YAML node containing the radiation input file + static AbsorberContainer CreateFrom(std::vector const& names, + std::string band_name, + YAML::Node const& rad); + + protected: + //! \brief Only create an absorber based on its name and class + //! + //! \param[in] name name of the absorber + //! \param[in] type class identifier of the absorber + static AbsorberPtr createAbsorberPartial(std::string name, std::string type); +}; + +#endif // SRC_OPACITY_ABSORBER_HPP_ diff --git a/src/opacity/absorber_ck.cpp b/src/opacity/absorber_ck.cpp new file mode 100644 index 00000000..158aca97 --- /dev/null +++ b/src/opacity/absorber_ck.cpp @@ -0,0 +1,170 @@ +// C/C++ +#include +#include +#include +#include +#include +#include +#include +#include +#include + +// canoe +#include +#include +#include + +// climath +#include + +// harp +#include + +// athena +#include + +// snap +#include + +// opacity +#include "absorber_ck.hpp" + +// netcdf +#ifdef NETCDFOUTPUT +extern "C" { +#include +} +#endif + +void AbsorberCK::LoadCoefficient(std::string fname, int) { +#ifdef NETCDFOUTPUT + int fileid, varid, err; + int dimids[10]; // Array to hold dimension IDs (assuming a maximum of 10 dimensions) + + size_t len_p, len_t, len_bin_centers, len_weights; + +std::map var_names = { + {"p", "p"}, // Pressure (logarithmic, in Pa) + {"t", "t"}, // Temperature (K) + {"bin_centers", "bin_centers"}, // Wavenumber bin centers (cm⁻¹) + {"bin_edges", "bin_edges"}, // Wavenumber bin edges (cm⁻¹) + {"weights", "weights"}, // Gaussian quadrature weights + {"gaussian_points", "gaussian_points"}, // Gaussian quadrature points + {"kcoeff", "kcoeff"} // Logarithmic absorption coefficients (log(m²/molecule)) +}; + + nc_open(fname.c_str(), NC_NOWRITE, &fileid); + + for (const auto& [role, var_name] : var_names) { + if (role == "kcoeff" || role == "bin_edges") continue; // Handle separately + + nc_inq_varid(fileid, var_name.c_str(), &varid); + nc_inq_vardimid(fileid, varid, dimids); + + // Assuming each variable is 1D, get the length of the first dimension + size_t dim_len; + nc_inq_dimlen(fileid, dimids[0], &dim_len); + + // Assign to len_ based on role + if (role == "p") { + len_p = dim_len; + len_[0] = len_p; + } + if (role == "t") { + len_t = dim_len; + len_[1] = len_t; + } + if (role == "bin_centers") { + len_bin_centers = dim_len; + len_[2] = len_bin_centers; + } + if (role == "weights") { + len_weights = dim_len; + len_[3] = len_weights; + } + } + + p_.resize(len_[0]); // Raw pressure + t_.resize(len_[1]); // Raw temperature + axis_.resize(len_[0] + len_[1]); // Pressure + Temperature + bin_centers_.resize(len_[2]); // Bin centers + weights_.resize(len_[3]); // Weights + bin_edges_.resize(len_[2] + 1); // Bin edges (len_bin_centers + 1) + kcoeff_.resize(len_[0] * len_[1] * len_[2] * len_[3]); // Absorption coefficients + + for (const auto& [role, var_name] : var_names) { + if (role == "kcoeff" || role == "bin_edges") continue; // Handle separately + + // Get the variable ID + nc_inq_varid(fileid, var_name.c_str(), &varid); + + // Read the data into the corresponding container + if (role == "p") { + nc_get_var_double(fileid, varid, axis_.data()); + nc_get_var_double(fileid, varid, p_.data()); + } + + if (role == "t") { + nc_get_var_double(fileid, varid, axis_.data() + len_p); + nc_get_var_double(fileid, varid, t_.data()); + } + if (role == "bin_centers") { + nc_get_var_double(fileid, varid, bin_centers_.data()); + } + if (role == "weights") { + nc_get_var_double(fileid, varid, weights_.data()); + } + } + + nc_inq_varid(fileid, var_names["bin_edges"].c_str(), &varid); + nc_get_var_double(fileid, varid, bin_edges_.data()); + + nc_inq_varid(fileid, var_names["kcoeff"].c_str(), &varid); + nc_get_var_double(fileid, varid, kcoeff_.data()); + + nc_close(fileid); + + // Optional: You can log some debug output to check that everything was loaded correctly + std::cout << "Loaded coefficient grid: " + << "Pressure levels = " << len_[0] << ", " + << "Temperature levels = " << len_[1] << ", " + << "Wavelength bins = " << len_[2] << ", " + << "Weights = " << len_[3] << std::endl; +#endif +} + +Real AbsorberCK::GetAttenuation(int m, AirParcel const& var) const { + int iw = m / len_[3]; // Wavelength bin index + int ig = m % len_[3]; // Weight index + Real val; + Real coord[2] = {log(var.w[IPR]), var.w[IDN]}; + + // Calculate the starting index for the 2D slice (p, t) at iw, ig + size_t base_index = iw * len_[3] * len_[0] * len_[1] + ig * len_[0] * len_[1]; + size_t len_2d[2] = {len_[0], len_[1]}; + + // Perform 2D interpolation using the flattened kcoeff_ + interpn(&val, coord, &kcoeff_[base_index], axis_.data(), len_2d, 2, 1); + + auto pthermo = Thermodynamics::GetInstance(); + Real dens = var.w[IPR] / (pthermo->GetRd() * var.w[IDN]); + + return exp(val) * dens; // ln(m^2/kg) -> 1/m +} + +void AbsorberCK::ModifySpectralGrid(std::vector& spec) const { + // Resize the spec grid to account for both wavelength bins and weights (Gauss points) + spec.resize(len_[2] * len_[3]); // len_[2] = nbin (wavelength bins), len_[3] = nwght (Gauss points) + + // Loop through both wavelength bins and weights + for (size_t iw = 0; iw < len_[2]; ++iw) { // len_[2] is the number of wavelength bins + for (size_t ig = 0; ig < len_[3]; ++ig) { // len_[3] is the number of weights (Gauss points) + size_t m = iw * len_[3] + ig; // Flattened index for the spectral grid + + // Use bin_edges for wav1 and wav2 + spec[m].wav1 = bin_edges_[iw]; // Lower boundary (wav1) + spec[m].wav2 = bin_edges_[iw + 1]; // Upper boundary (wav2) + spec[m].wght = weights_[ig]; // Gauss point weight (from `weights`) + } + } +} diff --git a/src/opacity/absorber_ck.hpp b/src/opacity/absorber_ck.hpp new file mode 100644 index 00000000..0e11ef53 --- /dev/null +++ b/src/opacity/absorber_ck.hpp @@ -0,0 +1,42 @@ +#ifndef SRC_OPACITY_ABSORBER_CK_HPP_ +#define SRC_OPACITY_ABSORBER_CK_HPP_ + +// C/C++ +#include +#include + +// opacity +#include "absorber.hpp" + +class AbsorberCK : public Absorber { + public: + AbsorberCK(std::string name) : Absorber(name) {} + virtual ~AbsorberCK() {} + + void LoadCoefficient(std::string fname, int bid) override; + Real GetAttenuation(int m, AirParcel const& var) const override; + void ModifySpectralGrid(std::vector& spec) const override; + + protected: + size_t len_[4]; + + std::vector axis_; + std::vector p_; + std::vector t_; + std::vector bin_centers_; + std::vector bin_edges_; + std::vector kcoeff_; + std::vector weights_; +}; + +class HeliosCK : public AbsorberCK { + public: + HeliosCK(std::string name); + virtual ~HeliosCK() {} + + void LoadCoefficient(std::string fname, int bid) override; + Real GetAttenuation(Real g1, Real g2, AirParcel const& var) const override; + void ModifySpectralGrid(std::vector& spec) const override; +}; + +#endif // SRC_OPACITY_ABSORBER_CK_HPP_ diff --git a/src/opacity/absorber_factory.cpp b/src/opacity/absorber_factory.cpp new file mode 100644 index 00000000..ee3a7504 --- /dev/null +++ b/src/opacity/absorber_factory.cpp @@ -0,0 +1,147 @@ +// athena +#include + +// external +#include + +// application +#include + +// canoe +#include +#include + +// opacity +#include "Giants/grey_gas.hpp" +#include "Giants/hydrogen_cia.hpp" +#include "Giants/microwave/mwr_absorbers.hpp" +#include "absorber_ck.hpp" +#include "hitran_absorber.hpp" +#include "nitrogen_cia.hpp" +#include "oxygen_cia.hpp" +#include "simple_cloud.hpp" + +AbsorberContainer AbsorberFactory::CreateFrom( + std::vector const& names, std::string band_name, + YAML::Node const& rad) { + AbsorberContainer absorbers; + + for (auto& name : names) { + AbsorberContainer atmp; + + for (auto& my : rad["opacity-sources"]) { + if (name == my["name"].as()) { + atmp.push_back(AbsorberFactory::CreateFrom(my, band_name)); + } + } + + if (atmp.empty()) { + throw NotFoundError("AbsorberFactory", "Opacity " + name); + } else if (atmp.size() > + 1) { // use band_name to uniquely identify the absorber + bool found = false; + for (auto& ab : atmp) { + if (ab->GetOpacityFile().find(band_name) != std::string::npos) { + absorbers.push_back(ab); + found = true; + break; + } + } + + if (!found) { + throw NotFoundError("AbsorberFactory", + "There are multiple absorbers with the same name " + + name + + ". Rename the `data` field in the absorber to " + "identify the band."); + } + } else { + absorbers.push_back(atmp[0]); + } + } + + return absorbers; +} + +AbsorberPtr AbsorberFactory::CreateFrom(YAML::Node const& my, + std::string band_name) { + if (!my["name"]) { + throw NotFoundError("AbsorberFactory", "'name' field in absorber"); + } + + if (!my["class"]) { + throw NotFoundError("AbsorberFactory", "'class' field in absorber " + + my["name"].as()); + } + + auto ab = createAbsorberPartial(my["name"].as(), + my["class"].as()); + + if (my["data"]) { + ab->SetOpacityFile(my["data"].as()); + } else { + ab->SetOpacityFile(""); + } + + if (my["dependent-species"]) { + auto species = my["dependent-species"].as>(); + ab->SetSpeciesIndex(species); + } + + if (my["model"]) { + ab->SetModel(my["model"].as()); + } + + if (my["parameters"]) { + ab->SetRealsFrom(my["parameters"]); + } + + return ab; +} + +namespace gp = GiantPlanets; + +AbsorberPtr AbsorberFactory::createAbsorberPartial(std::string name, + std::string type) { + AbsorberPtr ab; + + if (type == "XIZ-H2-H2-CIA") { + ab = std::make_shared(); + } else if (type == "XIZ-H2-He-CIA") { + ab = std::make_shared(); + } else if (type == "Hitran") { + ab = std::make_shared(name); + } else if (type == "HitranCK") { + ab = std::make_shared(name); + } else if (type == "SimpleGrey") { + ab = std::make_shared(name); + } else if (type == "FreedmanMean") { + ab = std::make_shared(name); + } else if (type == "JupGasv") { + ab = std::make_shared(name); + } else if (type == "JupGasir") { + ab = std::make_shared(name); + } else if (type == "AbsorberCK") { + ab = std::make_shared(name); + } else if (type == "HeliosCK") { + ab = std::make_shared(name); + } else if (type == "SimpleCloud") { + ab = std::make_shared(name); + } else if (type == "radio-NH3") { + ab = std::make_shared(); + } else if (type == "radio-H2O") { + ab = std::make_shared(); + } else if (type == "radio-H2S") { + ab = std::make_shared(); + } else if (type == "radio-PH3") { + ab = std::make_shared(); + } else if (type == "radio-CIA") { + ab = std::make_shared(); + } else if (type == "radio-Electron") { + ab = std::make_shared(); + } else { + throw NotFoundError("createAbsorberPartial", "Class = " + type); + } + + return ab; +} From 7ca54223db06a25934c8d54c943266d4e5cee6c1 Mon Sep 17 00:00:00 2001 From: Xi Zhang Date: Tue, 8 Apr 2025 08:43:22 -0700 Subject: [PATCH 02/10] wip --- README | 11 +- cmake/examples/ice_h2s.cmake | 20 ++ cmake/examples/rtcldexo.cmake | 10 +- examples/2024-XZhang-cloud-rt/CMakeLists.txt | 6 +- examples/2024-XZhang-cloud-rt/ice_h2s.cpp | 195 ++++++++++++ examples/2024-XZhang-cloud-rt/jupiter_crm.cpp | 279 +++++++++++++++++ examples/2024-XZhang-cloud-rt/rtcldu3d.cpp | 288 +++++++++++++++++ .../2024-XZhang-cloud-rt/svp_functions.hpp | 53 ++++ examples/2024-XZhang-minichem/jup_onetra.cpp | 290 ++++++++++++++++++ examples/CMakeLists.txt | 1 + 10 files changed, 1144 insertions(+), 9 deletions(-) create mode 100644 cmake/examples/ice_h2s.cmake create mode 100644 examples/2024-XZhang-cloud-rt/ice_h2s.cpp create mode 100644 examples/2024-XZhang-cloud-rt/jupiter_crm.cpp create mode 100644 examples/2024-XZhang-cloud-rt/rtcldu3d.cpp create mode 100644 examples/2024-XZhang-minichem/jup_onetra.cpp diff --git a/README b/README index f05caadb..d59c7683 100644 --- a/README +++ b/README @@ -4,6 +4,13 @@ when vapor pressure are too low (~1e-10 Pa) at low temperatuers, the saturationa when vapor pressures are too high (1e5 Pa) eg. for CO2, the saturation adjustment fails immediately (all nan) -exo can use cores 24/96/384/1536 -also permit 216/600/864/1176/ if specify NBLOCKS in the .cmake file +exo can use cores 24/96/384/1536: 6*4^n +also permit 216/600/864/1176/1536/1944/2400... (24*n^2) if specify NBLOCKS in the .cmake file +That is, if meshblock is x, the mesh should be x*4*n, and x*6*n, for example, if x is 10, mesh could be 40/60, or 240/360 etc. +meshblock = 5 is not very efficient, maybe meshblock is 8 or higher + +gammal in svp is delta[n] = (cp_ratio_mass[n] - cp_ratio_mass[i]) * mu_ratio[i] / (1. - 1. / gammad); + +cp_ratio_mass is rcp +mu_ratio is eps diff --git a/cmake/examples/ice_h2s.cmake b/cmake/examples/ice_h2s.cmake new file mode 100644 index 00000000..b0fd9320 --- /dev/null +++ b/cmake/examples/ice_h2s.cmake @@ -0,0 +1,20 @@ +# configure file for test jupiter crm + +macro(SET_IF_EMPTY _variable) + if("${${_variable}}" STREQUAL "") + set(${_variable} ${ARGN}) + endif() +endmacro() + +# athena variables +set_if_empty(NUMBER_GHOST_CELLS 3) + +# canoe configure +set(NVAPOR 1) +set(NCLOUD 2) +set(NPHASE_LEGACY 3) +set(PNETCDF ON) +set(MPI ON) +set(TASKLIST ImplicitHydroTasks) +#set_if_empty(RSOLVER hllc_transform) +set(RSOLVER lmars) diff --git a/cmake/examples/rtcldexo.cmake b/cmake/examples/rtcldexo.cmake index 5f9adc11..61423a5c 100644 --- a/cmake/examples/rtcldexo.cmake +++ b/cmake/examples/rtcldexo.cmake @@ -9,14 +9,14 @@ endmacro() # athena variables set_if_empty(NUMBER_GHOST_CELLS 3) -set(CUBED_SPHERE ON) -#set(NBLOCKS 600) +#set(CUBED_SPHERE ON) +#set(NBLOCKS 2400) #set(COORDINATE_SYSTEM "spherical_polar") -set(NVAPOR 1) -set(NCLOUD 2) +set(NVAPOR 2) +set(NCLOUD 4) set(NPHASE_LEGACY 3) set(NETCDF ON) -#set(PNETCDF OFF) +set(PNETCDF ON) set(MPI ON) set(DISORT ON) set(PYTHON_BINDINGS ON) diff --git a/examples/2024-XZhang-cloud-rt/CMakeLists.txt b/examples/2024-XZhang-cloud-rt/CMakeLists.txt index 2142c0c7..9844e3d5 100644 --- a/examples/2024-XZhang-cloud-rt/CMakeLists.txt +++ b/examples/2024-XZhang-cloud-rt/CMakeLists.txt @@ -1,11 +1,13 @@ #if (${NVAPOR} EQUAL 1 AND ${NCLOUD} EQUAL 2 AND ${NPHASE_LEGACY} EQUAL 3) if (${NVAPOR} EQUAL 1) #setup_problem(rt_silicate) - setup_problem(beta_rtsi) + setup_problem(ice_h2s) + #setup_problem(beta_rtsi) #setup_problem(rt_h2o) elseif (${NVAPOR} EQUAL 2) #setup_problem(rt_h2o_nh3) - setup_problem(jupiter_crm) + #setup_problem(jupiter_crm) + setup_problem(rtcldu3d) elseif (${NVAPOR} EQUAL 3) setup_problem(beta_rtjup3cld) else() diff --git a/examples/2024-XZhang-cloud-rt/ice_h2s.cpp b/examples/2024-XZhang-cloud-rt/ice_h2s.cpp new file mode 100644 index 00000000..8514ac13 --- /dev/null +++ b/examples/2024-XZhang-cloud-rt/ice_h2s.cpp @@ -0,0 +1,195 @@ +/* ------------------------------------------------------------------------------------- + * SNAP Example Program + * + * Contributer: + * Cheng Li, University of Michigan + * + * Year: 2023 + * Contact: chengcli@umich.edu + * Reference: Test Jupiter CRM + * ------------------------------------------------------------------------------------- + */ + +// athena +#include +#include +#include +#include +#include +#include +#include +#include +#include + +// application +#include +#include + +// canoe +#include +#include + +// climath +#include +#include + +// snap +#include + +// special includes +#include "svp_functions.hpp" + +Real grav, P0, T0, Tmin, prad, hrate; +int iH2S; + +void MeshBlock::InitUserMeshBlockData(ParameterInput *pin) { + AllocateUserOutputVariables(4 + NVAPOR); + SetUserOutputVariableName(0, "temp"); + SetUserOutputVariableName(1, "theta"); + SetUserOutputVariableName(2, "thetav"); + SetUserOutputVariableName(3, "mse"); + SetUserOutputVariableName(4, "rh_H2S"); +} + +void MeshBlock::UserWorkBeforeOutput(ParameterInput *pin) { + auto pthermo = Thermodynamics::GetInstance(); + + for (int k = ks; k <= ke; ++k) + for (int j = js; j <= je; ++j) + for (int i = is; i <= ie; ++i) { + user_out_var(0, k, j, i) = pthermo->GetTemp(this, k, j, i); + user_out_var(1, k, j, i) = pthermo->PotentialTemp(this, P0, k, j, i); + // theta_v + user_out_var(2, k, j, i) = + user_out_var(1, k, j, i) * pthermo->RovRd(this, k, j, i); + // mse + user_out_var(3, k, j, i) = + pthermo->MoistStaticEnergy(this, grav * pcoord->x1v(i), k, j, i); + // relative humidity + user_out_var(4, k, j, i) = pthermo->RelativeHumidity(this, 1, k, j, i); + } +} + +void Forcing(MeshBlock *pmb, Real const time, Real const dt, + AthenaArray const &w, AthenaArray const &r, + AthenaArray const &bcc, AthenaArray &du, + AthenaArray &s) { + int is = pmb->is, js = pmb->js, ks = pmb->ks; + int ie = pmb->ie, je = pmb->je, ke = pmb->ke; + auto pthermo = Thermodynamics::GetInstance(); + + for (int k = ks; k <= ke; ++k) + for (int j = js; j <= je; ++j) + for (int i = is; i <= ie; ++i) { + auto &&air = AirParcelHelper::gather_from_primitive(pmb, k, j, i); + + air.ToMoleFraction(); + Real cv = + get_cv_mass(air, 0, pthermo->GetRd(), pthermo->GetCvRatioMass()); + +// if (w(IPR, k, j, i) < prad) { +// du(IEN, k, j, i) += dt * hrate * w(IDN, k, j, i) * cv * +// (1. + 1.E-4 * sin(2. * M_PI * rand() / RAND_MAX)); +// } + } +} + +void Mesh::InitUserMeshData(ParameterInput *pin) { + grav = -pin->GetReal("hydro", "grav_acc1"); + + P0 = pin->GetReal("problem", "P0"); + T0 = pin->GetReal("problem", "T0"); + + Tmin = pin->GetReal("problem", "Tmin"); +// prad = pin->GetReal("problem", "prad"); +// hrate = pin->GetReal("problem", "hrate") / 86400.; + + // index + auto pindex = IndexMap::GetInstance(); + iH2S = pindex->GetVaporId("H2S"); + EnrollUserExplicitSourceFunction(Forcing); +} + +void MeshBlock::ProblemGenerator(ParameterInput *pin) { + srand(Globals::my_rank + time(0)); + + Application::Logger app("main"); + app->Log("ProblemGenerator: jupiter_crm"); + + auto pthermo = Thermodynamics::GetInstance(); + + // mesh limits + Real x1min = pmy_mesh->mesh_size.x1min; + Real x1max = pmy_mesh->mesh_size.x1max; + + // request temperature and pressure + app->Log("request T", T0); + app->Log("request P", P0); + + // thermodynamic constants + Real gamma = pin->GetReal("hydro", "gamma"); + Real Rd = pthermo->GetRd(); + Real cp = gamma / (gamma - 1.) * Rd; + + // set up an adiabatic atmosphere + int max_iter = 400, iter = 0; + Real Ttol = pin->GetOrAddReal("problem", "init_Ttol", 0.01); + + AirParcel air(AirParcel::Type::MoleFrac); + + // estimate surface temperature and pressure + Real Ts = T0 - grav / cp * x1min; + Real Ps = P0 * pow(Ts / T0, cp / Rd); + Real xH2S = pin->GetReal("problem", "qH2S.ppmv") / 1.E6; + + while (iter++ < max_iter) { + // read in vapors + air.w[iH2S] = xH2S; + air.w[IPR] = Ps; + air.w[IDN] = Ts; + + // stop at just above P0 + for (int i = is; i <= ie; ++i) { + pthermo->Extrapolate(&air, pcoord->dx1f(i), "pseudo", grav); + if (air.w[IPR] < P0) break; + } + + // make up for the difference + Ts += T0 - air.w[IDN]; + if (std::abs(T0 - air.w[IDN]) < Ttol) break; + + app->Log("Iteration #", iter); + app->Log("T", air.w[IDN]); + } + + if (iter > max_iter) { + throw RuntimeError("ProblemGenerator", "maximum iteration reached"); + } + + // construct atmosphere from bottom up + for (int k = ks; k <= ke; ++k) + for (int j = js; j <= je; ++j) { + air.SetZero(); + air.w[iH2S] = xH2S; + air.w[IPR] = Ps; + air.w[IDN] = Ts; + + // half a grid to cell center + pthermo->Extrapolate(&air, pcoord->dx1f(is) / 2., "reversible", grav); + + int i = is; + for (; i <= ie; ++i) { + if (air.w[IDN] < Tmin) break; + AirParcelHelper::distribute_to_conserved(this, k, j, i, air); + //pthermo->Extrapolate(&air, pcoord->dx1f(i), "pseudo", grav, 1.e-5); + pthermo->Extrapolate(&air, pcoord->dx1f(i), "pseudo", grav); + } + + // Replace adiabatic atmosphere with isothermal atmosphere if temperature + // is too low + for (; i <= ie; ++i) { + AirParcelHelper::distribute_to_conserved(this, k, j, i, air); + pthermo->Extrapolate(&air, pcoord->dx1f(i), "isothermal", grav); + } + } +} diff --git a/examples/2024-XZhang-cloud-rt/jupiter_crm.cpp b/examples/2024-XZhang-cloud-rt/jupiter_crm.cpp new file mode 100644 index 00000000..020a077e --- /dev/null +++ b/examples/2024-XZhang-cloud-rt/jupiter_crm.cpp @@ -0,0 +1,279 @@ +//======================================================================================== +// Athena++ astrophysical MHD code +// Copyright(C) 2014 James M. Stone and other code +// contributors Licensed under the 3-clause BSD License, see LICENSE file for +// details +//======================================================================================== +//! \file hs94.cpp +// \brief Problem generator for Held-Suarez-94 GCM bench mark. +// +// REFERENCE: I.M Held & M.J Suarez, "A Proposal for the Intercomparison of the +// Dynamical Cores of Atmospheric General Circulation Models" + +// C++ headers +#include +#include +#include +#include +#include + +// athena +#include +#include +#include +#include +#include +#include +#include +#include +#include + +// canoe +#include +#include +#include +#include + +// exo3 +#include +#include + +// snap +#include + +// application +#include +#include + +// climath +#include +#include + +// special includes +#include + + +#define _sqr(x) ((x) * (x)) +#define _qur(x) ((x) * (x) * (x) * (x)) + +using namespace std; + +Real grav, P0, T0, Tmin, prad, hrate, Omega, Rp; +int iH2O, iNH3; + +std::default_random_engine generator; +std::normal_distribution distribution(0.0, 1.0); + +namespace cs = CubedSphereUtility; + +void Forcing(MeshBlock *pmb, Real const time, Real const dt, + AthenaArray const &w, AthenaArray const &r, + AthenaArray const &bcc, AthenaArray &du, + AthenaArray &s) { + int is = pmb->is, js = pmb->js, ks = pmb->ks; + int ie = pmb->ie, je = pmb->je, ke = pmb->ke; + auto pthermo = Thermodynamics::GetInstance(); + + auto pexo3 = pmb->pimpl->pexo3; + Real om_earth = Omega; + + for (int k = pmb->ks; k <= pmb->ke; ++k) + for (int j = pmb->js; j <= pmb->je; ++j) + for (int i = pmb->is; i <= pmb->ie; ++i) { + } + + + /* Sponge Layer + for (int k = pmb->ks; k <= pmb->ke; ++k) { + for (int j = pmb->js; j <= pmb->je; ++j) { + for (int i = pmb->is; i <= pmb->ie; ++i) { + Real pres = w(IPR, k, j, i); + if (pres < piso) { // sponge layer at top + Real tau = sponge_tau * pow(pres / piso, 2); + u(IVX, k, j, i) = u(IVX, k, j, i) / (1 + dt / tau); + u(IVY, k, j, i) = u(IVY, k, j, i) / (1 + dt / tau); + u(IVZ, k, j, i) = u(IVZ, k, j, i) / (1 + dt / tau); + } + } + } + }*/ + for (int k = ks; k <= ke; ++k) + for (int j = js; j <= je; ++j) + for (int i = is; i <= ie; ++i) { + // coriolis force + Real lat, lon; + pexo3->GetLatLon(&lat, &lon, k, j, i); + // coriolis force + Real f = 2. * om_earth * sin(lat); + Real f2 = 2. * om_earth * cos(lat); + Real U, V; + + pexo3->GetUV(&U, &V, w(IVY, k, j, i), w(IVZ, k, j, i), k, j, i); + + Real m1 = w(IDN, k, j, i) * w(IVX, k, j, i); + Real m2 = w(IDN, k, j, i) * U; + Real m3 = w(IDN, k, j, i) * V; + + // du(IM1, k, j, i) += dt * f * m2; + Real ll_acc_U = f * m3; //- f2 * m1; + Real ll_acc_V = -f * m2; + Real acc1, acc2, acc3; + pexo3->GetVyVz(&acc2, &acc3, ll_acc_U, ll_acc_V, k, j, i); + pexo3->ContravariantVectorToCovariant(j, k, acc2, acc3, &acc2, &acc3); + du(IM2, k, j, i) += dt * acc2; + du(IM3, k, j, i) += dt * acc3; + } +} + +void MeshBlock::InitUserMeshBlockData(ParameterInput *pin) { + AllocateUserOutputVariables(8 + NVAPOR); + SetUserOutputVariableName(0, "temp"); + SetUserOutputVariableName(1, "theta"); + SetUserOutputVariableName(2, "thetav"); + SetUserOutputVariableName(3, "mse"); + SetUserOutputVariableName(4, "lat"); + SetUserOutputVariableName(5, "lon"); + SetUserOutputVariableName(6, "U"); + SetUserOutputVariableName(7, "V"); + for (int n = 1; n <= NVAPOR; ++n) { + std::string name = "rh" + std::to_string(n); + SetUserOutputVariableName(7 + n, name.c_str()); + } +} + +void MeshBlock::UserWorkBeforeOutput(ParameterInput *pin) { + auto pthermo = Thermodynamics::GetInstance(); + + for (int k = ks; k <= ke; ++k) + for (int j = js; j <= je; ++j) + for (int i = is; i <= ie; ++i) { + user_out_var(0, k, j, i) = pthermo->GetTemp(this, k, j, i); + user_out_var(1, k, j, i) = pthermo->PotentialTemp(this, P0, k, j, i); + // theta_v + user_out_var(2, k, j, i) = + user_out_var(1, k, j, i) * pthermo->RovRd(this, k, j, i); + // mse + user_out_var(3, k, j, i) = + pthermo->MoistStaticEnergy(this, grav * (pcoord->x1v(i)-Rp), k, j, i); + for (int n = 1; n <= NVAPOR; ++n) + user_out_var(7 + n, k, j, i) = + pthermo->RelativeHumidity(this, n, k, j, i); + Real x = tan(pcoord->x2v(j)); + + Real Vy = phydro->w(IVY, k, j, i); + Real Vz = phydro->w(IVZ, k, j, i); + Real lat, lon; + pimpl->pexo3->GetLatLon(&lat, &lon, k, j, i); + user_out_var(4, k, j, i) = lat / PI * 180.0; + user_out_var(5, k, j, i) = lon / PI * 180.0; + Real U, V; + pimpl->pexo3->GetUV(&U, &V, Vy, Vz, k, j, i); + user_out_var(6, k, j, i) = U; + user_out_var(7, k, j, i) = V; + } +} + +void Mesh::InitUserMeshData(ParameterInput *pin) { + grav = -pin->GetReal("hydro", "grav_acc1"); + + P0 = pin->GetReal("problem", "P0"); + T0 = pin->GetReal("problem", "T0"); + Omega = pin->GetReal("problem", "Omega"); + Rp = pin->GetReal("problem", "Rp"); + + Tmin = pin->GetReal("problem", "Tmin"); + prad = pin->GetReal("problem", "prad"); + hrate = pin->GetReal("problem", "hrate") / 86400.; + + // index + auto pindex = IndexMap::GetInstance(); + iH2O = pindex->GetVaporId("H2O"); + iNH3 = pindex->GetVaporId("NH3"); + EnrollUserExplicitSourceFunction(Forcing); +} + +void MeshBlock::ProblemGenerator(ParameterInput *pin) { + srand(Globals::my_rank + time(0)); + + Application::Logger app("main"); + app->Log("ProblemGenerator: jupiter_crm"); + + auto pthermo = Thermodynamics::GetInstance(); + + // mesh limits + Real x1min = pmy_mesh->mesh_size.x1min; + Real x1max = pmy_mesh->mesh_size.x1max; + + // request temperature and pressure + app->Log("request T", T0); + app->Log("request P", P0); + + // thermodynamic constants + Real gamma = pin->GetReal("hydro", "gamma"); + Real Rd = pthermo->GetRd(); + Real cp = gamma / (gamma - 1.) * Rd; + + // set up an adiabatic atmosphere + int max_iter = 400, iter = 0; + Real Ttol = pin->GetOrAddReal("problem", "init_Ttol", 0.01); + + AirParcel air(AirParcel::Type::MoleFrac); + + // estimate surface temperature and pressure + Real Ts = T0 - grav / cp * (x1min - Rp); + Real Ps = P0 * pow(Ts / T0, cp / Rd); + Real xH2O = pin->GetReal("problem", "qH2O.ppmv") / 1.E6; + Real xNH3 = pin->GetReal("problem", "qNH3.ppmv") / 1.E6; + + while (iter++ < max_iter) { + // read in vapors + air.w[iH2O] = xH2O; + air.w[iNH3] = xNH3; + air.w[IPR] = Ps; + air.w[IDN] = Ts; + + // stop at just above P0 + for (int i = is; i <= ie; ++i) { + pthermo->Extrapolate(&air, pcoord->dx1f(i), "pseudo", grav); + if (air.w[IPR] < P0) break; + } + + // make up for the difference + Ts += T0 - air.w[IDN]; + if (std::abs(T0 - air.w[IDN]) < Ttol) break; + + app->Log("Iteration #", iter); + app->Log("T", air.w[IDN]); + } + + if (iter > max_iter) { + throw RuntimeError("ProblemGenerator", "maximum iteration reached"); + } + + // construct atmosphere from bottom up + for (int k = ks; k <= ke; ++k) + for (int j = js; j <= je; ++j) { + air.SetZero(); + air.w[iH2O] = xH2O; + air.w[iNH3] = xNH3; + air.w[IPR] = Ps; + air.w[IDN] = Ts; + + // half a grid to cell center + pthermo->Extrapolate(&air, pcoord->dx1f(is) / 2., "reversible", grav); + + int i = is; + for (; i <= ie; ++i) { + if (air.w[IDN] < Tmin) break; + AirParcelHelper::distribute_to_conserved(this, k, j, i, air); + pthermo->Extrapolate(&air, pcoord->dx1f(i), "pseudo", grav, 1.e-5); + } + + // Replace adiabatic atmosphere with isothermal atmosphere if temperature + // is too low + for (; i <= ie; ++i) { + AirParcelHelper::distribute_to_conserved(this, k, j, i, air); + pthermo->Extrapolate(&air, pcoord->dx1f(i), "isothermal", grav); + } + } +} diff --git a/examples/2024-XZhang-cloud-rt/rtcldu3d.cpp b/examples/2024-XZhang-cloud-rt/rtcldu3d.cpp new file mode 100644 index 00000000..2c102735 --- /dev/null +++ b/examples/2024-XZhang-cloud-rt/rtcldu3d.cpp @@ -0,0 +1,288 @@ +// C++ headers +#include +#include +#include +#include +#include + +// athena +#include +#include +#include +#include +#include +#include +#include +#include + +// canoe +#include +#include +#include +#include + +// exo3 +#include +#include + +// snap +#include + +// application +#include +#include + +// climath +#include +#include + +// special includes +#include "svp_functions.hpp" + +// astro +#include + +// harp +#include "harp/radiation.hpp" + +#define _sqr(x) ((x) * (x)) +#define _qur(x) ((x) * (x) * (x) * (x)) + +using namespace std; + +Real grav, P0, T0, Tmin, Omega, radius; +Real Omega_f, center_deg, beta_f; +int iH2S, iCH4; + +std::default_random_engine generator; +std::normal_distribution distribution(0.0, 1.0); + +namespace cs = CubedSphereUtility; + +void MeshBlock::InitUserMeshBlockData(ParameterInput *pin) { +#ifdef CUBED_SPHERE + AllocateUserOutputVariables(11 + NVAPOR); +#else + AllocateUserOutputVariables(5 + NVAPOR); +#endif + + SetUserOutputVariableName(0, "temp"); + SetUserOutputVariableName(1, "theta"); + SetUserOutputVariableName(2, "thetav"); + SetUserOutputVariableName(3, "mse"); + SetUserOutputVariableName(4, "pres"); + for (int n = 1; n <= NVAPOR; ++n) { + std::string name = "rh" + std::to_string(n); + SetUserOutputVariableName(4 + n, name.c_str()); + } + +#ifdef CUBED_SPHERE + SetUserOutputVariableName(5 + NVAPOR, "lat"); + SetUserOutputVariableName(6 + NVAPOR, "lon"); + SetUserOutputVariableName(7 + NVAPOR, "vlat"); + SetUserOutputVariableName(8 + NVAPOR, "vlon"); + SetUserOutputVariableName(9 + NVAPOR, "w"); + SetUserOutputVariableName(10 + NVAPOR, "zenith"); +#endif +} + +void MeshBlock::UserWorkBeforeOutput(ParameterInput *pin) { + auto pthermo = Thermodynamics::GetInstance(); + + for (int k = ks; k <= ke; ++k) + for (int j = js; j <= je; ++j) + for (int i = is; i <= ie; ++i) { + user_out_var(0, k, j, i) = pthermo->GetTemp(this, k, j, i); + user_out_var(1, k, j, i) = pthermo->PotentialTemp(this, P0, k, j, i); + // theta_v + user_out_var(2, k, j, i) = + user_out_var(1, k, j, i) * pthermo->RovRd(this, k, j, i); + // mse + user_out_var(3, k, j, i) = pthermo->MoistStaticEnergy( + this, grav * (pcoord->x1v(i) - radius), k, j, i); + user_out_var(4, k, j, i) = phydro->w(IPR, k, j, i); + for (int n = 1; n <= NVAPOR; ++n) + user_out_var(4 + n, k, j, i) = + pthermo->RelativeHumidity(this, n, k, j, i); + } + +#ifdef CUBED_SPHERE + auto pexo3 = pimpl->pexo3; + Real lat, lon; + Real U, V; + Direction ray = pimpl->prad->GetRayInput(0); + Real zenith; + + for (int k = ks; k <= ke; ++k) + for (int j = js; j <= je; ++j) + for (int i = is; i <= ie; ++i) { + pexo3->GetLatLon(&lat, &lon, k, j, i); + pexo3->GetUV(&U, &V, phydro->w(IVY, k, j, i), phydro->w(IVZ, k, j, i), + k, j, i); + user_out_var(5 + NVAPOR, k, j, i) = lat; + user_out_var(6 + NVAPOR, k, j, i) = lon; + user_out_var(7 + NVAPOR, k, j, i) = U; + user_out_var(8 + NVAPOR, k, j, i) = V; + user_out_var(9 + NVAPOR, k, j, i) = phydro->w(IVX, k, j, i); + + ray = pimpl->planet->ParentZenithAngle(pmy_mesh->time, M_PI / 2. - lat, + lon); + zenith = std::acos(ray.mu) / M_PI * 180.0; + user_out_var(10 + NVAPOR, k, j, i) = zenith; + } +#endif +} + +void Forcing(MeshBlock *pmb, Real const time, Real const dt, + AthenaArray const &w, AthenaArray const &r, + AthenaArray const &bcc, AthenaArray &du, + AthenaArray &s) { + int is = pmb->is, js = pmb->js, ks = pmb->ks; + int ie = pmb->ie, je = pmb->je, ke = pmb->ke; + auto pthermo = Thermodynamics::GetInstance(); + +#ifndef CUBED_SPHERE + for (int k = pmb->ks; k <= pmb->ke; ++k) + for (int j = pmb->js; j <= pmb->je; ++j) + for (int i = pmb->is; i <= pmb->ie; ++i) { + Real dx2 = pmb->pcoord->dx2v(j); + Real dx3 = pmb->pcoord->dx3v(k); + // Coriolis force + Real x2 = pmb->pcoord->x2v(j); + // Load Latitude (constant f here) + Real theta0 = center_deg * M_PI / 180.; + // Coriolis force (shallow approximation) + Real local_f = 2. * sin(theta0) * Omega_f + x2 * beta_f; + + Real m2 = w(IDN, k, j, i) * w(IVY, k, j, i); + Real m3 = w(IDN, k, j, i) * w(IVZ, k, j, i); + du(IM2, k, j, i) -= dt * local_f * m3; + du(IM3, k, j, i) += dt * local_f * m2; + } +#endif + +#ifdef CUBED_SPHERE + auto pexo3 = pmb->pimpl->pexo3; + for (int k = pmb->ks; k <= pmb->ke; ++k) + for (int j = pmb->js; j <= pmb->je; ++j) + for (int i = pmb->is; i <= pmb->ie; ++i) { + Real lat, lon; + pexo3->GetLatLon(&lat, &lon, k, j, i); + + // coriolis force + Real f = 2. * Omega * sin(lat); + Real f2 = 2. * Omega * cos(lat); + Real U, V; + + pexo3->GetUV(&U, &V, w(IVY, k, j, i), w(IVZ, k, j, i), k, j, i); + + Real m1 = w(IDN, k, j, i) * w(IVX, k, j, i); + Real m2 = w(IDN, k, j, i) * U; + Real m3 = w(IDN, k, j, i) * V; + + Real ll_acc_U = f * m3; + Real ll_acc_V = -f * m2; + Real acc2, acc3; + pexo3->GetVyVz(&acc2, &acc3, ll_acc_U, ll_acc_V, k, j, i); + pexo3->ContravariantVectorToCovariant(j, k, acc2, acc3, &acc2, &acc3); + du(IM2, k, j, i) += dt * acc2; + du(IM3, k, j, i) += dt * acc3; + } +#endif +} + +void Mesh::InitUserMeshData(ParameterInput *pin) { + grav = -pin->GetReal("hydro", "grav_acc1"); + P0 = pin->GetReal("problem", "P0"); + T0 = pin->GetReal("problem", "T0"); + Omega = pin->GetReal("problem", "Omega"); + radius = pin->GetReal("problem", "radius"); + Tmin = pin->GetReal("problem", "Tmin"); + center_deg = pin->GetReal("problem", "center_deg"); + beta_f = pin->GetReal("problem", "beta_f"); + Omega_f = pin->GetReal("problem", "Omega_f"); + + // index + auto pindex = IndexMap::GetInstance(); + iH2S = pindex->GetVaporId("H2S"); + iCH4 = pindex->GetVaporId("CH4"); + EnrollUserExplicitSourceFunction(Forcing); +} + +void MeshBlock::ProblemGenerator(ParameterInput *pin) { + srand(Globals::my_rank + time(0)); + auto pthermo = Thermodynamics::GetInstance(); + AirParcel air(AirParcel::Type::MoleFrac); + + // mesh limits + Real x1min = pmy_mesh->mesh_size.x1min; + Real x1max = pmy_mesh->mesh_size.x1max; + + // thermodynamic constants + Real gammad = pin->GetReal("hydro", "gamma"); + Real Rd = pthermo->GetRd(); + Real cp = gammad / (gammad - 1.) * Rd; + + // set up an adiabatic atmosphere + int max_iter = 400, iter = 0; + Real Ttol = pin->GetOrAddReal("problem", "init_Ttol", 0.01); + + // estimate surface temperature and pressure + Real Ts = T0 - grav / cp * (x1min - radius); + Real Ps = P0 * pow(Ts / T0, cp / Rd); + Real xH2S = pin->GetReal("problem", "qH2S.ppmv") / 1.E6; + Real xCH4 = pin->GetReal("problem", "qCH4.ppmv") / 1.E6; + + while (iter++ < max_iter) { + // read in vapors + air.w[iH2S] = xH2S; + air.w[iCH4] = xCH4; + air.w[IPR] = Ps; + air.w[IDN] = Ts; + + // stop at just above P0 + for (int i = is; i <= ie; ++i) { + pthermo->Extrapolate(&air, pcoord->dx1f(i), "pseudo", grav); + //pthermo->Extrapolate(&air, pcoord->dx1f(i), "reversible", grav); + if (air.w[IPR] < P0) break; + } + + // make up for the difference + Ts += T0 - air.w[IDN]; + if (std::abs(T0 - air.w[IDN]) < Ttol) break; + } + + if (iter > max_iter) { + throw RuntimeError("ProblemGenerator", "maximum iteration reached"); + } + + // construct atmosphere from bottom up + for (int k = ks; k <= ke; ++k) + for (int j = js; j <= je; ++j) { + air.SetZero(); + air.w[iH2S] = xH2S; + air.w[iCH4] = xCH4; + air.w[IPR] = Ps; + air.w[IDN] = Ts; + + // half a grid to cell center + pthermo->Extrapolate(&air, pcoord->dx1f(is) / 2., "reversible", grav); + + int i = is; + for (; i <= ie; ++i) { + air.w[IVX] = 0.001 * (1. * rand() / RAND_MAX - 0.5); + if (air.w[IDN] < Tmin) break; + AirParcelHelper::distribute_to_conserved(this, k, j, i, air); + pthermo->Extrapolate(&air, pcoord->dx1f(i), "pseudo", grav); + //pthermo->Extrapolate(&air, pcoord->dx1f(i), "reversible", grav); + } + + // Replace adiabatic atmosphere with isothermal atmosphere if temperature + // is too low + for (; i <= ie; ++i) { + air.w[IVX] = 0.001 * (1. * rand() / RAND_MAX - 0.5); + AirParcelHelper::distribute_to_conserved(this, k, j, i, air); + pthermo->Extrapolate(&air, pcoord->dx1f(i), "isothermal", grav); + } + } +} diff --git a/examples/2024-XZhang-cloud-rt/svp_functions.hpp b/examples/2024-XZhang-cloud-rt/svp_functions.hpp index 3652e7da..e90a6c3d 100644 --- a/examples/2024-XZhang-cloud-rt/svp_functions.hpp +++ b/examples/2024-XZhang-cloud-rt/svp_functions.hpp @@ -137,10 +137,63 @@ void enroll_vapor_Function_fe(Thermodynamics::SVPFunc1Container &svp_func1, } } + +// CH4 svp +double sat_vapor_p_CH4(double T) { + double betal = 10.15, gammal = 1.87, tr = 90.67, + pr = 11690.; + return SatVaporPresIdeal(T / tr, pr, betal, gammal); +} + +void enroll_vapor_Function_CH4(Thermodynamics::SVPFunc1Container &svp_func1, + std::vector const &cloud_index_set) { + Application::Logger app("snap"); + app->Log("Enrolling CH4 vapor pressures"); + + auto pindex = IndexMap::GetInstance(); + if (!pindex->HasVapor("CH4")) return; + int iCH4 = pindex->GetVaporId("CH4"); + + svp_func1[iCH4][0] = [](AirParcel const &qfrac, int, int) { + return sat_vapor_p_CH4(qfrac.w[IDN]); + }; + + for (int n = 1; n < cloud_index_set[iCH4].size(); ++n) { + svp_func1[iCH4][n] = NullSatVaporPres1; + } +} + +// H2S svp +double sat_vapor_p_H2S(double T) { + double betal = 23.67, gammal = 3.96, tr = 187.7, + pr = 23200.; + return SatVaporPresIdeal(T / tr, pr, betal, gammal); +} + +void enroll_vapor_Function_H2S(Thermodynamics::SVPFunc1Container &svp_func1, + std::vector const &cloud_index_set) { + Application::Logger app("snap"); + app->Log("Enrolling H2S vapor pressures"); + + auto pindex = IndexMap::GetInstance(); + if (!pindex->HasVapor("H2S")) return; + int iH2S = pindex->GetVaporId("H2S"); + + svp_func1[iH2S][0] = [](AirParcel const &qfrac, int, int) { + return sat_vapor_p_H2S(qfrac.w[IDN]); + }; + + for (int n = 1; n < cloud_index_set[iH2S].size(); ++n) { + svp_func1[iH2S][n] = NullSatVaporPres1; + } +} + void Thermodynamics::enrollVaporFunctions() { enroll_vapor_Function_H2O(svp_func1_, cloud_index_set_); enroll_vapor_Function_NH3(svp_func1_, cloud_index_set_); enroll_vapor_Function_NH4SH(svp_func1_, cloud_index_set_); enroll_vapor_Function_mgsio3(svp_func1_, cloud_index_set_); enroll_vapor_Function_fe(svp_func1_, cloud_index_set_); + enroll_vapor_Function_H2S(svp_func1_, cloud_index_set_); + enroll_vapor_Function_CH4(svp_func1_, cloud_index_set_); } diff --git a/examples/2024-XZhang-minichem/jup_onetra.cpp b/examples/2024-XZhang-minichem/jup_onetra.cpp new file mode 100644 index 00000000..9cb472e5 --- /dev/null +++ b/examples/2024-XZhang-minichem/jup_onetra.cpp @@ -0,0 +1,290 @@ +/* ------------------------------------------------------------------------------------- + * SNAP Example Program + * + * Contributer: + * Cheng Li, University of Michigan + * + * Year: 2023 + * Contact: chengcli@umich.edu + * Reference: Test Jupiter CRM + * ------------------------------------------------------------------------------------- + */ +// C++ headers +#include +#include +#include +#include +#include + +// athena +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +// application +#include +#include + +// canoe +#include +#include + +// climath +#include +#include + +// snap +#include +#include + +// special includes +#include + +// minichem +//#include + +Real grav, P0, T0, Tmin, prad, hrate; +int iH2O; +//'OH','H2','H2O','H','CO','CO2','O','CH4','C2H2','NH3','N2','HCN', 'He' +std::vector vmass = {17.01, 2.02, 18.02, 1.01, 28.01, 44.01, 16., + 16.05, 26.04, 17.04, 28.02, 27.03, 4.}; +//MiniChem *mc; + +void MeshBlock::InitUserMeshBlockData(ParameterInput *pin) { + AllocateUserOutputVariables(5 + NVAPOR); + SetUserOutputVariableName(0, "temp"); + SetUserOutputVariableName(1, "theta"); + SetUserOutputVariableName(2, "thetav"); + SetUserOutputVariableName(3, "mse"); + SetUserOutputVariableName(4, "pres"); + for (int n = 1; n <= NVAPOR; ++n) { + std::string name = "rh" + std::to_string(n); + SetUserOutputVariableName(4 + n, name.c_str()); + } +} + +void MeshBlock::UserWorkBeforeOutput(ParameterInput *pin) { + auto pthermo = Thermodynamics::GetInstance(); + + for (int k = ks; k <= ke; ++k) + for (int j = js; j <= je; ++j) + for (int i = is; i <= ie; ++i) { + user_out_var(0, k, j, i) = pthermo->GetTemp(this, k, j, i); + user_out_var(1, k, j, i) = pthermo->PotentialTemp(this, P0, k, j, i); + // theta_v + user_out_var(2, k, j, i) = + user_out_var(1, k, j, i) * pthermo->RovRd(this, k, j, i); + // mse + user_out_var(3, k, j, i) = + pthermo->MoistStaticEnergy(this, grav * pcoord->x1v(i), k, j, i); + user_out_var(4, k, j, i) = phydro->w(IPR, k, j, i); + // relative humidity + for (int n = 1; n <= NVAPOR; ++n) + user_out_var(4 + n, k, j, i) = + pthermo->RelativeHumidity(this, n, k, j, i); + } +} + +void Forcing(MeshBlock *pmb, Real const time, Real const dt, + AthenaArray const &w, const AthenaArray &prim_scalar, + AthenaArray const &bcc, AthenaArray &du, + AthenaArray &cons_scalar) { + int is = pmb->is, js = pmb->js, ks = pmb->ks; + int ie = pmb->ie, je = pmb->je, ke = pmb->ke; + auto pthermo = Thermodynamics::GetInstance(); + auto phydro = pmb->phydro; + std::vector vmr(NCHEMISTRY); + std::vector vmr_sp(NCHEMISTRY - 1); + std::vector mmr(NCHEMISTRY); + Real sumVMR; + + for (int k = ks; k <= ke; ++k) + for (int j = js; j <= je; ++j) + for (int i = is; i <= ie; ++i) { + auto &&air = AirParcelHelper::gather_from_primitive(pmb, k, j, i); + + air.ToMoleFraction(); + Real cv = + get_cv_mass(air, 0, pthermo->GetRd(), pthermo->GetCvRatioMass()); + + if (w(IPR, k, j, i) < prad) { + du(IEN, k, j, i) += dt * hrate * w(IDN, k, j, i) * cv * + (1. + 1.E-4 * sin(2. * M_PI * rand() / RAND_MAX)); + } +/* + // minichem + Real temp = pthermo->GetTemp(w.at(k, j, i)); + Real pres = w(IPR, k, j, i); + + // change mmr to vmr and normalize + for (int n = 0; n < NCHEMISTRY; ++n) + vmr[n] = prim_scalar(NCLOUD + n, k, j, i) / vmass[n]; + sumVMR = std::accumulate(vmr.begin(), vmr.end(), static_cast(0)); + for (auto &value : vmr) value /= sumVMR; + + // call minichem for sp species + for (int n = 0; n < NCHEMISTRY - 1; ++n) vmr_sp[n] = vmr[n]; +// if (temp > 200.) mc->Run(temp, pres, dt, vmr_sp.data(), "NCHO"); + for (int n = 0; n < NCHEMISTRY - 1; ++n) vmr[n] = vmr_sp[n]; + + // change vmr to mmr and normalize + for (int n = 0; n < NCHEMISTRY; ++n) mmr[n] = vmr[n] * vmass[n]; + sumVMR = std::accumulate(mmr.begin(), mmr.end(), static_cast(0)); + for (auto &value : mmr) value /= sumVMR; + + for (int n = 0; n < NCHEMISTRY; ++n) { + cons_scalar(NCLOUD + n, k, j, i) += + phydro->w(IDN, k, j, i) * + (mmr[n] - prim_scalar(NCLOUD + n, k, j, i)); + } +*/ + + } +} + +void Mesh::InitUserMeshData(ParameterInput *pin) { + grav = -pin->GetReal("hydro", "grav_acc1"); + + P0 = pin->GetReal("problem", "P0"); + T0 = pin->GetReal("problem", "T0"); + + Tmin = pin->GetReal("problem", "Tmin"); + prad = pin->GetReal("problem", "prad"); + hrate = pin->GetReal("problem", "hrate") / 86400.; + + // index + auto pindex = IndexMap::GetInstance(); + iH2O = pindex->GetVaporId("H2O"); +// EnrollUserExplicitSourceFunction(Forcing); + + // minichem +/* mc = new MiniChem(); + mc->SetDataFile("chem_data/mini_chem_data_NCHO.txt"); + mc->SetSpeciesFile("chem_data/mini_chem_sp_NCHO.txt"); + mc->SetNetworkDir("chem_data/1x/"); + mc->SetMetallicityStr("1x"); + mc->Initialize(); +*/ +} + +void MeshBlock::ProblemGenerator(ParameterInput *pin) { + srand(Globals::my_rank + time(0)); + + Application::Logger app("main"); + app->Log("ProblemGenerator: jupiter_crm"); + + auto pthermo = Thermodynamics::GetInstance(); + + // mesh limits + Real x1min = pmy_mesh->mesh_size.x1min; + Real x1max = pmy_mesh->mesh_size.x1max; + + // request temperature and pressure + app->Log("request T", T0); + app->Log("request P", P0); + + // thermodynamic constants + Real gamma = pin->GetReal("hydro", "gamma"); + Real Rd = pthermo->GetRd(); + Real cp = gamma / (gamma - 1.) * Rd; + + // set up an adiabatic atmosphere + int max_iter = 400, iter = 0; + Real Ttol = pin->GetOrAddReal("problem", "init_Ttol", 0.01); + + AirParcel air(AirParcel::Type::MoleFrac); + + // estimate surface temperature and pressure + Real Ts = T0 - grav / cp * x1min; + Real Ps = P0 * pow(Ts / T0, cp / Rd); + Real xH2O = pin->GetReal("problem", "qH2O.ppmv") / 1.E6; + Real xtra = pin->GetReal("problem", "tracer.ppb") / 1.E9; + + while (iter++ < max_iter) { + // read in vapors + air.w[iH2O] = xH2O; + air.w[IPR] = Ps; + air.w[IDN] = Ts; + + // stop at just above P0 + for (int i = is; i <= ie; ++i) { + pthermo->Extrapolate(&air, pcoord->dx1f(i), "pseudo", grav); + if (air.w[IPR] < P0) break; + } + + // make up for the difference + Ts += T0 - air.w[IDN]; + if (std::abs(T0 - air.w[IDN]) < Ttol) break; + + app->Log("Iteration #", iter); + app->Log("T", air.w[IDN]); + } + + if (iter > max_iter) { + throw RuntimeError("ProblemGenerator", "maximum iteration reached"); + } + + // construct atmosphere from bottom up + for (int k = ks; k <= ke; ++k) + for (int j = js; j <= je; ++j) { + air.SetZero(); + air.w[iH2O] = xH2O; + air.w[IPR] = Ps; + air.w[IDN] = Ts; + + // half a grid to cell center + pthermo->Extrapolate(&air, pcoord->dx1f(is) / 2., "reversible", grav); + + int i = is; + for (; i <= ie; ++i) { + air.w[IVX] = 0.001 * (1. * rand() / RAND_MAX - 0.5); + if (air.w[IDN] < Tmin) break; + AirParcelHelper::distribute_to_primitive(this, k, j, i, air); + AirParcelHelper::distribute_to_conserved(this, k, j, i, air); + pthermo->Extrapolate(&air, pcoord->dx1f(i), "pseudo", grav, 1.e-5); + } + + // Replace adiabatic atmosphere with isothermal atmosphere if temperature + // is too low + for (; i <= ie; ++i) { + air.w[IVX] = 0.001 * (1. * rand() / RAND_MAX - 0.5); + AirParcelHelper::distribute_to_primitive(this, k, j, i, air); + AirParcelHelper::distribute_to_conserved(this, k, j, i, air); + pthermo->Extrapolate(&air, pcoord->dx1f(i), "isothermal", grav); + } + } + + // minichem initialize conserved variables + // interpolate from ce table + std::vector vmr_ic(NCHEMISTRY); + std::vector mmr(NCHEMISTRY); + std::string ic_file = "chem_data/IC/mini_chem_IC_FastChem_1x.txt"; + + for (int k = ks; k <= ke; k++) + for (int j = js; j <= je; j++) + for (int i = is; i <= ie; i++) { + double mu; + double T_in = pthermo->GetTemp(this, k, j, i); + double P_in = phydro->w(IPR, k, j, i); +// interp_ce_table(NCHEMISTRY, T_in, P_in, vmr_ic.data(), &mu, ic_file); + + // change vmr to mmr and normalize +// for (int n = 0; n < NCHEMISTRY; ++n) mmr[n] = vmr_ic[n] * vmass[n]; +// Real sumVMR = +// std::accumulate(mmr.begin(), mmr.end(), static_cast(0)); +// for (auto &value : mmr) value /= sumVMR; + if (P_in > 6.e7) { + for (int n = 0; n < NCHEMISTRY; ++n) { + mmr[n] = xtra; + pscalars->s(NCLOUD + n, k, j, i) += phydro->w(IDN, k, j, i) * mmr[n]; + } + } + } +} diff --git a/examples/CMakeLists.txt b/examples/CMakeLists.txt index 2430bf57..8346a142 100644 --- a/examples/CMakeLists.txt +++ b/examples/CMakeLists.txt @@ -42,6 +42,7 @@ endif() if (${TASK} STREQUAL "rtdry" OR ${TASK} STREQUAL "rtcld" OR + ${TASK} STREQUAL "ice_h2s" OR ${TASK} STREQUAL "rtcldexo") add_subdirectory(2024-XZhang-cloud-rt) endif() From 0626af61f3d031caaf044e64edd52b6b24e58203 Mon Sep 17 00:00:00 2001 From: mac/cli Date: Mon, 11 Aug 2025 12:10:28 -0400 Subject: [PATCH 03/10] minichem pgen --- examples/2024-XZhang-minichem/minichem.cpp | 357 +++++++++------------ src/CMakeLists.txt | 62 ++-- src/interface/eos.hpp | 15 + 3 files changed, 195 insertions(+), 239 deletions(-) diff --git a/examples/2024-XZhang-minichem/minichem.cpp b/examples/2024-XZhang-minichem/minichem.cpp index bf2683d7..28aa8049 100644 --- a/examples/2024-XZhang-minichem/minichem.cpp +++ b/examples/2024-XZhang-minichem/minichem.cpp @@ -1,14 +1,3 @@ -/* ------------------------------------------------------------------------------------- - * SNAP Example Program - * - * Contributer: - * Cheng Li, University of Michigan - * - * Year: 2023 - * Contact: chengcli@umich.edu - * Reference: Test Jupiter CRM - * ------------------------------------------------------------------------------------- - */ // C++ headers #include #include @@ -33,46 +22,45 @@ #include #include +// snap +#include + // canoe -#include -#include -#include +#include + #include +#include // climath #include #include -// snap -#include -#include - -// special includes -#include - // minichem #include Real grav, P0, T0, Tmin, prad, hrate, xH2O; Real Tbot, xtra, ptra, Ttop_tau, Tbot_tau, Trabot_tau; std::string met; -int iH2O; bool use_mini, use_mini_ic, use_tra_ic, fix_bot_tra; +Real p0 = 1.e5; + //'OH','H2','H2O','H','CO','CO2','O','CH4','C2H2','NH3','N2','HCN', 'He' std::vector vmass = {17.01, 2.02, 18.02, 1.01, 28.01, 44.01, 16., 16.05, 26.04, 17.04, 28.02, 27.03, 4.}; MiniChem *mc; void MeshBlock::InitUserMeshBlockData(ParameterInput *pin) { - if (NVAPOR > 0) { - AllocateUserOutputVariables(5 + NVAPOR); + int nvapor = pimpl->peos->pthermo->options.vapor_ids().size() - 1; + + if (nvapor > 0) { + AllocateUserOutputVariables(5 + nvapor); SetUserOutputVariableName(0, "temp"); SetUserOutputVariableName(1, "theta"); SetUserOutputVariableName(2, "pres"); SetUserOutputVariableName(3, "thetav"); SetUserOutputVariableName(4, "mse"); - for (int n = 1; n <= NVAPOR; ++n) { + for (int n = 1; n <= nvapor; ++n) { std::string name = "rh" + std::to_string(n); SetUserOutputVariableName(4 + n, name.c_str()); } @@ -85,16 +73,24 @@ void MeshBlock::InitUserMeshBlockData(ParameterInput *pin) { } void MeshBlock::UserWorkBeforeOutput(ParameterInput *pin) { - auto pthermo = Thermodynamics::GetInstance(); + auto w = phydro->w; + auto temp = get_temp(pimpl->peos, w); + auto pres = get_pres(w); + auto chi = (get_gammad() - 1.) / get_gammad(); + auto theta = temp * pow(p0 / pres, chi); + + auto temp_a = temp.accessor(); + auto theta_a = theta.accessor(); + auto pres_a = pres.accessor(); for (int k = ks; k <= ke; ++k) for (int j = js; j <= je; ++j) for (int i = is; i <= ie; ++i) { - user_out_var(0, k, j, i) = pthermo->GetTemp(this, k, j, i); - user_out_var(1, k, j, i) = pthermo->PotentialTemp(this, P0, k, j, i); - user_out_var(2, k, j, i) = phydro->w(IPR, k, j, i); + user_out_var(0, k, j, i) = temp_a[k][j][i]; + user_out_var(1, k, j, i) = theta_a[k][j][i]; + user_out_var(2, k, j, i) = pres_a[k][j][i]; - if (NVAPOR > 0) { + /*if (NVAPOR > 0) { // theta_v user_out_var(3, k, j, i) = user_out_var(1, k, j, i) * pthermo->RovRd(this, k, j, i); @@ -105,7 +101,7 @@ void MeshBlock::UserWorkBeforeOutput(ParameterInput *pin) { for (int n = 1; n <= NVAPOR; ++n) user_out_var(4 + n, k, j, i) = pthermo->RelativeHumidity(this, n, k, j, i); - } + }*/ } } @@ -115,82 +111,86 @@ void Forcing(MeshBlock *pmb, Real const time, Real const dt, AthenaArray &cons_scalar) { int is = pmb->is, js = pmb->js, ks = pmb->ks; int ie = pmb->ie, je = pmb->je, ke = pmb->ke; - auto pthermo = Thermodynamics::GetInstance(); auto phydro = pmb->phydro; + auto temp = get_temp(pmb->pimpl->peos, w); + auto cv = get_cv(pmb->pimpl->peos, w); + + auto temp_a = temp.accessor(); + auto cv_a = cv.accessor(); + for (int k = ks; k <= ke; ++k) for (int j = js; j <= je; ++j) for (int i = is; i <= ie; ++i) { - auto &&air = AirParcelHelper::gather_from_primitive(pmb, k, j, i); - - Real tem = pthermo->GetTemp(w.at(k,j,i)); - air.ToMoleFraction(); - Real cv = - get_cv_mass(air, 0, pthermo->GetRd(), pthermo->GetCvRatioMass()); - - if (w(IPR, k, j, i) < prad) { -// body cooling + auto tem = temp_a[k][j][i]; + auto cv = cv_a[k][j][i]; + if (w(IPR, k, j, i) < prad) { + // body cooling du(IEN, k, j, i) += dt * hrate * w(IDN, k, j, i) * cv * (1. + 1.E-4 * sin(2. * M_PI * rand() / RAND_MAX)); -// newtonian at the top - if (tem < Tmin) { - du(IEN, k, j, i) += dt * (Tmin - tem)/Ttop_tau * w(IDN, k, j, i) * cv; - } + // newtonian at the top + if (tem < Tmin) { + du(IEN, k, j, i) += + dt * (Tmin - tem) / Ttop_tau * w(IDN, k, j, i) * cv; + } } -// relax T at the bottom - if(i == is) du(IEN, k, j, i) += dt * (Tbot - tem)/Tbot_tau * w(IDN, k, j, i) * cv; - } - - if (NCHEMISTRY > 0) { + // relax T at the bottom + if (i == is) + du(IEN, k, j, i) += + dt * (Tbot - tem) / Tbot_tau * w(IDN, k, j, i) * cv; + } + + if (NCHEM > 0) { if (use_mini) { - std::vector vmr(NCHEMISTRY); - std::vector vmr_sp(NCHEMISTRY - 1); - std::vector mmr(NCHEMISTRY); - Real sumVMR; - for (int k = ks; k <= ke; ++k) - for (int j = js; j <= je; ++j) - for (int i = is; i <= ie; ++i) { - // minichem - Real temp = pthermo->GetTemp(w.at(k, j, i)); - Real pres = w(IPR, k, j, i); - - // change mmr to vmr and normalize - for (int n = 0; n < NCHEMISTRY; ++n) - vmr[n] = prim_scalar(NCLOUD + n, k, j, i) / vmass[n]; - sumVMR = - std::accumulate(vmr.begin(), vmr.end(), static_cast(0)); - for (auto &value : vmr) value /= sumVMR; - - // call minichem for sp species - for (int n = 0; n < NCHEMISTRY - 1; ++n) vmr_sp[n] = vmr[n]; - if (temp > 200.) mc->Run(temp, pres, dt, vmr_sp.data(), "NCHO"); - for (int n = 0; n < NCHEMISTRY - 1; ++n) vmr[n] = vmr_sp[n]; - - // change vmr to mmr and normalize - for (int n = 0; n < NCHEMISTRY; ++n) mmr[n] = vmr[n] * vmass[n]; - sumVMR = - std::accumulate(mmr.begin(), mmr.end(), static_cast(0)); - for (auto &value : mmr) value /= sumVMR; - - for (int n = 0; n < NCHEMISTRY; ++n) { - cons_scalar(NCLOUD + n, k, j, i) += - phydro->w(IDN, k, j, i) * - (mmr[n] - prim_scalar(NCLOUD + n, k, j, i)); + std::vector vmr(NCHEM); + std::vector vmr_sp(NCHEM - 1); + std::vector mmr(NCHEM); + Real sumVMR; + auto temp_a = temp.accessor(); + + for (int k = ks; k <= ke; ++k) + for (int j = js; j <= je; ++j) + for (int i = is; i <= ie; ++i) { + // minichem + Real temp = temp_a[k][j][i]; + Real pres = w(IPR, k, j, i); + + // change mmr to vmr and normalize + for (int n = 0; n < NCHEM; ++n) + vmr[n] = prim_scalar(n, k, j, i) / vmass[n]; + sumVMR = + std::accumulate(vmr.begin(), vmr.end(), static_cast(0)); + for (auto &value : vmr) value /= sumVMR; + + // call minichem for sp species + for (int n = 0; n < NCHEM - 1; ++n) vmr_sp[n] = vmr[n]; + if (temp > 200.) mc->Run(temp, pres, dt, vmr_sp.data(), "NCHO"); + for (int n = 0; n < NCHEM - 1; ++n) vmr[n] = vmr_sp[n]; + + // change vmr to mmr and normalize + for (int n = 0; n < NCHEM; ++n) mmr[n] = vmr[n] * vmass[n]; + sumVMR = + std::accumulate(mmr.begin(), mmr.end(), static_cast(0)); + for (auto &value : mmr) value /= sumVMR; + + for (int n = 0; n < NCHEM; ++n) { + cons_scalar(n, k, j, i) += + phydro->w(IDN, k, j, i) * (mmr[n] - prim_scalar(n, k, j, i)); + } } - } - } else if (fix_bot_tra) { - for (int k = ks; k <= ke; k++) - for (int j = js; j <= je; j++) - for (int i = is; i <= ie; i++) { - if (phydro->w(IPR, k, j, i) > ptra) { - for (int n = 0; n < NCHEMISTRY; ++n) { - cons_scalar(NCLOUD + n, k, j, i) += - dt/Trabot_tau * phydro->w(IDN, k, j, i) * - (xtra - prim_scalar(NCLOUD + n, k, j, i)); + } else if (fix_bot_tra) { + for (int k = ks; k <= ke; k++) + for (int j = js; j <= je; j++) + for (int i = is; i <= ie; i++) { + if (phydro->w(IPR, k, j, i) > ptra) { + for (int n = 0; n < NCHEM; ++n) { + cons_scalar(n, k, j, i) += dt / Trabot_tau * + phydro->w(IDN, k, j, i) * + (xtra - prim_scalar(n, k, j, i)); + } } } - } - } + } } } @@ -213,11 +213,6 @@ void Mesh::InitUserMeshData(ParameterInput *pin) { xtra = pin->GetReal("problem", "tracer.ppb") / 1.E9; ptra = pin->GetReal("problem", "tracer.pre"); - if (NVAPOR > 0) { - // index - auto pindex = IndexMap::GetInstance(); - iH2O = pindex->GetVaporId("H2O"); - } EnrollUserExplicitSourceFunction(Forcing); if (use_mini) { @@ -234,131 +229,77 @@ void Mesh::InitUserMeshData(ParameterInput *pin) { void MeshBlock::ProblemGenerator(ParameterInput *pin) { srand(Globals::my_rank + time(0)); - Application::Logger app("main"); - app->Log("ProblemGenerator: jupiter_crm"); - - auto pthermo = Thermodynamics::GetInstance(); - - // mesh limits - Real x1min = pmy_mesh->mesh_size.x1min; - Real x1max = pmy_mesh->mesh_size.x1max; - - // request temperature and pressure - app->Log("request T", T0); - app->Log("request P", P0); - // thermodynamic constants Real gamma = pin->GetReal("hydro", "gamma"); - Real Rd = pthermo->GetRd(); + auto mud = kintera::species_weights[0]; + auto Rd = kintera::constants::Rgas / mud; Real cp = gamma / (gamma - 1.) * Rd; // set up an adiabatic atmosphere int max_iter = 400, iter = 0; Real Ttol = pin->GetOrAddReal("problem", "init_Ttol", 0.01); - AirParcel air(AirParcel::Type::MoleFrac); - // estimate surface temperature and pressure + auto x1min = pcoord->x1f(is); Real Ts = T0 - grav / cp * x1min; Real Ps = P0 * pow(Ts / T0, cp / Rd); - if (NVAPOR > 0) xH2O = pin->GetReal("problem", "qH2O.ppmv") / 1.E6; - - while (iter++ < max_iter) { - // read in vapors - if (NVAPOR > 0) air.w[iH2O] = xH2O; - air.w[IPR] = Ps; - air.w[IDN] = Ts; - - // stop at just above P0 - for (int i = is; i <= ie; ++i) { - pthermo->Extrapolate(&air, pcoord->dx1f(i), "pseudo", grav); - if (air.w[IPR] < P0) break; - } - - // make up for the difference - Ts += T0 - air.w[IDN]; - if (std::abs(T0 - air.w[IDN]) < Ttol) break; - - app->Log("Iteration #", iter); - app->Log("T", air.w[IDN]); - } - - if (iter > max_iter) { - throw RuntimeError("ProblemGenerator", "maximum iteration reached"); - } - // construct atmosphere from bottom up + // Loop over the grids and set initial condition for (int k = ks; k <= ke; ++k) - for (int j = js; j <= je; ++j) { - air.SetZero(); - if (NVAPOR > 0) air.w[iH2O] = xH2O; - air.w[IPR] = Ps; - air.w[IDN] = Ts; - - // half a grid to cell center - pthermo->Extrapolate(&air, pcoord->dx1f(is) / 2., "reversible", grav); - - int i = is; - for (; i <= ie; ++i) { - air.w[IVX] = 0.001 * (1. * rand() / RAND_MAX - 0.5); - if (air.w[IDN] < Tmin) break; - AirParcelHelper::distribute_to_primitive(this, k, j, i, air); - AirParcelHelper::distribute_to_conserved(this, k, j, i, air); - //pthermo->Extrapolate(&air, pcoord->dx1f(i), "pseudo", grav, 1.e-5); - pthermo->Extrapolate(&air, pcoord->dx1f(i), "pseudo", grav, 0.); - } - - // Replace adiabatic atmosphere with isothermal atmosphere if temperature - // is too low - for (; i <= ie; ++i) { - air.w[IVX] = 0.001 * (1. * rand() / RAND_MAX - 0.5); - AirParcelHelper::distribute_to_primitive(this, k, j, i, air); - AirParcelHelper::distribute_to_conserved(this, k, j, i, air); - pthermo->Extrapolate(&air, pcoord->dx1f(i), "isothermal", grav); + for (int j = js; j <= je; ++j) + for (int i = is; i <= ie; ++i) { + Real x1 = pcoord->x1v(i); + Real x2 = pcoord->x2v(j); + Real x3 = pcoord->x3v(k); + Real temp = Ts - grav * x1 / cp; + phydro->w(IPR, k, j, i) = p0 * pow(temp / Ts, cp / Rd); + phydro->w(IDN, k, j, i) = phydro->w(IPR, k, j, i) / (Rd * temp); + phydro->w(IVX, k, j, i) = phydro->w(IVY, k, j, i) = 0.; } - } - Tbot = pthermo->GetTemp(this, ks, js, is);; - - if (NCHEMISTRY > 0) { - - if (use_mini_ic) { - // minichem initialize conserved variables - // interpolate from ce table - std::vector vmr_ic(NCHEMISTRY); - std::vector mmr(NCHEMISTRY); - std::string ic_file = "chem_data/IC/mini_chem_IC_FastChem_" + met + ".txt"; - - for (int k = ks; k <= ke; k++) - for (int j = js; j <= je; j++) - for (int i = is; i <= ie; i++) { - double mu; - double T_in = pthermo->GetTemp(this, k, j, i); - double P_in = phydro->w(IPR, k, j, i); - interp_ce_table(NCHEMISTRY, T_in, P_in, vmr_ic.data(), &mu, ic_file); - - // change vmr to mmr and normalize - for (int n = 0; n < NCHEMISTRY; ++n) mmr[n] = vmr_ic[n] * vmass[n]; - Real sumVMR = - std::accumulate(mmr.begin(), mmr.end(), static_cast(0)); - for (auto &value : mmr) value /= sumVMR; - - for (int n = 0; n < NCHEMISTRY; ++n) { - pscalars->s(NCLOUD + n, k, j, i) += - phydro->w(IDN, k, j, i) * mmr[n]; + auto temp = get_temp(pimpl->peos, phydro->w); + auto temp_a = temp.accessor(); + + if (NCHEM > 0) { + if (use_mini_ic) { + // minichem initialize conserved variables + // interpolate from ce table + std::vector vmr_ic(NCHEM); + std::vector mmr(NCHEM); + std::string ic_file = + "chem_data/IC/mini_chem_IC_FastChem_" + met + ".txt"; + + for (int k = ks; k <= ke; k++) + for (int j = js; j <= je; j++) + for (int i = is; i <= ie; i++) { + double mu; + double T_in = temp_a[k][j][i]; + double P_in = phydro->w(IPR, k, j, i); + interp_ce_table(NCHEM, T_in, P_in, vmr_ic.data(), &mu, ic_file); + + // change vmr to mmr and normalize + for (int n = 0; n < NCHEM; ++n) mmr[n] = vmr_ic[n] * vmass[n]; + Real sumVMR = + std::accumulate(mmr.begin(), mmr.end(), static_cast(0)); + for (auto &value : mmr) value /= sumVMR; + + for (int n = 0; n < NCHEM; ++n) { + pscalars->s(n, k, j, i) += phydro->w(IDN, k, j, i) * mmr[n]; + } } - } - } else if (use_tra_ic) { - for (int k = ks; k <= ke; k++) - for (int j = js; j <= je; j++) - for (int i = is; i <= ie; i++) { - if (phydro->w(IPR, k, j, i) > ptra) { - for (int n = 0; n < NCHEMISTRY; ++n) { - pscalars->s(NCLOUD + n, k, j, i) += - phydro->w(IDN, k, j, i) * xtra; + } else if (use_tra_ic) { + for (int k = ks; k <= ke; k++) + for (int j = js; j <= je; j++) + for (int i = is; i <= ie; i++) { + if (phydro->w(IPR, k, j, i) > ptra) { + for (int n = 0; n < NCHEM; ++n) { + pscalars->s(n, k, j, i) += phydro->w(IDN, k, j, i) * xtra; + } } } - } - } + } } + + peos->PrimitiveToConserved(phydro->w, pfield->bcc, phydro->u, pcoord, is, ie, + js, je, ks, ke); } diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 2b35c535..d736eb73 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -19,6 +19,7 @@ include_directories( SYSTEM ${HARP_INCLUDE_DIR} SYSTEM ${KINTERA_INCLUDE_DIR} SYSTEM ${SNAP_INCLUDE_DIR} + SYSTEM ${FMT_INCLUDE_DIR} ) set(namel canoe) @@ -27,7 +28,17 @@ string(TOUPPER ${namel} nameu) string(TOLOWER ${CMAKE_BUILD_TYPE} buildl) string(TOUPPER ${CMAKE_BUILD_TYPE} buildu) -file(GLOB src_files *.cpp) +file(GLOB src_files + *.cpp + climath/*.c + astro/*.cpp + astro/Jupiter/*.cpp + astro/Saturn/*.cpp + utils/*.cpp + tasklist/*.cpp + exo3/*.cpp + exchanger/*.cpp + ) add_library(${namel}_${buildl} OBJECT @@ -39,15 +50,27 @@ set_target_properties(${namel}_${buildl} COMPILE_FLAGS ${CMAKE_CXX_FLAGS_${buildu}} ) -add_subdirectory(climath) -add_subdirectory(astro) -add_subdirectory(utils) +target_link_libraries(${namel}_${buildl} + PUBLIC + ${NETCDF_LIBRARIES} + ${PNETCDF_LIBRARY} + ${TORCH_LIBRARY} + ${TORCH_CPU_LIBRARY} + ${C10_LIBRARY} + ${DISORT_LIBRARY} + ${HARP_LIBRARY} + ${KINTERA_LIBRARY} + ${SNAP_LIBRARY} + ${MPI_CXX_LIBRARIES} + ${GLOG_LIBRARIES} + ${OpenMP_libomp_LIBRARY} + ${EXTRA_LIBS} + fmt::fmt + yaml-cpp::yaml-cpp + ) + add_subdirectory(snap) add_subdirectory(outputs) -add_subdirectory(tasklist) -add_subdirectory(exo3) -add_subdirectory(exchanger) -add_subdirectory(special) if(UNIX AND NOT APPLE) set(EXTRA_LIBS dl stdc++fs) @@ -57,33 +80,10 @@ set(CANOE_LIBRARY_${buildu} "canoe_${buildl}" "athenapp_${buildl}" "application_${buildl}" - "climath_${buildl}" - "astro_${buildl}" - "utils_${buildl}" "snap_${buildl}" "outputs_${buildl}" - "tasklist_${buildl}" - "exo3_${buildl}" - "exchanger_${buildl}" $<$:minichem::minichem++> $<$:exofmsrt::exofmsrt++> - ${YAML_CPP_LIBRARIES} - ${CPPDISORT_LIBRARY_${buildu}} - ${PYTHON_LIBRARY_RELEASE} - ${NETCDF_LIBRARIES} - ${PNETCDF_LIBRARY} - ${MPI_CXX_LIBRARIES} - ${GLOG_LIBRARIES} - ${OpenMP_libomp_LIBRARY} - ${CANTERA_LIBRARY} - ${TORCH_LIBRARY} - ${TORCH_CPU_LIBRARY} - ${C10_LIBRARY} - ${DISORT_LIBRARY} - ${HARP_LIBRARY} - ${KINTERA_LIBRARY} - ${SNAP_LIBRARY} - ${EXTRA_LIBS} CACHE INTERNAL "canoe library") set(CANOE_INCLUDE_DIR diff --git a/src/interface/eos.hpp b/src/interface/eos.hpp index ae48da76..67b0ae85 100644 --- a/src/interface/eos.hpp +++ b/src/interface/eos.hpp @@ -4,6 +4,8 @@ #include // snap +#include + #include // canoe @@ -24,6 +26,19 @@ inline torch::Tensor get_temp(snap::IdealMoist peos, return peos->compute("W->T", {w}); } +inline torch::Tensor get_cv(snap::IdealMoist peos, + AthenaArray const& hydro_w) { + auto pthermo = peos->pthermo; + auto mud = kintera::species_weights[0]; + auto Rd = kintera::constants::Rgas / mud; + auto cvd = kintera::species_cref_R[0] * Rd; + int ny = pthermo->options.vapor_ids().size() + + pthermo->options.cloud_ids().size() - 1; + auto w = get_all(hydro_w); + auto fsig = peos->f_sig(w.narrow(0, snap::ICY, ny)); + return cvd * fsig; +} + //! \brief Effective polytropic index //! //! Eq.63 in Li2019 From 5cd6fc42d8d7cb1cd73fc85effd62e3faed991ae Mon Sep 17 00:00:00 2001 From: mac/cli Date: Tue, 12 Aug 2025 17:53:28 -0400 Subject: [PATCH 04/10] wip --- cmake/examples/minichem.cmake | 6 +- cmake/modules/FindKintera.cmake | 19 +---- examples/2024-XZhang-minichem/minichem.cpp | 89 ++++++++++++++++++++++ 3 files changed, 94 insertions(+), 20 deletions(-) diff --git a/cmake/examples/minichem.cmake b/cmake/examples/minichem.cmake index 93a3e5a8..f4895240 100644 --- a/cmake/examples/minichem.cmake +++ b/cmake/examples/minichem.cmake @@ -10,10 +10,8 @@ endmacro() set_if_empty(NUMBER_GHOST_CELLS 3) set(MINICHEM ON) -set(NCHEMISTRY 13) -#set(NVAPOR 1) -#set(NCLOUD 2) -#set(NPHASE_LEGACY 3) +set(NCHEM 13) +# set(NVAPOR 1) set(NCLOUD 2) set(NPHASE_LEGACY 3) set(NETCDF ON) set(PNETCDF ON) set(MPI ON) diff --git a/cmake/modules/FindKintera.cmake b/cmake/modules/FindKintera.cmake index 18efb724..5abc6ef1 100644 --- a/cmake/modules/FindKintera.cmake +++ b/cmake/modules/FindKintera.cmake @@ -4,8 +4,6 @@ # # KINTERA_LIBRARY: Link these libraries when using KINTERA # -# VAPORS_LIBRARY: Link these libraries when using KINTERA -# # KINTERA_FOUND: True if kintera found # # Normal usage would be: @@ -14,7 +12,7 @@ # # include_directories(${KINTERA_INCLUDE_DIR}) # -# target_link_libraries(${KINTERA_LIBRARY} ${VAPORS_LIBRARY}) +# target_link_libraries(${KINTERA_LIBRARY}) include(FindPackageHandleStandardArgs) @@ -102,22 +100,11 @@ find_library( $ENV{KINTERA_DIR}/lib $ENV{KINTERA_ROOT}/lib) -find_library( - VAPORS_LIBRARY vapors_release - HINTS ${kintera_lib_dir} - /opt/homebrew/lib - /usr/lib/x86_64-linux-gnu/ - KINTERA_DIR/lib - KINTERA_LIB - $ENV{KINTERA_LIB} - $ENV{KINTERA_DIR}/lib - $ENV{KINTERA_ROOT}/lib) - if(${CUDAToolKit_FOUND}) - set(kintera_required_vars VAPORS_LIBRARY KINTERA_LIBRARY KINTERA_CUDA_LIBRARY + set(kintera_required_vars KINTERA_LIBRARY KINTERA_CUDA_LIBRARY KINTERA_INCLUDE_DIR) else() - set(kintera_required_vars VAPORS_LIBRARY KINTERA_LIBRARY KINTERA_INCLUDE_DIR) + set(kintera_required_vars KINTERA_LIBRARY KINTERA_INCLUDE_DIR) endif() mark_as_advanced(${kintera_required_vars}) diff --git a/examples/2024-XZhang-minichem/minichem.cpp b/examples/2024-XZhang-minichem/minichem.cpp index 28aa8049..0653fe1b 100644 --- a/examples/2024-XZhang-minichem/minichem.cpp +++ b/examples/2024-XZhang-minichem/minichem.cpp @@ -5,6 +5,9 @@ #include #include +// yaml +#include + // athena #include #include @@ -24,6 +27,7 @@ // snap #include +#include // canoe #include @@ -226,6 +230,91 @@ void Mesh::InitUserMeshData(ParameterInput *pin) { } } +torch::Tensor setup_moist_adiabatic_profile(std::string infile) { + auto config = YAML::LoadFile(infile); + auto Ps = config["problem"]["Ps"].as(1.e5); + auto Ts = config["problem"]["Ts"].as(300.); + auto Tmin = config["problem"]["Tmin"].as(200.); + auto grav = -config["forcing"]["const-gravity"]["grav1"].as(); + + // initialize the block + auto block = snap::MeshBlock(snap::MeshBlockOptions::from_yaml(infile)); + + // useful modules + auto phydro = block->phydro; + auto pcoord = phydro->pcoord; + auto peos = phydro->peos; + auto m = block->named_modules()["hydro.eos.thermo"]; + auto thermo_y = std::dynamic_pointer_cast(m); + + // dimensions and indices + int nc3 = pcoord->x3v.size(0); + int nc2 = pcoord->x2v.size(0); + int nc1 = pcoord->x1v.size(0); + int ny = thermo_y->options.species().size() - 1; + int nvar = peos->nvar(); + + // construct an adiabatic atmosphere + kintera::ThermoX thermo_x(thermo_y->options); + + auto temp = Ts * torch::ones({nc3, nc2}, torch::kFloat64); + auto pres = Ps * torch::ones({nc3, nc2}, torch::kFloat64); + auto xfrac = torch::zeros({nc3, nc2, 1 + ny}, torch::kFloat64); + auto w = torch::zeros({nvar, nc3, nc2, nc1}, torch::kFloat64); + + // read in compositions + for (int i = 1; i <= ny; ++i) { + auto name = thermo_y->options.species()[i]; + auto xmixr = config["problem"]["x" + name].as(0.); + xfrac.select(2, i) = xmixr; + } + + // dry air mole fraction + xfrac.select(2, 0) = 1. - xfrac.narrow(-1, 1, ny).sum(-1); + + // adiabatic extrapolate half a grid to cell center + int is = pcoord->is(); + int ie = pcoord->ie(); + auto dz = pcoord->dx1f[is].item(); + thermo_x->extrapolate_ad(temp, pres, xfrac, grav, dz / 2.); + + int i = is; + int nvapor = thermo_x->options.vapor_ids().size(); + int ncloud = thermo_x->options.cloud_ids().size(); + for (; i <= ie; ++i) { + auto conc = thermo_x->compute("TPX->V", {temp, pres, xfrac}); + + w[IPR].select(2, i) = pres; + w[IDN].select(2, i) = thermo_x->compute("V->D", {conc}); + + auto result = thermo_x->compute("X->Y", {xfrac}); + w.narrow(0, snap::ICY, ny).select(3, i) = + thermo_x->compute("X->Y", {xfrac}); + + if ((temp < Tmin).any().item()) break; + dz = pcoord->dx1f[i].item(); + thermo_x->extrapolate_ad(temp, pres, xfrac, grav, dz); + } + + // isothermal extrapolation + for (; i <= ie; ++i) { + auto mu = (thermo_x->mu * xfrac).sum(-1); + dz = pcoord->dx1f[i].item(); + pres *= exp(-grav * mu * dz / (kintera::constants::Rgas * temp)); + auto conc = thermo_x->compute("TPX->V", {temp, pres, xfrac}); + w[IPR].select(2, i) = pres; + w[IDN].select(2, i) = thermo_x->compute("V->D", {conc}); + w.narrow(0, snap::ICY, ny).select(3, i) = + thermo_x->compute("X->Y", {xfrac}); + } + + // add noise + w[IVX] += 0.01 * torch::rand_like(w[IVX]); + w[IVY] += 0.01 * torch::rand_like(w[IVY]); + + return w; +} + void MeshBlock::ProblemGenerator(ParameterInput *pin) { srand(Globals::my_rank + time(0)); From 11d97214c75c19e393e65ed81a552baa7d24b965 Mon Sep 17 00:00:00 2001 From: mac/cli Date: Tue, 12 Aug 2025 18:14:04 -0400 Subject: [PATCH 05/10] wip --- cmake/modules/FindSnap.cmake | 20 +++----------------- 1 file changed, 3 insertions(+), 17 deletions(-) diff --git a/cmake/modules/FindSnap.cmake b/cmake/modules/FindSnap.cmake index e153d073..2afd8961 100644 --- a/cmake/modules/FindSnap.cmake +++ b/cmake/modules/FindSnap.cmake @@ -4,8 +4,6 @@ # # SNAP_LIBRARY: Link these libraries when using SNAP # -# BC_LIBRARY: Link these libraries when using SNAP -# # SNAP_FOUND: True if snap found # # Normal usage would be: @@ -14,7 +12,7 @@ # # include_directories(${SNAP_INCLUDE_DIR}) # -# target_link_libraries(${SNAP_LIBRARY} ${BC_LIBRARY}) +# target_link_libraries(${SNAP_LIBRARY}) include(FindPackageHandleStandardArgs) @@ -100,22 +98,10 @@ find_library( $ENV{SNAP_DIR}/lib $ENV{SNAP_ROOT}/lib) -find_library( - BC_LIBRARY bc_release - HINTS ${snap_lib_dir} - /opt/homebrew/lib - /usr/lib/x86_64-linux-gnu/ - SNAP_DIR/lib - SNAP_LIB - $ENV{SNAP_LIB} - $ENV{SNAP_DIR}/lib - $ENV{SNAP_ROOT}/lib) - if(${CUDAToolKit_FOUND}) - set(snap_required_vars BC_LIBRARY SNAP_LIBRARY SNAP_CUDA_LIBRARY - SNAP_INCLUDE_DIR) + set(snap_required_vars SNAP_LIBRARY SNAP_CUDA_LIBRARY SNAP_INCLUDE_DIR) else() - set(snap_required_vars BC_LIBRARY SNAP_LIBRARY SNAP_INCLUDE_DIR) + set(snap_required_vars SNAP_LIBRARY SNAP_INCLUDE_DIR) endif() mark_as_advanced(${snap_required_vars}) From e70061678250e474c7d7febe4d4aea251948a053 Mon Sep 17 00:00:00 2001 From: mac/cli Date: Tue, 12 Aug 2025 21:12:58 -0400 Subject: [PATCH 06/10] wip --- examples/2024-XZhang-minichem/jup-dry.yaml | 15 +++++++++++++++ examples/2024-XZhang-minichem/mini.inp | 2 ++ examples/2024-XZhang-minichem/minichem.cpp | 2 +- 3 files changed, 18 insertions(+), 1 deletion(-) create mode 100644 examples/2024-XZhang-minichem/jup-dry.yaml diff --git a/examples/2024-XZhang-minichem/jup-dry.yaml b/examples/2024-XZhang-minichem/jup-dry.yaml new file mode 100644 index 00000000..04f3b210 --- /dev/null +++ b/examples/2024-XZhang-minichem/jup-dry.yaml @@ -0,0 +1,15 @@ +reference-state: + Tref: 0. + Pref: 1.e5 + +species: + - name: dry + composition: {H: 1.5, He: 0.15} + cv_R: 2.5 + +geometry: + type: cartesian + +dynamics: + equation-of-state: + type: ideal-moist diff --git a/examples/2024-XZhang-minichem/mini.inp b/examples/2024-XZhang-minichem/mini.inp index 2e061720..9e5833a8 100644 --- a/examples/2024-XZhang-minichem/mini.inp +++ b/examples/2024-XZhang-minichem/mini.inp @@ -106,6 +106,8 @@ top_sponge_lyr.width = 20.E3 #relax_time = 1.E4 +config_file = jup-dry.yaml + qH2O.ppmv = 2500. #minichem diff --git a/examples/2024-XZhang-minichem/minichem.cpp b/examples/2024-XZhang-minichem/minichem.cpp index 0653fe1b..fdf8920f 100644 --- a/examples/2024-XZhang-minichem/minichem.cpp +++ b/examples/2024-XZhang-minichem/minichem.cpp @@ -55,7 +55,7 @@ std::vector vmass = {17.01, 2.02, 18.02, 1.01, 28.01, 44.01, 16., MiniChem *mc; void MeshBlock::InitUserMeshBlockData(ParameterInput *pin) { - int nvapor = pimpl->peos->pthermo->options.vapor_ids().size() - 1; + int nvapor = 0; if (nvapor > 0) { AllocateUserOutputVariables(5 + nvapor); From 3f7c1bf9e2b660f34cf0622925fe549c2d7abfc0 Mon Sep 17 00:00:00 2001 From: mac/cli Date: Tue, 12 Aug 2025 21:49:15 -0400 Subject: [PATCH 07/10] wip --- examples/2024-XZhang-minichem/minichem.cpp | 2 ++ 1 file changed, 2 insertions(+) diff --git a/examples/2024-XZhang-minichem/minichem.cpp b/examples/2024-XZhang-minichem/minichem.cpp index fdf8920f..67a2f60c 100644 --- a/examples/2024-XZhang-minichem/minichem.cpp +++ b/examples/2024-XZhang-minichem/minichem.cpp @@ -349,6 +349,8 @@ void MeshBlock::ProblemGenerator(ParameterInput *pin) { auto temp = get_temp(pimpl->peos, phydro->w); auto temp_a = temp.accessor(); + Tbot = temp_a[ks][js][is]; + if (NCHEM > 0) { if (use_mini_ic) { // minichem initialize conserved variables From a7006cc7310ba4a7adf3cecbc99b743b7205e6e1 Mon Sep 17 00:00:00 2001 From: mac/cli Date: Wed, 13 Aug 2025 14:02:49 -0400 Subject: [PATCH 08/10] wip --- examples/2024-XZhang-minichem/mini.inp | 15 ++------------- examples/2024-XZhang-minichem/minichem.cpp | 21 +++++++++++++-------- 2 files changed, 15 insertions(+), 21 deletions(-) diff --git a/examples/2024-XZhang-minichem/mini.inp b/examples/2024-XZhang-minichem/mini.inp index 9e5833a8..3b4f92ed 100644 --- a/examples/2024-XZhang-minichem/mini.inp +++ b/examples/2024-XZhang-minichem/mini.inp @@ -68,7 +68,8 @@ nx2 = 16 nx3 = 1 -gamma = 1.42 +Rd = 3936.04 +gamma = 1.4 grav_acc1 = -24.79 implicit_flag = 1 @@ -79,18 +80,6 @@ cloud = H2O(c), H2O(p) microphysics_config = water_ammonia.yaml - -Rd = 3777. # mu = 2.3175 g/mol -eps1 = 8.18 8.18 8.18 -beta1 = 0. 24.845 24.845 -rcp1 = 0.15 0.33 0.33 -Ttriple1 = 273.16 -Ptriple1 = 611.7 - -sa.relax = 1.0 -sa.max_iter = 5 -sa.ftol = 1.e-2 - #packages = bot_heating, top_cooling #bot_heating.flux = 8. diff --git a/examples/2024-XZhang-minichem/minichem.cpp b/examples/2024-XZhang-minichem/minichem.cpp index 67a2f60c..444731fb 100644 --- a/examples/2024-XZhang-minichem/minichem.cpp +++ b/examples/2024-XZhang-minichem/minichem.cpp @@ -47,8 +47,6 @@ Real Tbot, xtra, ptra, Ttop_tau, Tbot_tau, Trabot_tau; std::string met; bool use_mini, use_mini_ic, use_tra_ic, fix_bot_tra; -Real p0 = 1.e5; - //'OH','H2','H2O','H','CO','CO2','O','CH4','C2H2','NH3','N2','HCN', 'He' std::vector vmass = {17.01, 2.02, 18.02, 1.01, 28.01, 44.01, 16., 16.05, 26.04, 17.04, 28.02, 27.03, 4.}; @@ -81,7 +79,7 @@ void MeshBlock::UserWorkBeforeOutput(ParameterInput *pin) { auto temp = get_temp(pimpl->peos, w); auto pres = get_pres(w); auto chi = (get_gammad() - 1.) / get_gammad(); - auto theta = temp * pow(p0 / pres, chi); + auto theta = temp * pow(P0 / pres, chi); auto temp_a = temp.accessor(); auto theta_a = theta.accessor(); @@ -316,13 +314,16 @@ torch::Tensor setup_moist_adiabatic_profile(std::string infile) { } void MeshBlock::ProblemGenerator(ParameterInput *pin) { + // auto infile = pin->GetString("problem", "config_file"); + // auto w = setup_moist_adiabatic_profile(infile) + srand(Globals::my_rank + time(0)); // thermodynamic constants Real gamma = pin->GetReal("hydro", "gamma"); - auto mud = kintera::species_weights[0]; - auto Rd = kintera::constants::Rgas / mud; + Real Rd = pin->GetReal("hydro", "Rd"); Real cp = gamma / (gamma - 1.) * Rd; + std::cout << "Rd = " << Rd << ", cp = " << cp << std::endl; // set up an adiabatic atmosphere int max_iter = 400, iter = 0; @@ -330,6 +331,8 @@ void MeshBlock::ProblemGenerator(ParameterInput *pin) { // estimate surface temperature and pressure auto x1min = pcoord->x1f(is); + auto x1max = pcoord->x1f(ie + 1); + Real Ts = T0 - grav / cp * x1min; Real Ps = P0 * pow(Ts / T0, cp / Rd); @@ -340,10 +343,10 @@ void MeshBlock::ProblemGenerator(ParameterInput *pin) { Real x1 = pcoord->x1v(i); Real x2 = pcoord->x2v(j); Real x3 = pcoord->x3v(k); - Real temp = Ts - grav * x1 / cp; - phydro->w(IPR, k, j, i) = p0 * pow(temp / Ts, cp / Rd); + Real temp = T0 - grav * x1 / cp; + phydro->w(IPR, k, j, i) = P0 * pow(temp / T0, cp / Rd); phydro->w(IDN, k, j, i) = phydro->w(IPR, k, j, i) / (Rd * temp); - phydro->w(IVX, k, j, i) = phydro->w(IVY, k, j, i) = 0.; + phydro->w(IVX, k, j, i) = 0.001 * (1. * rand() / RAND_MAX - 0.5); } auto temp = get_temp(pimpl->peos, phydro->w); @@ -351,6 +354,8 @@ void MeshBlock::ProblemGenerator(ParameterInput *pin) { Tbot = temp_a[ks][js][is]; + std::cout << "ptra = " << ptra << " Pa" << std::endl; + if (NCHEM > 0) { if (use_mini_ic) { // minichem initialize conserved variables From 223b45f50f46b0836a589d25c750bf1e5b0f40c1 Mon Sep 17 00:00:00 2001 From: mac/cli Date: Wed, 13 Aug 2025 21:17:53 -0400 Subject: [PATCH 09/10] wip --- cmake/minichem.cmake | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/cmake/minichem.cmake b/cmake/minichem.cmake index 9850ecc2..c22f3fe4 100644 --- a/cmake/minichem.cmake +++ b/cmake/minichem.cmake @@ -7,4 +7,4 @@ set(patch_command set(PACKAGE_NAME minichem) set(REPO_URL "https://github.com/chengcli/mini_chem") set(REPO_TAG "3372a500f038f83228c9f8f944b3fb6b2dedc572") -add_package(${PACKAGE_NAME} ${REPO_URL} ${REPO_TAG} "${patch_commend}" OFF) +add_package(${PACKAGE_NAME} ${REPO_URL} ${REPO_TAG} "${patch_command}" OFF) From 356752062ad1f3e9cfc5f6157193d17cab1982da Mon Sep 17 00:00:00 2001 From: mac/cli Date: Wed, 13 Aug 2025 22:10:22 -0400 Subject: [PATCH 10/10] wip --- src/harp/radiation_band.cpp | 298 ---------- src/harp/radiation_band.hpp | 246 --------- src/harp/rt_solver_disort.cpp | 425 -------------- src/harp/rt_solver_toon.cpp | 797 --------------------------- src/harp/rt_solvers.hpp | 155 ------ src/harp/set_spectral_properties.cpp | 143 ----- src/opacity/absorber.hpp | 114 ---- src/opacity/absorber_ck.cpp | 170 ------ src/opacity/absorber_ck.hpp | 42 -- src/opacity/absorber_factory.cpp | 147 ----- 10 files changed, 2537 deletions(-) delete mode 100644 src/harp/radiation_band.cpp delete mode 100644 src/harp/radiation_band.hpp delete mode 100644 src/harp/rt_solver_disort.cpp delete mode 100644 src/harp/rt_solver_toon.cpp delete mode 100644 src/harp/rt_solvers.hpp delete mode 100644 src/harp/set_spectral_properties.cpp delete mode 100644 src/opacity/absorber.hpp delete mode 100644 src/opacity/absorber_ck.cpp delete mode 100644 src/opacity/absorber_ck.hpp delete mode 100644 src/opacity/absorber_factory.cpp diff --git a/src/harp/radiation_band.cpp b/src/harp/radiation_band.cpp deleted file mode 100644 index f1eddd6d..00000000 --- a/src/harp/radiation_band.cpp +++ /dev/null @@ -1,298 +0,0 @@ -// C/C++ -#include -#include -#include -#include - -// external -#include - -// athena -#include -#include -#include -#include - -// canoe -#include - -// climath -#include - -// application -#include -#include - -// utils -#include -#include -#include - -// opacity -#include - -// harp -#include "radiation.hpp" -#include "radiation_band.hpp" -#include "rt_solvers.hpp" - -RadiationBand::RadiationBand(std::string myname, YAML::Node const &rad, - bool load_opacity) - : NamedGroup(myname) { - Application::Logger app("harp"); - app->Log("Initialize RadiationBand " + myname); - - if (!rad[myname]) { - throw NotFoundError("RadiationBand", myname); - } - - auto my = rad[myname]; - - if (my["parameters"]) { - SetRealsFrom(my["parameters"]); - } - - pgrid_ = SpectralGridFactory::CreateFrom(my); - - wrange_ = RadiationHelper::parse_wave_range(my); - - if (my["outdir"]) { - if (!my["outdir"].IsSequence()) { - throw RuntimeError("RadiationBand", "outdir must be a sequence"); - } - - for (const auto &item : my["outdir"]) { - rayOutput_.push_back( - RadiationHelper::parse_radiation_direction(item.as())); - } - } - - // set absorbers - if (my["opacity"]) { - if (!my["opacity"].IsSequence()) { - throw RuntimeError("RadiationBand", "opacity must be a sequence"); - } - - auto names = my["opacity"].as>(); - absorbers_ = AbsorberFactory::CreateFrom(names, GetName(), rad); - - if (load_opacity) { - for (auto &ab : absorbers_) { - ab->LoadOpacity(RadiationBandsFactory::GetBandId(myname)); - // Correlated-k absorbers need to modify the spectral grid - ab->ModifySpectralGrid(pgrid_->spec); - } - } - } - - // set flags - if (my["flags"]) { - if (!my["flags"].IsSequence()) { - throw RuntimeError("RadiationBand", "flags must be a sequence"); - } - - auto flag_strs = my["flags"].as>(); - for (auto flag : flag_strs) { - SetFlag(RadiationHelper::parse_radiation_flags(flag)); - } - } - - // set rt solver - if (my["rt-solver"]) { - psolver_ = CreateRTSolverFrom(my["rt-solver"].as(), rad); - } else { - psolver_ = nullptr; - } - - char buf[80]; - snprintf(buf, sizeof(buf), "%.2f - %.2f", wrange_.first, wrange_.second); - app->Log("Spectral range", buf); - app->Log("Number of spectral bins", pgrid_->spec.size()); -} - -RadiationBand::~RadiationBand() { - Application::Logger app("harp"); - app->Log("Destroy RadiationBand " + GetName()); -} - -void RadiationBand::Resize(int nc1, int nc2, int nc3, int nstr, - MeshBlock const *pmb) { - // allocate memory for spectral properties - tem_.resize(nc1); - temf_.resize(nc1 + 1); - - tau_.NewAthenaArray(pgrid_->spec.size(), nc1); - tau_.ZeroClear(); - - ssa_.NewAthenaArray(pgrid_->spec.size(), nc1); - ssa_.ZeroClear(); - - pmom_.NewAthenaArray(pgrid_->spec.size(), nc1, nstr + 1); - pmom_.ZeroClear(); - - // spectral grids properties - toa_.NewAthenaArray(pgrid_->spec.size(), rayOutput_.size(), nc3, nc2); - toa_.ZeroClear(); - - flxup_.NewAthenaArray(pgrid_->spec.size(), nc3, nc2, nc1); - flxup_.ZeroClear(); - - flxdn_.NewAthenaArray(pgrid_->spec.size(), nc3, nc2, nc1); - flxdn_.ZeroClear(); - - // band properties - btau.NewAthenaArray(nc3, nc2, nc1); - bssa.NewAthenaArray(nc3, nc2, nc1); - bpmom.NewAthenaArray(nstr + 1, nc3, nc2, nc1); - - btoa.NewAthenaArray(rayOutput_.size(), nc3, nc2); - bflxup.NewAthenaArray(nc3, nc2, nc1 + 1); - bflxdn.NewAthenaArray(nc3, nc2, nc1 + 1); - - // exchange buffer - pexv = std::make_shared>(GetName()); - - int nlayers = GetNumLayers(); - int npmom = GetNumPhaseMoments(); - pexv->send_buffer[0].resize(temf_.size()); - pexv->send_buffer[1].resize(pgrid_->spec.size()*nlayers * (npmom + 3)); - - pexv->Regroup(pmb, X1DIR); - int nblocks = pexv->GetGroupSize(); - pexv->recv_buffer[0].resize(nblocks * pexv->send_buffer[0].size()); - pexv->recv_buffer[1].resize(nblocks * pexv->send_buffer[1].size()); - - if (psolver_ != nullptr) { - psolver_->Resize(nblocks * (nc1 - 2 * NGHOST), nstr); - } -} - -AbsorberPtr RadiationBand::GetAbsorberByName(std::string const &name) { - for (auto &absorber : absorbers_) { - if (absorber->GetName() == name) { - return absorber; - } - } - - throw NotFoundError("RadiationBand", "Absorber " + name); - - return nullptr; -} - -RadiationBand const *RadiationBand::CalBandFlux(MeshBlock const *pmb, int k, - int j) { - // reset flux of this column - for (int i = pmb->is; i <= pmb->ie + 1; ++i) { - bflxup(k, j, i) = 0.; - bflxdn(k, j, i) = 0.; - } - - psolver_->Prepare(pmb, k, j); - - psolver_->CalBandFlux(pmb, k, j); - - return this; -} - -RadiationBand const *RadiationBand::CalBandRadiance(MeshBlock const *pmb, int k, - int j) { - // reset radiance of this column - for (int n = 0; n < GetNumOutgoingRays(); ++n) { - btoa(n, k, j) = 0.; - } - - psolver_->Prepare(pmb, k, j); - psolver_->CalBandRadiance(pmb, k, j); - - return this; -} - -void RadiationBand::WriteAsciiHeader(OutputParameters const *pout) const { - if (!TestFlag(RadiationFlags::WriteBinRadiance)) return; - - char fname[80], number[6]; - snprintf(number, sizeof(number), "%05d", pout->file_number); - snprintf(fname, sizeof(fname), "%s.radiance.%s.txt", GetName().c_str(), - number); - FILE *pfile = fopen(fname, "w"); - - fprintf(pfile, "# Bin Radiances of Band %s: %.3g - %.3g\n", GetName().c_str(), - wrange_.first, wrange_.second); - fprintf(pfile, "# Ray output size: %lu\n", rayOutput_.size()); - - fprintf(pfile, "# Polar angles: "); - for (auto &ray : rayOutput_) { - fprintf(pfile, "%.3f", rad2deg(acos(ray.mu))); - } - fprintf(pfile, "\n"); - - fprintf(pfile, "# Azimuthal angles: "); - for (auto &ray : rayOutput_) { - fprintf(pfile, "%.3f", rad2deg(ray.phi)); - } - fprintf(pfile, "\n"); - - fprintf(pfile, "#%12s%12s", "wave1", "wave2"); - for (size_t j = 0; j < rayOutput_.size(); ++j) { - fprintf(pfile, "%12s%lu", "Radiance", j + 1); - } - fprintf(pfile, "%12s\n", "weight"); - - fclose(pfile); -} - -void RadiationBand::WriteAsciiData(OutputParameters const *pout) const { - if (!TestFlag(RadiationFlags::WriteBinRadiance)) return; - - char fname[80], number[6]; - snprintf(number, sizeof(number), "%05d", pout->file_number); - snprintf(fname, sizeof(fname), "%s.radiance.%s.txt", GetName().c_str(), - number); - FILE *pfile = fopen(fname, "w"); - - for (size_t i = 0; i < pgrid_->spec.size(); ++i) { - fprintf(pfile, "%13.3g%12.3g", pgrid_->spec[i].wav1, pgrid_->spec[i].wav2); - for (size_t j = 0; j < rayOutput_.size(); ++j) { - fprintf(pfile, "%12.3f", toa_(i, j)); - } - if (TestFlag(RadiationFlags::Normalize) && - (wrange_.first != wrange_.second)) { - fprintf(pfile, "%12.3g\n", - pgrid_->spec[i].wght / (wrange_.second - wrange_.first)); - } else { - fprintf(pfile, "%12.3g\n", pgrid_->spec[i].wght); - } - } - - fclose(pfile); -} - -std::string RadiationBand::ToString() const { - std::stringstream ss; - ss << "RadiationBand: " << GetName() << std::endl; - ss << "Absorbers: "; - for (auto &ab : absorbers_) { - ss << ab->GetName() << ", "; - } - ss << std::endl << "RT-Solver: " << psolver_->GetName(); - return ss.str(); -} - -std::shared_ptr RadiationBand::CreateRTSolverFrom( - std::string const &rt_name, YAML::Node const &rad) { - std::shared_ptr psolver; - - if (rt_name == "Lambert") { - psolver = std::make_shared(this, rad); -#ifdef RT_DISORT - } else if (rt_name == "Disort") { - psolver = std::make_shared(this, rad); - } else if (rt_name == "Toon") { - psolver = std::make_shared(this, rad); -#endif // RT_DISORT - } else { - throw NotFoundError("RadiationBand", rt_name); - } - - return psolver; -} diff --git a/src/harp/radiation_band.hpp b/src/harp/radiation_band.hpp deleted file mode 100644 index 76a34337..00000000 --- a/src/harp/radiation_band.hpp +++ /dev/null @@ -1,246 +0,0 @@ -#ifndef SRC_HARP_RADIATION_BAND_HPP_ -#define SRC_HARP_RADIATION_BAND_HPP_ - -// C/C++ headers -#include -#include -#include -#include - -// athena -#include // Real - -// canoe -#include -#include - -// exchanger -#include - -// harp -#include "spectral_grid.hpp" - -namespace YAML { -class Node; -} - -class OutputParameters; -class Absorber; -class SpectralGridBase; -struct Direction; - -class RadiationBand : public NamedGroup, - public FlagGroup, - public ParameterGroup, - public ASCIIOutputGroup, - public StringReprGroup { - public: // public access data - // implementation of RT Solver - class RTSolver; - class RTSolverLambert; - class RTSolverDisort; - class RTSolverToon; - - std::shared_ptr> pexv; - - //! band optical depth - AthenaArray btau; - - //! band single scattering albedo - AthenaArray bssa; - - //! band phase function moments - AthenaArray bpmom; - - //! band upward flux (shallow reference to fluxup) - AthenaArray bflxup; - - //! band downward flux (shallow reference to flxdn) - AthenaArray bflxdn; - - //! \brief band top-of-the-atmosphere radiance (shallow reference to radiance) - AthenaArray btoa; - - public: // constructor and destructor - RadiationBand(std::string name, YAML::Node const &rad, bool load_opacity); - virtual ~RadiationBand(); - - public: // member functions - //! \brief Allocate memory for radiation band - void Resize(int nc1, int nc2, int nc3, int nstr, MeshBlock const *pmb); - - //! \brief Create radiative transfer solver from YAML node - //! - //! \param[in] rad YAML node containing whole radiation configuration - std::shared_ptr CreateRTSolverFrom(std::string const &name, - YAML::Node const &rad); - - //! Get number of spectral grids - size_t GetNumSpecGrids() const { return pgrid_->spec.size(); } - - //! Get number of absorbers - size_t GetNumAbsorbers() const { return absorbers_.size(); } - - std::vector> const &Absorbers() const { - return absorbers_; - } - - //! Get number of phase function moments - size_t GetNumPhaseMoments() const { return pmom_.GetDim1() - 1; } - - //! Get number of phase function moments - size_t GetNumLayers() const { return tem_.size() - 2 * NGHOST; } - - //! Get range of wavenumbers - std::pair GetRange() const { return wrange_; } - - //! Get an individual absorber - std::shared_ptr GetAbsorber(int i) { return absorbers_[i]; } - - //! Get an individual absorber by name - std::shared_ptr GetAbsorberByName(std::string const &name); - - //! Get number of outgoing rays - size_t GetNumOutgoingRays() { return rayOutput_.size(); } - - //! Get cosine polar angle of the outgoing ray - Real GetCosinePolarAngle(int n) const { return rayOutput_[n].mu; } - - //! Get azimuthal angle of the outgoing ray - Real GetAzimuthalAngle(int n) const { return rayOutput_[n].phi; } - - public: // python bindings - Real const *_GetToaPtr() const { return toa_.data(); } - Real const *_GetTauPtr() const { return tau_.data(); } - - public: // inbound functions - //! \brief Set spectral properties for an air column - //! - //! \param[in] air mole fraction representation of air column - //! \param[in] pcoord coordinates - //! \param[in] gH0 grav * H0 - //! \param[in] k horizontal index - //! \param[in] j horizontal index - void SetSpectralProperties(AirColumn &air, Real const *x1f, Real gH0 = 0, - int k = 0, int j = 0); - - //! \brief Calculate band radiative fluxes - RadiationBand const *CalBandFlux(MeshBlock const *pmb, int k, int j); - - //! \brief Calculate band radiances - RadiationBand const *CalBandRadiance(MeshBlock const *pmb, int k, int j); - - //! \brief Set outgoing ray directions - //! - //! \param[in] rayOutput outgoing ray directions - void SetOutgoingRays(std::vector const &ray) { rayOutput_ = ray; } - - public: // ASCIIOutputGroup functions - void WriteAsciiHeader(OutputParameters const *) const override; - void WriteAsciiData(OutputParameters const *) const override; - - public: // StringReprGroup functions - std::string ToString() const override; - - public: // Exchanger functions - //! \brief Pack temperature at cell face into send buffer 0 - void packTemperature(); - - //! \brief Unpack temperature at cell face from receive buffer 0 - bool unpackTemperature(void *arg); - - //! \brief Pack data in all spectral grids into send buffer 1 - //! - //! \param[in] b spectral bin index - void packSpectralProperties(); - - //! \brief Unpack data from receive buffer 1 into spectral grid b - void unpackSpectralProperties(int b, void *arg); - - protected: - //! radiative transfer solver - std::shared_ptr psolver_; - - //! spectral grid - std::shared_ptr pgrid_; - - //! \brief identify wave- number (wavelength, frequency) range - std::pair wrange_; - - //! all absorbers - std::vector> absorbers_; - - //! outgoing rays - std::vector rayOutput_; - - //! \brief spectral bin optical depth - //! - //! This is a two-dimensional array. The first dimension is the size of the - //! spectral bin and the second dimension is the size of the vertical - //! dimension. - AthenaArray tau_; - - //! \brief spectral bin single scattering albedo - //! - //! This is a two-dimensional array. The first dimension is the size of the - //! spectral bin and the second dimension is the size of the vertical - //! dimension. - AthenaArray ssa_; - - //! \brief spectral bin top-of-the-atmosphere radiance - //! - //! This is a two-dimensional array. The first dimension is the size of the - //! spectral bin and the second dimension is the size of the vertical - //! dimension. - AthenaArray toa_; - - //! \brief spectral bin upward flux - //! - //! This is a two-dimensional array. The first dimension is the size of the - //! spectral bin and the second dimension is the size of the vertical - //! dimension. - AthenaArray flxup_; - - //! \brief spectral bin downward flux - //! - //! This is a two-dimensional array. The first dimension is the size of the - //! spectral bin and the second dimension is the size of the vertical - //! dimension. - AthenaArray flxdn_; - - //! \brief spectral bin phase function moments - //! - //! This is a three-dimensional array. The first dimension is the size of the - //! spectral bin. The second dimension is the size of the vertical dimension. - //! The third dimension is the size of the phase function moments (1 + npmom). - AthenaArray pmom_; - - //! temperature at cell center - std::vector tem_; - - //! temperature at cell boundary (face) - std::vector temf_; -}; - -using RadiationBandPtr = std::shared_ptr; -using RadiationBandContainer = std::vector; - -class RadiationBandsFactory { - public: - static RadiationBandContainer CreateFrom(ParameterInput *pin, - std::string key); - static RadiationBandContainer CreateFrom(std::string filename); - - static int GetBandId(std::string const &bname) { - if (band_id_.find(bname) == band_id_.end()) { - return -1; - } - return band_id_.at(bname); - } - - protected: - static std::map band_id_; - static int next_band_id_; -}; - -#endif // SRC_HARP_RADIATION_BAND_HPP_ diff --git a/src/harp/rt_solver_disort.cpp b/src/harp/rt_solver_disort.cpp deleted file mode 100644 index 6762d84e..00000000 --- a/src/harp/rt_solver_disort.cpp +++ /dev/null @@ -1,425 +0,0 @@ -//! \file rt_solver_disort.cpp -//! \brief Call DISORT to perform radiative transfer calculation - -// C/C++ -#include -#include - -// external -#include - -// athena -#include -#include - -// application -#include -#include - -// climath -#include - -// canoe -#include -#include - -// astro -#include - -// exo3 -#include -#include - -// harp -#include "radiation.hpp" -#include "rt_solvers.hpp" - -#ifdef RT_DISORT - -std::map to_map_bool(YAML::Node const &node) { - std::map flags; - - for (auto it = node.begin(); it != node.end(); ++it) { - flags[it->first.as()] = it->second.as(); - } - - return flags; -} - -RadiationBand::RTSolverDisort::RTSolverDisort(RadiationBand *pmy_band, - YAML::Node const &rad) - : RTSolver(pmy_band, "Disort") { - if (rad["Disort-flags"]) { - // SetFlags is a member function of DisortWrapper - SetFlags(to_map_bool(rad["Disort-flags"])); - } - - if (pmy_band->HasPar("btemp")) { - ds_.bc.btemp = pmy_band->GetPar("btemp"); - } - - if (pmy_band->HasPar("albedo")) { - ds_.bc.albedo = pmy_band->GetPar("albedo"); - } else { - // hard code hack to fix the surface albedo to 0, emissivity to 1, since it - // won't pickup from the yaml file - ds_.bc.albedo = 0.; - } - - if (pmy_band->HasPar("ttemp")) { - ds_.bc.ttemp = pmy_band->GetPar("ttemp"); - } - - if (pmy_band->HasPar("temis")) { - ds_.bc.temis = pmy_band->GetPar("temis"); - } else { - ds_.bc.temis = 0.; - } - - if (pmy_band->HasPar("fluor")) { - ds_.bc.fluor = pmy_band->GetPar("fluor"); - } else { - ds_.bc.fluor = 0.; - } - - if (pmy_band->HasPar("fisot")) { - ds_.bc.fisot = pmy_band->GetPar("fisot"); - } else { - ds_.bc.fisot = 0.; - } - - Application::Logger app("harp"); - - // override disort planck flag - if (pmy_band->TestFlag(RadiationFlags::ThermalEmission)) { - ds_.flag.planck = true; - app->Log("Planck function is enabled for band " + pmy_band->GetName()); - } else { - ds_.flag.planck = false; - app->Log("Planck function is disabled for band " + pmy_band->GetName()); - } - - SetHeader("Disort solving band " + pmy_band_->GetName()); -} - -//! \todo update based on band outdir -void RadiationBand::RTSolverDisort::Resize(int nlyr, int nstr) { - RadiationBand::RTSolver::Resize(nlyr, nstr); - Unseal(); - - auto &rayout = pmy_band_->rayOutput_; - auto &&uphi_umu = RadiationHelper::get_direction_grids(rayout); - - SetAtmosphereDimension(nlyr, nstr, nstr); - - dir_dim_[0] = uphi_umu.second.size(); // umu - dir_dim_[1] = uphi_umu.first.size(); // uphi - dir_axis_.resize(dir_dim_[0] + dir_dim_[1]); - - SetIntensityDimension(std::max(dir_dim_[1], 1lu), 1, - std::max(dir_dim_[0], 1lu)); - Seal(); -} - -//! \note Counting Disort Index -//! Example, il = 0, iu = 2, ds_.nlyr = 6, partition in to 3 blocks -//! face id -> 0 - 1 - 2 - 3 - 4 - 5 - 6 -//! cell id -> | 0 | 1 | 2 | 3 | 4 | 5 | -//! disort id -> 6 - 5 - 4 - 3 - 2 - 1 - 0 -//! blocks -> --------- * * -//! -> r = 0 * * * -//! -> --------- * -//! -> r = 1 * * -//! -> --------- -//! -> r = 2 -//! block r = 0 gets, 6 - 5 - 4 -//! block r = 1 gets, 4 - 3 - 2 -//! block r = 2 gets, 2 - 1 - 0 - -void RadiationBand::RTSolverDisort::Prepare(MeshBlock const *pmb, int k, - int j) { - auto &wmin = pmy_band_->wrange_.first; - auto &wmax = pmy_band_->wrange_.second; - - Real dist_au = 1.; - Direction ray = pmb->pimpl->prad->GetRayInput(0); - auto planet = pmb->pimpl->planet; - - if (planet && pmy_band_->TestFlag(RadiationFlags::TimeDependent)) { - Real time = pmb->pmy_mesh->time; - Real lat, lon; - - CubedSphereUtility::get_latlon_on_sphere(&lat, &lon, pmb, k, j, pmb->is); - - ray = planet->ParentZenithAngle(time, lat, lon); - dist_au = planet->ParentDistanceInAu(time); - } else { // constant zenith angle - if (pmy_band_->HasPar("umu0")) { - ray.mu = pmy_band_->GetPar("umu0"); - } - - if (pmy_band_->HasPar("phi0")) { - ray.phi = pmy_band_->GetPar("phi0"); - } - - if (pmy_band_->HasPar("dist_au")) { - dist_au = pmy_band_->GetPar("dist_au"); - } - } - - if (ds_.flag.ibcnd != 0) { - throw ValueError("RTSolverDisort::CalRadtranFlux", "ibcnd", ds_.flag.ibcnd, - 0); - } - - // pack temperature - if (ds_.flag.planck) { - pmy_band_->packTemperature(); - } - - // pack spectral properties - pmy_band_->packSpectralProperties(); - - ds_.bc.umu0 = ray.mu > 1.E-3 ? ray.mu : 1.E-3; - ds_.bc.phi0 = ray.phi; - - ds_.wvnmlo = wmin; - ds_.wvnmhi = wmax; - - if (pmy_band_->TestFlag(RadiationFlags::BroadBand)) { - // stellar source function overrides fbeam - if (pmy_band_->HasPar("S0")) { - ds_.bc.fbeam = pmy_band_->GetPar("S0"); - } else if (pmy_band_->HasPar("temp0")) { - Real temp0 = pmy_band_->GetPar("temp0"); - ds_.bc.fbeam = Constants::stefanBoltzmann * pow(temp0, 4); - } else if (planet && planet->HasParentFlux()) { - ds_.bc.fbeam = planet->ParentInsolationFlux(wmin, wmax, 1.); - } else { - ds_.bc.fbeam = 0.; - } - ds_.bc.fbeam /= dist_au * dist_au; - } - - pmb->pcoord->Face1Area(k, j, pmb->is, pmb->ie + 1, farea_); - pmb->pcoord->CellVolume(k, j, pmb->is, pmb->ie, vol_); - - auto &&uphi_umu = RadiationHelper::get_direction_grids(pmy_band_->rayOutput_); - auto &uphi = uphi_umu.first; - auto &umu = uphi_umu.second; - - if (umu.size() > 0) { - SetUserCosinePolarAngle(umu); - - for (int i = 0; i < umu.size(); ++i) { - dir_axis_[i] = umu[i]; - } - } - - if (uphi.size() > 0) { - SetUserAzimuthalAngle(uphi); - - for (int i = 0; i < uphi.size(); ++i) { - dir_axis_[ds_.numu + i] = uphi[i]; - } - } - - // set the surface temperature from the user input variables - // FIXME(cli) - // AthenaArray &ts = pmb->ruser_meshblock_data[1]; - // ds_.bc.btemp = ts(j); -} - -void RadiationBand::RTSolverDisort::CalBandFlux(MeshBlock const *pmb, int k, - int j) { - Real dist_au = 1.; - auto planet = pmb->pimpl->planet; - - if (planet && pmy_band_->TestFlag(RadiationFlags::TimeDependent)) { - dist_au = planet->ParentDistanceInAu(pmb->pmy_mesh->time); - } else if (pmy_band_->HasPar("dist_au")) { - dist_au = pmy_band_->GetPar("dist_au"); - } - - // bflxup and bflxdn has been reset in RadiationBand::CalBandFlux - - // loop over spectral grids in the band - bool override_with_stellar_spectra = false; - if (!pmy_band_->TestFlag(RadiationFlags::BroadBand) && - !pmy_band_->HasPar("S0") && !pmy_band_->HasPar("temp0") && planet && - planet->HasParentFlux()) { - override_with_stellar_spectra = true; - } - - pmy_band_->pexv->GatherAll(pmb); - if (ds_.flag.planck) { - pmy_band_->unpackTemperature(&ds_); - } - - int b = 0; - for (auto &spec : pmy_band_->pgrid_->spec) { - if (override_with_stellar_spectra) { - // stellar source function - ds_.bc.fbeam = - planet->ParentInsolationFlux(spec.wav1, spec.wav2, dist_au); - // planck source function - ds_.wvnmlo = spec.wav1; - ds_.wvnmhi = spec.wav2; - } - - // transfer spectral grid data - pmy_band_->unpackSpectralProperties(b, &ds_); - - if (ds_.nmom_nstr == 2) { - int nlay = ds_.nlyr; - int ierror[TWOSTR_NERR]; - std::vector gg(nlay, 0.); - for (int i = 0; i < nlay; ++i) gg[i] = ds_.pmom[i * (ds_.nmom_nstr + 1) + 1]; - c_twostr(&ds_, &ds_out_, 0, gg.data(), ierror, 1000.); - } else { - c_disort(&ds_, &ds_out_); - } - - // add spectral bin flux - addDisortFlux(pmb->pcoord, b++, k, j, pmb->is, pmb->ie + 1); - } -} - -void RadiationBand::RTSolverDisort::addDisortFlux(Coordinates const *pcoord, - int b, int k, int j, int il, - int iu) { - auto &bflxup = pmy_band_->bflxup; - auto &bflxdn = pmy_band_->bflxdn; - - auto &flxup = pmy_band_->flxup_; - auto &flxdn = pmy_band_->flxdn_; - auto const &spec = pmy_band_->pgrid_->spec; - - int rank_in_column = pmy_band_->pexv->GetRankInGroup(); - - /// accumulate flux from spectral bins - for (int i = il; i <= iu; ++i) { - int m = ds_.nlyr - (rank_in_column * (iu - il) + i - il); - //! \bug does not work for spherical geometry, need to scale area using - //! farea(il)/farea(i) - // flux up - flxup(b, k, j, i) = ds_out_.rad[m].flup; - - //! \bug does not work for spherical geomtry, need to scale area using - //! farea(il)/farea(i) - // flux down - flxdn(b, k, j, i) = ds_out_.rad[m].rfldir + ds_out_.rad[m].rfldn; - - bflxup(k, j, i) += spec[b].wght * flxup(b, k, j, i); - bflxdn(k, j, i) += spec[b].wght * flxdn(b, k, j, i); - } - - //! \note Spherical correction by XIZ - //! xiz 2022 flux scaling so that the heating rate is the same as the - //! plane-parallel scheme volheating scaling: first calculate flux divergence - //! from DISORT using Plane-parallel in a cell then mulitpled by the cell - //! volume divided by dx1f then solve for F using F1*S1-F2*S2 = volheating - //! the top fluxes are the still the same as the plane-paralell values - Real volh; - Real bflxup_iu = bflxup(k, j, iu); - Real bflxdn_iu = bflxdn(k, j, iu); - - for (int i = iu - 1; i >= il; --i) { - // upward - volh = (bflxup_iu - bflxup(k, j, i)) / pcoord->dx1f(i) * vol_(i); - bflxup_iu = bflxup(k, j, i); - bflxup(k, j, i) = (bflxup(k, j, i + 1) * farea_(i + 1) - volh) / farea_(i); - - // downward - volh = (bflxdn_iu - bflxdn(k, j, i)) / pcoord->dx1f(i) * vol_(i); - bflxdn_iu = bflxdn(k, j, i); - bflxdn(k, j, i) = (bflxdn(k, j, i + 1) * farea_(i + 1) - volh) / farea_(i); - } - - /* - for (int i = iu; i >= il; --i) { - std::cout << "i: " << iu-i+1 <<" flxup: " << bflxup(k, j, i) << " flxdn: " << bflxdn(k, j, i) << " fluxdiff: " << bflxup(k, j, i) - bflxdn(k, j, i) << std::endl; - } -*/ - -} - -void RadiationBand::RTSolverDisort::CalBandRadiance(MeshBlock const *pmb, int k, - int j) { - if (ds_.flag.onlyfl) { - throw RuntimeError("RTSolverDisort::CalBandRadiance", - "Radiance calculation disabled"); - } - - if (ds_.ntau != 1) { - throw RuntimeError("RTSolverDisort::CalBandRadiance", - "Only toa radiance (ds.ntau = 1) is supported"); - } - - int nrays = ds_.nphi * ds_.numu; - - if (nrays < pmy_band_->GetNumOutgoingRays()) { - throw RuntimeError("RTSolverDisort::CalBandRadiance", - "Number of outgoing rays more than DISORT can host"); - } - - // toa has been reset in RadiationBand::CalBandRadiance - - Real dist_au = 1.; - auto planet = pmb->pimpl->planet; - - if (planet && pmy_band_->TestFlag(RadiationFlags::TimeDependent)) { - dist_au = pmb->pimpl->planet->ParentDistanceInAu(pmb->pmy_mesh->time); - } else if (pmy_band_->HasPar("dist_au")) { - dist_au = pmy_band_->GetPar("dist_au"); - } - - pmy_band_->pexv->GatherAll(pmb); - if (ds_.flag.planck) { - pmy_band_->unpackTemperature(&ds_); - } - - // loop over spectral grids in the band - int b = 0; - for (auto &spec : pmy_band_->pgrid_->spec) { - // override source function for non-broadband radiation - if (!(pmy_band_->TestFlag(RadiationFlags::BroadBand))) { - // stellar source function - if (planet) { - ds_.bc.fbeam = - planet->ParentInsolationFlux(spec.wav1, spec.wav2, dist_au); - } else { - ds_.bc.fbeam = 0.0; - } - // planck source function - ds_.wvnmlo = spec.wav1; - ds_.wvnmhi = spec.wav2; - } - - // transfer spectral grid data - pmy_band_->unpackSpectralProperties(b, &ds_); - - // run disort - c_disort(&ds_, &ds_out_); - - // add spectral bin radiance - addDisortRadiance(b++, k, j); - } -} - -void RadiationBand::RTSolverDisort::addDisortRadiance(int b, int k, int j) { - auto &toa = pmy_band_->toa_; - auto &btoa = pmy_band_->btoa; - auto &spec = pmy_band_->pgrid_->spec; - auto &rayout = pmy_band_->rayOutput_; - - for (int n = 0; n < pmy_band_->GetNumOutgoingRays(); ++n) { - Real val; - Real coor[2] = {rayout[n].mu, rayout[n].phi}; - interpn(&val, coor, ds_out_.uu, dir_axis_.data(), dir_dim_, 2, 1); - toa(b, n, k, j) = val; - btoa(n, k, j) += spec[b].wght * toa(b, n, k, j); - } -} - -#endif // RT_DISORT diff --git a/src/harp/rt_solver_toon.cpp b/src/harp/rt_solver_toon.cpp deleted file mode 100644 index dd3c3235..00000000 --- a/src/harp/rt_solver_toon.cpp +++ /dev/null @@ -1,797 +0,0 @@ -// RT solvers based on Toon 1989 method by Xi Zhang -// Reference: Toon, O.B., 1989, JGR, 94, 16287-16301. - -#include -#include -#include -#include -#include - -// external -#include - -// athena -#include -#include - -// application -#include -#include - -// climath -#include - -// canoe -#include -#include - -// astro -#include - -// exo3 -#include -#include - -// harp -#include "radiation.hpp" -#include "rt_solvers.hpp" - -#ifdef RT_DISORT - -RadiationBand::RTSolverToon::RTSolverToon(RadiationBand *pmy_band, YAML::Node const &rad) - : RTSolver(pmy_band, "Toon") { - Application::Logger app("harp"); - app->Log("Toon solver initialized for band " + pmy_band_->GetName()); -} - - -//! \todo update based on band outdir -void RadiationBand::RTSolverToon::Resize(int nlyr, int nstr) { - RadiationBand::RTSolver::Resize(nlyr, nstr); - Unseal(); - SetAtmosphereDimension(nlyr, nstr, nstr); - Seal(); -} - -void RadiationBand::RTSolverToon::Prepare(MeshBlock const *pmb, int k, int j) { - auto &wmin = pmy_band_->wrange_.first; - auto &wmax = pmy_band_->wrange_.second; - - Real dist_au = 1.0; - Direction ray = pmb->pimpl->prad->GetRayInput(0); - auto planet = pmb->pimpl->planet; - - if (planet && pmy_band_->TestFlag(RadiationFlags::TimeDependent)) { - Real time = pmb->pmy_mesh->time; - Real lat, lon; - - CubedSphereUtility::get_latlon_on_sphere(&lat, &lon, pmb, k, j, pmb->is); - - ray = planet->ParentZenithAngle(time, lat, lon); - dist_au = planet->ParentDistanceInAu(time); - } else { - if (pmy_band_->HasPar("umu0")) { - ray.mu = pmy_band_->GetPar("umu0"); - } - - if (pmy_band_->HasPar("phi0")) { - ray.phi = pmy_band_->GetPar("phi0"); - } - - if (pmy_band_->HasPar("dist_au")) { - dist_au = pmy_band_->GetPar("dist_au"); - } - } - - // pack temperature - if (pmy_band_->TestFlag(RadiationFlags::ThermalEmission)) { - pmy_band_->packTemperature(); - } - - // pack spectral properties - pmy_band_->packSpectralProperties(); - ds_.bc.umu0 = ray.mu > 1.E-3 ? ray.mu : 1.E-3; - - if (pmy_band_->TestFlag(RadiationFlags::BroadBand)) { - // stellar source function overrides fbeam - if (pmy_band_->HasPar("S0")) { - ds_.bc.fbeam = pmy_band_->GetPar("S0"); - } else if (pmy_band_->HasPar("temp0")) { - Real temp0 = pmy_band_->GetPar("temp0"); - ds_.bc.fbeam = Constants::stefanBoltzmann * pow(temp0, 4); - } else if (planet && planet->HasParentFlux()) { - ds_.bc.fbeam = planet->ParentInsolationFlux(wmin, wmax, 1.); - } else { - ds_.bc.fbeam = 0.; - } - ds_.bc.fbeam /= dist_au * dist_au; - } - - pmb->pcoord->Face1Area(k, j, pmb->is, pmb->ie + 1, farea_); - pmb->pcoord->CellVolume(k, j, pmb->is, pmb->ie, vol_); -} - -void RadiationBand::RTSolverToon::CalBandFlux(MeshBlock const *pmb, int k, int j) { - Real dist_au = 1.0; - auto planet = pmb->pimpl->planet; - - if (planet && pmy_band_->TestFlag(RadiationFlags::TimeDependent)) { - dist_au = planet->ParentDistanceInAu(pmb->pmy_mesh->time); - } else if (pmy_band_->HasPar("dist_au")) { - dist_au = pmy_band_->GetPar("dist_au"); - } - - // loop over spectral grids in the band - bool override_with_stellar_spectra = false; - if (!pmy_band_->TestFlag(RadiationFlags::BroadBand) && - !pmy_band_->HasPar("S0") && !pmy_band_->HasPar("temp0") && planet && - planet->HasParentFlux()) { - override_with_stellar_spectra = true; - } - - pmy_band_->pexv->GatherAll(pmb); - if (pmy_band_->TestFlag(RadiationFlags::ThermalEmission)) { - pmy_band_->unpackTemperature(&ds_); - } - - int b = 0; - for (auto &spec : pmy_band_->pgrid_->spec) { - - if (override_with_stellar_spectra) { - // stellar source function - ds_.bc.fbeam = - planet->ParentInsolationFlux(spec.wav1, spec.wav2, dist_au); - } - - // Transfer spectral grid data - pmy_band_->unpackSpectralProperties(b, &ds_); - - // Determine flux for each spectral bin using Toon method - int nlay = ds_.nlyr; - Eigen::VectorXd w0(nlay), g(nlay); - Eigen::VectorXd tau_cum(nlay+1); - tau_cum.setZero(); - - for (int i = 0; i < nlay; ++i) { - tau_cum(i+1) = tau_cum(i) + ds_.dtauc[i]; - w0(i) = ds_.ssalb[i]; - g(i) = ds_.pmom[i * (ds_.nmom_nstr + 1) + 1]; - } - - Eigen::VectorXd flux_up(nlay), flux_down(nlay); - double mu0 = ds_.bc.umu0; - double surface_albedo = ds_.bc.albedo; - double Finc = ds_.bc.fbeam; - -// Call shortwave or longwave solver based on the type of radiation - if (pmy_band_->TestFlag(RadiationFlags::StellarBeam)) { - Eigen::VectorXd mu0_in(nlay+1); - mu0_in.setConstant(mu0); - toonShortwaveSolver(nlay, Finc, mu0_in, tau_cum, w0, g, surface_albedo, flux_up, flux_down); - - } else if (pmy_band_->TestFlag(RadiationFlags::ThermalEmission)) { - Eigen::VectorXd temp(nlay + 1); - Eigen::VectorXd be(nlay + 1); - for (int i = 0; i < nlay + 1; ++i) { - temp(i) = ds_.temper[i]; - be(i) = BB_integrate(ds_.temper[i], spec.wav1, spec.wav2); - } - double surface_emissivity = ds_.bc.temis; - toonLongwaveSolver(nlay, be, tau_cum, w0, g, surface_emissivity, flux_up, flux_down); - } - - // add spectral bin flux - addToonFlux(pmb->pcoord, b++, k, j, pmb->is, pmb->ie + 1, flux_up, flux_down); - } -} - -void RadiationBand::RTSolverToon::addToonFlux(Coordinates const *pcoord, int b, int k, int j, int il, int iu, - const Eigen::VectorXd &flux_up, const Eigen::VectorXd &flux_down) { - auto &bflxup = pmy_band_->bflxup; - auto &bflxdn = pmy_band_->bflxdn; - - auto &flxup = pmy_band_->flxup_; - auto &flxdn = pmy_band_->flxdn_; - auto const &spec = pmy_band_->pgrid_->spec; - - int rank_in_column = pmy_band_->pexv->GetRankInGroup(); - - // Accumulate flux from spectral bins - for (int i = il; i <= iu; ++i) { - int m = ds_.nlyr - (rank_in_column * (iu - il) + i - il); - // Flux up - flxup(b, k, j, i) = flux_up(m); - // Flux down - flxdn(b, k, j, i) = flux_down(m); - - bflxup(k, j, i) += spec[b].wght * flxup(b, k, j, i); - bflxdn(k, j, i) += spec[b].wght * flxdn(b, k, j, i); - } - - // Spherical correction - Real volh; - Real bflxup_iu = bflxup(k, j, iu); - Real bflxdn_iu = bflxdn(k, j, iu); - - for (int i = iu - 1; i >= il; --i) { - // Upward - volh = (bflxup_iu - bflxup(k, j, i)) / pcoord->dx1f(i) * vol_(i); - bflxup_iu = bflxup(k, j, i); - bflxup(k, j, i) = (bflxup(k, j, i + 1) * farea_(i + 1) - volh) / farea_(i); - - // Downward - volh = (bflxdn_iu - bflxdn(k, j, i)) / pcoord->dx1f(i) * vol_(i); - bflxdn_iu = bflxdn(k, j, i); - bflxdn(k, j, i) = (bflxdn(k, j, i + 1) * farea_(i + 1) - volh) / farea_(i); - } - - /* - for (int i = iu; i >= il; --i) { - std::cout << "i: " << iu-i+1 <<" flxup: " << bflxup(k, j, i) << " flxdn: " << bflxdn(k, j, i) << " fluxdiff: " << bflxup(k, j, i) - bflxdn(k, j, i) << std::endl; - } -*/ -} - - -// Toon 1989 shortwave solver, based on Elsie Lee's implementation, which was based on CHIMERA code by Mike Line. Reference: Toon, O.B., 1989, JGR, 94, 16287-16301. -void RadiationBand::RTSolverToon::toonShortwaveSolver(int nlay, double F0_in, - const Eigen::VectorXd& mu_in, - const Eigen::VectorXd& tau_in, - const Eigen::VectorXd& w_in, - const Eigen::VectorXd& g_in, - double w_surf_in, - Eigen::VectorXd& flx_up, - Eigen::VectorXd& flx_down) - { - int nlev = nlay + 1; - // Input validation - if(mu_in.size() != static_cast(nlev) || - tau_in.size() != static_cast(nlev) || - w_in.size() != static_cast(nlay) || - g_in.size() != static_cast(nlay)) - { - throw std::invalid_argument("Input vectors have incorrect sizes."); - } - - // Initialize output flux arrays - flx_down = Eigen::VectorXd::Zero(nlev); - flx_up = Eigen::VectorXd::Zero(nlev); - - // Constants - const double sqrt3 = std::sqrt(3.0); - const double sqrt3d2 = sqrt3 / 2.0; - const double bsurf = 0.0; - const double btop = 0.0; - - // Check if all albedos are effectively zero - bool all_w0_zero = (w_in.array() <= 1.0e-12).all(); - - if(all_w0_zero){ - // Direct beam only - double mu_top = mu_in(nlev - 1); - double mu_first = mu_in(0); - - if(mu_top > 0.0){ - if(std::abs(mu_top - mu_first) < 1e-12){ - // No zenith correction, use regular method - flx_down = F0_in * mu_top * (-tau_in.array() / mu_top).exp(); - } - else{ - // Zenith angle correction using cumulative transmission - Eigen::VectorXd cum_trans(nlev); - cum_trans(0) = tau_in(0) / mu_in(0); - for(int k = 1; k < nlev; ++k){ - cum_trans(k) = cum_trans(k - 1) + (tau_in(k) - tau_in(k - 1)) / mu_in(k); - } - flx_down = F0_in * mu_top * (-cum_trans.array()).exp(); - } - // Adjust the downward flux at the surface layer for surface albedo - flx_down(nlev - 1) *= (1.0 - w_surf_in); - } - // Upward flux remains zero - return; - } - - // Delta Eddington scaling - Eigen::VectorXd w0 = ((1.0 - g_in.array().square()) * w_in.array()) / (1.0 - w_in.array() * g_in.array().square()); - Eigen::VectorXd dtau = (1.0 - w_in.array() * g_in.array().square()) - .cwiseProduct( (tau_in.segment(1, nlay) - tau_in.head(nlay)).array() ) - .matrix(); - Eigen::VectorXd hg = g_in.array() / (1.0 + g_in.array()); - - // Initialize tau_total - Eigen::VectorXd tau_total(nlev); - tau_total(0) = 0.0; - for(int k = 0; k < nlay; ++k){ - tau_total(k + 1) = tau_total(k) + dtau(k); - } - - // Compute g1, g2, g3, g4 - Eigen::VectorXd g1 = sqrt3d2 * (2.0 - w0.array() * (1.0 + hg.array())); - Eigen::VectorXd g2 = (sqrt3d2 * w0.array()) * (1.0 - hg.array()); - // Prevent division by zero - for(int i = 0; i < nlay; ++i){ - if(std::abs(g2(i)) < 1.0e-10){ - g2(i) = 1.0e-10; - } - } - // Compute mu_zm at midpoints - Eigen::VectorXd mu_zm(nlay); - mu_zm = (mu_in.head(nlay) + mu_in.tail(nlay)) / 2.0; - Eigen::VectorXd g3 = (1.0 - sqrt3 * hg.array() * mu_zm.array()) / 2.0; - Eigen::VectorXd g4 = 1.0 - g3.array(); - - // Compute lam and gam - Eigen::VectorXd lam = (g1.array().square() - g2.array().square()).sqrt(); - Eigen::VectorXd gam = (g1.array() - lam.array()) / g2.array(); - - // Compute denom and handle denom == 0 - Eigen::VectorXd denom = lam.array().square() - (1.0 / (mu_in(nlev - 1) * mu_in(nlev - 1))); - for(int i = 0; i < nlay; ++i){ - if(std::abs(denom(i)) < 1e-10){ - denom(i) = 1.0e-10; - } - } - - // Compute Am and Ap - Eigen::VectorXd Am = F0_in * w0.array() * (g4.array() * (g1.array() + 1.0 / mu_in(nlev - 1)) + g2.array() * g3.array()) / denom.array(); - Eigen::VectorXd Ap = F0_in * w0.array() * (g3.array() * (g1.array() - 1.0 / mu_in(nlev - 1)) + g2.array() * g4.array()) / denom.array(); - - // Compute Cpm1 and Cmm1 at the top of the layer - Eigen::VectorXd Cpm1 = Ap.array() * (-tau_total.head(nlay).array() / mu_in(nlev - 1)).exp(); - Eigen::VectorXd Cmm1 = Am.array() * (-tau_total.head(nlay).array() / mu_in(nlev - 1)).exp(); - - // Compute Cp and Cm at the bottom of the layer - Eigen::VectorXd Cp = Ap.array() * (-tau_total.segment(1, nlay).array() / mu_in(nlev - 1)).exp(); - Eigen::VectorXd Cm = Am.array() * (-tau_total.segment(1, nlay).array() / mu_in(nlev - 1)).exp(); - - // Compute exponential terms, clamped to prevent overflow - Eigen::VectorXd exptrm = (lam.array() * dtau.array()).min(35.0); - Eigen::VectorXd Ep = exptrm.array().exp(); - Eigen::VectorXd Em = 1.0 / Ep.array(); - Eigen::VectorXd E1 = (Ep.array() + gam.array() * Em.array()).matrix(); - Eigen::VectorXd E2 = (Ep.array() - gam.array() * Em.array()).matrix(); - Eigen::VectorXd E3 = (gam.array() * Ep.array() + Em.array()).matrix(); - Eigen::VectorXd E4 = (gam.array() * Ep.array() - Em.array()).matrix(); - - // Initialize Af, Bf, Cf, Df - int l = 2 * nlay; - Eigen::VectorXd Af_vec = Eigen::VectorXd::Zero(l); - Eigen::VectorXd Bf_vec = Eigen::VectorXd::Zero(l); - Eigen::VectorXd Cf_vec = Eigen::VectorXd::Zero(l); - Eigen::VectorXd Df_vec = Eigen::VectorXd::Zero(l); - - // Boundary conditions at the top - Af_vec(0) = 0.0; - Bf_vec(0) = gam(0) + 1.0; - Cf_vec(0) = gam(0) - 1.0; - Df_vec(0) = btop - Cmm1(0); - for(int i =1, n=1; i < l -1; i +=2, ++n){ - if(n >= nlay){ - throw std::out_of_range("Index out of range in sw_Toon89 Af, Bf, Cf, Df population."); - } - Af_vec(i) = (E1(n-1) + E3(n-1)) * (gam(n) - 1.0); - Bf_vec(i) = (E2(n-1) + E4(n-1)) * (gam(n) - 1.0); - Cf_vec(i) = 2.0 * (1.0 - gam(n) * gam(n)); - Df_vec(i) = (gam(n) - 1.0) * (Cpm1(n) - Cp(n - 1)) + (1.0 - gam(n)) * (Cm(n - 1) - Cmm1(n)); - } - - // Populate Af, Bf, Cf, Df for even indices - // Start from n=1 to avoid negative indexing (Cp(n-1) when n=0) - for(int i = 2, n=1; i < l -1; i +=2, ++n){ - if(n >= nlay){ - throw std::out_of_range("Index out of range in sw_Toon89 Af, Bf, Cf, Df population."); - } - Af_vec(i) = 2.0 * (1.0 - gam(n) * gam(n)); - Bf_vec(i) = (E1(n-1) - E3(n-1)) * (1.0 + gam(n)); - Cf_vec(i) = (E1(n-1) + E3(n-1)) * (gam(n) - 1.0); - Df_vec(i) = E3(n-1) * (Cpm1(n) - Cp(n - 1)) + E1(n-1) * (Cm(n - 1) - Cmm1(n)); - } - - // Boundary conditions at l (last index) - Af_vec(l - 1) = E1(nlay - 1) - w_surf_in * E3(nlay - 1); - Bf_vec(l - 1) = E2(nlay - 1) - w_surf_in * E4(nlay - 1); - Cf_vec(l - 1) = 0.0; - Df_vec(l - 1) = bsurf - Cp(nlay - 1) + w_surf_in * Cm(nlay - 1); - - // Prepare a, b, c, d for the solver - Eigen::VectorXd a_tridiag = Af_vec.segment(1, l - 1); - Eigen::VectorXd b_tridiag = Bf_vec; - Eigen::VectorXd c_tridiag = Cf_vec.segment(0, l - 1); - Eigen::VectorXd d_tridiag = Df_vec; - - // Solve the tridiagonal system - Eigen::VectorXd xk = tridiagonal_solver(a_tridiag, b_tridiag, c_tridiag, d_tridiag); - - // Compute xk1 and xk2 from xk - Eigen::VectorXd xk1(nlay); - Eigen::VectorXd xk2(nlay); - for(int idx = 0; idx < nlay; ++idx){ - int two_n = 2 * idx; - if(two_n + 1 >= xk.size()){ - throw std::out_of_range("Index out of range when accessing xk."); - } - xk1(idx) = xk(two_n) + xk(two_n + 1); - xk2(idx) = xk(two_n) - xk(two_n + 1); - if(std::abs(xk2(idx) / xk(two_n)) < 1e-30){ - xk2(idx) = 0.0; - } - } - - // Populate flx_up and flx_down for layers 1 to nlay - flx_up.head(nlay) = (xk1.array() + gam.array() * xk2.array() + Cpm1.array()).matrix(); - flx_down.head(nlay) = (xk1.array() * gam.array() + xk2.array() + Cmm1.array()).matrix(); - - // Compute flx_up and flx_down at level nlev - flx_up(nlev - 1) = xk1(nlay - 1) * std::exp(1.0) + gam(nlay - 1) * xk2(nlay - 1) * std::exp(-1.0) + Cp(nlay - 1); - flx_down(nlev - 1) = xk1(nlay - 1) * std::exp(1.0) * gam(nlay - 1) + xk2(nlay - 1) * std::exp(-1.0) + Cm(nlay - 1); - - // Compute dir flux - Eigen::VectorXd dir = Eigen::VectorXd::Zero(nlev); - double mu_top_nonzero = mu_in(nlev - 1); - double mu_first_nonzero = mu_in(0); - if(std::abs(mu_top_nonzero - mu_first_nonzero) < 1e-12){ - // No zenith correction - dir = F0_in * mu_top_nonzero * (-tau_in.array() / mu_top_nonzero).exp(); - } - else{ - // Zenith angle correction - Eigen::VectorXd cum_trans(nlev); - cum_trans(0) = tau_total(0) / mu_in(0); - for(int k = 1; k < nlev; ++k){ - cum_trans(k) = cum_trans(k - 1) + (tau_total(k) - tau_total(k - 1)) / mu_in(k); - } - dir = F0_in * mu_top_nonzero * (-cum_trans.array()).exp(); - } - - // Adjust the downward flux at the surface layer for surface albedo - dir(nlev - 1) *= (1.0 - w_surf_in); - -//for(int i=0; i nlay){ - throw std::out_of_range("Index out of range in lw_Toon89 Af, Bf, Cf, Df population."); - } - Af_vec(i) = (E1(n-1) + E3(n-1)) * (gam(n) - 1.0); - Bf_vec(i) = (E2(n-1) + E4(n-1)) * (gam(n) - 1.0); - Cf_vec(i) = 2.0 * (1.0 - gam(n) * gam(n)); - Df_vec(i) = (gam(n) - 1.0) * (Cpm1(n) - Cp(n-1)) + (1.0 - gam(n)) * (Cm(n-1) - Cmm1(n)); - } - - // Populate Af, Bf, Cf, Df for even indices (i=2,4,...) - n = 1; // Reset n to 1 for even indices - for(int i=2; i < l -1; i +=2, ++n){ - if(n > nlay){ - throw std::out_of_range("Index out of range in lw_Toon89 Af, Bf, Cf, Df population."); - } - Af_vec(i) = 2.0 * (1.0 - gam(n-1) * gam(n-1)); - Bf_vec(i) = (E1(n-1) - E3(n-1)) * (1.0 + gam(n)); - Cf_vec(i) = (E1(n-1) + E3(n-1)) * (gam(n) - 1.0); - Df_vec(i) = E3(n-1) * (Cpm1(n) - Cp(n-1)) + E1(n-1) * (Cm(n-1) - Cmm1(n)); - } - - // Boundary conditions at the last index - Af_vec(l -1) = E1(nlay -1) - a_surf_in * E3(nlay -1); - Bf_vec(l -1) = E2(nlay -1) - a_surf_in * E4(nlay -1); - Cf_vec(l -1) = 0.0; - Df_vec(l -1) = Bsurf - Cp(nlay -1) + a_surf_in * Cm(nlay -1); - - // Prepare a, b, c, d for the tridiagonal solver - Eigen::VectorXd a_tridiag = Af_vec.segment(1, l -1); - Eigen::VectorXd b_tridiag = Bf_vec; - Eigen::VectorXd c_tridiag = Cf_vec.segment(0, l -1); - Eigen::VectorXd d_tridiag = Df_vec; - - // Solve the tridiagonal system - Eigen::VectorXd xkk = tridiagonal_solver(a_tridiag, b_tridiag, c_tridiag, d_tridiag); - - // Compute xk1 and xk2 from xkk - Eigen::VectorXd xk1 = Eigen::VectorXd::Zero(nlay); - Eigen::VectorXd xk2 = Eigen::VectorXd::Zero(nlay); - for(int idx=0; idx < nlay; ++idx){ - int two_n = 2 * idx; - if(two_n +1 >= xkk.size()){ - throw std::out_of_range("Index out of range when accessing xk."); - } - xk1(idx) = xkk(two_n) + xkk(two_n +1); - xk2(idx) = xkk(two_n) - xkk(two_n +1); - if(std::abs(xk2(idx) / xkk(two_n)) < 1e-30){ - xk2(idx) = 0.0; - } - } - - // Apply conditional "where (w0 <= 1e-4)" - // Using Eigen's select function for element-wise conditional assignments - Eigen::ArrayXd mask = (w0.array() <= 1e-4).cast(); - - Eigen::VectorXd g_var_final = mask * 0.0 + (1.0 - mask) * (twopi * w0.array() * xk1.array() * (1.0 + hg.array() * alp.array()) / (1.0 + alp.array())); - Eigen::VectorXd h_var_final = mask * 0.0 + (1.0 - mask) * (twopi * w0.array() * xk2.array() * (1.0 - hg.array() * alp.array()) / (1.0 + alp.array())); - Eigen::VectorXd xj_final = mask * 0.0 + (1.0 - mask) * (twopi * w0.array() * xk1.array() * (1.0 - hg.array() * alp.array()) / (1.0 + alp.array())); - Eigen::VectorXd xk_final = mask * 0.0 + (1.0 - mask) * (twopi * w0.array() * xk2.array() * (1.0 + hg.array() * alp.array()) / (1.0 + alp.array())); - Eigen::VectorXd alpha1_vec_final = mask * (twopi * B0.array()) + (1.0 - mask) * (twopi * (B0.array() + B1.array() * (ubari * w0.array() * hg.array() / (1.0 - w0.array() * hg.array())))); - Eigen::VectorXd alpha2_vec_final = twopi * B1.array(); - Eigen::VectorXd sigma1_final = mask * (twopi * B0.array()) + (1.0 - mask) * (twopi * (B0.array() - B1.array() * (ubari * w0.array() * hg.array() / (1.0 - w0.array() * hg.array())))); - Eigen::VectorXd sigma2_final = alpha2_vec_final; - - // Compute obj, epp, em1, obj2, epp2 - Eigen::VectorXd obj_vec = (lam.array() * dtau.array()).min(35.0); - Eigen::VectorXd epp_vec_final = obj_vec.array().exp(); - Eigen::VectorXd em1_final = 1.0 / epp_vec_final.array(); - Eigen::VectorXd obj2_vec = (0.5 * lam.array() * dtau.array()).min(35.0); - Eigen::VectorXd epp2_vec_final = obj2_vec.array().exp(); - - // Initialize and resize flux arrays to ensure correct dimensions - flx_up = Eigen::VectorXd::Zero(nlev); - flx_down = Eigen::VectorXd::Zero(nlev); - - // Initialize lw_up_g and lw_down_g - Eigen::VectorXd lw_up_g = Eigen::VectorXd::Zero(nlev); - Eigen::VectorXd lw_down_g = Eigen::VectorXd::Zero(nlev); - - // Loop over m=1 to nmu (0 to nmu-1 in C++) - for(int m=0; m < nmu; ++m){ - // Compute em2 and em3 - Eigen::VectorXd em2_vec = (-dtau.array() / uarr(m)).exp(); - Eigen::VectorXd em3_vec = em1_final.array() * em2_vec.array(); - - // Downward loop - lw_down_g(0) = twopi * (1.0 - std::exp(-tautop / uarr(m))) * be(0); - for(int k=0; k < nlay; ++k){ - lw_down_g(k+1) = lw_down_g(k) * em2_vec(k) + - (xj_final(k) / (lam(k) * uarr(m) + 1.0)) * (epp_vec_final(k) - em2_vec(k)) + - (xk_final(k) / (lam(k) * uarr(m) - 1.0)) * (em2_vec(k) - em1_final(k)) + - sigma1_final(k) * (1.0 - em2_vec(k)) + - sigma2_final(k) * (uarr(m) * em2_vec(k) + dtau(k) - uarr(m)); - } - - // Upward loop - lw_up_g(nlev -1) = twopi * (Bsurf + B1(nlay -1) * uarr(m)); - for(int k = nlay -1; k >=0; --k){ - lw_up_g(k) = lw_up_g(k+1) * em2_vec(k) + - (g_var_final(k) / (lam(k) * uarr(m) - 1.0)) * (epp_vec_final(k) * em2_vec(k) - 1.0) + - (h_var_final(k) / (lam(k) * uarr(m) + 1.0)) * (1.0 - em3_vec(k)) + - alpha1_vec_final(k) * (1.0 - em2_vec(k)) + - alpha2_vec_final(k) * (uarr(m) - (dtau(k) + uarr(m)) * em2_vec(k)); - } - - // Sum up flux arrays with Gaussian quadrature weights - flx_down += (lw_down_g.array() * wuarr(m)).matrix(); - flx_up += (lw_up_g.array() * wuarr(m)).matrix(); - } -//for(int i=0; i (m) * v[i]; - exm *= ex; - d[i] += exm * (6.0 + mv * (6.0 + mv * (3.0 + mv))) / (m * m); - } - d[i] *= conc; - } - } - - double ans; - if (smallv == 2) { - // Both wavenumbers are small - ans = p[1] - p[0]; - } else if (smallv == 1) { - // One wavenumber is small, the other is large - ans = 1.0 - p[0] - d[1]; - } else { - // Both wavenumbers are large - ans = d[0] - d[1]; - } - - ans *= sigdpi * T * T * T * T; - - if (ans == 0.0) { - std::cerr << "BB_integrate: Warning - result is zero; possible underflow" << std::endl; - } - - return ans; -} - - -//Tridiagonal Solver using the Thomas Algorithm -inline Eigen::VectorXd RadiationBand::RTSolverToon::tridiagonal_solver( - const Eigen::VectorXd& a, const Eigen::VectorXd& b, - const Eigen::VectorXd& c, const Eigen::VectorXd& d) -{ - int l = b.size(); - if(a.size() != static_cast(l - 1) || c.size() != static_cast(l - 1) || d.size() != static_cast(l)){ - throw std::invalid_argument("Incorrect vector sizes for tridiagonal_solver."); - } - - Eigen::VectorXd c_prime(l - 1); - Eigen::VectorXd d_prime(l); - - // Forward sweep - c_prime(0) = c(0) / b(0); - d_prime(0) = d(0) / b(0); - for(int i = 1; i < l - 1; ++i){ - double denom = b(i) - a(i - 1) * c_prime(i - 1); - if(std::abs(denom) < 1e-12){ - throw std::runtime_error("Tridiagonal solver failed: near-zero denominator."); - } - c_prime(i) = c(i) / denom; - d_prime(i) = (d(i) - a(i - 1) * d_prime(i - 1)) / denom; - } - - // Last equation - double denom_last = b(l - 1) - a(l - 2) * c_prime(l - 2); - if(std::abs(denom_last) < 1e-12){ - throw std::runtime_error("Tridiagonal solver failed: near-zero denominator at last equation."); - } - d_prime(l - 1) = (d(l - 1) - a(l - 2) * d_prime(l - 2)) / denom_last; - - // Back substitution - Eigen::VectorXd x(l); - x(l - 1) = d_prime(l - 1); - for(int i = l - 2; i >= 0; --i){ - x(i) = d_prime(i) - c_prime(i) * x(i + 1); - } - - return x; -} - -#endif - diff --git a/src/harp/rt_solvers.hpp b/src/harp/rt_solvers.hpp deleted file mode 100644 index 5f069d39..00000000 --- a/src/harp/rt_solvers.hpp +++ /dev/null @@ -1,155 +0,0 @@ -#ifndef SRC_HARP_RT_SOLVERS_HPP_ -#define SRC_HARP_RT_SOLVERS_HPP_ - -// C/C++ -#include -#include // Needed for Eigen::VectorXd -#include -#include -#include -#include -#include -#include - -// athena -#include - -// application -#include -#include - -// canoe -#include -#include - -// cppdisort -#ifdef RT_DISORT -#include -#endif - -// harp -#include "radiation_band.hpp" - -class RadiationBand::RTSolver : public NamedGroup { - public: // constructor and destructor - RTSolver(RadiationBand *pmy_band, std::string name) - : NamedGroup(name), pmy_band_(pmy_band) { - Application::Logger app("harp"); - app->Log("Initialize RTSolver " + GetName()); - } - - virtual ~RTSolver() { - Application::Logger app("harp"); - app->Log("Destroy RTSolver " + GetName()); - } - - public: // member functions - //! \brief Prepare and seal the solver for the current column - virtual void Prepare(MeshBlock const *pmb, int k, int j) {} - - //! \brief Allocate memory for radiation solver - virtual void Resize(int nlyr, int nstr) { - farea_.DeleteAthenaArray(); - vol_.DeleteAthenaArray(); - - farea_.NewAthenaArray(nlyr + 2 * NGHOST); - vol_.NewAthenaArray(nlyr + 2 * NGHOST); - } - - public: // inbound functions - virtual void CalBandFlux(MeshBlock const *pmb, int k, int j) { - throw NotImplementedError("CalBandFlux not implemented."); - } - virtual void CalBandRadiance(MeshBlock const *pmb, int k, int j) { - throw NotImplementedError("CalBandRadiance not implemented."); - } - - protected: - RadiationBand *pmy_band_; - AthenaArray farea_, vol_; -}; - -class RadiationBand::RTSolverLambert : public RadiationBand::RTSolver { - public: // constructor and destructor - RTSolverLambert(RadiationBand *pmy_band, YAML::Node const &rad) - : RTSolver(pmy_band, "Lambert") {} - ~RTSolverLambert() {} - - public: // inbound functions - void CalBandRadiance(MeshBlock const *pmb, int k, int j) override; -}; - -#ifdef RT_DISORT -class RadiationBand::RTSolverDisort : public RadiationBand::RTSolver, - protected DisortWrapper { - public: // constructor and destructor - RTSolverDisort(RadiationBand *pmy_band, YAML::Node const &rad); - ~RTSolverDisort() {} - - public: // member functions - void Prepare(MeshBlock const *pmb, int k, int j) override; - void Resize(int nlyr, int nstr) override; - - public: // inbound functions - void CalBandFlux(MeshBlock const *pmb, int k, int j) override; - void CalBandRadiance(MeshBlock const *pmb, int k, int j) override; - - protected: - size_t dir_dim_[2]; - std::vector dir_axis_; - - void addDisortFlux(Coordinates const *pcoord, int n, int k, int j, int il, - int iu); - - void addDisortRadiance(int n, int k, int j); -}; - - -class RadiationBand::RTSolverToon : public RadiationBand::RTSolver, - protected DisortWrapper { - public: // Constructor and Destructor - RTSolverToon(RadiationBand *pmy_band, YAML::Node const &rad); - ~RTSolverToon() {} - - public: // Member Functions - void Prepare(MeshBlock const *pmb, int k, int j) override; - void Resize(int nlyr, int nstr) override; - - public: // Inbound Functions - void CalBandFlux(MeshBlock const *pmb, int k, int j) override; - - protected: - void addToonFlux(Coordinates const *pcoord, int b, int k, int j, int il, int iu, - const Eigen::VectorXd &flux_up, const Eigen::VectorXd &flux_down); - - void toonShortwaveSolver(int nlay, double F0_in, - const Eigen::VectorXd& mu_in, - const Eigen::VectorXd& tau_in, - const Eigen::VectorXd& w_in, - const Eigen::VectorXd& g_in, - double w_surf_in, - Eigen::VectorXd& flx_up, - Eigen::VectorXd& flx_down); - - void toonLongwaveSolver(int nlay, - const Eigen::VectorXd& be, - const Eigen::VectorXd& tau_in, - const Eigen::VectorXd& w_in, - const Eigen::VectorXd& g_in, - double a_surf_in, - Eigen::VectorXd& flx_up, - Eigen::VectorXd& flx_down); - - double BB_integrate(double T, double wn1, double wn2); - - private: - // Private Tridiagonal Solver using the Thomas Algorithm - inline Eigen::VectorXd tridiagonal_solver(const Eigen::VectorXd& a, // Sub-diagonal (size n-1) - const Eigen::VectorXd& b, // Main diagonal (size n) - const Eigen::VectorXd& c, // Super-diagonal (size n-1) - const Eigen::VectorXd& d); // Right-hand side (size n) -}; - -#endif - -#endif // SRC_HARP_RT_SOLVERS_HPP_ diff --git a/src/harp/set_spectral_properties.cpp b/src/harp/set_spectral_properties.cpp deleted file mode 100644 index b5b2b51e..00000000 --- a/src/harp/set_spectral_properties.cpp +++ /dev/null @@ -1,143 +0,0 @@ -// cnaoe -#include - -// athena -#include -#include -#include -#include -#include - -// canoe -#include -#include -#include - -// snap -#include -#include - -// opacity -#include -#include - -// harp -#include "radiation.hpp" -#include "radiation_band.hpp" - -// setting optical properties -void RadiationBand::SetSpectralProperties(AirColumn& ac, Real const* x1f, - Real gH0, int k, int j) { - int nspec = GetNumSpecGrids(); - int npmom = GetNumPhaseMoments(); - - // set tau, ssalb, pmom, etc... - tau_.ZeroClear(); - ssa_.ZeroClear(); - pmom_.ZeroClear(); - - std::vector mypmom(1 + npmom); - - for (int i = 0; i < ac.size(); ++i) { - auto& air = ac[i]; - air.ToMoleFraction(); - tem_[i] = air.w[IDN]; - - for (auto& a : absorbers_) { - for (int m = 0; m < nspec; ++m) { - auto& spec = pgrid_->spec[m]; - Real kcoeff; // Variable to store attenuation coefficient - if (auto* ck = dynamic_cast(a.get())) { - kcoeff = ck->GetAttenuation(m, air); - } else { - kcoeff = a->GetAttenuation(spec.wav1, spec.wav2, air); - } - Real dssalb = - a->GetSingleScatteringAlbedo(spec.wav1, spec.wav2, air) * kcoeff; - // tau - tau_(m, i) += kcoeff; - - // ssalb - ssa_(m, i) += dssalb; - // pmom - a->GetPhaseMomentum(mypmom.data(), spec.wav1, spec.wav2, air, npmom); - for (int p = 0; p <= npmom; ++p) pmom_(m, i, p) += mypmom[p] * dssalb; - - } - } - } - - // set temperature at cell interface - int il = NGHOST, iu = ac.size() - 1 - NGHOST; - temf_[il] = (3. * tem_[il] - tem_[il + 1]) / 2.; - temf_[il + 1] = (tem_[il] + tem_[il + 1]) / 2.; - for (int i = il + 2; i <= iu - 1; ++i) - temf_[i] = interp_cp4(tem_[i - 2], tem_[i - 1], tem_[i], tem_[i + 1]); - temf_[iu] = (tem_[iu] + tem_[iu - 1]) / 2.; - // temf_[iu + 1] = (3. * tem_[iu] - tem_[iu - 1]) / 2.; - temf_[iu + 1] = tem_[iu]; // isothermal top boundary - - for (int i = 0; i < il; ++i) temf_[i] = tem_[il]; - for (int i = iu + 2; i < ac.size(); ++i) temf_[i] = tem_[iu + 1]; - - bool error = false; - for (int i = 0; i < ac.size(); ++i) { - if (temf_[i] < 0.) { - temf_[i] = tem_[i]; - // error = true; - } - } - for (int i = il; i <= iu; ++i) { - if (tem_[i] < 0.) error = true; - } - if (error) { - for (int i = il; i <= iu; ++i) { - std::cout << "--- temf[" << i << "] = " << temf_[i] << std::endl; - std::cout << "tem[" << i << "] = " << tem_[i] << std::endl; - } - std::cout << "--- temf[" << iu + 1 << "] = " << temf_[iu + 1] << std::endl; - throw std::runtime_error("Negative temperature at cell interface"); - } - - // absorption coefficients -> optical thickness - for (int m = 0; m < nspec; ++m) { - for (int i = 0; i < ac.size(); ++i) { - if (tau_(m, i) > 1e-10 && ssa_(m, i) > 1e-10) { // has scattering - for (int p = 0; p <= npmom; ++p) pmom_(m, i, p) /= ssa_(m, i); - ssa_(m, i) /= tau_(m, i); - } else { - ssa_(m, i) = 0.; - pmom_(m, i, 0) = 1.; - for (int p = 1; p <= npmom; ++p) pmom_(m, i, p) = 0.; - } - -#ifdef HYDROSTATIC - auto pthermo = Thermodynamics::GetInstance(); - Real Rgas = get_rovrd(ac[i], pthermo->GetMuRatio()) * pthermo->GetRd(); - // TODO(cli) check this - // \delta z = \delt Z * (R T)/(g H0) - tau_(m, i) *= (x1f[i + 1] - x1f[i]) * (Rgas * tem_[i]) / (gH0); -#else - tau_(m, i) *= x1f[i + 1] - x1f[i]; -#endif - } - } - - // aggregated band properties - for (int i = 0; i < ac.size(); ++i) { - btau(k, j, i) = 0; - bssa(k, j, i) = 0; - for (int p = 0; p <= npmom; ++p) bpmom(p, k, j, i) = 0.; - - for (int m = 0; m < nspec; ++m) { - btau(k, j, i) += tau_(m, i); - bssa(k, j, i) += ssa_(m, i) * tau_(m, i); - for (int p = 0; p <= npmom; ++p) - bpmom(p, k, j, i) += pmom_(m, i, p) * ssa_(m, i); - } - - for (int p = 0; p <= npmom; ++p) bpmom(p, k, j, i) /= bssa(k, j, i); - bssa(k, j, i) /= btau(k, j, i); - btau(k, j, i) /= nspec; - } -} diff --git a/src/opacity/absorber.hpp b/src/opacity/absorber.hpp deleted file mode 100644 index 91067da0..00000000 --- a/src/opacity/absorber.hpp +++ /dev/null @@ -1,114 +0,0 @@ -#ifndef SRC_OPACITY_ABSORBER_HPP_ -#define SRC_OPACITY_ABSORBER_HPP_ - -// C/C++ -#include -#include -#include -#include - -// athena -#include - -// canoe -#include - -class AirParcel; -class SpectralBin; -namespace YAML { -class Node; -} - -//! \brief base class of all absorbers -class Absorber : public NamedGroup, - public ParameterGroup, - public SpeciesIndexGroup, - public StringReprGroup { - public: // constructor and destructor - Absorber(std::string name); - virtual ~Absorber(); - - public: // member functions - //! Set absorption model - void SetModel(std::string name) { model_name_ = name; } - - //! Combines SetOpacityFile() and LoadOpacity() - void LoadOpacityFromFile(std::string filename); - - //! Set opacity filename to internal variable, does not load opacity - void SetOpacityFile(std::string filename); - - //! Get opacity filename - std::string GetOpacityFile() const { return opacity_filename_; } - - //! Load opacity from internal variable - void LoadOpacity(int bid); - - //! Load absorption coefficient from file - virtual void LoadCoefficient(std::string fname, int bid) {} - - //! Get attenuation coefficient [1/m] - virtual Real GetAttenuation(Real wave1, Real wave2, - AirParcel const& var) const { - return 0.; - } - - //! Get attenuation coefficient for CK - virtual Real GetAttenuation(int m, AirParcel const& var) const { - return 0.; - } - - //! Get single scattering albedo [1] - virtual Real GetSingleScatteringAlbedo(Real wave1, Real wave2, - AirParcel const& var) const { - return 0.; - } - - //! Get phase function [1] - virtual void GetPhaseMomentum(Real* pp, Real wave1, Real wave2, - AirParcel const& var, int np) const {} - - virtual void CheckFail() const {} - - virtual void ModifySpectralGrid(std::vector& spec) const {} - - public: // StringRepr - std::string ToString() const override; - - protected: - //! absorption model model - std::string model_name_; - - //! opacity filename - std::string opacity_filename_; -}; - -using AbsorberPtr = std::shared_ptr; -using AbsorberContainer = std::vector; - -class AbsorberFactory { - public: - //! \brief Create an absorber from YAML node - //! - //! \param[in] my YAML node containing the current absorber - //! \param[in] band_name name of the radiation band - static AbsorberPtr CreateFrom(YAML::Node const& my, std::string band_name); - - //! \brief Create absorbers from YAML node - //! - //! \param[in] names names of absorbers - //! \param[in] band_name name of the radiation band - //! \param[in] rad YAML node containing the radiation input file - static AbsorberContainer CreateFrom(std::vector const& names, - std::string band_name, - YAML::Node const& rad); - - protected: - //! \brief Only create an absorber based on its name and class - //! - //! \param[in] name name of the absorber - //! \param[in] type class identifier of the absorber - static AbsorberPtr createAbsorberPartial(std::string name, std::string type); -}; - -#endif // SRC_OPACITY_ABSORBER_HPP_ diff --git a/src/opacity/absorber_ck.cpp b/src/opacity/absorber_ck.cpp deleted file mode 100644 index 158aca97..00000000 --- a/src/opacity/absorber_ck.cpp +++ /dev/null @@ -1,170 +0,0 @@ -// C/C++ -#include -#include -#include -#include -#include -#include -#include -#include -#include - -// canoe -#include -#include -#include - -// climath -#include - -// harp -#include - -// athena -#include - -// snap -#include - -// opacity -#include "absorber_ck.hpp" - -// netcdf -#ifdef NETCDFOUTPUT -extern "C" { -#include -} -#endif - -void AbsorberCK::LoadCoefficient(std::string fname, int) { -#ifdef NETCDFOUTPUT - int fileid, varid, err; - int dimids[10]; // Array to hold dimension IDs (assuming a maximum of 10 dimensions) - - size_t len_p, len_t, len_bin_centers, len_weights; - -std::map var_names = { - {"p", "p"}, // Pressure (logarithmic, in Pa) - {"t", "t"}, // Temperature (K) - {"bin_centers", "bin_centers"}, // Wavenumber bin centers (cm⁻¹) - {"bin_edges", "bin_edges"}, // Wavenumber bin edges (cm⁻¹) - {"weights", "weights"}, // Gaussian quadrature weights - {"gaussian_points", "gaussian_points"}, // Gaussian quadrature points - {"kcoeff", "kcoeff"} // Logarithmic absorption coefficients (log(m²/molecule)) -}; - - nc_open(fname.c_str(), NC_NOWRITE, &fileid); - - for (const auto& [role, var_name] : var_names) { - if (role == "kcoeff" || role == "bin_edges") continue; // Handle separately - - nc_inq_varid(fileid, var_name.c_str(), &varid); - nc_inq_vardimid(fileid, varid, dimids); - - // Assuming each variable is 1D, get the length of the first dimension - size_t dim_len; - nc_inq_dimlen(fileid, dimids[0], &dim_len); - - // Assign to len_ based on role - if (role == "p") { - len_p = dim_len; - len_[0] = len_p; - } - if (role == "t") { - len_t = dim_len; - len_[1] = len_t; - } - if (role == "bin_centers") { - len_bin_centers = dim_len; - len_[2] = len_bin_centers; - } - if (role == "weights") { - len_weights = dim_len; - len_[3] = len_weights; - } - } - - p_.resize(len_[0]); // Raw pressure - t_.resize(len_[1]); // Raw temperature - axis_.resize(len_[0] + len_[1]); // Pressure + Temperature - bin_centers_.resize(len_[2]); // Bin centers - weights_.resize(len_[3]); // Weights - bin_edges_.resize(len_[2] + 1); // Bin edges (len_bin_centers + 1) - kcoeff_.resize(len_[0] * len_[1] * len_[2] * len_[3]); // Absorption coefficients - - for (const auto& [role, var_name] : var_names) { - if (role == "kcoeff" || role == "bin_edges") continue; // Handle separately - - // Get the variable ID - nc_inq_varid(fileid, var_name.c_str(), &varid); - - // Read the data into the corresponding container - if (role == "p") { - nc_get_var_double(fileid, varid, axis_.data()); - nc_get_var_double(fileid, varid, p_.data()); - } - - if (role == "t") { - nc_get_var_double(fileid, varid, axis_.data() + len_p); - nc_get_var_double(fileid, varid, t_.data()); - } - if (role == "bin_centers") { - nc_get_var_double(fileid, varid, bin_centers_.data()); - } - if (role == "weights") { - nc_get_var_double(fileid, varid, weights_.data()); - } - } - - nc_inq_varid(fileid, var_names["bin_edges"].c_str(), &varid); - nc_get_var_double(fileid, varid, bin_edges_.data()); - - nc_inq_varid(fileid, var_names["kcoeff"].c_str(), &varid); - nc_get_var_double(fileid, varid, kcoeff_.data()); - - nc_close(fileid); - - // Optional: You can log some debug output to check that everything was loaded correctly - std::cout << "Loaded coefficient grid: " - << "Pressure levels = " << len_[0] << ", " - << "Temperature levels = " << len_[1] << ", " - << "Wavelength bins = " << len_[2] << ", " - << "Weights = " << len_[3] << std::endl; -#endif -} - -Real AbsorberCK::GetAttenuation(int m, AirParcel const& var) const { - int iw = m / len_[3]; // Wavelength bin index - int ig = m % len_[3]; // Weight index - Real val; - Real coord[2] = {log(var.w[IPR]), var.w[IDN]}; - - // Calculate the starting index for the 2D slice (p, t) at iw, ig - size_t base_index = iw * len_[3] * len_[0] * len_[1] + ig * len_[0] * len_[1]; - size_t len_2d[2] = {len_[0], len_[1]}; - - // Perform 2D interpolation using the flattened kcoeff_ - interpn(&val, coord, &kcoeff_[base_index], axis_.data(), len_2d, 2, 1); - - auto pthermo = Thermodynamics::GetInstance(); - Real dens = var.w[IPR] / (pthermo->GetRd() * var.w[IDN]); - - return exp(val) * dens; // ln(m^2/kg) -> 1/m -} - -void AbsorberCK::ModifySpectralGrid(std::vector& spec) const { - // Resize the spec grid to account for both wavelength bins and weights (Gauss points) - spec.resize(len_[2] * len_[3]); // len_[2] = nbin (wavelength bins), len_[3] = nwght (Gauss points) - - // Loop through both wavelength bins and weights - for (size_t iw = 0; iw < len_[2]; ++iw) { // len_[2] is the number of wavelength bins - for (size_t ig = 0; ig < len_[3]; ++ig) { // len_[3] is the number of weights (Gauss points) - size_t m = iw * len_[3] + ig; // Flattened index for the spectral grid - - // Use bin_edges for wav1 and wav2 - spec[m].wav1 = bin_edges_[iw]; // Lower boundary (wav1) - spec[m].wav2 = bin_edges_[iw + 1]; // Upper boundary (wav2) - spec[m].wght = weights_[ig]; // Gauss point weight (from `weights`) - } - } -} diff --git a/src/opacity/absorber_ck.hpp b/src/opacity/absorber_ck.hpp deleted file mode 100644 index 0e11ef53..00000000 --- a/src/opacity/absorber_ck.hpp +++ /dev/null @@ -1,42 +0,0 @@ -#ifndef SRC_OPACITY_ABSORBER_CK_HPP_ -#define SRC_OPACITY_ABSORBER_CK_HPP_ - -// C/C++ -#include -#include - -// opacity -#include "absorber.hpp" - -class AbsorberCK : public Absorber { - public: - AbsorberCK(std::string name) : Absorber(name) {} - virtual ~AbsorberCK() {} - - void LoadCoefficient(std::string fname, int bid) override; - Real GetAttenuation(int m, AirParcel const& var) const override; - void ModifySpectralGrid(std::vector& spec) const override; - - protected: - size_t len_[4]; - - std::vector axis_; - std::vector p_; - std::vector t_; - std::vector bin_centers_; - std::vector bin_edges_; - std::vector kcoeff_; - std::vector weights_; -}; - -class HeliosCK : public AbsorberCK { - public: - HeliosCK(std::string name); - virtual ~HeliosCK() {} - - void LoadCoefficient(std::string fname, int bid) override; - Real GetAttenuation(Real g1, Real g2, AirParcel const& var) const override; - void ModifySpectralGrid(std::vector& spec) const override; -}; - -#endif // SRC_OPACITY_ABSORBER_CK_HPP_ diff --git a/src/opacity/absorber_factory.cpp b/src/opacity/absorber_factory.cpp deleted file mode 100644 index ee3a7504..00000000 --- a/src/opacity/absorber_factory.cpp +++ /dev/null @@ -1,147 +0,0 @@ -// athena -#include - -// external -#include - -// application -#include - -// canoe -#include -#include - -// opacity -#include "Giants/grey_gas.hpp" -#include "Giants/hydrogen_cia.hpp" -#include "Giants/microwave/mwr_absorbers.hpp" -#include "absorber_ck.hpp" -#include "hitran_absorber.hpp" -#include "nitrogen_cia.hpp" -#include "oxygen_cia.hpp" -#include "simple_cloud.hpp" - -AbsorberContainer AbsorberFactory::CreateFrom( - std::vector const& names, std::string band_name, - YAML::Node const& rad) { - AbsorberContainer absorbers; - - for (auto& name : names) { - AbsorberContainer atmp; - - for (auto& my : rad["opacity-sources"]) { - if (name == my["name"].as()) { - atmp.push_back(AbsorberFactory::CreateFrom(my, band_name)); - } - } - - if (atmp.empty()) { - throw NotFoundError("AbsorberFactory", "Opacity " + name); - } else if (atmp.size() > - 1) { // use band_name to uniquely identify the absorber - bool found = false; - for (auto& ab : atmp) { - if (ab->GetOpacityFile().find(band_name) != std::string::npos) { - absorbers.push_back(ab); - found = true; - break; - } - } - - if (!found) { - throw NotFoundError("AbsorberFactory", - "There are multiple absorbers with the same name " + - name + - ". Rename the `data` field in the absorber to " - "identify the band."); - } - } else { - absorbers.push_back(atmp[0]); - } - } - - return absorbers; -} - -AbsorberPtr AbsorberFactory::CreateFrom(YAML::Node const& my, - std::string band_name) { - if (!my["name"]) { - throw NotFoundError("AbsorberFactory", "'name' field in absorber"); - } - - if (!my["class"]) { - throw NotFoundError("AbsorberFactory", "'class' field in absorber " + - my["name"].as()); - } - - auto ab = createAbsorberPartial(my["name"].as(), - my["class"].as()); - - if (my["data"]) { - ab->SetOpacityFile(my["data"].as()); - } else { - ab->SetOpacityFile(""); - } - - if (my["dependent-species"]) { - auto species = my["dependent-species"].as>(); - ab->SetSpeciesIndex(species); - } - - if (my["model"]) { - ab->SetModel(my["model"].as()); - } - - if (my["parameters"]) { - ab->SetRealsFrom(my["parameters"]); - } - - return ab; -} - -namespace gp = GiantPlanets; - -AbsorberPtr AbsorberFactory::createAbsorberPartial(std::string name, - std::string type) { - AbsorberPtr ab; - - if (type == "XIZ-H2-H2-CIA") { - ab = std::make_shared(); - } else if (type == "XIZ-H2-He-CIA") { - ab = std::make_shared(); - } else if (type == "Hitran") { - ab = std::make_shared(name); - } else if (type == "HitranCK") { - ab = std::make_shared(name); - } else if (type == "SimpleGrey") { - ab = std::make_shared(name); - } else if (type == "FreedmanMean") { - ab = std::make_shared(name); - } else if (type == "JupGasv") { - ab = std::make_shared(name); - } else if (type == "JupGasir") { - ab = std::make_shared(name); - } else if (type == "AbsorberCK") { - ab = std::make_shared(name); - } else if (type == "HeliosCK") { - ab = std::make_shared(name); - } else if (type == "SimpleCloud") { - ab = std::make_shared(name); - } else if (type == "radio-NH3") { - ab = std::make_shared(); - } else if (type == "radio-H2O") { - ab = std::make_shared(); - } else if (type == "radio-H2S") { - ab = std::make_shared(); - } else if (type == "radio-PH3") { - ab = std::make_shared(); - } else if (type == "radio-CIA") { - ab = std::make_shared(); - } else if (type == "radio-Electron") { - ab = std::make_shared(); - } else { - throw NotFoundError("createAbsorberPartial", "Class = " + type); - } - - return ab; -}