MatrixManager.h

Go to the documentation of this file.

#ifndef MATRIX_MANAGER_H
#define MATRIX_MANAGER_H
 
#include <algorithm>
#include <memory>
#include <string>
#include <vector>
 
#include "GenericMatrix.h"
#include "MatrixFactory.h"
#include "SimulationClassBase.h"
 
template <typename T>
class MatrixManager
{
    private:
        std::vector<std::pair<int, int>> m_non_zero_elements;
        SimulationClassBase<T>* m_sim;
 
        void flushMatrices();
 
    public:
        // It only makes sense for them to be public as I need to get references to these objects.
        std::shared_ptr<GenericMatrix<T>> lhs_matrix{};
        std::valarray<T> rhs_vector{};
 
        MatrixManager(int size, SimulationClassBase<T>* sim,
                      std::shared_ptr<MatrixFactory<T>> matrix_factory);
 
        void assembleGlobalLse(std::shared_ptr<GenericMatrix<T>> capacitance_matrix,
                               std::shared_ptr<GenericMatrix<T>> stiffness_matrix,
                               std::valarray<T>& forcing_vector);
};
 
template <typename T>
MatrixManager<T>::MatrixManager(const int size, SimulationClassBase<T>* sim,
                                std::shared_ptr<MatrixFactory<T>> matrix_factory)
    : m_sim{sim}, lhs_matrix(matrix_factory->createShared())
{
    // CSR build-up needs to be postponed until elements could be registered.
    if (lhs_matrix->getMatrixType() == MatrixType::Dense)
    {
        lhs_matrix->buildAtRuntime(size, size);
    }
    rhs_vector.resize(size);
}
 
template <typename T>
void MatrixManager<T>::assembleGlobalLse(std::shared_ptr<GenericMatrix<T>> capacitance_matrix,
                                         std::shared_ptr<GenericMatrix<T>> stiffness_matrix,
                                         std::valarray<T>& forcing_vector)
{
    if (m_non_zero_elements.empty())
    {
        stiffness_matrix->updateNonZeroElements();
        m_non_zero_elements = stiffness_matrix->getNonZeroElements();
        if ((lhs_matrix->getMatrixType() == MatrixType::Csr))
        {
            for (const auto& pos : m_non_zero_elements)
            {
                lhs_matrix->registerElement(0.0, pos.first, pos.second);
            }
            int dim{(stiffness_matrix->size()).first};
            lhs_matrix->buildAtRuntime(dim, dim);
        }
    }
    T delta_t{m_sim->m_simulation_config.getDoubleParameters("time_delta")};
    const std::valarray<T>& temperature{m_sim->m_field_map.at("Temperature")};
    // Old dense matrix assembly only left for comparisons to MoCSI v1.0.0. Deprecated and highly
    // inefficient!
    if (lhs_matrix->getMatrixType() == MatrixType::Dense)
    {
        flushMatrices();
        for (const auto& [row, col] : m_non_zero_elements)
        {
            (*lhs_matrix)(row, col) =
                (*capacitance_matrix)(row, col) + (*stiffness_matrix)(row, col) * delta_t;
            if (row == col)
            {
                rhs_vector[row] += forcing_vector[row] * delta_t;
            }
            rhs_vector[row] += (*capacitance_matrix)(row, col) * temperature[col];
        }
    }
    else
    {
        // lhs_matrix = capacitance_matrix + stiffness_matrix * delta_t
        lhs_matrix->overwriteMatrixValues(
            capacitance_matrix->vectorElemAdd(stiffness_matrix->scalarMultiply(delta_t)));
        // rhs_vector = capacitance_matrix * temperature(vector) + forcing_vector * delta_t
        rhs_vector =
            capacitance_matrix->matrixVectorMultiplication(temperature) + forcing_vector * delta_t;
    }
}
 
template <typename T>
void MatrixManager<T>::flushMatrices()
{
    lhs_matrix->zeroMatrix();
    rhs_vector *= 0.0;
}
 
#endif