polytropeOperator.h

Go to the documentation of this file.

/* ***********************************************************************
//
//   Copyright (C) 2025 -- The 4D-STAR Collaboration
//   File Author: Emily Boudreaux
//   Last Modified: April 21, 2025
//
//   4DSSE is free software; you can use it and/or modify
//   it under the terms and restrictions the GNU General Library Public
//   License version 3 (GPLv3) as published by the Free Software Foundation.
//
//   4DSSE is distributed in the hope that it will be useful,
//   but WITHOUT ANY WARRANTY; without even the implied warranty of
//   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
//   See the GNU Library General Public License for more details.
//
//   You should have received a copy of the GNU Library General Public License
//   along with this software; if not, write to the Free Software
//   Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
//
// *********************************************************************** */
#pragma once
 
#include "mfem.hpp"
#include "4DSTARTypes.h"
#include <memory>
 
#include "probe.h"
 
namespace serif::polytrope {

class SchurCompliment final : public mfem::Operator {
public:
    SchurCompliment(
        const mfem::SparseMatrix &QOp,
        const mfem::SparseMatrix &DOp,
        const mfem::SparseMatrix &MOp,
        const mfem::Solver &GradInvOp
    );

    SchurCompliment(
        const mfem::SparseMatrix &QOp,
        const mfem::SparseMatrix &DOp,
        const mfem::SparseMatrix &MOp
    );

    ~SchurCompliment() override = default;

    void Mult(const mfem::Vector &x, mfem::Vector &y) const override;

    void SetOperator(const mfem::SparseMatrix &QOp, const mfem::SparseMatrix &DOp, const mfem::SparseMatrix &MOp, const mfem::Solver &GradInvOp);

    void updateInverseNonlinearJacobian(const mfem::Solver &gradInv);
 
private:
    void updateConstantTerms(const mfem::SparseMatrix &QOp, const mfem::SparseMatrix &DOp, const mfem::SparseMatrix &MOp);
 
private:
    // Pointers to external operators (not owned by this class)
    const mfem::SparseMatrix* m_QOp = nullptr;         
    const mfem::SparseMatrix* m_DOp = nullptr;         
    const mfem::SparseMatrix* m_MOp = nullptr;         
    const mfem::Solver* m_GradInvOp = nullptr;         
 
    // Dimensions
    int m_nPhi = 0;                                    
    int m_nTheta = 0;                                  
 
    // Owned resources
    mutable std::unique_ptr<mfem::SparseMatrix> m_matrixForm; 
};

class GMRESInverter final : public mfem::Operator {
public:
    explicit GMRESInverter(const SchurCompliment& op);

    ~GMRESInverter() override = default;

    void Mult(const mfem::Vector &x, mfem::Vector &y) const override;
 
private:
    const SchurCompliment& m_op; 
    mfem::GMRESSolver m_solver;   
};

class PolytropeOperator final : public mfem::Operator {
public:
    PolytropeOperator(
        std::unique_ptr<mfem::MixedBilinearForm> M,
        std::unique_ptr<mfem::MixedBilinearForm> Q,
        std::unique_ptr<mfem::BilinearForm> D,
        std::unique_ptr<mfem::BilinearForm> S,
        std::unique_ptr<mfem::NonlinearForm> f,
        const mfem::Array<int> &blockOffsets);

    ~PolytropeOperator() override = default;

    void Mult(const mfem::Vector &xFree, mfem::Vector &yFree) const override;

    mfem::Operator& GetGradient(const mfem::Vector &xFree) const override;

    void finalize(const mfem::Vector &initTheta);

    bool isFinalized() const { return m_isFinalized; }

    void set_essential_true_dofs(const serif::types::MFEMArrayPair& theta_ess_tdofs, const serif::types::MFEMArrayPair& phi_ess_tdofs);

    void set_essential_true_dofs(const serif::types::MFEMArrayPairSet& ess_tdof_pair_set);
 

    [[nodiscard]] const mfem::Vector &reconstruct_full_state_vector(const mfem::Vector &reducedState) const;

    [[nodiscard]] const mfem::BlockVector reconstruct_full_block_state_vector(const mfem::Vector &reducedState) const;

    void populate_free_dof_array();


    const mfem::BlockOperator &get_jacobian_operator() const;

    mfem::BlockDiagonalPreconditioner &get_preconditioner() const;
 


    const mfem::Array<int>& get_free_dofs() const { return m_freeDofs; } 

    int get_reduced_system_size() const;

    serif::types::MFEMArrayPairSet get_essential_true_dofs() const;

    const mfem::Array<int>& get_block_offsets() const { return m_blockOffsets; }

    const mfem::Array<int>& get_reduced_block_offsets() const {return m_reducedBlockOffsets; }
 
private:
    // --- Logging ---
    serif::probe::LogManager& m_logManager = serif::probe::LogManager::getInstance(); 
    quill::Logger* m_logger = m_logManager.getLogger("log");            
 
    // --- Input Bilinear/Nonlinear Forms (owned) ---
    std::unique_ptr<mfem::MixedBilinearForm> m_M; 
    std::unique_ptr<mfem::MixedBilinearForm> m_Q; 
    std::unique_ptr<mfem::BilinearForm> m_D;      
    std::unique_ptr<mfem::BilinearForm> m_S;      
    std::unique_ptr<mfem::NonlinearForm> m_f;     
 
    // --- Full Matrix Representations (owned, derived from forms) ---
    std::unique_ptr<mfem::SparseMatrix> m_Mmat;   
    std::unique_ptr<mfem::SparseMatrix> m_Qmat;   
    std::unique_ptr<mfem::SparseMatrix> m_Dmat;   
    std::unique_ptr<mfem::SparseMatrix> m_Smat;
 
    // --- Reduced Matrix Representations (owned, after eliminating essential DOFs) ---
    std::unique_ptr<mfem::SparseMatrix> m_MReduced; 
    std::unique_ptr<mfem::SparseMatrix> m_QReduced; 
    std::unique_ptr<mfem::SparseMatrix> m_DReduced; 
    std::unique_ptr<mfem::SparseMatrix> m_SReduced; 
    mutable std::unique_ptr<mfem::SparseMatrix> m_gradReduced; 
 
    // --- Scaled Reduced Matrix Representations (owned, after eliminating essential DOFs and scaling by stabilization coefficients) ---
    std::unique_ptr<mfem::SparseMatrix> m_MScaledReduced; 
    std::unique_ptr<mfem::SparseMatrix> m_QScaledReduced; 
    std::unique_ptr<mfem::SparseMatrix> m_DScaledReduced; 
    std::unique_ptr<mfem::SparseMatrix> m_ScaledSReduced; 
 
 
    // --- Stabilization Coefficient --- (Perhapses this should be a user parameter...) // TODO: Sort out why this is negative (sign error?)
    static constexpr double m_stabilizationCoefficient = -2.0;       
    static constexpr double m_IncrementedStabilizationCoefficient = 1.0 + m_stabilizationCoefficient; 
 
    // --- State Vectors and DOF Management ---
    mutable mfem::Vector m_state;                 
    mfem::Array<int> m_freeDofs;                  
 
    // --- Block Offsets ---
    const mfem::Array<int> m_blockOffsets;        
    mfem::Array<int> m_reducedBlockOffsets;       
 
    // --- Essential Boundary Conditions ---
    serif::types::MFEMArrayPair m_theta_ess_tdofs;         
    serif::types::MFEMArrayPair m_phi_ess_tdofs;           
 
    // --- Jacobian and Preconditioner Components (owned, mutable) ---
    std::unique_ptr<mfem::ScaledOperator> m_negQ_mat; 
    mutable std::unique_ptr<mfem::BlockOperator> m_jacobian; 
    mutable std::unique_ptr<SchurCompliment> m_schurCompliment; 
    mutable std::unique_ptr<GMRESInverter> m_invSchurCompliment; 
    mutable std::unique_ptr<mfem::Solver> m_invNonlinearJacobian; 
    mutable std::unique_ptr<mfem::BlockDiagonalPreconditioner> m_schurPreconditioner; 
 
    // --- State Flags ---
    bool m_isFinalized = false;                   
 
private:
    void construct_matrix_representations();

    void construct_reduced_block_offsets();

    void construct_jacobian_constant_terms();

    void scatter_boundary_conditions();

    void update_inverse_nonlinear_jacobian(const mfem::Operator &grad) const;

    void update_inverse_schur_compliment() const;

    void update_preconditioner(const mfem::Operator &grad) const;
};
 
} // namespace serif::polytrope