b2free_vibration_solver.H

//------------------------------------------------------------------------
// b2free_vibration_solver.H --
//
//
// written by Mathias Doreille
//            Neda Ebrahimi Pour <neda.ebrahimipour@dlr.de>
//            Thomas Blome <thomas.blome@dlr.de>
//
// (c) 2004-2012,2016 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 B2FREE_VIBRATION_SOLVER_H_
#define B2FREE_VIBRATION_SOLVER_H_
 
#include <filesystem>
 
#include "b2ppconfig.h"
#include "model/b2boundary_condition.H"
#include "model/b2domain.H"
#include "model/b2element.H"
#include "model/b2initial_condition.H"
#include "model/b2model.H"
#include "model/b2solution.H"
#include "model/b2solver.H"
#include "utils/b2exception.H"
#include "utils/b2logging.H"
#include "utils/b2object.H"
 
namespace b2000::solver {
 
template <typename T, typename MATRIX_FORMAT>
struct EigenDecompDense {
    void resolve(
          b2linalg::Matrix<T, typename MATRIX_FORMAT::sparse_inversible>& K_augm,
          b2linalg::Matrix<T, typename MATRIX_FORMAT::sparse_inversible>& M_augm,
          const std::string& which, b2linalg::Vector<T, b2linalg::Vdense>& eigenvalue,
          b2linalg::Matrix<T, b2linalg::Mrectangle>& eigenvector) {
        UnimplementedError() << THROW;
    }
};
 
template <typename MATRIX_FORMAT>
struct EigenDecompDense<double, MATRIX_FORMAT> {
    void resolve(
          b2linalg::Matrix<double, typename MATRIX_FORMAT::sparse_inversible>& K_augm,
          b2linalg::Matrix<double, typename MATRIX_FORMAT::sparse_inversible>& M_augm,
          const std::string& which, b2linalg::Vector<double, b2linalg::Vdense>& eigenvalue,
          b2linalg::Matrix<double, b2linalg::Mrectangle>& eigenvector) {
        assert(eigenvalue.size() <= K_augm.size1());
        if (which != "SM" && which != "LM") { UnimplementedError() << THROW; }
 
        b2linalg::Matrix<double, b2linalg::Mrectangle> K(K_augm.size1(), K_augm.size2());
        for (size_t i = 0; i != K.size1(); ++i) {
            for (size_t j = 0; j != K.size2(); ++j) { K(i, j) = K_augm(i, j); }
        }
 
        b2linalg::Matrix<double, b2linalg::Mrectangle> M(M_augm.size1(), M_augm.size2());
        for (size_t i = 0; i != M.size1(); ++i) {
            for (size_t j = 0; j != M.size2(); ++j) { M(i, j) = M_augm(i, j); }
        }
 
        const int n = int(K.size1());
        b2linalg::Vector<double, b2linalg::Vdense> eigvalue_all(n);
        eigvalue_all.set_zero();
 
        std::vector<double> work(2 * (1 + 6 * n + 2 * n * n));
        std::vector<int> iwork(2 * (3 + 5 * n));
        int info = 0;
        lapack::sygvd(
              1, 'V', 'U', n, &K(0, 0), n, &M(0, 0), n, &eigvalue_all[0], &work[0],
              int(work.size()), &iwork[0], int(iwork.size()), info);
 
        assert(info >= 0);
        if (info > 0 && info <= n) {
            Exception() << "Eigendecomposition with lapack routine dsygvd "
                           "failed, "
                        << info
                        << " off-diagonal elements of an "
                           "intermediate tridiagonal form did not converge to zero."
                        << THROW;
        }
        if (info > n) {
            Exception() << "Eigendecomposition with lapack routine dsygvd "
                           "failed, the "
                        << info - n
                        << ". leading minor of the "
                           "("
                        << M.size1() << "x" << M.size2()
                        << ") mass matrix is "
                           "not positive-definite."
                        << THROW;
        }
 
        const size_t o = (which == "SM" ? 0 : eigvalue_all.size() - eigenvalue.size());
        eigenvector.resize(K.size1(), eigenvalue.size());
        for (size_t i = 0; i != eigenvalue.size(); ++i) {
            eigenvalue[i] = eigvalue_all[o + i];
            eigenvector[i] = K[o + i];
        }
    }
};
 
template <typename T, typename MATRIX_FORMAT>
void eigen_decomposition_dense(
      b2linalg::Matrix<T, typename MATRIX_FORMAT::sparse_inversible>& K_augm,
      b2linalg::Matrix<T, typename MATRIX_FORMAT::sparse_inversible>& M_augm,
      const std::string& which, b2linalg::Vector<T, b2linalg::Vdense>& eigenvalue,
      b2linalg::Matrix<T, b2linalg::Mrectangle>& eigenvector) {
    EigenDecompDense<T, MATRIX_FORMAT> ed;
    ed.resolve(K_augm, M_augm, which, eigenvalue, eigenvector);
}
 
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>
class FreeVibrationSolver : public Solver {
public:
    using eig_matrix_type = eigen_matrix_type<T, MATRIX_FORMAT>;
    using ebc_t = TypedEssentialBoundaryCondition<T>;
    using mrbc_t = TypedModelReductionBoundaryCondition<T>;
    using type_t = ObjectTypeComplete<FreeVibrationSolver, Solver::type_t>;
 
    void solve() override {
        while (solve_iteration()) { ; }
    }
 
    bool solve_iteration() override {
        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_) {
            if (case_->get_stage_size() != 1.0) {
                logging::get_logger("solver.free_vibration")
                      << logging::warning << "The stage size of the case " << case_->get_id()
                      << " is equal to " << case_->get_stage_size()
                      << " but the free vibration solver only handle a case with"
                         " one stage of size equal to one."
                         " The specified stage size is ignored."
                      << logging::LOGFLUSH;
            }
            if (case_->get_number_of_stage() != 1) {
                logging::get_logger("solver.free_vibration")
                      << logging::warning << "The number of stage of the case" << case_->get_id()
                      << " is " << case_->get_stage_size()
                      << " but the free vibration solver only handle a case with"
                         " one stage of size equal to one. The extra stages are"
                         " ignored."
                      << logging::LOGFLUSH;
            }
            solve_on_case(*case_);
            return true;
        }
        return false;
    }
 
    static type_t type;
 
protected:
    void solve_on_case(Case& case_);
    // Number of eigen modes
    size_t nb_eigenmodes_;
    // Contains eigen value for each mode. Size: nModes
    b2linalg::Vector<T, b2linalg::Vdense> eigenvalues_;
    // Contains eigen vector for each mode. Size: nDofs, nModes
    b2linalg::Matrix<T, b2linalg::Mrectangle> eigenvectors_;
};
 
template <typename T, typename MATRIX_FORMAT>
typename FreeVibrationSolver<T, MATRIX_FORMAT>::type_t FreeVibrationSolver<T, MATRIX_FORMAT>::type(
      "FreeVibrationSolver", type_name<T, MATRIX_FORMAT>(),
      StringList("FREE_VIBRATION", "VIBRATION"), module, &Solver::type);
 
}  // namespace b2000::solver
 
/* Perform a free vibration analysis on a dynamic linear problem using the
    Lagrange Multiplier Adjunction method. This version accept only
    one essential boundary condition. The non-linear part of the
    essential boundary condition are ignored.
 
Solve the free vibration problem:
  K * x + M * x'' = 0
with the essential boundary conditions:
  x = R * x_r
  L * x = 0
 
Where:
  - K (matrix) is the derivative of the value relatively to the degree of freedom.
  - M (matrix) is the derivative of the value relatively to the acceleration.
  - x (vector) is the degree of freedom.
  - x'' (vector) is the acceleration (the double derivative of the
    degree of freedom relatively to the time).
  - x_r (vector) is the degree of freedom in the reduction model.
  - f (vector) is the natural boundary conditions.
  - R (matrix) define the linear model reduction. r (vector) is ignored.
  - L (matrix) define the linear essential boundary conditions. l
    (vector) is ignored.
 
Using the model reduction equation, the system to solve becomes:
  trans(R) * K * R * x_r + trans(R) * M * R * x_r'' = 0
Where:
  - x_r'' (vector) is the acceleration in the reduced system (the
    double derivative of x_r relatively to the time).
 
with the boundary conditions:
  L * R * x_r = 0
 
Using the Lagrange Multiplier Adjunction method for imposing the
displacement boundary conditions, the generalised eigenproblem of the
free vibration equilibrium are:
 
  | trans(R) * K * R       trans(L * R) |   | x_r |        | trans(R) * M * R      0 |   | x_r |
  |                                     | * |     | =  e * |                         | * |     |
  | L * R                             0 |   | m   |        | 0                     0 |   | m   |
 
The eigenvector x is then obtained from the eigenvector x_r using the model
reduction equation.
 
*/
template <typename T, typename MATRIX_FORMAT>
void b2000::solver::FreeVibrationSolver<T, MATRIX_FORMAT>::solve_on_case(Case& case_) {
    logging::Logger& logger = logging::get_logger("solver.free_vibration");
 
    logger << logging::info << "Start the free vibration solver" << logging::LOGFLUSH;
 
    model_->set_case(case_);
 
    Domain& domain = model_->get_domain();
    SolutionFreeVibration<T>* solution =
          &dynamic_cast<SolutionFreeVibration<T>&>(model_->get_solution());
    if (solution == 0) {
        TypeError() << "The solver only accepts a solution of type"
                    << " SolutionFreeVibration<T>" << THROW;
    }
 
    T sigma;
    T sigma_max;
    char which[3];
    {
        int n;
        solution->which_modes(n, which, sigma, sigma_max);
        assert(n > 0);
        nb_eigenmodes_ = size_t(n);
    }
    const std::string constraint_string = case_.get_string("CONSTRAINT_METHOD", "REDUCTION");
 
    const std::string matrices_dir = case_.get_string("MATRICES_DIR", "");
    const bool export_sys_matrices = (matrices_dir != "");
 
    const size_t number_of_dof = domain.get_number_of_dof();
 
    b2linalg::Vector<T, b2linalg::Vdense> dof;
    double time = 0;
    int stage;
    std::string domain_state_id;
    const bool nonline = !model_->get_initial_condition().is_zero();
    if (nonline) {
        dof.resize(number_of_dof);
        dynamic_cast<TypedInitialCondition<T>&>(model_->get_initial_condition())
              .get_static_initial_condition_value(dof, time, stage, domain_state_id);
        if (stage != 0) { case_.set_stage(stage); }
        if (!domain_state_id.empty()) { model_->get_domain().load_state(domain_state_id); }
    }
    EquilibriumSolution esf(false);
 
    // 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::Matrix<T, b2linalg::Mcompressed_col> trans_R;
    mrbc->get_linear_value(b2linalg::Vector<T, b2linalg::Vdense_ref>::null, trans_R);
 
    // 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;
    if (nonline) {
        dbc->get_nonlinear_value(
              dof, time, esf, l, trans_L, b2linalg::Vector<T, b2linalg::Vdense_ref>::null, 0);
    } else {
        dbc->get_linear_value(l, trans_L);
    }
 
    if (!l.empty()) {
        // scale the ebc constraint
        b2linalg::Vector<double> trans_L_scale(trans_L.size2());
        trans_L.scale_col_to_norminf(trans_L_scale);
        for (size_t i = 0; i != l.size(); ++i) { l[i] *= trans_L_scale[i]; }
        trans_L.remove_zero();
 
        // procedure to remove dependency on the Lagrange multiplier
        if (case_.get_bool("REMOVE_DEPENDENT_CONSTRAINTS", true)
            && (constraint_string == "LAGRANGE" || constraint_string == "AUGMENTED_LAGRANGE")) {
            b2linalg::Index P;
            trans_L.get_dep_columns(P, 3e-13, 1.5);
            if (!P.empty()) {
                l.remove(P);
                trans_L.remove_col(P);
                logger << logging::info << "Remove " << P.size() << " dependent constraints."
                       << logging::LOGFLUSH;
            }
        }
    }
 
    // Assembly the matrices K and M by computing
    // (trans(R) * K * R) = K_augm = sum (trans(R) * K_element * R)
    // (trans(R) * M * R) = M_augm = sum (trans(R) * M_element * R)
    logger << logging::info << "Element matrix assembly" << logging::LOGFLUSH;
    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_);
    b2linalg::Matrix<T, typename MATRIX_FORMAT::sparse_inversible> M_augm(
          K_augm_size, domain.get_dof_connectivity_type(), case_);
 
    auto Kgg = export_sys_matrices
                     ? b2linalg::Matrix<T, typename MATRIX_FORMAT::sparse_inversible>(
                           number_of_dof, domain.get_dof_connectivity_type(), case_)
                     : b2linalg::Matrix<T, typename MATRIX_FORMAT::sparse_inversible>();
    auto Mgg = export_sys_matrices
                     ? b2linalg::Matrix<T, typename MATRIX_FORMAT::sparse_inversible>(
                           number_of_dof, domain.get_dof_connectivity_type(), case_)
                     : b2linalg::Matrix<T, typename MATRIX_FORMAT::sparse_inversible>();
 
    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);
        d_value_d_dof_flags[0] = d_value_d_dof_flags[2] = true;
        std::vector<b2linalg::Matrix<T, typename MATRIX_FORMAT::dense> > d_value_d_dof(3);
        const 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_, !nonline,
                      true,  // equilibrium_solution
                      time,
                      0,  // delta_time
                      nonline ? b2linalg::Matrix<T, b2linalg::Mrectangle_constref>(dof)
                              : 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_augm += trans_R * StructuredMatrix(number_of_dof, index, d_value_d_dof[0])
                      * transposed(trans_R);
            M_augm += trans_R * StructuredMatrix(number_of_dof, index, d_value_d_dof[2])
                      * transposed(trans_R);
 
            if (export_sys_matrices) {
                Kgg += StructuredMatrix(number_of_dof, index, d_value_d_dof[0]);
                Mgg += 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[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;
 
    if (trans_L.size2() > 0) {
        b2linalg::Matrix<T, b2linalg::Mcompressed_col> trans_LR;
        trans_LR = trans_R * trans_L;
        if (constraint_string == "PENALTY" || constraint_string == "AUGMENTED_LAGRANGE") {
            const double penalty_factor =
                  case_.get_double("CONSTRAINT_PENALTY", constraint_string == "PENALTY" ? 1e10 : 1);
            K_augm += (trans_LR * transposed(trans_LR)) * penalty_factor;
        }
        if (constraint_string == "LAGRANGE" || constraint_string == "AUGMENTED_LAGRANGE") {
            const double lagrange_mult_scale = case_.get_double("LAGRANGE_MULT_SCALE_PENALTY", 1.0);
            trans_LR *= lagrange_mult_scale;
            if (!MATRIX_FORMAT::symmetric) {
                K_augm(
                      b2linalg::Interval(0, trans_R.size1()),
                      b2linalg::Interval(trans_R.size1(), K_augm_size)) += trans_LR;
            }
            b2linalg::Matrix<T, b2linalg::Mcompressed_col> LR = transposed(trans_LR);
            K_augm(
                  b2linalg::Interval(trans_R.size1(), K_augm_size),
                  b2linalg::Interval(0, trans_R.size1())) += LR;
        }
    }
 
    if (logger.is_enabled_for(logging::data)) {
        logger << logging::data << "d_value_d_dof "
               << logging::DBName("D_VALUE_D_DOF.GLOB", 0, 0, case_.get_id()) << K_augm
               << logging::LOGFLUSH;
        logger << logging::data << "d_value_d_dofdotdot "
               << logging::DBName("D_VALUE_D_ACCELERATION.GLOB", 0, 0, case_.get_id()) << M_augm
               << logging::LOGFLUSH;
    }
 
    if (export_sys_matrices) {
        // This section exports the system matrices to the
        // directory given in the MATRICES_DIR option as csv-files.
        // The resulting tables can, for example, be read with np.loadtxt() in
        // Python.
 
        std::filesystem::path out_path(matrices_dir);
        if (std::filesystem::create_directories(out_path)) {
            logger << logging::info << "Created new directory " << matrices_dir
                   << logging::LOGFLUSH;
        }
        // give the user some feedback where to look for the output
        logger << logging::info << "Writing system matrices as CSV file(s) to " << out_path
               << logging::LOGFLUSH;
 
        // K_augm output
        logger << logging::info << " - writing "
               << "Kaa"
               << " with (" << K_augm.size1() << " x " << K_augm.size2() << ") DoF"
               << logging::LOGFLUSH;
        export_csv(K_augm, "Kaa", out_path / "Kaa.csv");
 
        // M_augm output
        logger << logging::info << " - writing "
               << "Maa"
               << " with (" << M_augm.size1() << " x " << M_augm.size2() << ") DoF"
               << logging::LOGFLUSH;
        export_csv(M_augm, "Maa", out_path / "Maa.csv");
 
        // R output
        logger << logging::info << " - writing "
               << "Rtrans"
               << " with (" << trans_R.size1() << " x " << trans_R.size2() << ") DoF"
               << logging::LOGFLUSH;
        export_csv(trans_R, "Rtrans", out_path / "Rtrans.csv");
 
        // L output
        if (trans_L.size2() > 0) {
            logger << logging::info << " - writing "
                   << "Ltrans"
                   << " with (" << trans_L.size1() << " x " << trans_L.size2() << ") DoF"
                   << logging::LOGFLUSH;
            export_csv(trans_L, "Ltrans", out_path / "Ltrans.csv");
        } else {
            logger << logging::info << " - skipping "
                   << "Ltrans"
                   << " with (" << trans_L.size1() << " x " << trans_L.size2() << ") DoF"
                   << logging::LOGFLUSH;
        }
 
        // Kgg output
        logger << logging::info << " - writing "
               << "Kgg"
               << " with (" << Kgg.size1() << " x " << Kgg.size2() << ") DoF" << logging::LOGFLUSH;
        export_csv(Kgg, "Kgg", out_path / "Kgg.csv");
 
        // Mgg output
        logger << logging::info << " - writing "
               << "Mgg"
               << " with (" << Mgg.size1() << " x " << Mgg.size2() << ") DoF" << logging::LOGFLUSH;
        export_csv(Mgg, "Mgg", out_path / "Mgg.csv");
    }
 
    logger << logging::info << "Eigenvalue problem resolution" << logging::LOGFLUSH;
 
    if (nb_eigenmodes_ > K_augm.size1()) {
        Exception() << "The number of requested free-vibration modes "
                       "(NMODES="
                    << nb_eigenmodes_
                    << ") is greater than the size of the "
                       "reduced system of equations ("
                    << K_augm.size1() << ")." << THROW;
    }
    eigenvalues_.resize(nb_eigenmodes_);
    b2linalg::Matrix<T, b2linalg::Mrectangle> eigenvector_red;
    if (nb_eigenmodes_ == K_augm.size1() && K_augm.size1() < 100 && typeid(T) == typeid(double)
        && eig_matrix_type::A == 'P' && eig_matrix_type::B == 'S' && sigma == T(0)
        && constraint_string != "LAGRANGE" && constraint_string != "AUGMENTED_LAGRANGE") {
        eigen_decomposition_dense<T, MATRIX_FORMAT>(
              K_augm, M_augm, which, eigenvalues_, eigenvector_red);
    } else {
        eigenvector_red = eigenvector(
              K_augm, eig_matrix_type::A, M_augm, eig_matrix_type::B, eigenvalues_, which, sigma);
    }
 
    eigenvectors_.resize(number_of_dof, nb_eigenmodes_);
    eigenvectors_ = transposed(trans_R)
                    * eigenvector_red(
                          b2linalg::Interval(0, trans_R.size1()),
                          b2linalg::Interval(0, eigenvector_red.size2()));
 
    solution->set_modes(eigenvalues_(b2linalg::Interval(0, eigenvectors_.size2())), eigenvectors_);
    RTable attributes;
    solution->terminate_case(true, attributes);
 
    case_.warn_on_non_used_key();
 
    logger << logging::info << "End of free vibration solver" << logging::LOGFLUSH;
}
 
#endif  // B2FREE_VIBRATION_SOLVER_H_