Merge pull request #284 from Magritte-code/283-port-cooling-computati…

…on-to-c++

Added cooling rate computation to C++ magritte.

.github/workflows/build-and-test-on-pr.yml

@@ -109,7 +109,7 @@ jobs:
           activate-environment: magritte
           environment-file: dependencies/conda_env.yml
           auto-activate-base: false
-          miniforge-variant: Mambaforge
+          miniforge-version: latest
           channels: conda-forge, defaults
       - name: wait 5 s if attempt failed
         if: steps.setupcondatry1.outcome == 'failure'
@@ -122,7 +122,7 @@ jobs:
           activate-environment: magritte
           environment-file: dependencies/conda_env.yml
           auto-activate-base: false
-          miniforge-variant: Mambaforge
+          miniforge-version: latest
           channels: conda-forge, defaults
       - name: wait 5 s if attempt failed
         if: steps.setupcondatry2.outcome == 'failure'
@@ -135,7 +135,7 @@ jobs:
           activate-environment: magritte
           environment-file: dependencies/conda_env.yml
           auto-activate-base: false
-          miniforge-variant: Mambaforge
+          miniforge-version: latest
           channels: conda-forge, defaults
       # end retrying conda
 

.github/workflows/build-and-test.yml

@@ -130,7 +130,7 @@ jobs:
           activate-environment: magritte
           environment-file: dependencies/conda_env.yml
           auto-activate-base: false
-          miniforge-variant: Mambaforge
+          miniforge-version: latest
           channels: conda-forge, defaults
       - name: wait 5 s if attempt failed
         if: steps.setupcondatry1.outcome == 'failure'
@@ -143,7 +143,7 @@ jobs:
           activate-environment: magritte
           environment-file: dependencies/conda_env.yml
           auto-activate-base: false
-          miniforge-variant: Mambaforge
+          miniforge-version: latest
           channels: conda-forge, defaults
       - name: wait 5 s if attempt failed
         if: steps.setupcondatry2.outcome == 'failure'
@@ -156,7 +156,7 @@ jobs:
           activate-environment: magritte
           environment-file: dependencies/conda_env.yml
           auto-activate-base: false
-          miniforge-variant: Mambaforge
+          miniforge-version: latest
           channels: conda-forge, defaults
       # end retrying conda
 

src/bindings/pybindings.cpp

@@ -258,6 +258,9 @@ PYBIND11_MODULE(core, module) {
             "Compute an image of the optical depth for the model along the "
             "given ray direction, using the new imager, specifying the ray "
             "direction and image resolution.")
+        .def("compute_line_cooling_rates", &Model::compute_line_cooling_rates,
+            "Compute the line cooling rates for the model. Requires the level populations to be in "
+            "NLTE.")
         .def("set_eta_and_chi", &Model::set_eta_and_chi,
             "Set latest emissivity and opacity for the model in the eta and chi "
             "variables respectively.")
@@ -561,6 +564,7 @@ PYBIND11_MODULE(core, module) {
         .def_readonly("RT", &LineProducingSpecies::RT)
         .def_readonly("LambdaStar", &LineProducingSpecies::LambdaStar)
         .def_readonly("LambdaTest", &LineProducingSpecies::LambdaTest)
+        .def_readonly("line_cooling_rate", &LineProducingSpecies::line_cooling_rate)
         // functions
         .def("read", &LineProducingSpecies::read, "Read object from file.")
         .def("write", &LineProducingSpecies::write, "Write object to file.")

src/model/lines/lineProducingSpecies/lineProducingSpecies.hpp

@@ -43,8 +43,8 @@ struct LineProducingSpecies {
 
     // For ng-acceleration purposes, the level populations must be stored as accurately as possible
 
-    VectorXld population; ///< level population (most recent)
-    Real1 population_tot; ///< total level population (sum over levels)
+    VectorXld population;        ///< level population (most recent)
+    Vector<Real> population_tot; ///< total level population (sum over levels)
 
     vector<VectorXld> populations; ///< list of populations in previous iterations
     vector<VectorXld> residuals;   ///< list of residuals in the populations
@@ -62,6 +62,8 @@ struct LineProducingSpecies {
     SparseMatrix<long double> LambdaTest;
     SparseMatrix<long double> LambdaStar;
 
+    Real1 line_cooling_rate; ///< cooling rate due to line emission [W/m^3]
+
     LineProducingSpecies(std::shared_ptr<Parameters> params) :
         parameters(params), quadrature(params), lambda(params){};
 
@@ -95,6 +97,9 @@ struct LineProducingSpecies {
     inline void update_using_acceleration_trial(const Size order);
     inline void correct_negative_populations(
         const Double2& abundance, const Vector<Real>& temperature);
+
+    inline void compute_line_cooling_rates(
+        const Double2& abundance, const Vector<Real>& temperature);
 };
 
 #include "lineProducingSpecies.tpp"

src/model/lines/lineProducingSpecies/lineProducingSpecies.tpp

@@ -800,3 +800,37 @@ inline void LineProducingSpecies::correct_negative_populations(
                   << std::endl;
     }
 }
+
+inline void LineProducingSpecies ::compute_line_cooling_rates(
+    const Double2& abundance, const Vector<Real>& temperature) {
+    line_cooling_rate.resize(parameters->npoints());
+    // zero the vector manually, as it might already contain data
+    for (Size p = 0; p < parameters->npoints(); p++) {
+        line_cooling_rate[p] = 0.0;
+    }
+    threaded_for(p, parameters->npoints(), {
+        // Collisional transitions
+        for (CollisionPartner& colpar : linedata.colpar) {
+            Real abn = abundance[p][colpar.num_col_partner];
+            Real tmp = temperature[p];
+
+            colpar.adjust_abundance_for_ortho_or_para(tmp, abn);
+            colpar.interpolate_collision_coefficients(tmp);
+
+            for (Size k = 0; k < colpar.ncol; k++) {
+                const long double v_IJ = colpar.Cd_intpld()[k] * abn;
+                const long double v_JI = colpar.Ce_intpld()[k] * abn;
+
+                const Size I = colpar.icol[k];
+                const Size J = colpar.jcol[k];
+
+                const long double total_rate_IJ = v_IJ * population(index(p, I));
+                const long double total_rate_JI = v_JI * population(index(p, J));
+
+                const Real energy_diff = linedata.energy[I] - linedata.energy[J];
+
+                line_cooling_rate[p] += (total_rate_JI - total_rate_IJ) * energy_diff;
+            }
+        }
+    })
+}

src/model/lines/lines.cpp

@@ -149,6 +149,12 @@ void Lines ::iteration_using_statistical_equilibrium_sparse(
     // gather_emissivities_and_opacities ();
 }
 
+void Lines ::compute_line_cooling_rates(const Double2& abundance, const Vector<Real>& temperature) {
+    for (LineProducingSpecies& lspec : lineProducingSpecies) {
+        lspec.compute_line_cooling_rates(abundance, temperature);
+    }
+}
+
 // DEPRECATED: Now try to keep emissivities and opacities local.
 //
 // void Lines :: gather_emissivities_and_opacities ()

src/model/lines/lines.hpp

@@ -43,6 +43,8 @@ struct Lines {
 
     void trial_iteration_using_adaptive_Ng_acceleration(const Real pop_prec, const Size order);
 
+    void compute_line_cooling_rates(const Double2& abundance, const Vector<Real>& temperature);
+
     inline Size index(const Size p, const Size line_index) const;
     inline Size line_index(const Size l, const Size k) const;
     inline Size index(const Size p, const Size l, const Size k) const;

src/model/model.cpp

@@ -1133,6 +1133,12 @@ int Model ::set_eta_and_chi(const Size rr) {
     return (0);
 }
 
+int Model::compute_line_cooling_rates() {
+    std::cout << "Computing line cooling rates based on current level populations." << std::endl;
+    lines.compute_line_cooling_rates(chemistry.species.abundance, thermodynamics.temperature.gas);
+    return (0);
+}
+
 int Model ::set_boundary_condition() {
     Solver solver;
     solver.set_boundary_condition(*this);

src/model/model.hpp

@@ -88,6 +88,8 @@ struct Model {
     int compute_image_optical_depth_new(const Size ray_nr); // most similar function formulation
                                                             // to old imager
 
+    int compute_line_cooling_rates();
+
     Double1 error_max;
     Double1 error_mean;