Revise pyharp API

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)

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<std::string, YAML::Node>())
+      .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_<Absorber, AbsorberPtr>(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<double>(elem); });
              return ab.GetSingleScatteringAlbedo(wave1, wave2, air);
-           })
-
-      .def("__str__", [](Absorber &ab) {
-        std::stringstream ss;
-        ss << "Absorber: " << ab.GetName();
-        return ss.str();
-      });
+           });
 }

src/harp/radiation.cpp

@@ -1,4 +1,5 @@
 // C/C++ headers
+#include <algorithm>
 #include <fstream>
 #include <sstream>
 #include <stdexcept>
@@ -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<Real>, std::vector<Real>> get_direction_grids(
+    std::vector<Direction> const &dirs) {
+  std::vector<Real> uphi;
+  std::vector<Real> 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.;
@@ -241,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") {

src/harp/radiation.hpp

@@ -10,6 +10,9 @@
 // athena
 #include <athena/athena.hpp>
 
+// canoe
+#include <common.hpp>  // Direction
+
 // harp
 #include "radiation_band.hpp"
 
@@ -88,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;
 
@@ -101,7 +101,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<Real>, std::vector<Real>> get_direction_grids(
+    std::vector<Direction> 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<Direction> 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);

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,18 @@ 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();
+  ss << "RT-Solver: " << psolver_->GetName() << std::endl;
+  return ss.str();
+}
+
 std::shared_ptr<RadiationBand::RTSolver> RadiationBand::CreateRTSolverFrom(
     std::string const &rt_name, YAML::Node const &rad) {
   std::shared_ptr<RTSolver> psolver;

src/harp/radiation_band.hpp

@@ -33,6 +33,7 @@ class RadiationBand : public NamedGroup,
                       public FlagGroup,
                       public ParameterGroup,
                       public ASCIIOutputGroup,
+                      public StringReprGroup,
                       public LinearExchanger<RadiationBand> {
  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<RTSolver> psolver_;