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"
 
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;
 
    private:
        std::string m_init_str{"InitChain"};
        std::string m_preTs_str{"PreTimeStepChain"};
        std::string m_postNonLinIter_str{"PostNonLinIterChain"};
        std::string m_postTs_str{"PostTimeStepChain"};
        std::string m_output_str{"OutputChain"};
};
 
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();
    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();
 
    this->m_field_map.insert(
        {m_dependent_variable_name,
         std::valarray<double>(
             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"};
    std::vector<std::shared_ptr<GenericManagingModule<T>>> all_modules{
        ModuleFactory<T>::createAllManagingModules(this)};
    std::vector<std::shared_ptr<GenericSubmodule<T>>> all_submodules{
        ModuleFactory<T>::createAllSubmodules(this)};
 
    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;
            }
            /*this->printField("HeatConductivity");
            this->printField("HeatCapacity");
            this->printField("Density");
            this->printField("QSol");
            this->print_parameter("HeliocentricDistance");
            this->printField("SolarVector");*/
 
            // 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);
 
            if ((abs(dependent_variable_previous_iteration
                     - this->m_field_map[m_dependent_variable_name]))
                        .max()
                    < epsilon
                || inner_iterations == 1)
            {
                break;
            }
            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);
 
    this->printField(m_dependent_variable_name);
 
    SnapshotCreator{this, all_modules, all_submodules};
    // this->printField("HeatConductivity");
    // this->printField("HeatCapacity");
    // this->printField("Density");
 
    return 0;
}
 
#endif