b2frequency_dependent_free_vibration_solver.H

//------------------------------------------------------------------------
// b2frequency_dependent_free_vibration_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 B2FREQUENCY_DEPENDENT_FREE_VIBRATION_SOLVER_H_
#define B2FREQUENCY_DEPENDENT_FREE_VIBRATION_SOLVER_H_
 
#include "b2ppconfig.h"
#include "model/b2boundary_condition.H"
#include "model/b2domain.H"
#include "model/b2element.H"
#include "model/b2model.H"
#include "model/b2solution.H"
#include "solvers/b2free_vibration_solver.H"
#include "utils/b2constants.H"
#include "utils/b2exception.H"
#include "utils/b2logging.H"
#include "utils/b2object.H"
 
namespace b2000::solver {
 
template <typename T, typename MATRIX_FORMAT>
class FrequencyDependentFreeVibrationSolver : 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<FrequencyDependentFreeVibrationSolver, 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_);
    logging::Logger& logger{logging::get_logger("solver.frequency_dependent_free_vibration")};
    // Number of eigen modes
    int 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_;
 
    void compute_K_and_M(
          b2linalg::Matrix<T, typename MATRIX_FORMAT::sparse_inversible>& K_augm,
          b2linalg::Matrix<T, typename MATRIX_FORMAT::sparse_inversible>& M_augm,
          b2linalg::Matrix<T, b2linalg::Mcompressed_col>& trans_R,
          b2linalg::Matrix<T, b2linalg::Mcompressed_col>& trans_L) {
        Domain& domain = model_->get_domain();
        size_t number_of_dof = domain.get_number_of_dof();
 
        K_augm.set_zero_same_struct();
        M_augm.set_zero_same_struct();
 
        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);
        bool log_el_flag = logger.is_enabled_for(logging::data);
        for (Element* e = i->next(); e != nullptr; e = i->next()) {
            d_value_d_dof[0].resize(index.size(), index.size());
            d_value_d_dof[2].resize(index.size(), index.size());
            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,  // dof
                      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 (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;
            }
        }
 
        K_augm(
              b2linalg::Interval(trans_R.size1(), K_augm.size1()),
              b2linalg::Interval(0, trans_R.size1())) +=
              b2linalg::transposed(trans_L) * b2linalg::transposed(trans_R);
    }
};
 
template <typename T, typename MATRIX_FORMAT>
typename FrequencyDependentFreeVibrationSolver<T, MATRIX_FORMAT>::type_t
      FrequencyDependentFreeVibrationSolver<T, MATRIX_FORMAT>::type(
            "FrequencyDependentFreeVibrationSolver", type_name<T, MATRIX_FORMAT>(),
            StringList("FREQUENCY_DEPENDENT_FREE_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 accepts only
    one essential boundary condition. The non-linear part of the
    essential boundary condition are ignored.
 
Solve the free vibration frequency dependent 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 become:
  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 eigen-problem 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 eigen-vector x is then obtained from the eigen-vector x_r using the model
reduction equation.
 
 
 
K and M are frequency dependent. In the first step, an approximation
of all the eigenvalues are computed with the initial frequency equal
to the shift value. In the second step each exact frequency dependent
eigenvalues are computed with an iterative approach.
 
*/
template <typename T, typename MATRIX_FORMAT>
void b2000::solver::FrequencyDependentFreeVibrationSolver<T, MATRIX_FORMAT>::solve_on_case(
      Case& case_) {
    logger << logging::info << "Start the frequency dependent 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 solutions of type "
                       "SolutionFreeVibration<T>"
                    << THROW;
    }
 
    T sigma;
    T sigma_max;
    char which[3];
    solution->which_modes(nb_eigenmodes_, which, sigma, sigma_max);
 
    double tol_exact_value = case_.get_double("TOL", 0.0000001);
    double tol_predictor_value = case_.get_double("TOL_PREDICTOR", 0.01);
 
    // Get the number of dof of the model
    size_t number_of_dof = domain.get_number_of_dof();
 
    // 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);
 
    // 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);
 
    size_t K_augm_size = trans_R.size1() + 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_);
 
    logger << logging::info
           << "Compute an approximation of all the required modes with spectral shift value = "
           << sqrt(sigma) / T(2 * M_PIl) << "." << logging::LOGFLUSH;
 
    // First compute an approximation of all the eigenvalues.
    // use the value that come form the database material.
    // case_.set("FREQ_INIT", sqrt(sigma) / T(2 * M_PIl));
    // domain.set_case(case_);
 
    compute_K_and_M(K_augm, M_augm, trans_R, trans_L);
 
    eigenvalues_.resize(nb_eigenmodes_);
    eigenvectors_.resize(number_of_dof, nb_eigenmodes_);
    b2linalg::Matrix<T, b2linalg::Mrectangle> eigenvectors_red;
    eigenvectors_red = eigenvector(
          K_augm, eig_matrix_type::A, M_augm, eig_matrix_type::B, eigenvalues_, which, sigma);
 
    // Iterate on the approximated eigenvalues.
    int mode = 0;
    while (mode != nb_eigenmodes_) {
        T local_sigma = eigenvalues_[mode];
 
        logger << logging::info << "Start to compute the exact (frequency dependent) mode " << mode
               << " with the spectral shift value " << sqrt(local_sigma) / T(2 * M_PIl) << "."
               << logging::LOGFLUSH;
 
        T previous_eigenvalue;
        b2linalg::Vector<T, b2linalg::Vdense> current_eigenvalue(1);
        current_eigenvalue[0] = eigenvalues_[mode];
        b2linalg::Matrix<T, b2linalg::Mrectangle> eigenvector_red;
 
        b2linalg::Vector<T, b2linalg::Vdense> restart_eigenvector;
        restart_eigenvector = eigenvectors_red[mode];
 
        int ncv = 2;
 
        // Iterate at one frequency until the frequency residue (the residue
        // between the previous current_eigenvalue and the new current_eigenvalue
        // after an update of the frequency of the domain) become small.
        do {
            logger << logging::info << "Update the domain frequency to "
                   << sqrt(current_eigenvalue[0]) / T(2 * M_PIl) << logging::LOGFLUSH;
 
            // Update the frequency of the domain.
            case_.set("FREQ_INIT", sqrt(current_eigenvalue[0]) / T(2 * M_PIl));
            domain.set_case(case_);
 
            // Compute the stiffness and mass matrix for the new frequency.
            compute_K_and_M(K_augm, M_augm, trans_R, trans_L);
 
            previous_eigenvalue = current_eigenvalue[0];
 
            while (1) {
                try {
                    eigenvector_red = eigenvector(
                          K_augm, eig_matrix_type::A, M_augm, eig_matrix_type::B,
                          current_eigenvalue, "SM", local_sigma, ncv, -1, -1, restart_eigenvector);
                    restart_eigenvector = eigenvector_red[0];
                    local_sigma = current_eigenvalue[0];
                    break;
                } catch (Exception& e) {
                    ncv = std::min(32, ncv * 2);
                    logger << logging::info << "The current_eigenvalue solver cannot converge ("
                           << e
                           << "). Retry by increasing the number of Lanczos basis vectors (ncv="
                           << ncv << "." << logging::LOGFLUSH;
                }
            }
 
            logger << logging::info << "Frequency residue = "
                   << std::abs(current_eigenvalue[0] - previous_eigenvalue)
                            / std::abs(current_eigenvalue[0])
                   << logging::LOGFLUSH;
        } while (std::abs(current_eigenvalue[0]) > tol_exact_value
                 && (std::abs(current_eigenvalue[0] - previous_eigenvalue)
                           / std::abs(current_eigenvalue[0])
                     > tol_exact_value));
 
        if (std::abs(eigenvalues_[mode] - current_eigenvalue[0]) / std::abs(current_eigenvalue[0])
            > tol_predictor_value) {
            logger << logging::info
                   << "The approximation of all the required modes is too far from the exact "
                      "computed mode, update the approximation with the domain frequency equal to "
                   << sqrt(current_eigenvalue[0]) / T(2 * M_PIl) << logging::LOGFLUSH;
 
            eigenvectors_red = eigenvector(
                  K_augm, eig_matrix_type::A, M_augm, eig_matrix_type::B, eigenvalues_, which,
                  sigma);
 
            if (std::abs(eigenvalues_[mode] - current_eigenvalue[0])
                      / std::abs(current_eigenvalue[0])
                > tol_exact_value) {
                logger << logging::info
                       << "The exact computed mode is a jump to an other mode, restart the "
                          "computation from the mode "
                       << mode << "." << logging::LOGFLUSH;
                continue;
            }
        }
 
        logger << logging::info << "End of computation of the mode " << mode
               << ". Eigenfrequency = " << sqrt(current_eigenvalue[0]) / T(2 * M_PIl) << "."
               << logging::LOGFLUSH;
 
        // Convert the modes from the reduction model
        eigenvectors_[mode] = transposed(trans_R) * eigenvector_red[0];
 
        // Store the computed mode.
        solution->set_mode(mode, current_eigenvalue[0], eigenvectors_[mode]);
 
        // go to the next mode
        ++mode;
    }
 
    RTable attributes;
    solution->terminate_case(true, attributes);
    case_.warn_on_non_used_key();
 
    logger << logging::info << "End of frequency dependent free vibration solver"
           << logging::LOGFLUSH;
}
 
#endif  // B2FREQUENCY_DEPENDENT_FREE_VIBRATION_SOLVER_H_