SolarFluxVariablePotter2023Radiosity.h

Go to the documentation of this file.

#ifndef SOLAR_FLUX_VARIABLE_POTTER_2023_RADIOSITY_H
#define SOLAR_FLUX_VARIABLE_POTTER_2023_RADIOSITY_H
 
#define VERSION "6"
 
#ifdef VTK_PRESENT
#include <cmath>
#include <iostream>
#include <memory>
#include <numbers>
#include <string>
#include <valarray>
 
#ifdef MPI_PRESENT
#include <mpi.h>
 
#include "MpiCommunicationPatterns.h"
#endif
 
#include "../../src/CsrMatrix.h"
#include "../../src/GenericSubmodule.h"
#include "../../src/ModuleFactory.h"
#include "../../src/PhysicalConstants.h"
#include "../../src/VtkRayCasting.h"
 
template <typename T>
class SolarFluxVariablePotter2023Radiosity : public GenericSubmodule<T>
{
    private:
        static bool m_registered;
        std::string m_ini_filepath{"solar_flux/SolarFluxVariablePotter2023Radiosity.ini"};
 
        std::unique_ptr<VtkRayCasting> m_ray_caster;
        int m_num_facets;
        int m_numerical_layers;
        int m_check_timestep{0};
        std::string m_shape_model_path;
        // // std::vector<std::vector<double>> m_viewfactor_matrix;
        std::vector<double> direct_solar_flux;
        CsrSparseMatrix<double> m_viewfactor_matrix;
        int m_number_non_zero_elements{};
        std::valarray<T> m_flux_sol;
        std::valarray<T> m_albedo;
        std::valarray<T> m_emissivity;
        std::valarray<T> m_flux_radiative;
        std::valarray<T> m_flux_reflected;
        std::valarray<T> m_flux_infrared;
        std::valarray<T> m_flux_reflected_temp;
        std::valarray<T> m_flux_infrared_temp;
        std::valarray<T> m_heat_conduction_flux;
        int id_start{0};
        int id_end;
        T m_stefan_boltzmann_const;
        T m_solar_const;
        void calculateViewfactorMatrix();
        void calculateFluxTotal(const std::valarray<T>& sol_vec, T initial_flux);
        void radiosityCalcFirstTimeStep(const std::vector<T>& direct_solar_flux);
        void radiosityCalc(const std::vector<T>& direct_solar_flux);
 
#ifdef MPI_PRESENT
        std::valarray<T> m_global_temperature_array;
        void calculateFluxTotalMPI(const std::valarray<T>& sol_vec, T initial_flux);
#endif
 
    public:
        SolarFluxVariablePotter2023Radiosity(SimulationClassBase<T>* sim);
        bool setup(std::vector<std::shared_ptr<GenericSubmodule<T>>> all_submodules)
            override;  // loads module into necessary module chains (auto load for calculation
                       // chains)
        bool exec(std::string_view param) override;  // main functionality of the module
 
        bool init() override;
        bool preTimeStep() override;
        bool postTimeStep() override { return true; };
        bool output() override { return true; };
 
        void setFieldPtr(std::shared_ptr<std::valarray<T>> field_ptr) override
        {
            this->module_field = field_ptr;
        }
 
        static std::shared_ptr<GenericSubmodule<T>> createMethode(SimulationClassBase<T>* sim)
        {
            return std::make_shared<SolarFluxVariablePotter2023Radiosity<T>>(sim);
        }
        static std::string getName() { return "SolarFluxVariablePotter2023Radiosity"; }
        std::string_view getNameLocal() const override
        {
            return "SolarFluxVariablePotter2023Radiosity";
        };
        std::vector<std::string> getDependencies() const override
        {
            return this->ini_file_data.getStringVectorParameters("dependencies");
        };
};
 
// ================= Implementation =================
 
template <typename T>
SolarFluxVariablePotter2023Radiosity<T>::SolarFluxVariablePotter2023Radiosity(
    SimulationClassBase<T>* sim)
    : GenericSubmodule<T>(sim)
{
    try
    {
        std::string ini_folder_path{
            this->sim->m_simulation_config.getStringParameters("ini_folder_path")};
        this->ini_file_data.loadUserInput(ini_folder_path + m_ini_filepath);
        this->m_generic_submodules = this->ini_file_data.getStringVectorParameters("submodules");
    }
    catch (const BadInput& e)
    {
        std::cerr << e.what() << '\n';
    }
}
 
template <typename T>
bool SolarFluxVariablePotter2023Radiosity<T>::setup(
    std::vector<std::shared_ptr<GenericSubmodule<T>>> all_submodules)
{
    m_shape_model_path = this->sim->m_simulation_config.getStringParameters("shape_model_path");
    bool calculate_viewfactor_matrix{true};
#ifdef MPI_PRESENT
    id_start = this->sim->m_simulation_config.getIntParameters("global_starting_facet_number");
    id_end = (this->sim->m_simulation_config.getIntParameters("global_starting_facet_number")
              + this->sim->m_simulation_config.getIntParameters("number_of_facets"));
    if (this->sim->world_size > 1)
    {
        if (this->sim->my_rank == 0)
        {
            // Calculate the Viewfactormatrix only once and broadcast it to all threads
            // to minimize VTK calls. (Will later be updated to using shared memory)
            m_ray_caster =
                std::make_unique<VtkRayCasting>(m_shape_model_path.c_str(), id_start,
                                                id_end - id_start, calculate_viewfactor_matrix);
            m_num_facets = m_ray_caster->getFacetNumber();
            m_number_non_zero_elements = m_ray_caster->getMatrixElementNumber();
        }
        else
        {
            calculate_viewfactor_matrix = false;
            m_ray_caster =
                std::make_unique<VtkRayCasting>(m_shape_model_path.c_str(), id_start,
                                                id_end - id_start, calculate_viewfactor_matrix);
        }
    }
    else
    {
        m_ray_caster = std::make_unique<VtkRayCasting>(m_shape_model_path.c_str(),
                                                       calculate_viewfactor_matrix);
        m_num_facets = m_ray_caster->getFacetNumber();
    }
#else
    m_ray_caster =
        std::make_unique<VtkRayCasting>(m_shape_model_path.c_str(), calculate_viewfactor_matrix);
    m_num_facets = m_ray_caster->getFacetNumber();
#endif
    m_numerical_layers = this->sim->m_simulation_config.getIntParameters("numerical_layers");
#ifdef MPI_PRESENT
    int number_of_facets_sim{
        this->sim->m_simulation_config.getIntParameters("global_number_of_facets")};
    m_global_temperature_array.resize(number_of_facets_sim * m_numerical_layers);
    if (this->sim->world_size > 1 && this->sim->my_rank > 0)
    {
        m_num_facets = number_of_facets_sim;
    }
    MPI_Bcast(&m_number_non_zero_elements, 1, MPI_INT, 0, MPI_COMM_WORLD);
#else
    int number_of_facets_sim{this->sim->m_simulation_config.getIntParameters("number_of_facets")};
    id_end = m_num_facets;
#endif
    if (m_num_facets != number_of_facets_sim)
    {
        throw std::runtime_error(
            "Number of facets does not match the number of facets on the shape model!\n The "
            "supplied number of facets is "
            + std::to_string(this->sim->m_simulation_config.getIntParameters("number_of_facets"))
            + ", while the number of facets on the shape model is " + std::to_string(m_num_facets)
            + ".\n");
    }
    direct_solar_flux = std::vector<double>(m_num_facets);
    m_flux_sol.resize(m_num_facets);
    m_albedo.resize(m_num_facets);
    m_emissivity.resize(m_num_facets);
    m_flux_radiative.resize(m_num_facets);
    m_flux_reflected.resize(m_num_facets);
    m_flux_infrared.resize(m_num_facets);
    m_flux_reflected_temp.resize(m_num_facets);
    m_flux_infrared_temp.resize(m_num_facets);
    m_heat_conduction_flux.resize(m_num_facets);
    m_solar_const = PhysicalConstants::solar_const;
    m_stefan_boltzmann_const = PhysicalConstants::stefan_boltzmann_const;
    return true;
}
 
template <typename T>
bool SolarFluxVariablePotter2023Radiosity<T>::exec(std::string_view param)
{
    if (param == "InitChain")
    {
        return init();
    }
    if (param == "PreTimeStepChain")
    {
        return preTimeStep();
    }
    if (param == "PostTimeStepChain")
    {
        return postTimeStep();
    }
    if (param == "OutputChain")
    {
        return output();
    }
    return false;
}
 
template <typename T>
void SolarFluxVariablePotter2023Radiosity<T>::radiosityCalcFirstTimeStep(
    const std::vector<T>& direct_solar_flux)
{
    for (int i{id_start}; i < id_end; i++)
    {
        m_flux_sol[i] = (1 - m_albedo[i]) * direct_solar_flux[i] + m_heat_conduction_flux[i];
        m_flux_infrared[i] = 0;
        m_flux_reflected[i] = 0;
    }
}
 
template <typename T>
void SolarFluxVariablePotter2023Radiosity<T>::radiosityCalc(const std::vector<T>& direct_solar_flux)
{
    for (int i{0}; i < m_num_facets; i++)
    {
        m_flux_reflected_temp[i] = direct_solar_flux[i] + m_flux_reflected[i];
        m_flux_infrared_temp[i] = m_flux_radiative[i] + (1 - m_emissivity[i]) * m_flux_infrared[i];
        m_flux_reflected[i] = 0;
        m_flux_infrared[i] = 0;
    }
    m_flux_reflected =
        m_viewfactor_matrix.matrixVectorMultiplication(m_flux_reflected_temp) * m_albedo;
    m_flux_infrared = m_viewfactor_matrix.matrixVectorMultiplication(m_flux_infrared_temp);
    for (int i{id_start}; i < id_end; i++)
    {
        m_flux_sol[i] = (1 - m_albedo[i]) * (direct_solar_flux[i] + m_flux_reflected[i])
                        + m_emissivity[i] * m_flux_infrared[i];
    }
}
 
template <typename T>
void SolarFluxVariablePotter2023Radiosity<T>::calculateFluxTotal(const std::valarray<T>& sol_vec,
                                                                 T initial_flux)
{
    double sol_vec_c[3] = {sol_vec[0], sol_vec[1], sol_vec[2]};
    direct_solar_flux = m_ray_caster->getSolarIrradiation(sol_vec_c, initial_flux);
    const std::valarray<T>& temperature_field{this->sim->getField("Temperature")};
    const std::valarray<T>& heat_conductivity_field{this->sim->getField("HeatConductivity")};
    T top_cell_length{this->sim->getFieldValue("CellLength", 0)};
    for (int i{0}; i < m_num_facets; i++)
    {
        m_flux_radiative[i] = m_emissivity[i] * m_stefan_boltzmann_const
                              * pow(temperature_field[i * m_numerical_layers], 4);
    }
    if (this->sim->time_step == 1)
    {
        for (int i{0}; i < m_num_facets; i++)
        {
            m_heat_conduction_flux[i] = heat_conductivity_field[i * m_numerical_layers]
                                        * (temperature_field[i * m_numerical_layers]
                                           - temperature_field[i * m_numerical_layers + 1])
                                        / top_cell_length;
        }
        radiosityCalcFirstTimeStep(direct_solar_flux);
    }
    else
    {
        radiosityCalc(direct_solar_flux);
    }
}
 
#ifdef MPI_PRESENT
template <typename T>
void SolarFluxVariablePotter2023Radiosity<T>::calculateFluxTotalMPI(const std::valarray<T>& sol_vec,
                                                                    T initial_flux)
{
    double sol_vec_c[3] = {sol_vec[0], sol_vec[1], sol_vec[2]};
    /*if (this->sim->my_rank == 0)
    {
        direct_solar_flux = m_ray_caster->getSolarIrradiation(sol_vec_c, initial_flux);
    }
    // This Bcast will be replaced by shared memory.
    MPI_Bcast(&direct_solar_flux[0], direct_solar_flux.size(), MPI_DOUBLE, 0, MPI_COMM_WORLD);*/
    direct_solar_flux = m_ray_caster->getSolarIrradiation(sol_vec_c, initial_flux);
    commPatternAllDataToAllInPlace(direct_solar_flux, id_end - id_start, this->sim);
    const std::valarray<T>& temperature_field{this->sim->getField("Temperature")};
    T top_cell_length{this->sim->getFieldValue("CellLength", 0)};
    m_flux_radiative = 0;
    for (int i{id_start}; i < id_end; i++)
    {
        m_flux_radiative[i] = m_emissivity[i] * m_stefan_boltzmann_const
                              * pow(temperature_field[(i - id_start) * m_numerical_layers], 4);
    }
    commPatternAllDataToAllInPlace(m_flux_radiative, id_end - id_start, this->sim);
    if (this->sim->time_step == 1)
    {
        std::valarray<T> heat_conductivity_field(m_numerical_layers * m_num_facets);
        commPatternAllDataToAll(heat_conductivity_field, this->sim->getField("HeatConductivity"),
                                this->sim);
        for (int i{0}; i < m_num_facets; i++)
        {
            m_heat_conduction_flux[i] = heat_conductivity_field[i * m_numerical_layers]
                                        * (m_global_temperature_array[i * m_numerical_layers]
                                           - m_global_temperature_array[i * m_numerical_layers + 1])
                                        / top_cell_length;
        }
        radiosityCalcFirstTimeStep(direct_solar_flux);
    }
    else
    {
        radiosityCalc(direct_solar_flux);
        commPatternAllDataToAllInPlace(m_flux_reflected, id_end - id_start, this->sim);
        commPatternAllDataToAllInPlace(m_flux_infrared, id_end - id_start, this->sim);
    }
}
 
#endif  // MPI_PRESENT
 
template <typename T>
void SolarFluxVariablePotter2023Radiosity<T>::calculateViewfactorMatrix()
{
    m_viewfactor_matrix = m_ray_caster->getViewfactorMatrix();
}
 
template <typename T>
bool SolarFluxVariablePotter2023Radiosity<T>::init()
{
#ifdef MPI_PRESENT
    if (this->sim->world_size > 1)
    {
        std::valarray<T> temp_viewfactor_matrix_values;
        std::valarray<int> temp_viewfactor_column_indices;
        std::valarray<int> temp_viewfactor_row_pointers;
        int number_of_facets_sim{
            this->sim->m_simulation_config.getIntParameters("global_number_of_facets")};
        if (this->sim->my_rank == 0)
        {
            calculateViewfactorMatrix();
            temp_viewfactor_matrix_values = m_viewfactor_matrix.getMatrixElements();
            temp_viewfactor_column_indices = m_viewfactor_matrix.getColumnIndices();
            temp_viewfactor_row_pointers = m_viewfactor_matrix.getRowPointers();
        }
        if (this->sim->world_size > 1 && this->sim->my_rank > 0)
        {
            temp_viewfactor_matrix_values.resize(m_number_non_zero_elements);
            temp_viewfactor_column_indices.resize(m_number_non_zero_elements);
            // One more, since row pointers always needs a end of matrix marker.
            temp_viewfactor_row_pointers.resize(number_of_facets_sim + 1);
        }
        // Broadcast the matrix values and element positions from rank 0 to all
        MPI_Request request_handles[3];
        MPI_Ibcast(&temp_viewfactor_matrix_values[0], temp_viewfactor_matrix_values.size(),
                   MPI_DOUBLE, 0, MPI_COMM_WORLD, &request_handles[0]);
        MPI_Ibcast(&temp_viewfactor_column_indices[0], temp_viewfactor_column_indices.size(),
                   MPI_INT, 0, MPI_COMM_WORLD, &request_handles[1]);
        MPI_Ibcast(&temp_viewfactor_row_pointers[0], temp_viewfactor_row_pointers.size(), MPI_INT,
                   0, MPI_COMM_WORLD, &request_handles[2]);
        MPI_Waitall(3, request_handles, MPI_STATUSES_IGNORE);
        // Use the creatiom routine from filled matrices in rank > 0 to get everyone to have the
        // viewfactor matrix.
        if (this->sim->world_size > 1 && this->sim->my_rank > 0)
        {
            m_viewfactor_matrix = CsrSparseMatrix<T>(
                temp_viewfactor_matrix_values, temp_viewfactor_column_indices,
                temp_viewfactor_row_pointers, number_of_facets_sim, number_of_facets_sim);
        }
        commPatternAllDataToAll(m_emissivity, this->sim->getField("Emissivity"), this->sim);
        commPatternAllDataToAll(m_albedo, this->sim->getField("Albedo"), this->sim);
    }
    else
    {
        calculateViewfactorMatrix();
        m_emissivity = this->sim->getField("Emissivity");
        m_albedo = this->sim->getField("Albedo");
    }
#else
    calculateViewfactorMatrix();
    m_albedo = this->sim->getField("Albedo");
    m_emissivity = this->sim->getField("Emissivity");
#endif
    return true;
}
 
template <typename T>
bool SolarFluxVariablePotter2023Radiosity<T>::preTimeStep()
{
    if (this->sim->time_step > m_check_timestep)
    {
        m_check_timestep = this->sim->time_step;
        const std::valarray<T>& sol_vec{this->sim->getField("SolarVector")};
        const std::valarray<T>& heliocentric_distance{this->sim->getField("HeliocentricDistance")};
        T initial_flux{m_solar_const / std::pow(heliocentric_distance[0], 2)};
 
#ifdef MPI_PRESENT
        if (this->sim->world_size > 1)
        {
            calculateFluxTotalMPI(sol_vec, initial_flux);
        }
        else
        {
            calculateFluxTotal(sol_vec, initial_flux);
        }
#else
        calculateFluxTotal(sol_vec, initial_flux);
#endif
    }
    (*this->module_field) = m_flux_sol[std::slice(id_start, id_end - id_start, 1)];
    return true;
}
 
template <typename T>
bool SolarFluxVariablePotter2023Radiosity<T>::m_registered =
    ModuleFactory<T>::registerModule(getName(), createMethode);
 
#endif  // VTK_PRESENT
 
#endif  // SOLAR_FLUX_VARIABLE_POTTER_2023_RADIOSITY_H