b2static_nonlinear_solver.H

//------------------------------------------------------------------------
// b2static_nonlinear_solver.H --
//
//
// written by Mathias Doreille
//            Neda Ebrahimi Pour <neda.ebrahimipour@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 B2STATIC_NONLINEAR_SOLVER_H_
#define B2STATIC_NONLINEAR_SOLVER_H_
 
#include <memory>
 
#include "b2ppconfig.h"
#include "model/b2boundary_condition.H"
#include "model/b2domain.H"
#include "model/b2initial_condition.H"
#include "model/b2model.H"
#include "model/b2solution.H"
#include "model/b2solver.H"
#include "solvers/b2static_nonlinear_utile.H"
#include "utils/b2logging.H"
#include "utils/b2object.H"
 
// This is part of the Monolithic solver and should be removed here:
namespace {
const int MAX_STEPS_DEFAULT = 99999;
}
 
namespace b2000::solver {
 
template <typename T, typename MATRIX_FORMAT>
class StaticNonLinearSolver : public Solver {
public:
    using mrbc_t = TypedModelReductionBoundaryCondition<T>;
    using type_t = ObjectTypeComplete<StaticNonLinearSolver, Solver::type_t>;
 
    void solve() override {
        while (solve_iteration()) { ; }
    }
 
    bool solve_iteration() override;
 
    static type_t type;
 
protected:
    void solve_init();
 
    Domain* domain{};
    SolutionStaticNonlinear<T>* solution{};
 
    logging::Logger& logger{logging::get_logger("solver.static_nonlinear")};
 
    std::unique_ptr<StaticResidueFunction<T, MATRIX_FORMAT>> residue_function{};
    std::unique_ptr<IncrementConstraintStrategy<T, MATRIX_FORMAT>> constraint{};
    std::unique_ptr<CorrectionStrategy<T, MATRIX_FORMAT>> correction{};
 
    int increment_iteration_number{};
    int first_increment_iteration_number_of_stage{};
    double time{};
    // Contains dof solution. Size: ndofs
    b2linalg::Vector<T, b2linalg::Vdense> d_;
    // Contains reduced dof solution
    b2linalg::Vector<T, b2linalg::Vdense> d_red_;
    b2linalg::Matrix<T, typename MATRIX_FORMAT::sparse_inversible_ext> K;
    b2linalg::Vector<T, b2linalg::Vdense> q;
 
    int max_number_of_iteration_by_stage{};
};
 
template <typename T, typename MATRIX_FORMAT>
typename StaticNonLinearSolver<T, MATRIX_FORMAT>::type_t
      StaticNonLinearSolver<T, MATRIX_FORMAT>::type(
            "StaticNonLinearSolver", type_name<T, MATRIX_FORMAT>(),
            StringList("NONLINEAR", "COLLAPSE"), module, &Solver::type);
 
template <typename T, typename MATRIX_FORMAT>
void b2000::solver::StaticNonLinearSolver<T, MATRIX_FORMAT>::solve_init() {
    if (model_ == nullptr) { Exception() << THROW; }
 
    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 static nonlinear solver must have a case that contains at least one "
                       "stage but the case "
                    << case_->get_id() << " does not have a stage." << THROW;
    }
 
    logger << logging::info << "Start static nonlinear solver for case " << case_->get_id()
           << logging::LOGFLUSH;
 
    model_->set_case(*case_);
    max_number_of_iteration_by_stage = case_->get_int("MAX_STEPS", MAX_STEPS_DEFAULT);
    domain = &model_->get_domain();
    solution = dynamic_cast<SolutionStaticNonlinear<T>*>(&model_->get_solution());
    if (solution == 0) {
        TypeError() << "The solver only accepts solutions of type "
                    << typeid(SolutionStaticNonlinear<T>) << THROW;
    }
 
    std::string control_type_name = case_->get_string("INCREMENT_CONTROL_TYPE", "LOAD");
    if (control_type_name != "") { control_type_name += "_CONTROL_CONSTRAINT_STRATEGY"; }
    control_type_name = type_name<T, MATRIX_FORMAT>(control_type_name);
    typename IncrementConstraintStrategy<T, MATRIX_FORMAT>::type_t* control_type =
          IncrementConstraintStrategy<T, MATRIX_FORMAT>::type.get_subtype(
                control_type_name, solver::module);
    constraint.reset(control_type->new_object());
    constraint->set_max_time(case_->get_stage_size());
    constraint->init(*model_, *case_);
 
    std::string residue_function_name = case_->get_string("RESIDUE_FUNCTION_TYPE", "DEFAULT");
    if (residue_function_name != "") { residue_function_name += "_STATIC_RESIDUE_FUNCTION"; }
    residue_function_name = type_name<T, MATRIX_FORMAT>(residue_function_name);
    typename StaticResidueFunction<T, MATRIX_FORMAT>::type_t* residue_function_type =
          StaticResidueFunction<T, MATRIX_FORMAT>::type.get_subtype(
                residue_function_name, solver::module);
    residue_function.reset(residue_function_type->new_object());
    residue_function->init(*model_, *case_);
 
    std::string correction_type_name = case_->get_string("CORRECTION_TYPE", "NEWTON");
    if (!correction_type_name.empty()) { correction_type_name += "_METHOD"; }
    correction_type_name = type_name<T, MATRIX_FORMAT>(correction_type_name);
    typename CorrectionStrategy<T, MATRIX_FORMAT>::type_t* correction_type =
          CorrectionStrategy<T, MATRIX_FORMAT>::type.get_subtype(
                correction_type_name, solver::module);
    correction.reset(correction_type->new_object());
    correction->init(*model_, *case_);
 
    increment_iteration_number = 0;
    first_increment_iteration_number_of_stage = 0;
    size_t s = residue_function->get_number_of_dof();
 
    d_red_.resize(s);
    d_red_.set_zero();
    K.set_zero();
    K.resize(s, 1, domain->get_dof_connectivity_type(), *case_);
    q.resize(s);
    {
        std::string domain_state_id;
        int stage;
 
        b2linalg::Vector<T, b2linalg::Vdense> dof_init(domain->get_number_of_dof());
        dynamic_cast<TypedInitialCondition<T>&>(model_->get_initial_condition())
              .get_static_initial_condition_value(dof_init, time, stage, domain_state_id);
 
        dynamic_cast<mrbc_t&>(
              model_->get_model_reduction_boundary_condition(mrbc_t::type.get_subtype(
                    "TypedModelReductionBoundaryConditionListComponent" + type_name<T>())))
              .get_nonlinear_inverse_value(
                    b2linalg::Matrix<T, b2linalg::Mrectangle_constref>(dof_init), time,
                    b2linalg::Matrix<T, b2linalg::Mrectangle_ref>(d_red_));
        case_->set_stage(stage);
        if (!domain_state_id.empty()) { domain->load_state(domain_state_id); }
    }
}
 
template <typename T, typename MATRIX_FORMAT>
bool b2000::solver::StaticNonLinearSolver<T, MATRIX_FORMAT>::solve_iteration() {
    if (case_ == nullptr) {
        solve_init();
        if (case_ == nullptr) { return false; }
    }
 
    ++increment_iteration_number;
    std::string domain_state_id;
    typename IncrementConstraintStrategy<T, MATRIX_FORMAT>::Result result;
    const double old_time = time;
    try {
        if (!case_->get_bool("GRADIENTS_ONLY", false)) {
            EquilibriumSolution equilibrium_solution(0, 1);
            logger << logging::info << "Start increment " << increment_iteration_number
                   << " of stage " << case_->get_stage_id() << " for time " << time
                   << " and time step " << constraint->get_delta_time() << "." << logging::LOGFLUSH;
            while (1) {
                // prediction step
                logger << logging::debug << "Start prediction of increment "
                       << increment_iteration_number << " of stage " << case_->get_stage_id()
                       << " for time " << time << "." << logging::LOGFLUSH;
                equilibrium_solution.distance_to_iterative_convergence = -1;
                result = constraint->set_increment(
                      *residue_function, K, q, d_red_, time, equilibrium_solution);
 
                if (result != IncrementConstraintStrategy<T, MATRIX_FORMAT>::in_stage) { break; }
 
                // correction step
                logger << logging::debug << "Start correction of increment "
                       << increment_iteration_number << " of stage " << case_->get_stage_id()
                       << " for time " << time << "." << logging::LOGFLUSH;
 
                try {
                    result = correction->compute_correction(
                          *residue_function, K, q, *constraint, d_red_, time, equilibrium_solution);
                } catch (b2linalg::SingularMatrixError& e) {
                    result = constraint->set_correction_result(d_red_, time, false, 0);
                }
 
                if (result == IncrementConstraintStrategy<T, MATRIX_FORMAT>::in_stage_and_retry) {
                    logger << logging::debug
                           << "Correction has diverged, retrying the prediction-correction step "
                              "with a different predictor."
                           << logging::LOGFLUSH;
                    equilibrium_solution.input_level = 0;
                } else {
                    break;
                }
            }
 
            if (result == IncrementConstraintStrategy<T, MATRIX_FORMAT>::end_of_stage_with_error) {
                solution->terminate_stage(false);
                logger << logging::error << "Static nonlinear solver terminates with error(s)."
                       << logging::LOGFLUSH;
                Exception() << THROW;
            }
 
            // save the solution for the current stage and increment
            d_.resize(domain->get_number_of_dof());
            residue_function->get_solution(d_red_, time, true, d_);
            solution->new_cycle(
                  result
                  == IncrementConstraintStrategy<T, MATRIX_FORMAT>::end_of_stage_with_success);
            solution->set_dof(
                  time, d_, residue_function->get_nbc_sol(),
                  residue_function->get_residuum_sol(d_red_, time));
            domain_state_id = solution->get_domain_state_id();
        } else {
            solution->new_cycle(true);
            d_.resize(domain->get_number_of_dof());
            solution->get_dof(time, d_);
            if (std::abs(time - case_->get_stage_size()) < 1e-8) {
                result = IncrementConstraintStrategy<T, MATRIX_FORMAT>::end_of_stage_with_success;
            } else {
                result = IncrementConstraintStrategy<T, MATRIX_FORMAT>::in_stage;
            }
        }
 
        {
            typename SolutionStaticNonlinear<T>::gradient_container gc =
                  solution->get_gradient_container();
            domain->set_dof(d_);
            EquilibriumSolution es(1, 0);
            residue_function->get_residue(
                  d_red_, time, time - old_time, es,
                  b2linalg::Vector<T, b2linalg::Vdense_ref>::null,
                  b2linalg::Matrix<T, typename MATRIX_FORMAT::sparse_inversible_ext>::null,
                  b2linalg::Vector<T, b2linalg::Vdense_ref>::null, 0, 0, gc.get());
        }
        if (!domain_state_id.empty()) { domain->save_state(domain_state_id); }
        RTable attributes;
        solution->terminate_cycle(time, time - old_time, attributes);
    } catch (...) {
        solution->terminate_stage(false);
        RTable attributes;
        solution->terminate_case(false, attributes);
        throw;
    }
 
    if (result == IncrementConstraintStrategy<T, MATRIX_FORMAT>::end_of_stage_with_success) {
        solution->terminate_stage(true);
        if (case_->next_stage()) {
            model_->set_case(*case_);
            time = 0;
            constraint->set_max_time(case_->get_stage_size());
            constraint->init(*model_, *case_);
            residue_function->init(*model_, *case_);
            correction->init(*model_, *case_);
            first_increment_iteration_number_of_stage = increment_iteration_number;
            max_number_of_iteration_by_stage = case_->get_int("MAX_STEPS", MAX_STEPS_DEFAULT);
            return true;
        }
        RTable attributes;
        solution->terminate_case(true, attributes);
        logger << logging::info << "End static nonlinear analysis for case " << case_->get_id()
               << " with success." << logging::LOGFLUSH;
        case_->warn_on_non_used_key();
        case_ = nullptr;
        return true;
    }
 
    if (increment_iteration_number - first_increment_iteration_number_of_stage
        > max_number_of_iteration_by_stage) {
        solution->terminate_stage(false);
        RTable attributes;
        solution->terminate_case(false, attributes);
        logger << logging::error
               << "Static nonlinear solver terminates with error(s). The number "
                  "of iterations for stage "
               << case_->get_stage_id() << " exceeds the maximum specified." << logging::LOGFLUSH;
        Exception() << THROW;
    }
 
    return true;
}
 
}  // namespace b2000::solver
 
#endif  // B2STATIC_NONLINEAR_SOLVER_H_