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 "StaticModuleContainer.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"};
 
        std::vector<std::shared_ptr<GenericManagingModule<T>>> m_all_modules{};
        std::vector<std::shared_ptr<GenericSubmodule<T>>> m_all_submodules{};
 
        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::shared_ptr<ShapeBase<T>> simulation_shape;
    try
    {
        shape_creator = std::make_unique<ShapeFactory<T>>(
            this, this->m_simulation_config.getBoolParameters("use_shape_file"),
            this->m_simulation_config.getStringParameters("spatial_dimension"));
        simulation_shape = shape_creator->getShapePtr();
    }
    catch (const BadInput& e)
    {
        // Map access errors.
        std::cerr << "[SimulationClassDefault]: Failed creating ShapeBase object from ShapeFactory "
                  << e.what() << '\n';
        throw BadInput(
            static_cast<std::string>(
                "[SimulationClassDefault]: Failed creating ShapeBase object from ShapeFactory ")
            + static_cast<std::string>(e.what()));
    }
    catch (const std::runtime_error& e)
    {
        // ShapeFactory internal errors.
        std::cerr << "[SimulationClassDefault]: Failed creating ShapeBase object from ShapeFactory "
                  << e.what() << '\n';
        throw std::runtime_error(
            static_cast<std::string>(
                "[SimulationClassDefault]: Failed creating ShapeBase object from ShapeFactory ")
            + static_cast<std::string>(e.what()));
    }
 
#ifdef MPI_PRESENT
    domainDecompositionOneDim(this);
#endif
 
    std::unique_ptr<GridFactory<T>> grid_creator;
    std::shared_ptr<GridBase<T>> simulation_grid;
    try
    {
        grid_creator = std::make_unique<GridFactory<T>>(
            this, simulation_shape,
            this->m_simulation_config.getStringParameters("spatial_dimension"));
        simulation_grid = grid_creator->getGridPtr();
    }
    catch (const BadInput& e)
    {
        // Map access errors or GridBase internal errors.
        std::cerr << "[SimulationClassDefault]: Failed creating GridBase object from GridFactory "
                  << e.what() << '\n';
        throw BadInput(
            static_cast<std::string>(
                "[SimulationClassDefault]: Failed creating GridBase object from GridFactory ")
            + static_cast<std::string>(e.what()));
    }
    catch (const std::runtime_error& e)
    {
        // GridFactory internal errors.
        std::cerr << "[SimulationClassDefault]: Failed creating GridBase object from GridFactory "
                  << e.what() << '\n';
        throw std::runtime_error(
            static_cast<std::string>(
                "[SimulationClassDefault]: Failed creating GridBase object from GridFactory ")
            + static_cast<std::string>(e.what()));
    }
 
#ifdef MPI_PRESENT
    // After (potentially) updating the numerical layers, the final splitting point of larger to
    // smaller processes can be calculated.
    try
    {
        this->global_array_split = this->number_of_larger_procs
                                   * (this->min_number_facets_per_proc + 1)
                                   * this->m_simulation_config.getIntParameters("numerical_layers");
    }
    catch (const BadInput& e)
    {
        std::cerr << "[SimulationClassDefault]: Failed creating GridBase object from GridFactory "
                  << e.what() << '\n';
        throw BadInput(
            static_cast<std::string>(
                "[SimulationClassDefault]: Failed creating GridBase object from GridFactory ")
            + static_cast<std::string>(e.what()));
    }
#endif
 
    try
    {
        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"))});
    }
    catch (const BadInput& e)
    {
        std::cerr << "[SimulationClassDefault]: Failed to create main independent variable field "
                     "(usually temperature) "
                  << e.what() << '\n';
        throw BadInput(static_cast<std::string>("[SimulationClassDefault]: Failed to create main "
                                                "independent variable field (usually temperature) ")
                       + static_cast<std::string>(e.what()));
    }
 
    std::shared_ptr<MatrixFactory<T>> matrix_factory;
    std::shared_ptr<ElementManager<T>> element_manager;
    std::shared_ptr<MatrixManager<T>> matrix_manager;
    try
    {
        matrix_factory = std::make_shared<MatrixFactory<T>>("Csr");
 
        element_manager = std::make_shared<ElementManager<T>>(
            simulation_grid->getGridSize(), simulation_grid->getCells(), matrix_factory);
 
        matrix_manager = std::make_shared<MatrixManager<T>>(simulation_grid->getGridSize(), this,
                                                            matrix_factory);
    }
    catch (const BadInput& e)
    {
        // Map access errors or MatrixManager internal errors.
        std::cerr << "[SimulationClassDefault]: Failed matrix build step " << e.what() << '\n';
        throw BadInput(
            static_cast<std::string>("[SimulationClassDefault]: Failed matrix build step ")
            + static_cast<std::string>(e.what()));
    }
    catch (const std::runtime_error& e)
    {
        // ElementManager internal errors.
        std::cerr << "[SimulationClassDefault]: Failed matrix build step at ElementManager "
                  << e.what() << '\n';
        throw std::runtime_error(
            static_cast<std::string>(
                "[SimulationClassDefault]: Failed matrix build step at ElementManager ")
            + static_cast<std::string>(e.what()));
    }
    catch (const std::invalid_argument& e)
    {
        // MatrixFactory internal errors.
        std::cerr << "[SimulationClassDefault]: Failed matrix build step at MatrixFactory "
                  << e.what() << '\n';
        throw std::runtime_error(
            static_cast<std::string>(
                "[SimulationClassDefault]: Failed matrix build step at MatrixFactory ")
            + static_cast<std::string>(e.what()));
    }
 
    ChainManager<GenericManagingModule<T>> global_chain{"global"};
 
    try
    {
        //  Create the modules and allow static access
        StaticModuleContainer<T> module_container{this};
        m_all_modules = StaticModuleContainer<T>::getAllSelectedManagingModules();
        m_all_submodules = StaticModuleContainer<T>::getAllSelectedSubmodules();
    }
    catch (const std::runtime_error& e)
    {
        std::cerr << "[SimulationClassDefault]: Failed at module or submodule creation step "
                  << e.what() << '\n';
        throw std::runtime_error(
            static_cast<std::string>(
                "[SimulationClassDefault]: Failed at module or submodule creation step ")
            + static_cast<std::string>(e.what()));
    }
 
    try
    {
        snapshot_data.overwriteIniFiles(this, m_all_modules, m_all_submodules);
    }
    catch (const std::runtime_error& e)
    {
        // This branch gets triggered when a module is in the snapshot but it was not loaded in this
        // simulation/ModuleFactory.cpp.
        std::cerr << "[SimulationClassDefault]: Failed at loading snapshot ini values due do "
                     "unloaded module that was requested in snapshot "
                  << e.what() << '\n';
        throw std::runtime_error(
            static_cast<std::string>(
                "[SimulationClassDefault]: Failed at loading snapshot ini values ")
            + static_cast<std::string>(e.what()));
    }
    catch (const BadInput& e)
    {
        // This branch should only be triggered when some sort of corruption occured to the snapshot
        // data.
        std::cerr << "[SimulationClassDefault]: Failed at loading snapshot ini values " << e.what()
                  << '\n';
        throw std::runtime_error(
            static_cast<std::string>(
                "[SimulationClassDefault]: Failed at loading snapshot ini values ")
            + static_cast<std::string>(e.what()));
    }
 
    try
    {
        // If a user supplied managing module list is present, use this
        global_chain.autoLoader(
            m_all_modules, this->m_simulation_config.getStringVectorParameters("managing_modules"));
    }
    catch (const DependencyError<GenericManagingModule<T>>& e)
    {
        std::cerr
            << "[SimulationClassDefault]: Failed to resolve dependencies for selected modules "
            << e.what() << '\n';
        throw std::runtime_error(
            static_cast<std::string>("[SimulationClassDefault]: (DependencyError) Failed to "
                                     "resolve dependencies for selected modules ")
            + static_cast<std::string>(e.what()));
    }
    catch (const BadInput& e)
    {
        std::cerr
            << "[SimulationClassDefault]: Warning - Failed to load managing modules from config: "
            << e.what() << "\n"
            << "Falling back to loading all modules.\n";
        try
        {
            global_chain.autoLoader(m_all_modules);
        }
        catch (const DependencyError<GenericManagingModule<T>>& e)
        {
            std::cerr << "[SimulationClassDefault]: Failed to resolve dependencies for all modules "
                      << e.what() << '\n';
            throw std::runtime_error(
                static_cast<std::string>("[SimulationClassDefault]: (DependencyError) Failed to "
                                         "resolve dependencies for all modules ")
                + static_cast<std::string>(e.what()));
        }
    }
    try
    {
        runSetup<T>(m_all_modules, m_all_submodules, global_chain.getManagingModuleNames());
    }
    catch (const BadInput& e)
    {
        // runSetup itself does not throw an error but has a high chance to trigger a BadInput
        // within the setup functions of the modules.
        std::cerr << "[SimulationClassDefault]: Failed at setting up the modules " << e.what()
                  << '\n';
        throw std::runtime_error(
            static_cast<std::string>("[SimulationClassDefault]: Failed at setting up the modules ")
            + static_cast<std::string>(e.what()));
    }
 
    try
    {
        snapshot_data.overwriteSimFields(this);
    }
    catch (const BadInput& e)
    {
        // This branch should only be triggered when some sort of corruption occured to the snapshot
        // data.
        std::cerr << "[SimulationClassDefault]: Failed at loading snapshot field values "
                  << e.what() << '\n';
        throw std::runtime_error(
            static_cast<std::string>(
                "[SimulationClassDefault]: Failed at loading snapshot field values ")
            + static_cast<std::string>(e.what()));
    }
 
    // 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;
    try
    {
        solver = std::make_shared<TridiagonalMatrixSolver<T>>(this, matrix_manager);
    }
    catch (const MatrixOutOfBoundsError& e)
    {
        // If there somehow was a mismatch between elements and matrix size this would trigger.
        std::cerr << "[SimulationClassDefault]: Failed at creating Solver " << e.what() << '\n';
        throw std::runtime_error(
            static_cast<std::string>(
                "[SimulationClassDefault]: (MatrixOutOfBoundsError) Failed at creating Solver ")
            + static_cast<std::string>(e.what()));
    }
    catch (const StaticSparseMatrixAccessError& e)
    {
        // This should only trigger if there is a mismatch of solver type and grid type.
        std::cerr << "[SimulationClassDefault]: Failed at creating Solver " << e.what() << '\n';
        throw std::runtime_error(
            static_cast<std::string>("[SimulationClassDefault]: (StaticSparseMatrixAccessError) "
                                     "Failed at creating Solver ")
            + static_cast<std::string>(e.what()));
    }
 
    int timesteps, inner_iterations;
    T epsilon_abs, epsilon_rel, delta_t;
#ifdef MPI_PRESENT
    T max_var_previous, max_var_local;
    MPI_Request request_handles[2];
#endif
    try
    {
        timesteps = static_cast<int>(this->m_simulation_config.getDoubleParameters("time_max")
                                     / this->m_simulation_config.getDoubleParameters("time_delta"));
        inner_iterations = this->m_simulation_config.getIntParameters("inner_iterations");
        epsilon_abs = this->m_simulation_config.getDoubleParameters("epsilon_inner_iterations");
        epsilon_rel =
            this->m_simulation_config.getDoubleParameters("relative_epsilon_inner_iterations");
        delta_t = this->m_simulation_config.getDoubleParameters("time_delta");
    }
    catch (const BadInput& e)
    {
        std::cerr << "[SimulationClassDefault]" << e.what() << '\n';
        throw BadInput(static_cast<std::string>("[SimulationClassDefault]")
                       + static_cast<std::string>(e.what()));
    }
 
    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_Iallreduce(&this->local_convergence_diff, &this->global_convergence_diff, 1,
                           MPI_DOUBLE, MPI_MAX, MPI_COMM_WORLD, &request_handles[0]);
            max_var_previous = dependent_variable_previous_iteration.max();
            max_var_local = (this->m_field_map[m_dependent_variable_name]).max();
            if (max_var_previous > max_var_local)
            {
                max_var_local = max_var_previous;
            }
            MPI_Iallreduce(&max_var_local, &this->max_var_value, 1, MPI_DOUBLE, MPI_MAX,
                           MPI_COMM_WORLD, &request_handles[1]);
            MPI_Waitall(2, request_handles, MPI_STATUSES_IGNORE);
            if (this->global_convergence_diff < epsilon_abs || inner_iterations == 1)
            {
                break;
            }
            else if (this->global_convergence_diff / this->max_var_value < epsilon_rel)
            {
                break;
            }
#else
            this->global_convergence_diff = abs(dependent_variable_previous_iteration
                                                - this->m_field_map[m_dependent_variable_name])
                                                .max();
            if (this->global_convergence_diff < epsilon_abs || inner_iterations == 1)
            {
                break;
            }
            // Relative criterion: normalise the inter-iteration difference by the
            // larger of the previous and current field maxima. This matches the
            // MPI path, where max_var_value is reduced with MPI_MAX.
            else if (this->global_convergence_diff
                         / std::max(dependent_variable_previous_iteration.max(),
                                    (this->m_field_map[m_dependent_variable_name]).max())
                     < epsilon_rel)
            {
                break;
            }
#endif
            if (iter < inner_iterations - 1)
            {
                dependent_variable_previous_iteration =
                    this->m_field_map[m_dependent_variable_name];
            }
            else
            {
                std::cerr << "[SimulationClassDefault]: Warning - 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, m_all_modules, m_all_submodules};
    // this->printField("HeatConductivity");
    // this->printField("HeatCapacity");
    // this->printField("Density");
 
    return 0;
}
 
#endif