b2monolithic_solver.H

//------------------------------------------------------------------------
// b2monolithic_solver.H --
//
//
// written by Thomas Blome <thomas.blome@dlr.de>
//            Neda Ebrahimi Pour <neda.ebrahimipour@dlr.de>
//
// (c) 2004-2012 SMR Engineering & Development SA
//               2502 Bienne, Switzerland
//
// (c) 2024 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 B2MONOLITHIC_SOLVER_H_
#define B2MONOLITHIC_SOLVER_H_
 
#include "b2typed_solver.H"
#include "time_integration/b2newmark.H"
#include "utils/linear_algebra/b2linear_algebra_utils.H"
 
namespace b2000::solver {
 
template <typename T, typename MATRIX_FORMAT>
class MonolithicSolver final : public TypedSolver<T, MATRIX_FORMAT> {
public:
    using type_t =
          ObjectTypeComplete<MonolithicSolver, typename TypedSolver<T, MATRIX_FORMAT>::type_t>;
 
    void solve() override;
 
    static type_t type;
 
protected:
    void InitializeSolverOnCase() override;
 
    bool solve_iteration() override;
};
 
template <typename T, typename MATRIX_FORMAT>
void MonolithicSolver<T, MATRIX_FORMAT>::InitializeSolverOnCase() {
    TypedSolver<T, MATRIX_FORMAT>::InitializeSolverOnCase();
 
    if (this->is_dynamic_ || this->is_damped_) {
        std::string domain_state_id;  
        dynamic_cast<TypedInitialCondition<T>&>(this->model_->get_initial_condition())
              .get_dynamic_initial_condition_value(
                    this->dof_, this->time_, this->current_stage_, domain_state_id);
 
        int n = this->is_dynamic_ ? 3 : 2;
        this->dof_.resize(this->dof_.size1(), n);
        this->case_->set_stage(this->current_stage_);  
 
        // if (!domain_state_id.empty()) { domain.load_state(domain_state_id); }   //! Wofür??
        // Löschen?
 
        if (this->is_linear_) {  
            this->time_integrator_.reset(NewmarkIntegrator<T>::type.new_object());
            this->time_integrator_->InitializeTimeIntegrator(this->dof_, this->step_size_);
        } else {  
            // time_integrator_.reset(NordsieckIntegrator<T, MATRIX_FORMAT>::type.new_object());
            // time_integrator_->InitializeTimeIntegrator(...)
        }
    }
 
    // this->essential_boundary_conditions_->UpdateEssentialBoundaryConditions(
    //       this->time_, this->fixed_dof_, this->dof_, this->dU_);
 
    // this->natural_boundary_conditions_->UpdateNaturalBoundaryConditions(
    //       this->time_, this->fixed_dof_, this->F_eff_);
}
 
template <typename T, typename MATRIX_FORMAT>
void MonolithicSolver<T, MATRIX_FORMAT>::solve() {
    CaseList& case_list = this->model_->get_case_list();
    this->case_iterator_ = case_list.get_case_iterator();
 
    while ((this->case_ = this->case_iterator_->next())) {
        // solve_on_case(*case_);
 
        // if (case_ == nullptr) { solve_init(); }
 
        logging::Logger& logger = logging::get_logger("solver.monolithic");
        logger << logging::info << "Starting the monolithic solver for case "
               << this->case_->get_id() << logging::LOGFLUSH;
        logger << logging::info << "Stage is: " << this->case_->get_stage_id() << logging::LOGFLUSH;
 
        InitializeSolverOnCase();
 
        while (this->time_ < this->time_max_) {
            if (this->is_dynamic_ || this->is_damped_) { this->solution_->new_cycle(); }
 
            this->iteration_count_ = 0;
 
            this->SetSolutionFieldsToZero();
 
            this->essential_boundary_conditions_->UpdateEssentialBoundaryConditions(
                  this->time_, this->fixed_dof_, this->dof_, this->dU_);
 
            this->natural_boundary_conditions_->UpdateNaturalBoundaryConditions(
                  this->time_, this->fixed_dof_, this->F_eff_);
 
            this->time_ += this->step_size_;
 
            logger << logging::info << "Solve for time " << this->time_ << logging::LOGFLUSH;
 
            if (this->is_damped_ || this->is_dynamic_) {
                this->time_integrator_->InitializeFieldsForThisTimeStep(this->dof_, this->dU_);
            }
 
            // clang-format off
            while (solve_iteration());
            // clang-format on
 
            if (!this->is_dynamic_ && !this->is_damped_) {
                this->solution_->set_dof(this->dof_[0]);
            } else {
                this->solution_->set_solution(
                      this->dof_, this->time_, b2linalg::Vector<T, b2linalg::Vdense>::null,
                      b2linalg::Vector<T, b2linalg::Vdense>::null);
            }
 
            if (this->case_->get_int("GRADIENTS", 0)) {
                this->MainLoopOverElements(
                      TypedSolver<T, MATRIX_FORMAT>::b2Task::compute_gradients);
            }
 
            this->solution_->set_natural_boundary_condition(
                  this->natural_boundary_conditions_->GetPrescribedNodalForces());
 
            this->solution_->terminate_cycle(this->time_, this->step_size_, RTable{});
        }
 
        this->solution_->terminate_stage();
        this->solution_->terminate_case(true, RTable{});
 
        logger << logging::info << "Finished the monolithic solver for case "
               << this->case_->get_id() << logging::LOGFLUSH;
    }
}
 
template <typename T, typename MATRIX_FORMAT>
bool MonolithicSolver<T, MATRIX_FORMAT>::solve_iteration() {
    bool result{};  
 
 
    this->MainLoopOverElements(TypedSolver<T, MATRIX_FORMAT>::b2Task::assemble_effective_system);
 
    try {
        // MainLoopOverElements(b2Task::compute_error);
        this->SolveEffectiveSystem();
    } catch (...) {
        // SolutionMonolithic<T>& solution =
        //      dynamic_cast<SolutionMonolithic<T>&>(model_->get_solution());
        // solution.terminate_stage(false);
        // RTable attributes;
        // solution.terminate_case(false, attributes);
        throw;
    }
 
    auto UpdateSolutionfields{[&](size_t i) {
        if (this->fixed_dof_[i] >= 0) {  
            this->dU_(i, 1) = this->F_eff_[this->fixed_dof_[i]];
            this->dU_(i, 0) += this->dU_(i, 1);
            this->dof_(i, 0) += this->dU_(i, 1);
        } else if (this->iteration_count_ == 0) {  
            this->dof_(i, 0) += this->dU_(i, 0);
            this->dU_(i, 0) = 0;
        }
    }};
 
#ifdef HAVE_LIB_TBB
    tbb::parallel_for(static_cast<size_t>(0), this->dof_.size1(), UpdateSolutionfields);
#else
    for (size_t i{}; i < this->dof_.size1(); ++i) { UpdateSolutionfields(i); }
#endif
 
    if (this->is_damped_ || this->is_dynamic_) {
        this->time_integrator_->UpdateFields(this->dof_, this->dU_);
    }
 
    // if (!result) {  //! < Wofür ist das???
    //     RTable attributes;
    //     // dynamic_cast<SolutionMonolithic<T>&>(model_->get_solution())
    //     //      .terminate_case(true, attributes);
    //     case_ = nullptr;
    //     solve_init();
    //     if (case_ == nullptr) { return false; }
    // }
 
    if (this->is_linear_) { return false; }
 
    this->iteration_count_++;
 
    return true;
}
 
template <typename T, typename MATRIX_FORMAT>
typename MonolithicSolver<T, MATRIX_FORMAT>::type_t MonolithicSolver<T, MATRIX_FORMAT>::type(
      "MonolithicSolver", type_name<T, MATRIX_FORMAT>(),
      StringList(
            "MLINEAR", "MSTATIC_LINEAR_SOLUTION", "MDYNAMIC_LINEAR", "MTRANSIENT_LINEAR",
            "MNONLINEAR", "MSTATIC_NONLINEAR_SOLUTION", "MDYNAMIC_NONLINEAR",
            "MTRANSIENT_NONLINEAR"),  // Add also "MCOLLAPSE"? (TB)
      b2000::solver::module, &TypedSolver<T, MATRIX_FORMAT>::type);
 
}  // namespace b2000::solver
 
#endif  // B2MONOLITHIC_SOLVER_H_