From 6cf8b84684758412f8d99a4178187e0ce583dc1b Mon Sep 17 00:00:00 2001 From: stormy/cli Date: Fri, 10 Nov 2023 17:14:46 -0500 Subject: [PATCH 1/2] wip --- cmake/disort.cmake | 2 +- python/pyharp.cpp | 32 +++----- src/harp/radiation.cpp | 39 +++++++++- src/harp/radiation.hpp | 17 +++++ src/harp/radiation_band.cpp | 28 +++++-- src/harp/radiation_band.hpp | 9 ++- src/harp/rt_solver_disort.cpp | 110 +++++++++++++++++++++------- src/harp/rt_solver_disort_flags.cpp | 12 --- src/harp/rt_solvers.hpp | 7 +- src/opacity/absorber.cpp | 9 +++ src/opacity/absorber.hpp | 4 + src/virtual_groups.hpp | 12 +++ tools/example_earth_rt.yaml | 6 +- 13 files changed, 203 insertions(+), 84 deletions(-) diff --git a/cmake/disort.cmake b/cmake/disort.cmake index 5925b644..aa52264c 100644 --- a/cmake/disort.cmake +++ b/cmake/disort.cmake @@ -18,7 +18,7 @@ if(DISORT) FetchContent_Declare( pydisort DOWNLOAD_EXTRACT_TIMESTAMP TRUE - URL https://${ACCOUNT}:${TOKEN}@github.com/zoeyzyhu/pydisort/archive/refs/tags/v0.5.2.tar.gz + URL https://${ACCOUNT}:${TOKEN}@github.com/zoeyzyhu/pydisort/archive/refs/tags/v0.6.tar.gz ) FetchContent_MakeAvailable(pydisort) diff --git a/python/pyharp.cpp b/python/pyharp.cpp index 08aa6882..51987e73 100644 --- a/python/pyharp.cpp +++ b/python/pyharp.cpp @@ -44,7 +44,8 @@ PYBIND11_MODULE(pyharp, m) { .def("get_num_bands", &Radiation::GetNumBands) .def("get_band", &Radiation::GetBand) - .def("cal_radiative_flux", &Radiation::CalRadiativeFlux) + + .def("cal_flux", &Radiation::CalRadiativeFlux) .def("cal_radiance", &Radiation::CalRadiance); // RadiationBand @@ -56,17 +57,20 @@ PYBIND11_MODULE(pyharp, m) { .def_readonly("bflxdn", &RadiationBand::bflxdn) .def(py::init()) + .def("__str__", &RadiationBand::ToString) .def("resize", &RadiationBand::Resize, py::arg("nc1"), py::arg("nc2") = 1, py::arg("nc3") = 1, py::arg("nstr") = 4) - .def("resize_solver", &RadiationBand::ResizeSolver, py::arg("nlyr"), - py::arg("nstr") = 4, py::arg("nuphi") = 1, py::arg("numu") = 1) .def("get_num_spec_grids", &RadiationBand::GetNumSpecGrids) .def("get_num_absorbers", &RadiationBand::GetNumAbsorbers) .def("absorbers", &RadiationBand::Absorbers) .def("get_absorber", &RadiationBand::GetAbsorber) .def("get_absorber_by_name", &RadiationBand::GetAbsorberByName) .def("get_range", &RadiationBand::GetRange) + + .def("cal_flux", &RadiationBand::CalBandFlux, py::arg("pmb") = nullptr, + py::arg("k") = 0, py::arg("j") = 0, py::arg("is") = 0, + py::arg("ie") = 0) .def("cal_radiance", &RadiationBand::CalBandRadiance, py::arg("pmb") = nullptr, py::arg("k") = 0, py::arg("j") = 0) @@ -143,21 +147,11 @@ PYBIND11_MODULE(pyharp, m) { } band.SetSpectralProperties(ac, x1f.data()); - }) - - .def("__str__", [](RadiationBand &band) { - std::stringstream ss; - ss << "RadiationBand: " << band.GetName() << std::endl; - ss << "Absorbers: "; - for (int n = 0; n < band.GetNumAbsorbers() - 1; ++n) { - ss << band.GetAbsorber(n)->GetName() << ", "; - } - ss << band.GetAbsorber(band.GetNumAbsorbers() - 1)->GetName(); - return ss.str(); - }); + }); // Absorber py::class_(m, "absorber") + .def("__str__", &Absorber::ToString) .def("load_opacity_from_file", &Absorber::LoadOpacityFromFile) .def("get_name", &Absorber::GetName) @@ -178,11 +172,5 @@ PYBIND11_MODULE(pyharp, m) { lst.begin(), lst.end(), air.w, [](const py::handle &elem) { return py::cast(elem); }); return ab.GetSingleScatteringAlbedo(wave1, wave2, air); - }) - - .def("__str__", [](Absorber &ab) { - std::stringstream ss; - ss << "Absorber: " << ab.GetName(); - return ss.str(); - }); + }); } diff --git a/src/harp/radiation.cpp b/src/harp/radiation.cpp index 2226d6f9..1f0a3ec1 100644 --- a/src/harp/radiation.cpp +++ b/src/harp/radiation.cpp @@ -1,4 +1,5 @@ // C/C++ headers +#include #include #include #include @@ -56,8 +57,6 @@ Radiation::Radiation(MeshBlock *pmb, ParameterInput *pin) { // radiation configuration int nstr = pin->GetOrAddInteger("radiation", "nstr", 8); - int nuphi = pin->GetOrAddInteger("radiation", "nuphi", 1); - int numu = pin->GetOrAddInteger("radiation", "numu", 1); for (auto &p : bands_) { // outgoing radiation direction (mu,phi) in degrees @@ -71,7 +70,6 @@ Radiation::Radiation(MeshBlock *pmb, ParameterInput *pin) { // allocate memory p->Resize(ncells1, ncells2, ncells3, nstr); - p->ResizeSolver(ncells1 - 2 * NGHOST, nstr, nuphi, numu); } // output radiance @@ -201,6 +199,41 @@ size_t Radiation::LoadRestartData(char *psrc) { } namespace RadiationHelper { +bool real_close(Real num1, Real num2, Real tolerance) { + return std::fabs(num1 - num2) <= tolerance; +} + +std::pair, std::vector> get_direction_grids( + std::vector const &dirs) { + std::vector uphi; + std::vector umu; + + for (auto &dir : dirs) { + // find phi + bool found = false; + for (auto &phi : uphi) + if (real_close(phi, dir.phi, 1.e-3)) { + found = true; + break; + } + if (!found) uphi.push_back(dir.phi); + + // find mu + found = false; + for (auto &mu : umu) + if (real_close(mu, dir.mu, 1.e-3)) { + found = true; + break; + } + if (!found) umu.push_back(dir.mu); + } + + std::sort(uphi.begin(), uphi.end()); + std::sort(umu.begin(), umu.end()); + + return std::make_pair(uphi, umu); +} + Direction parse_radiation_direction(std::string_view str) { Direction ray; ray.phi = 0.; diff --git a/src/harp/radiation.hpp b/src/harp/radiation.hpp index 22c1dd10..f96d6354 100644 --- a/src/harp/radiation.hpp +++ b/src/harp/radiation.hpp @@ -10,6 +10,9 @@ // athena #include +// canoe +#include // Direction + // harp #include "radiation_band.hpp" @@ -101,7 +104,21 @@ const uint64_t WriteBinRadiance = 1LL << 8; } // namespace RadiationFlags namespace RadiationHelper { +//! \brief Get the number of grids in the outgoing ray directions +std::pair, std::vector> get_direction_grids( + std::vector const &dirs); + +//! \brief Parse radiation direction string +//! +//! First number is the polar angle (degrees), second number is the azimuthal +//! angle (degrees) \param[in] str radiation direction string, e.g. (45, 30) +//! \return radiation direction Direction parse_radiation_direction(std::string_view str); + +//! \brief Parse radiation directions string, sperated by comma +//! +//! Example input string: "(45, 30), (45, 60)" +//! \param[in] str radiation directions string std::vector parse_radiation_directions(std::string str); uint64_t parse_radiation_flags(std::string str); void get_phase_momentum(Real *pmom, int iphas, Real gg, int npmom); diff --git a/src/harp/radiation_band.cpp b/src/harp/radiation_band.cpp index a74031e6..41424d2e 100644 --- a/src/harp/radiation_band.cpp +++ b/src/harp/radiation_band.cpp @@ -115,16 +115,17 @@ void RadiationBand::Resize(int nc1, int nc2, int nc3, int nstr) { ssa_.NewAthenaArray(pgrid_->spec.size(), nc1); ssa_.ZeroClear(); - toa_.NewAthenaArray(pgrid_->spec.size(), rayOutput_.size()); - toa_.ZeroClear(); - pmom_.NewAthenaArray(pgrid_->spec.size(), nc1, nstr + 1); pmom_.ZeroClear(); - flxup_.NewAthenaArray(pgrid_->spec.size(), nc1); + // 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(), nc1); + flxdn_.NewAthenaArray(pgrid_->spec.size(), nc3, nc2, nc1); flxdn_.ZeroClear(); // band properties @@ -135,10 +136,10 @@ void RadiationBand::Resize(int nc1, int nc2, int nc3, int nstr) { btoa.NewAthenaArray(rayOutput_.size(), nc3, nc2); bflxup.NewAthenaArray(nc3, nc2, nc1 + 1); bflxdn.NewAthenaArray(nc3, nc2, nc1 + 1); -} -void RadiationBand::ResizeSolver(int nlyr, int nstr, int nuphi, int numu) { - if (psolver_ != nullptr) psolver_->Resize(nlyr, nstr, nuphi, numu); + if (psolver_ != nullptr) { + psolver_->Resize(nc1 - 2 * NGHOST, nstr); + } } AbsorberPtr RadiationBand::GetAbsorberByName(std::string const &name) { @@ -241,6 +242,17 @@ void RadiationBand::WriteAsciiData(OutputParameters const *pout) const { 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 << absorbers_.back()->GetName(); + return ss.str(); +} + std::shared_ptr RadiationBand::CreateRTSolverFrom( std::string const &rt_name, YAML::Node const &rad) { std::shared_ptr psolver; diff --git a/src/harp/radiation_band.hpp b/src/harp/radiation_band.hpp index e938e901..71c80936 100644 --- a/src/harp/radiation_band.hpp +++ b/src/harp/radiation_band.hpp @@ -33,6 +33,7 @@ class RadiationBand : public NamedGroup, public FlagGroup, public ParameterGroup, public ASCIIOutputGroup, + public StringReprGroup, public LinearExchanger { public: // public access data // implementation of RT Solver @@ -64,10 +65,7 @@ class RadiationBand : public NamedGroup, public: // member functions //! \brief Allocate memory for radiation band - void Resize(int nc1, int nc2 = 1, int nc3 = 1, int nstr = 8); - - //! \brief Allocate memory for radiation solver - void ResizeSolver(int nlyr, int nstr = 8, int nuphi = 1, int numu = 1); + void Resize(int nc1, int nc2, int nc3, int nstr); //! \brief Create radiative transfer solver from YAML node //! @@ -158,6 +156,9 @@ class RadiationBand : public NamedGroup, void Transfer(MeshBlock const *pmb, int n) override; + public: // StringReprGroup functions + std::string ToString() const override; + protected: //! radiative transfer solver std::shared_ptr psolver_; diff --git a/src/harp/rt_solver_disort.cpp b/src/harp/rt_solver_disort.cpp index 1fb19ab5..85d28549 100644 --- a/src/harp/rt_solver_disort.cpp +++ b/src/harp/rt_solver_disort.cpp @@ -16,6 +16,9 @@ #include #include +// climath +#include + // canoe #include #include @@ -29,10 +32,22 @@ #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"), DisortWrapper() { - if (rad["Disort"]) setFlagsFromNode(rad["Disort"]); + : RTSolver(pmy_band, "Disort") { + if (rad["Disort-flags"]) { + SetFlags(to_map_bool(rad["Disort-flags"])); + } if (pmy_band->TestFlag(RadiationFlags::Planck)) { ds_.flag.planck = 1; @@ -42,20 +57,20 @@ RadiationBand::RTSolverDisort::RTSolverDisort(RadiationBand *pmy_band, } //! \todo update based on band outdir -void RadiationBand::RTSolverDisort::Resize(int nlyr, int nstr, int nuphi, - int numu) { +void RadiationBand::RTSolverDisort::Resize(int nlyr, int nstr) { + Unseal(); + + auto &rayout = pmy_band_->rayOutput_; + auto &&uphi_umu = RadiationHelper::get_direction_grids(rayout); + SetAtmosphereDimension(nlyr, nstr, nstr, nstr); - //! \todo revise this - SetIntensityDimension(nuphi, 1, pmy_band_->GetNumOutgoingRays()); - Finalize(); - Real utau = 0.; - Real uphi = 0.; - Real umu = 1.; + 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]); - SetUserAzimuthalAngle(&uphi, 1); - SetUserOpticalDepth(&utau, 1); - SetUserCosinePolarAngle(&pmy_band_->rayOutput_[0].mu, 1); + SetIntensityDimension(dir_dim_[1], 1, dir_dim_[0]); + Seal(); } //! \note Counting Disort Index @@ -133,11 +148,45 @@ void RadiationBand::RTSolverDisort::Prepare(MeshBlock const *pmb, int k, pmb->pcoord->CellVolume(k, j, pmb->is, pmb->ie, vol_); } - ds_.bc.fluor = 0.; - ds_.bc.albedo = 0.; + //! \todo update this + ds_.bc.fluor = fluor; + ds_.bc.fisot = fisot; + ds_.bc.albedo = albedo; ds_.bc.temis = 1.; - Finalize(); + auto &&uphi_umu = RadiationHelper::get_direction_grids(pmy_band_->rayOutput_); + auto &uphi = uphi_umu.first; + auto &umu = uphi_umu.second; + + if (umu.size() <= ds_.numu) { + SetUserCosinePolarAngle(umu); + + for (int i = 0; i < umu.size(); ++i) { + dir_axis_[i] = umu[i]; + } + + for (int i = umu.size(); i < ds_.numu; ++i) { + dir_axis_[i] = 1.; + } + } else { + throw RuntimeError("RTSolverDisort::Prepare", + "Number of polar angles in Disort is too small"); + } + + if (uphi.size() <= ds_.nphi) { + SetUserAzimuthalAngle(uphi); + + for (int i = 0; i < uphi.size(); ++i) { + dir_axis_[ds_.numu + i] = uphi[i]; + } + + for (int i = uphi.size(); i < ds_.nphi; ++i) { + dir_axis_[ds_.numu + i] = 2. * M_PI; + } + } else { + throw RuntimeError("RTSolverDisort::Prepare", + "Number of azimuthal angles in Disort is too small"); + } } void RadiationBand::RTSolverDisort::CalBandFlux(MeshBlock const *pmb, int k, @@ -195,15 +244,15 @@ void RadiationBand::RTSolverDisort::addDisortFlux(Coordinates const *pcoord, //! \bug does not work for spherical geometry, need to scale area using //! farea(il)/farea(i) // flux up - flxup(n, i) = ds_out_.rad[m].flup; + flxup(n, 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(n, i) = ds_out_.rad[m].rfldir + ds_out_.rad[m].rfldn; + flxdn(n, k, j, i) = ds_out_.rad[m].rfldir + ds_out_.rad[m].rfldn; - bflxup(k, j, i) += spec[n].wght * flxup(n, i); - bflxdn(k, j, i) += spec[n].wght * flxdn(n, i); + bflxup(k, j, i) += spec[n].wght * flxup(n, k, j, i); + bflxdn(k, j, i) += spec[n].wght * flxdn(n, k, j, i); } //! \note Spherical correction by XIZ @@ -233,7 +282,7 @@ void RadiationBand::RTSolverDisort::CalBandRadiance(MeshBlock const *pmb, int k, int j) { if (ds_.flag.onlyfl) { throw RuntimeError("RTSolverDisort::CalBandRadiance", - "Only flux calculation is requested"); + "Radiance calculation disabled"); } if (ds_.ntau != 1) { @@ -241,11 +290,11 @@ void RadiationBand::RTSolverDisort::CalBandRadiance(MeshBlock const *pmb, int k, "Only toa radiance (ds.ntau = 1) is supported"); } - int nrays = ds_.nphi * ds_.ntau * ds_.numu; + int nrays = ds_.nphi * ds_.numu; - if (nrays != pmy_band_->GetNumOutgoingRays()) { + if (nrays < pmy_band_->GetNumOutgoingRays()) { throw RuntimeError("RTSolverDisort::CalBandRadiance", - "Number of outgoing rays does not match"); + "Number of outgoing rays more than DISORT can host"); } Real dist_au; @@ -294,12 +343,17 @@ void RadiationBand::RTSolverDisort::CalBandRadiance(MeshBlock const *pmb, int k, } 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; - - for (int c = 0; c < ds_.nphi * ds_.ntau * ds_.numu; ++c) { - pmy_band_->toa_(b, 0) = ds_out_.uu[c]; - btoa(c, k, j) += spec[b].wght * ds_out_.uu[c]; + 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); } } diff --git a/src/harp/rt_solver_disort_flags.cpp b/src/harp/rt_solver_disort_flags.cpp index 4ec3360f..aad44261 100644 --- a/src/harp/rt_solver_disort_flags.cpp +++ b/src/harp/rt_solver_disort_flags.cpp @@ -16,18 +16,6 @@ void RadiationBand::RTSolverDisort::setFlagsFromNode(YAML::Node const& my) { ds_.flag.ibcnd = 0; } - if (my["accur"]) { - ds_.accur = my["accur"].as(); - } else { - ds_.accur = 1.e-3; - } - - if (my["fisot"]) { - ds_.bc.fisot = my["fisot"].as(); - } else { - ds_.bc.fisot = 0.; - } - if (my["usrtau"]) { ds_.flag.usrtau = my["usrtau"].as(); } else { diff --git a/src/harp/rt_solvers.hpp b/src/harp/rt_solvers.hpp index 4994ab36..e9440e73 100644 --- a/src/harp/rt_solvers.hpp +++ b/src/harp/rt_solvers.hpp @@ -43,7 +43,7 @@ class RadiationBand::RTSolver : public NamedGroup { virtual void Prepare(MeshBlock const *pmb, int k, int j) {} //! \brief Allocate memory for radiation solver - virtual void Resize(int nlyr, int nstr, int nuphi, int numu) {} + virtual void Resize(int nlyr, int nstr) {} public: // inbound functions virtual void CalBandFlux(MeshBlock const *pmb, int k, int j, int il, int iu) { @@ -78,13 +78,16 @@ class RadiationBand::RTSolverDisort : public RadiationBand::RTSolver, public: // member functions void Prepare(MeshBlock const *pmb, int k, int j) override; - void Resize(int nlyr, int nstr, int nuphi, int numu) override; + void Resize(int nlyr, int nstr) override; public: // inbound functions void CalBandFlux(MeshBlock const *pmb, int k, int j, int il, int iu) override; void CalBandRadiance(MeshBlock const *pmb, int k, int j) override; protected: + size_t dir_dim_[2]; + std::vector dir_axis_; + void setFlagsFromNode(YAML::Node const &flags); void addDisortFlux(Coordinates const *pcoord, int n, int k, int j, int il, diff --git a/src/opacity/absorber.cpp b/src/opacity/absorber.cpp index 6fdd682e..f9a502f4 100644 --- a/src/opacity/absorber.cpp +++ b/src/opacity/absorber.cpp @@ -1,4 +1,5 @@ // C/C++ +#include #include // application @@ -43,3 +44,11 @@ void Absorber::LoadOpacity() { log->Warn(ss.str()); } } + +std::string Absorber::ToString() const { + std::stringstream ss; + ss << "Absorber: " << GetName(); + ss << "Opacity file: " << opacity_filename_; + + return ss.str(); +} diff --git a/src/opacity/absorber.hpp b/src/opacity/absorber.hpp index 7a7b6d30..abef676d 100644 --- a/src/opacity/absorber.hpp +++ b/src/opacity/absorber.hpp @@ -22,6 +22,7 @@ class AirParcel; class Absorber : public NamedGroup, public ParameterGroup, public SpeciesIndexGroup, + public StringReprGroup, public CheckGroup { public: // constructor and destructor Absorber(std::string name); @@ -59,6 +60,9 @@ class Absorber : public NamedGroup, virtual void GetPhaseMomentum(Real* pp, Real wave1, Real wave2, AirParcel const& var, int np) const {} + public: // StringRepr + std::string ToString() const override; + protected: //! absorption model model std::string model_name_; diff --git a/src/virtual_groups.hpp b/src/virtual_groups.hpp index b9faf4cd..291a279e 100644 --- a/src/virtual_groups.hpp +++ b/src/virtual_groups.hpp @@ -44,6 +44,12 @@ class FlagGroup { uint64_t myflags_; }; +class StringReprGroup { + public: + virtual ~StringReprGroup() {} + virtual std::string ToString() const = 0; +}; + class ParameterGroup { public: virtual ~ParameterGroup() {} @@ -54,6 +60,12 @@ class ParameterGroup { } } + void SetIntsFrom(YAML::Node const &node) { + for (auto it = node.begin(); it != node.end(); ++it) { + params_int_[it->first.as()] = it->second.as(); + } + } + //! Set real parameter void SetPar(std::string const &name, Real value) { params_real_[name] = value; diff --git a/tools/example_earth_rt.yaml b/tools/example_earth_rt.yaml index cde021b2..aeb03391 100644 --- a/tools/example_earth_rt.yaml +++ b/tools/example_earth_rt.yaml @@ -23,7 +23,5 @@ B1: outdir: ["(45,)"] flags: [static, planck] -Disort: - prnt: [0,0,0,0,0] - quiet: true - accur: 1.e-6 +Disort-flags: + quiet: false From 10eca6f0e917895917ec87fa13d05b11c39d57cf Mon Sep 17 00:00:00 2001 From: mac/cli Date: Fri, 10 Nov 2023 19:10:48 -0500 Subject: [PATCH 2/2] wip --- src/harp/radiation.cpp | 6 ------ src/harp/radiation.hpp | 3 --- src/harp/radiation_band.cpp | 1 + src/harp/rt_solver_disort.cpp | 5 +---- tools/example_earth.yaml | 9 ++++++++- tools/example_earth_rt.yaml | 19 +++++++++++++++---- tools/run_disort_rt.py | 7 +------ 7 files changed, 26 insertions(+), 24 deletions(-) diff --git a/src/harp/radiation.cpp b/src/harp/radiation.cpp index 1f0a3ec1..73e3ec08 100644 --- a/src/harp/radiation.cpp +++ b/src/harp/radiation.cpp @@ -274,14 +274,8 @@ uint64_t parse_radiation_flags(std::string str) { flags |= RadiationFlags::LineByLine; } else if (dstr[i] == "ck") { flags |= RadiationFlags::CorrelatedK; - } else if (dstr[i] == "planck") { - flags |= RadiationFlags::Planck; } else if (dstr[i] == "star") { flags |= RadiationFlags::Star; - } else if (dstr[i] == "spher") { - flags |= RadiationFlags::Sphere; - } else if (dstr[i] == "only") { - flags |= RadiationFlags::FluxOnly; } else if (dstr[i] == "normalize") { flags |= RadiationFlags::Normalize; } else if (dstr[i] == "write_bin_radiance") { diff --git a/src/harp/radiation.hpp b/src/harp/radiation.hpp index f96d6354..aef6aba4 100644 --- a/src/harp/radiation.hpp +++ b/src/harp/radiation.hpp @@ -91,10 +91,7 @@ const uint64_t None = 0LL; const uint64_t Dynamic = 1LL << 0; const uint64_t LineByLine = 1LL << 1; const uint64_t CorrelatedK = 1LL << 2; -const uint64_t Planck = 1LL << 3; const uint64_t Star = 1LL << 4; -const uint64_t Sphere = 1LL << 5; -const uint64_t FluxOnly = 1LL << 6; const uint64_t Normalize = 1LL << 7; const uint64_t WriteBinRadiance = 1LL << 8; diff --git a/src/harp/radiation_band.cpp b/src/harp/radiation_band.cpp index 41424d2e..280d9e4a 100644 --- a/src/harp/radiation_band.cpp +++ b/src/harp/radiation_band.cpp @@ -250,6 +250,7 @@ std::string RadiationBand::ToString() const { ss << ab->GetName() << ", "; } ss << absorbers_.back()->GetName(); + ss << "RT-Solver: " << psolver_->GetName() << std::endl; return ss.str(); } diff --git a/src/harp/rt_solver_disort.cpp b/src/harp/rt_solver_disort.cpp index 85d28549..539486cf 100644 --- a/src/harp/rt_solver_disort.cpp +++ b/src/harp/rt_solver_disort.cpp @@ -49,10 +49,6 @@ RadiationBand::RTSolverDisort::RTSolverDisort(RadiationBand *pmy_band, SetFlags(to_map_bool(rad["Disort-flags"])); } - if (pmy_band->TestFlag(RadiationFlags::Planck)) { - ds_.flag.planck = 1; - } - SetHeader("Disort solver"); } @@ -169,6 +165,7 @@ void RadiationBand::RTSolverDisort::Prepare(MeshBlock const *pmb, int k, dir_axis_[i] = 1.; } } else { + std::cout << "umu.size() = " << umu.size() << std::endl; throw RuntimeError("RTSolverDisort::Prepare", "Number of polar angles in Disort is too small"); } diff --git a/tools/example_earth.yaml b/tools/example_earth.yaml index b608b33a..8f32d506 100644 --- a/tools/example_earth.yaml +++ b/tools/example_earth.yaml @@ -8,7 +8,7 @@ opacity-sources: - name: O3 class: Hitran -bands: [B1] +bands: [B1, B2] B1: units: cm-1 @@ -16,3 +16,10 @@ B1: wavenumber-range: [600., 700.] resolution: 0.01 opacity: [H2O, CO2] + +B2: + units: cm-1 + grid-type: regular + wavenumber-range: [200., 600.] + resolution: 0.01 + opacity: [H2O, CO2] diff --git a/tools/example_earth_rt.yaml b/tools/example_earth_rt.yaml index aeb03391..718b1f21 100644 --- a/tools/example_earth_rt.yaml +++ b/tools/example_earth_rt.yaml @@ -11,7 +11,7 @@ opacity-sources: class: Hitran dependent-species: [tracer.O3] -bands: [B1] +bands: [B1, B2] B1: units: cm-1 @@ -20,8 +20,19 @@ B1: resolution: 0.01 opacity: [H2O, CO2] rt-solver: Disort - outdir: ["(45,)"] - flags: [static, planck] + outdir: ["(0,)","(45,)"] + flags: [static] Disort-flags: - quiet: false + ibcnd: false + usrtau: true + usrang: true + lamber: true + planck: true + onlyfl: false + spher: false + intensity_correction: true + old_intensity_correction: false + general_source: false + output_uum: false + quiet: true diff --git a/tools/run_disort_rt.py b/tools/run_disort_rt.py index da4a909d..a0bdced2 100755 --- a/tools/run_disort_rt.py +++ b/tools/run_disort_rt.py @@ -2,8 +2,7 @@ from canoe import * from pyharp import radiation_band, subscribe_species from utilities import load_file -from collections import OrderedDict -from numpy import * +from numpy import linspace, ones, exp # create atmosphere dictionary @@ -25,9 +24,6 @@ def create_atmosphere(nlyr: int) -> dict: atm["CO2"] = 400.0 * ones(nlyr) # ppmv atm["H2O"] = 4000.0 * exp(-atm["HGT"] / (Hscale_km / 5.0)) - atm["O2"] = (1e6 - atm["CO2"] - atm["H2O"]) * 0.21 - atm["N2"] = (1e6 - atm["CO2"] - atm["H2O"]) * 0.78 - atm["Ar"] = (1e6 - atm["CO2"] - atm["H2O"]) * 0.01 return atm @@ -50,7 +46,6 @@ def create_atmosphere(nlyr: int) -> dict: num_layers = 100 band.resize(num_layers) - band.resize_solver(num_layers) atm = create_atmosphere(num_layers) band.set_spectral_properties(atm)