b2model_reduction_solver.H

//------------------------------------------------------------------------
// b2model_reduction_solver.H --
//
// written by doreille
//            Thomas Blome <thomas.blome@dlr.de>
//
// Copyright (c) 2010 SMR Engineering & Development SA
//                    2502 Bienne, Switzerland
//
// Copyright (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 B2MODEL_REDUCTION_SOLVER_H_
#define B2MODEL_REDUCTION_SOLVER_H_
 
#include <string.h>
 
#include "b2solver_utils.H"
#include "model/b2boundary_condition.H"
#include "model/b2domain.H"
#include "model/b2element.H"
#include "model/b2model.H"
#include "model/b2node.H"
#include "model/b2solution.H"
#include "model/b2solver.H"
#include "utils/b2linear_algebra.H"
#include "utils/b2logging.H"
#include "utils/b2object.H"
 
namespace b2000::solver {
 
template <typename T, typename MATRIX_FORMAT>
class ModelReductionSolver : public Solver {
public:
    using ebc_t = TypedEssentialBoundaryCondition<T>;
    using mrbc_t = TypedModelReductionBoundaryCondition<T>;
    using type_t = ObjectTypeComplete<ModelReductionSolver, Solver::type_t>;
 
    void solve() override {
        while (solve_iteration()) { ; }
    }
 
    bool solve_iteration() override;
 
    static type_t type;
 
private:
    SolutionStaticNonlinear<T>* solution;
};
 
template <typename T, typename MATRIX_FORMAT>
typename ModelReductionSolver<T, MATRIX_FORMAT>::type_t
      ModelReductionSolver<T, MATRIX_FORMAT>::type(
            "ModelReductionSolver", type_name<T, MATRIX_FORMAT>(), StringList("MODEL_REDUCTION"),
            module, &Solver::type);
 
template <typename T, typename MATRIX_FORMAT>
struct eigen_matrix_type {};
 
template <typename T>
struct eigen_matrix_type<std::complex<T>, b2linalg::Munsymmetric> {
    static const char A = 'N';
    static const char B = 'N';
};
 
template <typename T>
struct eigen_matrix_type<T, b2linalg::Msymmetric> {
    static const char A = 'P';
    static const char B = 'S';
};
 
template <typename T, typename MATRIX_FORMAT>
bool b2000::solver::ModelReductionSolver<T, MATRIX_FORMAT>::solve_iteration() {
    logging::Logger& logger = logging::get_logger("solver.model_reduction");
 
    if (case_iterator_.get() == nullptr) {
        case_iterator_ = model_->get_case_list().get_case_iterator();
    }
 
    case_ = case_iterator_->next();
    if (case_ == nullptr) {
        // all the cases are treated
        return false;
    }
    logger.LOG(logging::info, "Start model reduction solver");
    model_->set_case(*case_);
 
    Domain& domain = model_->get_domain();
    const size_t number_of_dof = domain.get_number_of_dof();
 
    SolutionModelReduction<T, typename MATRIX_FORMAT::dense>* solution =
          &dynamic_cast<SolutionModelReduction<T, typename MATRIX_FORMAT::dense>&>(
                model_->get_solution());
    if (solution == 0) {
        TypeError() << "The solver only accepts solutions of type SolutionModelReduction<T>"
                    << THROW;
    }
 
    const std::string constraint_string = case_->get_string("CONSTRAINT_METHOD", "REDUCTION");
 
    const std::string reduction_method = case_->get_string("REDUCTION_METHOD", "CRAIG_BAMPTON");
 
    // get the dof interface
    b2linalg::Index interface_dof;
    solution->get_interface_dof(interface_dof);
 
    // Assemble the matrix K, D and M
    b2linalg::Matrix<T, typename MATRIX_FORMAT::sparse_inversible> K(number_of_dof);
    b2linalg::Matrix<T, typename MATRIX_FORMAT::sparse_inversible> D(number_of_dof);
    b2linalg::Matrix<T, typename MATRIX_FORMAT::sparse_inversible> M(number_of_dof);
    {
        logger << logging::info << "Element matrix assembly" << logging::LOGFLUSH;
        domain.set_dof(b2linalg::Vector<T, b2linalg::Vdense_constref>::null);
        Domain::element_iterator i = domain.get_element_iterator();
        b2linalg::Index index;
        std::vector<bool> d_value_d_dof_flags(3, true);
        std::vector<b2linalg::Matrix<T, typename MATRIX_FORMAT::dense> > d_value_d_dof(3);
        bool log_el_flag = logger.is_enabled_for(logging::data);
        for (Element* e = i->next(); e != nullptr; e = i->next()) {
            if (auto te = dynamic_cast<TypedElement<T>*>(e); te != nullptr) {
                te->get_value(
                      *model_,
                      true,  // linear
                      true,  // equilibrium_solution
                      0,     // time
                      0,     // delta_time
                      b2linalg::Matrix<T, b2linalg::Mrectangle_constref>::null,
                      0,  // gradient_container
                      0,  // solver_hints
                      index, b2linalg::Vector<T, b2linalg::Vdense>::null, d_value_d_dof_flags,
                      d_value_d_dof, b2linalg::Vector<T, b2linalg::Vdense>::null);
            }
            K += b2linalg::StructuredMatrix(number_of_dof, index, d_value_d_dof[0]);
            D += b2linalg::StructuredMatrix(number_of_dof, index, d_value_d_dof[1]);
            M += b2linalg::StructuredMatrix(number_of_dof, index, d_value_d_dof[2]);
            if (log_el_flag) {
                b2linalg::Matrix<T> tmp_element(d_value_d_dof[0]);
                logger << logging::data << "elementary second variation "
                       << logging::DBName("ELSV", 0, e->get_object_name()) << tmp_element
                       << " of element id " << e->get_object_name() << logging::LOGFLUSH;
                tmp_element = d_value_d_dof[1];
                logger << logging::data << "elementary mass matrix "
                       << logging::DBName("ELSV", 1, e->get_object_name()) << tmp_element
                       << " of element id " << e->get_object_name() << logging::LOGFLUSH;
                tmp_element = d_value_d_dof[2];
                logger << logging::data << "elementary mass matrix "
                       << logging::DBName("ELSV", 2, e->get_object_name()) << tmp_element
                       << " of element id " << e->get_object_name() << logging::LOGFLUSH;
            }
        }
        logger << logging::info << "Element matrix assembly finished" << logging::LOGFLUSH;
    }
 
    // Get the displacement boundary condition (L * x + l = 0)
    ebc_t* dbc =
          &dynamic_cast<ebc_t&>(model_->get_essential_boundary_condition(ebc_t::type.get_subtype(
                "TypedEssentialBoundaryConditionListComponent" + type_name<T>())));
    b2linalg::Vector<T, b2linalg::Vdense> l(dbc->get_size(true));
    b2linalg::Matrix<T, b2linalg::Mcompressed_col> trans_L;
    dbc->get_linear_value(l, trans_L);
 
    // Get the model reduction (x = R * x_r + r)
    mrbc_t* mrbc = &dynamic_cast<mrbc_t&>(
          model_->get_model_reduction_boundary_condition(mrbc_t::type.get_subtype(
                "TypedModelReductionBoundaryConditionListComponent" + type_name<T>())));
    b2linalg::Vector<T, b2linalg::Vdense> r(number_of_dof);
    b2linalg::Matrix<T, b2linalg::Mcompressed_col> trans_R;
    mrbc->get_linear_value(r, trans_R);
 
    if (trans_L.size2() > 0) {
        if (constraint_string == "PENALTY" || constraint_string == "AUGMENTED_LAGRANGE") {
            const double penalty_factor = case_->get_double(
                  "CONSTRAINT_PENALTY", constraint_string == "PENALTY" ? 1e10 : 1);
            K += (trans_L * transposed(trans_L)) * penalty_factor;
        }
        if (constraint_string == "LAGRANGE" || constraint_string == "AUGMENTED_LAGRANGE") {
            UnimplementedError() << THROW;
        }
    }
 
    b2linalg::Matrix<T, b2linalg::Mrectangle> B;
 
    if (reduction_method == "CRAIG_BAMPTON" || reduction_method == "GUYAN") {
        // renumber the reduction to put the interface dof at the beginning
        {
            b2linalg::Index row_to_remove;
            if (!trans_R.get_nonzero_unit_row_for_column(interface_dof, row_to_remove)) {
                Exception() << "Some boundaries conditions are set on the interface dof." << THROW;
            }
            trans_R.remove_row(row_to_remove);
        }
 
        // Compute the interface constrain stiffness
        const size_t K_augm_size =
              trans_R.size1() + (constraint_string == "PENALTY" ? 0 : trans_L.size2());
        b2linalg::Matrix<T, typename MATRIX_FORMAT::sparse_inversible> K_augm(
              K_augm_size, domain.get_dof_connectivity_type(), *case_);
        K_augm += trans_R * K * b2linalg::transposed(trans_R);
        if (constraint_string == "LAGRANGE" || constraint_string == "AUGMENTED_LAGRANGE") {
            K_augm(
                  b2linalg::Interval(trans_R.size1(), K_augm_size),
                  b2linalg::Interval(0, trans_R.size1())) +=
                  b2linalg::transposed(trans_L) * b2linalg::transposed(trans_R);
        }
 
        // The reduction base
        const size_t number_of_mode =
              reduction_method == "GUYAN"
                    ? 0
                    : std::min(
                          size_t(case_->get_int("NMODES", interface_dof.size())), K_augm_size - 1);
        B.resize(number_of_dof, interface_dof.size() + number_of_mode);
 
        // compute the Guyan reduction
        {
            logger << logging::info << "Compute the " << interface_dof.size() << " Guyan modes"
                   << logging::LOGFLUSH;
            b2linalg::Matrix<T, b2linalg::Mrectangle> Mtmp;
            b2linalg::Matrix<T, b2linalg::Mcompressed_col> K1(
                  b2linalg::Matrix<T, b2linalg::Msym_compressed_col_st_constref>(K), interface_dof);
            Mtmp = trans_R * K1;
            Mtmp = inverse(K_augm) * Mtmp;
            b2linalg::Matrix<T, b2linalg::Mrectangle> Mtmp1;
            Mtmp1 = transposed(trans_R) * Mtmp;
            B(b2linalg::Interval(0, interface_dof.size())) = -Mtmp1;
            for (size_t i = 0; i != interface_dof.size(); ++i) { B(interface_dof[i], i) = 1; }
        }
 
        // compute the fixed interface dynamic modes
        if (number_of_mode > 0) {
            logger << logging::info << "Compute the " << number_of_mode
                   << " fixed interface dynamic modes" << logging::LOGFLUSH;
            b2linalg::Matrix<T, typename MATRIX_FORMAT::sparse_inversible> M_augm(
                  K_augm_size, domain.get_dof_connectivity_type(), *case_);
            M_augm += trans_R * M * b2linalg::transposed(trans_R);
 
            T sigma = case_->get_double("SHIFT", 0);
            std::string which = case_->get_string("WHICH_EIGENVALUES", "SM");
 
            typedef eigen_matrix_type<T, MATRIX_FORMAT> eig_matrix_type;
            b2linalg::Vector<T, b2linalg::Vdense> eigvalue(number_of_mode);
            b2linalg::Matrix<T, b2linalg::Mrectangle> eigvector_red = eigenvector(
                  K_augm, eig_matrix_type::A, M_augm, eig_matrix_type::B, eigvalue, which.c_str(),
                  sigma);
 
            b2linalg::Matrix<T, b2linalg::Mrectangle> eigvector =
                  transposed(trans_R)
                  * eigvector_red(
                        b2linalg::Interval(0, trans_R.size1()),
                        b2linalg::Interval(0, eigvector_red.size2()));
            B(b2linalg::Interval(interface_dof.size(), B.size2())) += eigvector;
        }
 
    } else if (reduction_method == "MODE_DISPLACEMENT" || reduction_method == "MODE_ACCELERATION") {
        const bool mode_acc = reduction_method == "MODE_ACCELERATION";
 
        const size_t K_augm_size =
              trans_R.size1() + (constraint_string == "PENALTY" ? 0 : trans_L.size2());
        b2linalg::Matrix<T, typename MATRIX_FORMAT::sparse_inversible> K_augm(
              K_augm_size, domain.get_dof_connectivity_type(), *case_);
        K_augm += trans_R * K * b2linalg::transposed(trans_R);
        if (constraint_string == "LAGRANGE" || constraint_string == "AUGMENTED_LAGRANGE") {
            K_augm(
                  b2linalg::Interval(trans_R.size1(), K_augm_size),
                  b2linalg::Interval(0, trans_R.size1())) +=
                  b2linalg::transposed(trans_L) * b2linalg::transposed(trans_R);
        }
 
        const size_t number_of_mode =
              std::min(size_t(case_->get_int("NMODES", interface_dof.size())), K_augm_size - 1);
        B.resize(number_of_dof, number_of_mode + (mode_acc ? interface_dof.size() : 0));
 
        // compute the mode displacement vector
        {
            b2linalg::Matrix<T, typename MATRIX_FORMAT::sparse_inversible> M_augm(
                  K_augm_size, domain.get_dof_connectivity_type(), *case_);
            M_augm += trans_R * M * b2linalg::transposed(trans_R);
 
            typedef eigen_matrix_type<T, MATRIX_FORMAT> eig_matrix_type;
 
            std::string which_eigenvalues = case_->get_string("WHICH_EIGENVALUES", "SM");
            b2linalg::Vector<T, b2linalg::Vdense> eigvalue(number_of_mode);
            b2linalg::Matrix<T, b2linalg::Mrectangle> eigvector_red = b2linalg::eigenvector(
                  K_augm, eig_matrix_type::A, M_augm, eig_matrix_type::B, eigvalue,
                  which_eigenvalues.c_str(), case_->get_double("SHIFT", 0.0));
 
            B(b2linalg::Interval(0, number_of_mode)) =
                  b2linalg::transposed(trans_R)
                  * eigvector_red(
                        b2linalg::Interval(0, trans_R.size1()),
                        b2linalg::Interval(0, eigvector_red.size2()));
        }
 
        // compute the static correction
        if (mode_acc) {
            b2linalg::Matrix<T, b2linalg::Mrectangle_ref> f =
                  B(b2linalg::Interval(number_of_mode, B.size2()));
            for (size_t i = 0; i != interface_dof.size(); ++i) { f(interface_dof[i], i) = 1; }
            b2linalg::Matrix<T, b2linalg::Mrectangle> f_r(K_augm_size, interface_dof.size());
            f_r(b2linalg::Interval(0, trans_R.size1()), b2linalg::Interval(0, f_r.size2())) =
                  -trans_R * f;
            f_r = b2linalg::inverse(K_augm) * f_r;
            f = b2linalg::transposed(trans_R)
                * f_r(b2linalg::Interval(0, trans_R.size1()), b2linalg::Interval(0, f_r.size2()));
        }
 
        // compute a "better" base vectors of the subspace generated by B
        const size_t rank = b2linalg::Matrix<T, b2linalg::Mrectangle_increment_ref>(B)
                                  .diagonalise_base_vector_on_observation_matrix(interface_dof);
        B.resize(B.size1(), rank);
 
    } else {
        Exception() << "unknown reduction method " << reduction_method << "." << THROW;
    }
 
    // Now B contain the reduction basis.
    // Compute the reduction at the interface dof
    // and store it in the solution object
    {
        logger << logging::info << "Compute the reduced stiffness and mass matrix"
               << logging::LOGFLUSH;
        std::vector<b2linalg::Matrix<T, typename MATRIX_FORMAT::dense> > RedSys(3);
        RedSys[0] = b2linalg::transposed(B) * K * B;
        if (D.get_nb_nonzero()) {
            RedSys[1] = b2linalg::transposed(B) * D * B;
        } else {
            RedSys[1].resize(B.size2(), B.size2());
        }
        RedSys[2] = b2linalg::transposed(B) * M * B;
        solution->set_reduction(interface_dof.size(), B, RedSys);
    }
 
    RTable attributes;
    solution->terminate_case(true, attributes);
 
    logger.LOG(logging::info, "End of model reduction solver");
    case_->warn_on_non_used_key();
    return true;
}
 
}  // namespace b2000::solver
 
#endif /* B2MODEL_REDUCTION_SOLVER_H_ */