FETI_operations.h

Go to the documentation of this file.

/*
 * FETI_operations.h
 *
 *  Created on: Apr 23, 2017
 *      Author: Thiago Milanetto Schlittler
 */

#ifndef FETI_OPERATIONS_H_
#define FETI_OPERATIONS_H_

#include "carl_headers.h"
#include "common_enums.h"
#include "PETSC_matrix_operations.h"

namespace carl
{

class   FETI_Operations
{
protected:

    //  --- Miscellaneous declarations
    libMesh::Parallel::Communicator& m_comm;    
    bool        m_bScratchFolderSet;            
    bool        m_bCouplingFolderSet;           
    std::string m_scratch_folder_path;          
    std::string m_coupling_folder_path;         
    IterationStatus m_FETI_solver_status;       
    int         m_kkk;                          

    //  --- Coupling matrix and preconditioner declarations

    // Coupling matrices
    Mat m_C_R_micro;    
    Mat m_C_R_BIG;      
    Mat m_C_RR;         

    // Matrix dimensions
    /*
     *                   M     x N
     * m_C_R_micro     : n_med x n_sys_micro
     * m_C_R_BIG       : n_med x n_sys_BIG
     * m_C_RR          : n_med x n_med
     *
     */
    PetscInt    m_C_R_micro_M, m_C_R_micro_N, m_C_R_micro_M_local, m_C_R_micro_N_local;
    PetscInt    m_C_R_BIG_M, m_C_R_BIG_N, m_C_R_BIG_M_local, m_C_R_BIG_N_local;
    PetscInt    m_C_RR_M, m_C_RR_M_local;

    // Objects used by the preconditioner
    KSP         m_coupling_precond_solver;  
    Vec         m_coupling_jacobi_precond_vec;  

    // Boolean flags
    bool        m_bC_R_micro_MatrixSet; 
    bool        m_bC_R_BIG_MatrixSet;   
    bool        m_bC_RR_MatrixSet;      
    bool        m_bCouplingMatricesSet; 

    bool        m_bmicro_sizes_set;     
    bool        m_bBIG_sizes_set;       
    bool        m_bR_sizes_set;         

    BaseCGPrecondType   m_precond_type; 

    //  --- Null space / rigid body modes declarations
    RBModesSystem m_RB_modes_system;    

    bool        m_bCreatedPrecondSolver;    
    bool        m_bCreatedPrecondJacobiVec; 

    // Vectors
    /*
     *      Note: using Vec arrays does not pose a memory problem, since a PETSc `Vec` is essentially a 
     *  pointer in the stack. The `Vec`'s contents arre allocated in the heap when calling functions such as
     *  `VecCreate` and such.
     *  
     */ 
    PetscInt    m_null_nb_vecs; 
    enum        { maxVecLimit = 6 };    
    Vec         m_null_vecs[maxVecLimit];           
    Vec         m_null_coupled_vecs[maxVecLimit];   

    // Null vector dimensions
    PetscInt    m_null_vecs_N, m_null_vecs_N_local;

    // Matrices
    Mat m_RITRI_mat;        
    Mat m_inv_RITRI_mat;    

    // Boolean flags
    bool        m_bUsingNullVecs;               
    bool        m_bNullVecsSet;                 
    bool        m_bNullVecsDimensionsSet;       
    bool        m_binvRITRIMatSet;              

    //  --- FETI operations declarations
    // Vectors
    Vec     m_u_0_micro;                
    Vec     m_u_0_BIG;                  

    Vec     m_ext_solver_sol_micro;     
    Vec     m_ext_solver_sol_BIG;       

    Vec     m_current_residual;         
    Vec     m_current_z;                
    Vec     m_current_p;                

    Vec     m_current_phi;              
    Vec     m_current_rb_correction;    

    Vec     m_previous_residual;        
    Vec     m_previous_phi;             
    Vec *   m_previous_p_ptr;           
    Vec *   m_previous_q_ptr;           

    // Real values
    double m_gamma;                     
    double m_rho_0;                     
    double m_current_rho;               
    double m_current_RB_mode_corr;      
    double m_previous_rho;              
    double m_previous_RB_mode_corr;     

    // Boolean flags
    bool    m_bSet_u_0;                     
    bool    m_bSet_ext_solver_sol;          

    bool    m_bSet_current_residual;        
    bool    m_bSet_current_z;               
    bool    m_bSet_current_p;               
    bool    m_bSet_current_phi;             
    bool    m_bSet_current_RB_correction;   

    bool    m_bSet_previous_residual;   
    bool    m_bSet_previous_phi;        
    bool    m_bSet_previous_p_ptr;      
    bool    m_bSet_previous_q_ptr;      

    //  --- Vectors containing the scalar data
    std::vector<double> m_p_dot_q;      

    // Boolean flags
    bool    m_bReadPreviousScalar;      
    bool    m_bCalculatedScalar;        

    // --- Convergence parameters
    double m_abs_residual_conv;         
    double m_rel_residual_conv;         
    double m_rb_modes_conv;             
    double m_rel_residual_div;          
    int m_max_iter_div;                 

    // Boolean flags
    bool    m_bConvResidualAbs;         
    bool    m_bConvResidualRel;         
    bool    m_bConvRBCorrRel;           
    bool    m_bDivResidualRel;          
    bool    m_bDivResidualNeg;          
    bool    m_bDivIter;                 
    bool    m_bConv;                    
    bool    m_bDiv;                     

    //  --- Converged coupled solution vectors
    Vec     m_coupled_sol_micro;    
    Vec     m_coupled_sol_BIG;      

    bool    m_bCoupled_sols_set;    

    //  --- Protected methods
    // Preconditioner methods
    void set_inverse_precond_solver();  

    void set_jacobi_precond_vector();   

    void read_jacobi_precond_vector();  

    void apply_inverse_coupling_precond(Vec vec_in, Vec vec_out);   

    void apply_jacobi_coupling_precond(Vec vec_in, Vec vec_out);    

    void apply_precond(Vec vec_in, Vec vec_out);    

    // Projection methods
    void apply_RB_projection(Vec vec_in, Vec vec_out);  

    void export_ext_solver_rhs(Vec vec_in);

    void clear_PETSc();

    FETI_Operations();

public:
    FETI_Operations(libMesh::Parallel::Communicator& comm, const std::string& scratch_folder_path , const std::string& coupling_folder_path) :
        m_comm { comm },
        m_bScratchFolderSet { true },
        m_bCouplingFolderSet { true },
        m_scratch_folder_path { scratch_folder_path },
        m_coupling_folder_path { coupling_folder_path },
        m_FETI_solver_status { IterationStatus::ITERATING },
        m_kkk { 0 },
        m_C_R_micro_M { -1},
        m_C_R_micro_N { -1},
        m_C_R_micro_M_local { -1},
        m_C_R_micro_N_local { -1},
        m_C_R_BIG_M { -1},
        m_C_R_BIG_N { -1},
        m_C_R_BIG_M_local { -1},
        m_C_R_BIG_N_local { -1},
        m_C_RR_M { -1},
        m_C_RR_M_local { -1},
        m_bC_R_micro_MatrixSet  { false },
        m_bC_R_BIG_MatrixSet  { false },
        m_bC_RR_MatrixSet  { false },
        m_bCouplingMatricesSet  { false },
        m_bmicro_sizes_set { false },
        m_bBIG_sizes_set { false },
        m_bR_sizes_set { false },
        m_precond_type { BaseCGPrecondType::NO_PRECONDITIONER },
        m_RB_modes_system { RBModesSystem::MICRO },
        m_bCreatedPrecondSolver { false },
        m_bCreatedPrecondJacobiVec { false },
        m_null_nb_vecs { -1 },
        m_null_vecs_N { -1 },
        m_null_vecs_N_local { -1 },
        m_bUsingNullVecs { false },
        m_bNullVecsSet { false },
        m_bNullVecsDimensionsSet { false },
        m_binvRITRIMatSet { false }, 
        m_gamma { -1 },
        m_rho_0 { -1 },
        m_current_rho { -1 },
        m_current_RB_mode_corr { -1 },
        m_previous_rho { -1 },
        m_previous_RB_mode_corr { -1 },
        m_bSet_u_0 { false },
        m_bSet_ext_solver_sol { false },
        m_bSet_current_residual { false },
        m_bSet_current_z { false },
        m_bSet_current_p { false },
        m_bSet_current_phi { false },
        m_bSet_current_RB_correction { false },
        m_bSet_previous_residual { false },
        m_bSet_previous_phi { false },
        m_bSet_previous_p_ptr { false },
        m_bSet_previous_q_ptr { false },
        m_bReadPreviousScalar { false },
        m_bCalculatedScalar { false },
        m_abs_residual_conv { -1 },
        m_rel_residual_conv { -1 },
        m_rb_modes_conv { -1 },
        m_rel_residual_div { -1 },
        m_max_iter_div { -1 },
        m_bConvResidualAbs { false },
        m_bConvResidualRel { false },
        m_bConvRBCorrRel { false },
        m_bDivResidualRel { false },
        m_bDivResidualNeg { false },
        m_bDivIter { false },
        m_bConv { false },
        m_bDiv { false },
        m_bCoupled_sols_set {false}
    {
    };

    FETI_Operations(libMesh::Parallel::Communicator& comm, const std::string& scratch_folder_path) :
        m_comm { comm },
        m_bScratchFolderSet { true },
        m_bCouplingFolderSet { false },
        m_scratch_folder_path { scratch_folder_path },
        m_coupling_folder_path { "" },
        m_FETI_solver_status { IterationStatus::ITERATING },
        m_kkk { 0 },
        m_C_R_micro_M { -1},
        m_C_R_micro_N { -1},
        m_C_R_micro_M_local { -1},
        m_C_R_micro_N_local { -1},
        m_C_R_BIG_M { -1},
        m_C_R_BIG_N { -1},
        m_C_R_BIG_M_local { -1},
        m_C_R_BIG_N_local { -1},
        m_C_RR_M { -1},
        m_C_RR_M_local { -1},
        m_bC_R_micro_MatrixSet  { false },
        m_bC_R_BIG_MatrixSet  { false },
        m_bC_RR_MatrixSet  { false },
        m_bCouplingMatricesSet  { false },
        m_bmicro_sizes_set { false },
        m_bBIG_sizes_set { false },
        m_bR_sizes_set { false },
        m_precond_type { BaseCGPrecondType::NO_PRECONDITIONER },
        m_RB_modes_system { RBModesSystem::MICRO },
        m_bCreatedPrecondSolver { false },
        m_bCreatedPrecondJacobiVec { false },
        m_null_nb_vecs { -1 },
        m_null_vecs_N { -1 },
        m_null_vecs_N_local { -1 },
        m_bUsingNullVecs { false },
        m_bNullVecsSet { false },
        m_bNullVecsDimensionsSet { false },
        m_binvRITRIMatSet { false }, 
        m_gamma { -1 },
        m_rho_0 { -1 },
        m_current_rho { -1 },
        m_current_RB_mode_corr { -1 },
        m_previous_rho { -1 },
        m_previous_RB_mode_corr { -1 },
        m_bSet_u_0 { false },
        m_bSet_ext_solver_sol { false },
        m_bSet_current_residual { false },
        m_bSet_current_z { false },
        m_bSet_current_p { false },
        m_bSet_current_phi { false },
        m_bSet_current_RB_correction { false },
        m_bSet_previous_residual { false },
        m_bSet_previous_phi { false },
        m_bSet_previous_p_ptr { false },
        m_bSet_previous_q_ptr { false },
        m_bReadPreviousScalar { false },
        m_bCalculatedScalar { false },
        m_abs_residual_conv { -1 },
        m_rel_residual_conv { -1 },
        m_rb_modes_conv { -1 },
        m_rel_residual_div { -1 },
        m_max_iter_div { -1 },
        m_bConvResidualAbs { false },
        m_bConvResidualRel { false },
        m_bConvRBCorrRel { false },
        m_bDivResidualRel { false },
        m_bDivResidualNeg { false },
        m_bDivIter { false },
        m_bConv { false },
        m_bDiv { false },
        m_bCoupled_sols_set {false}
    {
    };

    ~FETI_Operations()
    {
        this->clear_PETSc();
    };

    //  --- Coupling matrix and preconditioner methods
    void set_coupling_matrix_R_micro();

    void set_coupling_matrix_R_BIG();

    void set_coupling_matrix_RR();

    void read_coupling_matrices();

    void set_preconditioner(BaseCGPrecondType CG_precond_type, bool bInitialSet = true);

    //  --- Null space / rigid body modes methods
    void using_rb_modes(bool bUseRigidBodyModes);

    void set_null_space(const std::string& input_filename_base, int nb_of_vecs);

    void read_null_space_vecs(const std::string& RB_vectors_base, int nb_of_rb_vectors);

    void read_null_space_inv_RITRI_mat();

    void calculate_null_space_phi_0(const std::string& force_path);

    //  --- Read methods
    void read_decoupled_solutions();

    void read_ext_solver_output();

    void read_rb_corr();

    void read_previous_phi();

    void read_previous_r();

    void read_all_previous_p();

    void read_all_previous_q();

    void read_scalar_data();

    void read_vector_data();

    //  --- FETI steps methods
    void calculate_initial_p();

    void calculate_initial_r();

    void calculate_p();

    void calculate_q();

    void calculate_r();

    void calculate_z();

    void calculate_phi();
    
    void calculate_rb_correction();

    void calculate_scalar_data();

    void calculate_coupled_solution();

    IterationStatus check_convergence(double rel_residual_conv, double abs_residual_conv, int max_iter_div, double rel_residual_div, double rb_modes_conv = -1);

    void increase_iter_counter();

    //  --- Write methods
    void export_ext_solver_rhs_Ct_p();

    void export_ext_solver_rhs_initial();

    void export_ext_solver_rhs_Ct_phi();

    void export_rb_correction_vector();

    void export_p();

    void export_q();

    void export_r();

    void export_phi();

    void export_inital_vecs();

    void export_initial_scalar_data();

    void export_scalar_data();

    // Export the iteration vectors
    void export_iter_vecs();

    // Export the final coupled solution
    void export_coupled_solution(std::string output_base);

    // Print on `std::cout` the current values of the convergence parameters - and if we converged
    void print_previous_iters_conv(int nb_of_iters = 5);
};
}

#endif /* FETI_OPERATIONS_H_ */