polySolver.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 <memory>
#include <utility>
 
#include "integrators.h"
#include "4DSTARTypes.h"
#include "polytropeOperator.h"
#include "config.h"
#include "meshIO.h"
#include "probe.h"
#include "quill/Logger.h"
 
namespace serif {
namespace polytrope {

namespace laneEmden {
    double a (const int k, const double n);

    double c(const int m, const double n);

    double thetaSeriesExpansion(const double xi, const double n, const int order);
} // namespace laneEmden

struct solverBundle {
    mfem::GMRESSolver solver; 
    mfem::NewtonSolver newton; 
};

struct formBundle {
    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;
};

class PolySolver final{
public: // Public methods
    PolySolver(const double n, const double order);

    ~PolySolver();

    void solve() const;

    double getN() const { return m_polytropicIndex; }

    double getOrder() const { return m_feOrder; }

    mfem::Mesh& getMesh() const { return m_mesh; }

    mfem::GridFunction& getTheta() const { return *m_theta; }

    mfem::GridFunction& getPhi() const { return *m_phi; }
 
private: // Private Attributes
    // --- Configuration and Logging ---
    serif::config::Config& m_config;                                 
    serif::probe::LogManager& m_logManager;                  
    quill::Logger* m_logger;                          
 
    // --- Physical and Discretization Parameters ---
    double m_polytropicIndex;                         
    double m_feOrder;                                 
 
    // --- MFEM Core Objects ---
    mfem::Mesh& m_mesh;                               
    std::unique_ptr<mfem::H1_FECollection> m_fecH1;   
    std::unique_ptr<mfem::RT_FECollection> m_fecRT;   
 
    std::unique_ptr<mfem::FiniteElementSpace> m_feTheta; 
    std::unique_ptr<mfem::FiniteElementSpace> m_fePhi;   
 
    // --- Solution Vectors ---
    std::unique_ptr<mfem::GridFunction> m_theta;      
    std::unique_ptr<mfem::GridFunction> m_phi;        
 
    // --- Operator and Solver Components ---
    std::unique_ptr<PolytropeOperator> m_polytropOperator; 
    std::unique_ptr<mfem::OperatorJacobiSmoother> m_prec;  
 
private: // Private methods
    PolySolver(mfem::Mesh& mesh, double n, double order);

    static mfem::Mesh& prepareMesh(double n);

    void assembleBlockSystem();

    mfem::Array<int> computeBlockOffsets() const;

    std::unique_ptr<formBundle> buildIndividualForms(const mfem::Array<int>& blockOffsets);

    static void assembleAndFinalizeForm(auto &f);

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

    std::pair<mfem::Array<int>, mfem::Array<int>> findCenterElement() const;

    void setInitialGuess() const;

    void saveAndViewSolution(const mfem::BlockVector& state_vector) const;

    solverBundle setupNewtonSolver() const;

    void setOperatorEssentialTrueDofs() const;

    void LoadSolverUserParams(double &newtonRelTol, double &newtonAbsTol, int &newtonMaxIter, int &newtonPrintLevel,
                              double &gmresRelTol, double &gmresAbsTol, int &gmresMaxIter, int &gmresPrintLevel) const;

    static void GetDofCoordinates(const mfem::FiniteElementSpace &fes, const std::string& filename);
 
};
 
} // namespace polytrope
} // namespace serif