b2dynamic_nonlinear_solver.H

//------------------------------------------------------------------------
// b2dynamic_nonlinear_solver.H --
//
//
// written by Mathias Doreille
//            Neda Ebrahimi Pour <neda.ebrahimipour@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 B2DYNAMIC_NONLINEAR_SOLVER_H_
#define B2DYNAMIC_NONLINEAR_SOLVER_H_
 
#include "b2dynamic_nonlinear_utile.H"
#include "b2ppconfig.h"
#include "model/b2boundary_condition.H"
#include "model/b2domain.H"
#include "model/b2initial_condition.H"
#include "model/b2solution.H"
#include "model/b2solver.H"
#include "utils/b2logging.H"
#include "utils/b2object.H"
 
namespace b2000::solver {
 
template <typename T, typename MATRIX_FORMAT>
class DynamicNonLinearSolver : public Solver {
public:
    using mrbc_t = TypedModelReductionBoundaryCondition<T>;
    using type_t = ObjectTypeComplete<DynamicNonLinearSolver, Solver::type_t>;
 
    DynamicNonLinearSolver()
        : residue_function(0),
          implicit_parameter(nullptr),
          explicit_parameter(nullptr),
          bdf_solver(0) {}
 
    ~DynamicNonLinearSolver() override {
        if (residue_function) { residue_function->decref_or_free(); }
        if (implicit_parameter) { implicit_parameter->decref_or_free(); }
        if (explicit_parameter) { explicit_parameter->decref_or_free(); }
        if (bdf_solver) { bdf_solver->decref_or_free(); }
    }
 
    void solve() override {
        while (solve_iteration()) { ; }
    }
 
    bool solve_iteration() override {
        if (case_ == nullptr) { solve_init(); }
        bool res;
        try {
            res = bdf_solver->next();
        } catch (...) {
            SolutionDynamicNonlinear<T>& solution =
                  dynamic_cast<SolutionDynamicNonlinear<T>&>(model_->get_solution());
            solution.terminate_stage(false);
            RTable attributes;
            solution.terminate_case(false, attributes);
            throw;
        }
        if (!res) {
            RTable attributes;
            dynamic_cast<SolutionDynamicNonlinear<T>&>(model_->get_solution())
                  .terminate_case(true, attributes);
            case_ = nullptr;
            solve_init();
            if (case_ == nullptr) { return false; }
        }
        return true;
    }
 
    static type_t type;
 
protected:
    void 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 dynamic nonlinear solver must have a case that contains at least "
                           "one stage but the case "
                        << case_->get_id() << " does not have a stage." << THROW;
        }
 
        model_->set_case(*case_);
 
        if (residue_function) { residue_function->decref_or_free(); }
        if (implicit_parameter) { implicit_parameter->decref_or_free(); }
        if (explicit_parameter) { explicit_parameter->decref_or_free(); }
        if (bdf_solver) { bdf_solver->decref_or_free(); }
 
        b2linalg::Matrix<T, b2linalg::Mrectangle> dof;
        double time;
        int stage;
        std::string domain_state_id;
 
        dynamic_cast<TypedInitialCondition<T>&>(model_->get_initial_condition())
              .get_dynamic_initial_condition_value(dof, time, stage, domain_state_id);
        case_->set_stage(stage);
        if (!domain_state_id.empty()) { model_->get_domain().load_state(domain_state_id); }
 
        std::string residue_function_name = case_->get_string("RESIDUE_FUNCTION_TYPE", "ORDER_N");
        if (residue_function_name != "") { residue_function_name += "_DYNAMIC_RESIDUE_FUNCTION"; }
        residue_function_name += type_name<T, MATRIX_FORMAT>();
        typename DynamicResidueFunction<T, MATRIX_FORMAT>::type_t* residue_function_type =
              DynamicResidueFunction<T, MATRIX_FORMAT>::type.get_subtype(
                    residue_function_name, solver::module);
        residue_function = residue_function_type->new_object(allocator);
        residue_function->set_case(case_);
        residue_function->init(*model_, *case_);
 
        std::string implicit_parameter_name =
              case_->get_string("IMPLICIT_NORDSIECK_PARAMETER", "STANDARD_BDF");
        if (implicit_parameter_name != "") { implicit_parameter_name += "_PARAMETER"; }
        NordsieckParameter::type_t* implicit_parameter_type =
              NordsieckParameter::type.get_subtype(implicit_parameter_name, solver::module);
        implicit_parameter = implicit_parameter_type->new_object(allocator);
        implicit_parameter->set_attribute(*case_);
 
        std::string explicit_parameter_name =
              case_->get_string("EXPLICIT_NORDSIECK_PARAMETER", "ADAMS_BASHFORTH");
        if (explicit_parameter_name != "") { explicit_parameter_name += "_PARAMETER"; }
        NordsieckParameter::type_t* explicit_parameter_type =
              NordsieckParameter::type.get_subtype(explicit_parameter_name, solver::module);
        explicit_parameter = explicit_parameter_type->new_object(allocator);
        explicit_parameter->set_attribute(*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);
 
        bdf_solver = NordsieckSolver<T, MATRIX_FORMAT>::type.new_object(allocator);
        bdf_solver->init(
              *model_, residue_function, implicit_parameter, explicit_parameter, dof_red, time,
              case_);
    }
 
    CaseList::case_iterator case_iterator;
 
    Allocator allocator;
 
    DynamicResidueFunction<T, MATRIX_FORMAT>* residue_function;
    NordsieckParameter* implicit_parameter;
    NordsieckParameter* explicit_parameter;
    NordsieckSolver<T, MATRIX_FORMAT>* bdf_solver;
};
 
template <typename T, typename MATRIX_FORMAT>
typename DynamicNonLinearSolver<T, MATRIX_FORMAT>::type_t
      DynamicNonLinearSolver<T, MATRIX_FORMAT>::type(
            "DynamicNonLinearSolver", type_name<T, MATRIX_FORMAT>(),
            StringList("DYNAMIC_NONLINEAR", "TRANSIENT_NONLINEAR"), b2000::solver::module,
            &Solver::type);
 
}  // namespace b2000::solver
 
#endif  // B2DYNAMIC_NONLINEAR_SOLVER_H_