SimulationClassDefault.h

Go to the documentation of this file.

#ifndef SIMULATION_CLASS_DEFAULT_H
#define SIMULATION_CLASS_DEFAULT_H
 
#include <iostream>
#include <map>
#include <memory>
#include <string>
#include <valarray>
#include <vector>
 
#include "ChainManager.h"
#include "CsvParser.h"
#include "ElementManager.h"
#include "GenericSolver.h"
#include "GridFactory.h"
#include "InputManager.h"
#include "MatrixFactory.h"
#include "MatrixManager.h"
#include "ModuleFactory.h"
#include "SetupRunner.h"
#include "ShapeFactory.h"
#include "SimulationClassBase.h"
#include "SnapshotCreator.h"
#include "SnapshotLoader.h"
#include "TridiagonalMatrixSolver.h"
 
#ifdef MPI_PRESENT
#include <mpi.h>
 
#include "MpiCommunicationPatterns.h"
#include "MpiDomainDecomposition.h"
#endif
 
template <typename T>
class SimulationClassDefault : public SimulationClassBase<T>
{
    public:
        std::string m_dependent_variable_name{"Temperature"};
        std::string m_dependent_variable_initial_state{"initial_temperature"};
 
        SimulationClassDefault();
        SimulationClassDefault(int argc, const char* argv[]);
        SimulationClassDefault(const std::string& user_ini);
 
        int run() override;
 
        ~SimulationClassDefault() override = default;
};
 
template <typename T>
SimulationClassDefault<T>::SimulationClassDefault()
{
}
 
template <typename T>
SimulationClassDefault<T>::SimulationClassDefault(const int argc, const char* argv[])
    : SimulationClassBase<T>(argc, argv)
{
}
 
template <typename T>
SimulationClassDefault<T>::SimulationClassDefault(const std::string& user_ini)
    : SimulationClassBase<T>(user_ini)
{
}
 
template <typename T>
int SimulationClassDefault<T>::run()
{
    SnapshotLoader<T> snapshot_data{this};
    snapshot_data.overwriteSimConfig(this);
 
    // Loading of the shape model so that all the mesh info is available at grid creation.
    std::unique_ptr<ShapeFactory<T>> shape_creator{std::make_unique<ShapeFactory<T>>(
        this, this->m_simulation_config.getBoolParameters("use_shape_file"),
        this->m_simulation_config.getStringParameters("spatial_dimension"))};
    std::shared_ptr<ShapeBase<T>> simulation_shape = shape_creator->getShapePtr();
 
#ifdef MPI_PRESENT
    domainDecompositionOneDim(this);
#endif
 
    std::unique_ptr<GridFactory<T>> grid_creator{std::make_unique<GridFactory<T>>(
        this, simulation_shape,
        this->m_simulation_config.getStringParameters("spatial_dimension"))};
    std::shared_ptr<GridBase<T>> simulation_grid = grid_creator->getGridPtr();
 
#ifdef MPI_PRESENT
    // After (potentially) updating the numerical layers, the final splitting point of larger to
    // smaller processes can be calculated.
    this->global_array_split = this->number_of_larger_procs * (this->min_number_facets_per_proc + 1)
                               * this->m_simulation_config.getIntParameters("numerical_layers");
#endif
 
    this->m_field_map.insert(
        {m_dependent_variable_name,
         std::valarray<T>(
             this->m_simulation_config.getDoubleParameters(m_dependent_variable_initial_state),
             this->m_simulation_config.getIntParameters("numerical_layers")
                 * this->m_simulation_config.getIntParameters("number_of_facets"))});
 
    std::shared_ptr<MatrixFactory<T>> matrix_factory{std::make_shared<MatrixFactory<T>>("Csr")};
 
    std::shared_ptr<ElementManager<T>> element_manager{std::make_shared<ElementManager<T>>(
        simulation_grid->getGridSize(), simulation_grid->getCells(), matrix_factory)};
 
    std::shared_ptr<MatrixManager<T>> matrix_manager{
        std::make_shared<MatrixManager<T>>(simulation_grid->getGridSize(), this, matrix_factory)};
 
    ChainManager<GenericManagingModule<T>> global_chain{"global"};
    // Load all managing modules if no list is specified or only the listed ones if it is.
    std::vector<std::shared_ptr<GenericManagingModule<T>>> all_modules;
    try
    {
        all_modules = ModuleFactory<T>::createAllSelectedManagingModules(
            this, this->m_simulation_config.getStringVectorParameters("managing_modules"));
    }
    catch (const BadInput&)
    {
        all_modules = ModuleFactory<T>::createAllManagingModules(this);
    }
    std::vector<std::shared_ptr<GenericSubmodule<T>>> all_submodules{
        ModuleFactory<T>::createAllSelectedSubmodules(this, all_modules)};
 
    snapshot_data.overwriteIniFiles(this, all_modules, all_submodules);
 
    try
    {
        // If a user supplied managing module list is present, use this
        global_chain.autoLoader(
            all_modules, this->m_simulation_config.getStringVectorParameters("managing_modules"));
    }
    catch (const std::exception& e)
    {
        std::cerr << "[Warning] Failed to load managing modules from config: " << e.what() << "\n"
                  << "Falling back to loading all modules.\n";
        global_chain.autoLoader(all_modules);
    }
    runSetup<T>(all_modules, all_submodules, global_chain.getManagingModuleNames());
 
    snapshot_data.overwriteSimFields(this);
 
    // Insertion point "Init"
    this->setInitChainStr();
    global_chain.runChain(this->m_current_position);
 
    // this->printField("HeliocentricDistance");
 
    // Solver reads matrix size from sim class, so has to be created after matrix size is set.
    std::shared_ptr<SolverInterface<T>> solver{
        std::make_shared<TridiagonalMatrixSolver<T>>(this, matrix_manager)};
 
    // Hardcoded timesteps here, of course taken from ini in real class
    int timesteps{static_cast<int>(this->m_simulation_config.getDoubleParameters("time_max")
                                   / this->m_simulation_config.getDoubleParameters("time_delta"))};
    int inner_iterations{this->m_simulation_config.getIntParameters("inner_iterations")};
    T epsilon{this->m_simulation_config.getDoubleParameters("epsilon_inner_iterations")};
    T delta_t{this->m_simulation_config.getDoubleParameters("time_delta")};
 
    while (this->time_step * delta_t < this->m_simulation_config.getDoubleParameters("time_max"))
    {
        std::valarray<T> dependent_variable_previous_timestep{
            this->m_field_map[m_dependent_variable_name]};
        std::valarray<T> dependent_variable_previous_iteration{
            this->m_field_map[m_dependent_variable_name]};
        this->elapsed_time += delta_t;
        this->time_step++;
 
        // this->m_field_map["SolarVector"][0] = cos(this->time_step * 0.36 * std::numbers::pi /
        // 180); this->m_field_map["SolarVector"][2] = sin(this->time_step * 0.36 * std::numbers::pi
        // / 180);
        for (int iter{0}; iter < inner_iterations; iter++)
        {
            // Insertion point "PreTimeStep"
            this->setPreTimeStepChainStr();
            global_chain.runChain(this->m_current_position);
 
            if (iter < inner_iterations - 1 && iter > 0)
            {
                // Change temperature field after the module run, so they use the updated
                // temperatures!
                this->m_field_map[m_dependent_variable_name] = dependent_variable_previous_timestep;
            }
 
            // Do something with the temperature field
            element_manager->assembleFemMatrices();
            matrix_manager->assembleGlobalLse(element_manager->m_capacitance_matrix,
                                              element_manager->m_stiffness_matrix,
                                              element_manager->m_forcing_vector);
            solver->solveMatrixSystem(m_dependent_variable_name);
 
            // Insertion point "PostNonLinIter"
            this->setPostNonLinIterChainStr();
            global_chain.runChain(this->m_current_position);
 
#ifdef MPI_PRESENT
            this->local_convergence_diff = (abs(dependent_variable_previous_iteration
                                                - this->m_field_map[m_dependent_variable_name]))
                                               .max();
            MPI_Allreduce(&this->local_convergence_diff, &this->global_convergence_diff, 1,
                          MPI_DOUBLE, MPI_MAX, MPI_COMM_WORLD);
            if (this->global_convergence_diff < epsilon || inner_iterations == 1)
            {
                break;
            }
#else
            if ((abs(dependent_variable_previous_iteration
                     - this->m_field_map[m_dependent_variable_name]))
                        .max()
                    < epsilon
                || inner_iterations == 1)
            {
                break;
            }
#endif
            if (iter < inner_iterations - 1)
            {
                dependent_variable_previous_iteration =
                    this->m_field_map[m_dependent_variable_name];
            }
            else
            {
                std::cerr << "Reached specified inner iterations limit without fullfilling the "
                             "requested convergence criterion.\n";
            }
        }
        // this->time_step++;
        // Insertion point "PostTimeStep"
        this->setPostTimeStepChainStr();
        global_chain.runChain(this->m_current_position);
 
        delta_t = this->m_simulation_config.getDoubleParameters("time_delta");
 
        // this->printField(m_dependent_variable_name);
    }
 
    // Insertion point "Output"
    this->setOutputChainStr();
    global_chain.runChain(this->m_current_position);
 
#ifdef MPI_PRESENT
    if (this->world_size > 1)
    {
        std::valarray<T> global_temperature(
            this->m_simulation_config.getIntParameters("numerical_layers")
            * this->m_simulation_config.getIntParameters("global_number_of_facets"));
        commPatternAllDataToAll(global_temperature, this->getField("Temperature"), this);
        if (this->my_rank == 0)
        {
            std::cout << "Full temperature:\n";
            for (const auto& val : global_temperature)
            {
                std::cout << val << " ";
            }
            std::cout << "\n";
        }
    }
    else
    {
        this->printField(m_dependent_variable_name);
    }
#else
    this->printField(m_dependent_variable_name);
#endif
 
    SnapshotCreator{this, all_modules, all_submodules};
    // this->printField("HeatConductivity");
    // this->printField("HeatCapacity");
    // this->printField("Density");
 
    return 0;
}
 
#endif