b2typed_solver.H

//------------------------------------------------------------------------
// b2typed_solver.H --
//
//
// written by Neda Ebrahimi Pour <neda.ebrahimipour@dlr.de>
//            Harald Klimach <harald.klimach@dlr.de>
//            Thomas Blome <thomas.blome@dlr.de>
//
// (c) 2004-2012 SMR Engineering & Development SA
//               2502 Bienne, Switzerland
//
// (c) 2023 Deutsches Zentrum für Luft- und Raumfahrt (DLR) e.V.
//          Linder Höhe, 51147 Köln
//
// All Rights Reserved.  Proprietary source code.  The contents of
// this file may not be disclosed to third parties, copied or
// duplicated in any form, in whole or in part, without the prior
// written permission of SMR.
//------------------------------------------------------------------------
 
#ifndef B2TYPED_SOLVER_H_
#define B2TYPED_SOLVER_H_
 
#include <atomic>
#include <memory>
#include <type_traits>
 
#include "b2ppconfig.h"
// #include "model/b2boundary_condition.H"
#include "io/b2000_db/b2000_db_v3/b2solution_database.H"
#include "model/b2solution.H"
#include "model/b2solver.H"
#include "model_imp/b2monolithic_boundary_conditions.H"
#include "time_integration/b2time_integrator.H"
#include "utils/b2logging.H"
#include "utils/b2object.H"
#include "utils/b2timing.H"
 
#ifdef HAVE_LIB_TBB
#include "utils/b2threading.H"
#endif
 
namespace b2000::solver {
 
template <typename T, typename MATRIX_FORMAT>
class TypedSolver : public Solver {
public:
    // using mrbc_t = TypedModelReductionBoundaryCondition<T>;
    using EssentialBoundaryConditionPointer =
          std::unique_ptr<MonolithicEssentialBoundaryConditions<T>>;
    using NaturalBoundaryConditionPointer =
          std::unique_ptr<MonolithicNaturalBoundaryConditions<T, MATRIX_FORMAT>>;
    using type_t = ObjectTypeIncomplete<TypedSolver, Solver::type_t>;
 
    using type_scalar = T;
    using type_matrix = MATRIX_FORMAT;
 
    void solve() override = 0;
 
    const auto& GetDof() const { return dof_; }
 
    const auto& GetKeff() const { return K_eff_; }
 
    const auto& GetFeff() const { return F_eff_; }
 
    auto GetTimeIntegrator() const { return time_integrator_; }
 
    static type_t type;
 
protected:
    enum class b2Task {
        assemble_mass_matrix,
        assemble_effective_system,
        assemble_effective_matrix,
        assemble_effective_vector,
        compute_gradients,
        compute_J_integral,
        update_internal_variables,
        compute_error,
    };
 
    virtual void InitializeSolverOnCase() override;
 
    bool solve_iteration() override = 0;
 
    void MainLoopOverElements(const b2Task task);
 
    void SetSolutionFieldsToZero();
 
    void SolveEffectiveSystem();
 
 
    b2linalg::Matrix<T, b2linalg::Mrectangle> dof_{};
 
    b2linalg::Matrix<T, b2linalg::Mrectangle> dU_{};
 
    b2linalg::Matrix<T, typename MATRIX_FORMAT::sparse_inversible>
          K_eff_{};                                  
    b2linalg::Vector<T, b2linalg::Vdense> F_eff_{};  
 
    EssentialBoundaryConditionPointer essential_boundary_conditions_{};
    NaturalBoundaryConditionPointer natural_boundary_conditions_{};
 
    std::shared_ptr<TimeIntegrator<T>> time_integrator_{};
    std::shared_ptr<TimeIntegrator<T>> corrector_{};
 
    std::vector<int> fixed_dof_{};
 
    T reduced_scalar_value{};
 
    SolutionMonolithic<T>* solution_{nullptr};
 
private:
    enum class NewtonMethod { conventional, modified, delayed_modified, modified_postponed };
 
    void IterateOverElements(auto& elements, auto func, auto timer, std::string_view log);
 
    void SeriallyIterateOverElements(auto& elements, auto func, auto timer, std::string_view log);
 
    void IterateOverElementsParallelReduction(
          auto& elements, auto func, auto reduction, T identity, auto timer, std::string_view log);
 
    void AllocateAndAssembleEffectiveSystem(std::vector<std::reference_wrapper<Element>>& elements);
 
    void AssembleEffectiveSystem(TypedElement<T>& element);
 
    void AssembleEffectiveMatrix(TypedElement<T>& element);
 
    void AssembleEffectiveVector(TypedElement<T>& element);
 
    void ComputeGradients(TypedElement<T>& element, GradientContainer* const gradient_container);
 
    void ComputeError(
          std::vector<std::reference_wrapper<Element>>& elements, auto timer, std::string_view log);
};
 
template <typename T, typename MATRIX_FORMAT>
void TypedSolver<T, MATRIX_FORMAT>::InitializeSolverOnCase() {
    // if (case_iterator_.get() == nullptr)
    //     case_iterator_ = model_->get_case_list().get_case_iterator();
 
    // case_ = case_iterator_->next();
    // if (case_ == nullptr) return;
 
    // if (case_->get_number_of_stage() == 0)
    //     Exception() << "The monolithic solver must have a case that contains
    //     at least "
    //                    "one stage but the case "
    //                 << case_->get_id() << " does not have a stage." << THROW;
 
    Solver::InitializeSolverOnCase();
 
    solution_ = &dynamic_cast<SolutionMonolithic<T>&>(model_->get_solution());
 
    // const AttributeList& attribute = *case_;
    // std::string analysis = case_->get_string("ANALYSIS", "LINEAR");
 
    essential_boundary_conditions_ = std::make_unique<MonolithicEssentialBoundaryConditions<T>>();
    natural_boundary_conditions_ =
          std::make_unique<MonolithicNaturalBoundaryConditions<T, MATRIX_FORMAT>>();
    // model_->InitializeMonolithicEssentialBoundaryConditions(*case_);
 
    number_of_fixed_dof_ =
          essential_boundary_conditions_->InitializeEssentialBoundaryConditions(*model_);
    natural_boundary_conditions_->InitializeNaturalBoundaryConditions(*model_, *case_);
 
    // mrbc_t& mrbc = dynamic_cast<mrbc_t&>(
    //       model_->get_model_reduction_boundary_condition(mrbc_t::type.get_subtype(
    //             "TypedModelReductionBoundaryConditionListComponent" +
    //             type_name<T>())));
 
    // b2linalg::Matrix<T, b2linalg::Mrectangle> dof_red(
    //       mrbc.get_size(false, b2linalg::Vector<T, b2linalg::Vdense>::null,
    //       time_), dof_.size2());
    // mrbc.get_nonlinear_inverse_value(dof_, time_, dof_red);
 
    // mrbc.InitializeEssentialBoundaryConditions(fixed_dof_, dof_);
 
    effective_num_dof_ = total_num_dof_ - number_of_fixed_dof_;
 
    fixed_dof_.resize(total_num_dof_);
    std::iota(fixed_dof_.begin(), fixed_dof_.end(), 0);
 
    size_t columns{1};  
    // if (is_dynamic_) {
    //     columns = 3;  //!< Displacements & velocities and accelerations.
    // } else if (is_damped_) {
    //     columns = 2;  //!< Displacements & velocities.
    // }
 
    F_eff_.reset(effective_num_dof_);
 
    dof_.resize(total_num_dof_, columns);
    dU_.resize(total_num_dof_, 2);
 
    dU_.set_zero();
    dof_.set_zero();
    K_eff_.set_zero();
}
 
template <typename T, typename MATRIX_FORMAT>
void TypedSolver<T, MATRIX_FORMAT>::IterateOverElements(
      auto& elements, auto func, auto timer, std::string_view log) {
    logging::Logger& logger = logging::get_logger("solver.monolithic");
 
    logger << logging::info << "Starting to " << log << logging::LOGFLUSH;
 
    (*timer).second.start();
#ifdef HAVE_LIB_TBB
    // tbb::combinable<T> private_value;
    // T reduced_scalar_value{};
 
    // Optionally, specify a partitioner as last argument
    tbb::parallel_for(static_cast<size_t>(0), elements.size(), func);
#else
    // std::unique_ptr<Domain::ElementIterator> i =
    // domain.get_element_iterator(); Func b2task(*this, elements, task);
 
    // for (Element* element = i->next(); element != 0; element = i->next()) {
    for (size_t i{}; i < elements.size(); ++i) { func(i); }
#endif
    (*timer).second.stop();
 
    logger << logging::info << "Finished to " << log << logging::LOGFLUSH;
}
 
template <typename T, typename MATRIX_FORMAT>
void TypedSolver<T, MATRIX_FORMAT>::SeriallyIterateOverElements(
      auto& elements, auto func, auto timer, std::string_view log) {
    logging::Logger& logger = logging::get_logger("solver.monolithic");
 
    logger << logging::info << "Serially Starting to " << log << logging::LOGFLUSH;
 
    (*timer).second.start();
    for (size_t i{}; i < elements.size(); ++i) { func(i); }
    (*timer).second.stop();
 
    logger << logging::info << "Serially Finished to " << log << logging::LOGFLUSH;
}
 
#ifdef HAVE_LIB_TBB
template <typename T, typename MATRIX_FORMAT>
void TypedSolver<T, MATRIX_FORMAT>::IterateOverElementsParallelReduction(
      auto& elements, auto func, auto reduction, T identity, auto timer, std::string_view log) {
    logging::Logger& logger = logging::get_logger("solver.monolithic");
 
    logger << logging::info << "Starting to " << log << logging::LOGFLUSH;
 
    (*timer).second.start();
    reduced_scalar_value = tbb::parallel_reduce(
          tbb::blocked_range<size_t>(0, elements.size()), identity, func, reduction);
    (*timer).second.stop();
 
    logger << logging::info << "Finished to " << log << logging::LOGFLUSH;
}
#endif
 
template <typename T, typename MATRIX_FORMAT>
void TypedSolver<T, MATRIX_FORMAT>::MainLoopOverElements(const b2Task task) {
    // domain.set_dof(dof);
    // b2linalg::Index index;
    //
    // Domain& domain = model_->get_domain();
    auto& elements = model_->get_domain().GetElementContainer();
 
    // using test = void (TypedSolver<T,
    // MATRIX_FORMAT>::*methodPtr)(TypedElement<T>&);
 
    // std::function<void(TypedSolver<T, MATRIX_FORMAT>, TypedElement<T> &
    // element)> f1{
    //       &TypedSolver<T, MATRIX_FORMAT>::AssembleEffectiveMatrix};
 
    // void (TypedSolver<T, MATRIX_FORMAT>::*methodPtr)(TypedElement<T>&){
    //       &TypedSolver<T, MATRIX_FORMAT>::AssembleEffectiveMatrix};
 
    // b2000::b2dbv3::Domain& domain_test =
    // dynamic_cast<b2dbv3::Domain&>(model->get_domain());
    // b2000::Domain& b2_domain = dynamic_cast<b2000::Domain&>(domain);
 
    switch (task) {
        case (b2Task::assemble_mass_matrix): {
            auto time_assemble_mass_mat = obtain_timer("solver_assemble_mass_mat");
            break;
        }
            // clang-format off
        [[likely]] case (b2Task::assemble_effective_system):
            if (K_eff_.size1() == 0) {
                AllocateAndAssembleEffectiveSystem(elements);
            } else {
            IterateOverElements(
                elements,
                [this, &elements](const size_t i) {
                    AssembleEffectiveSystem(dynamic_cast<TypedElement<T>&>(elements[i].get()));
                },
                obtain_timer("solver_assemble_eff_sys"), "assemble the effective system");
            }
            
            break;
            // clang-format on
 
        case (b2Task::assemble_effective_matrix):
            IterateOverElements(
                  elements,
                  [this, &elements](const size_t i) {
                      AssembleEffectiveMatrix(dynamic_cast<TypedElement<T>&>(elements[i].get()));
                  },
                  obtain_timer("solver_assemble_eff_mat"), "assemble the effective matrix");
 
            break;
 
        case (b2Task::assemble_effective_vector):
            IterateOverElements(
                  elements,
                  [this, &elements](const size_t i) {
                      AssembleEffectiveVector(dynamic_cast<TypedElement<T>&>(elements[i].get()));
                  },
                  obtain_timer("solver_assemble_eff_vec"), "assemble the effective vector");
 
            break;
 
        case (b2Task::compute_gradients):  // post-processing
        {
            auto container = solution_->get_gradient_container();
            GradientContainer* gradients = container.get();
 
            IterateOverElements(
                  elements,
                  [this, &elements, gradients](const size_t i) {
                      ComputeGradients(
                            dynamic_cast<TypedElement<T>&>(elements[i].get()), gradients);
                  },
                  obtain_timer("compute_gradients"), "compute the gradients");
 
            break;
        }
 
        case (b2Task::compute_J_integral): {
            // logger << logging::info << "Starting to compute the J integral
            // for the solver"
            //        << logging::LOGFLUSH;
 
            auto time_comp_j_integral = obtain_timer("solver_comp_j_integral");
            break;
        }
 
        case (b2Task::update_internal_variables): {
            // logger << logging::info << "Starting to update the internal
            // variables for the solver"
            //        << logging::LOGFLUSH;
 
            auto time_update_inter_var = obtain_timer("solver_update_inter_var");
            break;
        }
 
        case (b2Task::compute_error):
            ComputeError(elements, obtain_timer("solver_comp_err"), "update the error");
            break;
 
        default:
            // Hier Error oder Exception;
            break;
    }
}
 
template <typename T, typename MATRIX_FORMAT>
void TypedSolver<T, MATRIX_FORMAT>::AllocateAndAssembleEffectiveSystem(
      std::vector<std::reference_wrapper<Element>>& elements) {
    K_eff_.resize(effective_num_dof_, model_->get_domain().get_dof_connectivity_type(), *case_);
 
#ifdef HAVE_LIB_TBB
    using my_mutex_t = tbb::spin_mutex;
    std::vector<my_mutex_t> row_blocker(effective_num_dof_);
#endif
    std::vector<std::map<size_t, T>> row_contributions(effective_num_dof_);
 
    logging::Logger& logger = logging::get_logger("solver.monolithic");
    auto timer = obtain_timer("solver_allocate_eff_sys");
 
    auto ElementLoop = [&](size_t i) {
        auto& element{dynamic_cast<TypedElement<T>&>(elements[i].get())};
 
        b2linalg::Matrix<T, typename MATRIX_FORMAT::dense> k_eff;
        b2linalg::Vector<T, b2linalg::Vdense> f_eff;
        b2linalg::Index index;  
 
        element.AssembleElementEffectiveSystem(*this, k_eff, f_eff, index);
 
        f_eff = -f_eff;
 
        size_t N{index.size()};
        bool swapped{};
        int global_index_i{};
        int global_index_j{};
        for (size_t i{}; i < N; ++i) {
            global_index_i = fixed_dof_[index[i]];
 
            if (global_index_i < 0) { continue; }
 
            // size_t j{};
            // if constexpr (std::is_same_v<MATRIX_FORMAT, b2linalg::Msymmetric>) {
            //     j = i;  //!< Fill only "half" of the matrix in case of symmetry.
            // }
            for (size_t j{}; j < N; ++j) {
                global_index_j = fixed_dof_[index[j]];
 
                if (global_index_j < 0) {  
                    if (k_eff(i, j) != T{}) { f_eff[i] -= k_eff(i, j) * dU_(index[j], 0); }
                } else {  
                    if constexpr (std::is_same_v<MATRIX_FORMAT, b2linalg::Msymmetric>) {
                        if (j < i) { continue; }  
 
                        if (global_index_j < global_index_i) {
                            std::swap(global_index_i, global_index_j);
                            swapped = true;
                        }
                    }
 
#ifdef HAVE_LIB_TBB
                    my_mutex_t::scoped_lock row_lock{row_blocker[global_index_i]};
#endif
                    row_contributions[global_index_i][global_index_j] += k_eff(i, j);
 
                    if constexpr (std::is_same_v<MATRIX_FORMAT, b2linalg::Msymmetric>) {
                        if (swapped) {
                            std::swap(global_index_i, global_index_j);
                            swapped = false;
                        }
                    }
                }
            }
 
#ifdef HAVE_LIB_TBB
            std::atomic_ref<T> atomic_T{F_eff_[global_index_i]};
            atomic_T += f_eff[i];
#else
            F_eff_[global_index_i] += f_eff[i];
#endif
        }
    };
 
    logger << logging::info << "Starting to allocate the system matrix" << logging::LOGFLUSH;
    (*timer).second.start();
#ifdef HAVE_LIB_TBB
    tbb::parallel_for(static_cast<size_t>(0), elements.size(), ElementLoop);
#else
    for (size_t i{}; i != elements.size(); ++i) { ElementLoop(i); }
#endif
    (*timer).second.stop();
    logger << logging::info << "Finished to allocate the system matrix" << logging::LOGFLUSH;
 
    timer = obtain_timer("solver_assemble_eff_sys");
 
    auto RowLoop = [&](size_t row) { K_eff_.InitializeRow(row, row_contributions[row]); };
 
    logger << logging::info << "Starting to assemble the system matrix" << logging::LOGFLUSH;
    (*timer).second.start();
#ifdef HAVE_LIB_TBB
    tbb::parallel_for(static_cast<size_t>(0), effective_num_dof_, RowLoop);
#else
    for (size_t row{}; row < effective_num_dof_; ++row) { RowLoop(row); }
#endif
    (*timer).second.stop();
    logger << logging::info << "Finished to assemble the system matrix" << logging::LOGFLUSH;
}
 
template <typename T, typename MATRIX_FORMAT>
void TypedSolver<T, MATRIX_FORMAT>::AssembleEffectiveSystem(TypedElement<T>& element) {
    // domain.set_dof(dof);
 
    b2linalg::Matrix<T, typename MATRIX_FORMAT::dense> k_eff;
    b2linalg::Vector<T, b2linalg::Vdense> f_eff;
    b2linalg::Index index;  
 
    element.AssembleElementEffectiveSystem(*this, k_eff, f_eff, index);
 
    f_eff = -f_eff;
 
    size_t N{index.size()};
    int global_index_i{};
    int global_index_j{};
    for (size_t i{}; i != N; ++i) {
        global_index_i = fixed_dof_[index[i]];
        if (global_index_i < 0) {
            continue;  
        }
 
        for (size_t j{}; j < N; ++j) {
            global_index_j = fixed_dof_[index[j]];
 
            if (global_index_j < 0) {  
                if (k_eff(i, j) != T{}) { f_eff[i] -= k_eff(i, j) * dU_(index[j], 0); }
            } else {  
                if constexpr (std::is_same_v<MATRIX_FORMAT, b2linalg::Msymmetric>) {
                    if (j < i) { continue; }  
                }
 
#ifdef HAVE_LIB_TBB
                std::atomic_ref<T> atomic_U{K_eff_(global_index_i, global_index_j)};
                atomic_U += k_eff(i, j);
#else
                K_eff_(global_index_i, global_index_j) += k_eff(i, j);
#endif
            }
        }
 
#ifdef HAVE_LIB_TBB
        std::atomic_ref<T> atomic_T{F_eff_[global_index_i]};
        atomic_T += f_eff[i];
#else
        F_eff_[global_index_i] += f_eff[i];
#endif
    }
}
 
template <typename T, typename MATRIX_FORMAT>
void TypedSolver<T, MATRIX_FORMAT>::AssembleEffectiveMatrix(TypedElement<T>& element) {
    // domain.set_dof(dof);
    // b2linalg::Index index;
 
    b2linalg::Matrix<T, typename MATRIX_FORMAT::dense>
          k_eff;            
    b2linalg::Index index;  
 
    k_eff = element.AssembleElementEffectiveMatrix(*this, index);
 
    // K_eff_ += StructuredMatrix(total_num_dof_, index, k_eff);  //!< Atomic
    // operation.
}
 
template <typename T, typename MATRIX_FORMAT>
void TypedSolver<T, MATRIX_FORMAT>::AssembleEffectiveVector(TypedElement<T>& element) {
    // domain.set_dof(dof);
    // b2linalg::Index index;
 
    // Element array to calculate and to store into global array
    b2linalg::Vector<T, b2linalg::Vdense> f_eff;  // or Vdense_ref ??? (TB)
    b2linalg::Index index;  
 
    f_eff = element.AssembleElementEffectiveVector(*this, index);
 
    // This method requires some form of atomic operation for TBB, as this will
    // be a critical section.
}
 
template <typename T, typename MATRIX_FORMAT>
inline void TypedSolver<T, MATRIX_FORMAT>::ComputeGradients(
      TypedElement<T>& element, GradientContainer* const gradient_container) {
    gradient_container->set_current_element(element);
 
    element.ComputeElementGradient(*this, gradient_container);
}
 
template <typename T, typename MATRIX_FORMAT>
inline void TypedSolver<T, MATRIX_FORMAT>::ComputeError(
      std::vector<std::reference_wrapper<Element>>& elements, auto timer, std::string_view log) {
#ifdef HAVE_LIB_TBB
    IterateOverElementsParallelReduction(
          elements,
          [this, &elements](const tbb::blocked_range<size_t>& r, T init) -> T {
              for (size_t i{r.begin()}; i != r.end(); ++i) {
                  init = std::max(
                        init, (dynamic_cast<TypedElement<T>&>(elements[i].get()))
                                    .ComputeElementError(*this));
              }
 
              return init;
          },
          [](T x, T y) { return std::max(x, y); },  
          T{},                                      
          timer, log);
#else
    IterateOverElements(
          elements,
          [this, &elements](const size_t i) {
              reduced_scalar_value = std::max(
                    reduced_scalar_value,
                    dynamic_cast<TypedElement<T>&>(elements[i].get()).ComputeElementError(*this));
          },
          timer, log);
#endif
}
 
template <typename T, typename MATRIX_FORMAT>
inline void TypedSolver<T, MATRIX_FORMAT>::SetSolutionFieldsToZero() {
    auto SetdUToZero{[&](size_t i) { dU_(i, 0) = 0; }};
 
#ifdef HAVE_LIB_TBB
    tbb::parallel_for(static_cast<size_t>(0), dof_.size1(), SetdUToZero);
#else
    for (size_t i{}; i < dof_.size1(); ++i) { SetdUToZero(i); }
#endif
 
    K_eff_.set_zero_same_struct();
    F_eff_.set_zero();
}
 
template <typename T, typename MATRIX_FORMAT>
void TypedSolver<T, MATRIX_FORMAT>::SolveEffectiveSystem() {
    b2linalg::Matrix<T, b2linalg::Mrectangle> nks_r;
    b2linalg::SingularMatrixError ee;
    bool is_singular = false;
 
    const int nbnsv = case_->get_int("COMPUTE_NULL_SPACE_VECTOR", 0);
 
    K_eff_.remove_zero();  
 
    try {
        logging::Logger& logger = logging::get_logger("solver.monolithic");
        logger << logging::info << "Starting to solve the effective system" << logging::LOGFLUSH;
 
        F_eff_ = inverse(K_eff_, nks_r, nbnsv) * F_eff_;
 
        logger << logging::info << "Finished to solve the effective system" << logging::LOGFLUSH;
    } catch (b2linalg::SingularMatrixError& e) {
        ee = e;
        is_singular = true;
 
        if (!e.singular_dofs.empty()) {
            if (nks_r.size2() == 0) {
                nks_r.resize(effective_num_dof_, 1);
                nks_r.set_zero();
            }
            for (auto i : e.singular_dofs) {
                assert(i < nks_r.size1());
                for (size_t j{}; j != nks_r.size2(); ++j) { nks_r(i, j) = T(1); }
            }
        }
    }
}
 
template <typename T, typename MATRIX_FORMAT>
typename TypedSolver<T, MATRIX_FORMAT>::type_t TypedSolver<T, MATRIX_FORMAT>::type(
      "TypedSolver", type_name<T, MATRIX_FORMAT>(), StringList("TypedSolver"),
      b2000::solver::module, &Solver::type);
 
}  // namespace b2000::solver
 
#endif  // B2TYPED_SOLVER_H_