tboudreaux/SERiF
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

feat(poly): moved to a block form for poly

Browse Source 
essential dofs can be applied to both theta and phi (grad theta) if we move to a block form. I have done this derivation and made that change so that we can properly apply the central boundary condition to the slope
This commit is contained in:
Emily Boudreaux
2025-04-02 14:57:37 -04:00
parent 407eef4e48
commit e3afe90f37
12 changed files with 640 additions and 731 deletions
Download Patch File Download Diff File Expand all files Collapse all files

300
src/poly/solver/private/polySolver.cpp
View File

@@ -23,12 +23,16 @@
#include <memory>
#include <string>
#include <stdexcept>
#include <utility>
#include "polySolver.h"
#include "polyMFEMUtils.h"
#include "integrators.h"
#include "polyCoeff.h"
#include "probe.h"
#include "config.h"
#include "probe.h"
#include "resourceManager.h"
#include "resourceManagerTypes.h"
#include "operator.h"
#include "quill/LogMacros.h"
@@ -63,116 +67,152 @@ namespace laneEmden {
}
}
PolySolver::PolySolver(double n, double order, mfem::Mesh& mesh_)
: logger(logManager.getLogger("log")),
n(n),
order(order),
mesh(mesh_),
feCollection(std::make_unique<mfem::H1_FECollection>(order, mesh.SpaceDimension())),
feSpace(std::make_unique<mfem::FiniteElementSpace>(&mesh, feCollection.get())),
compositeIntegrator(std::make_unique<polyMFEMUtils::CompositeNonlinearIntegrator>()),
nonlinearForm(std::make_unique<mfem::NonlinearForm>(feSpace.get())),
u(std::make_unique<mfem::GridFunction>(feSpace.get())) {
PolySolver::PolySolver(double n, double order) {
// --- Check the polytropic index ---
if (n > 4.99 || n < 0.0) {
LOG_ERROR(logger, "The polytropic index n must be less than 5.0 and greater than 0.0. Currently it is {}", n);
throw std::runtime_error("The polytropic index n must be less than 5.0 and greater than 0.0. Currently it is " + std::to_string(n));
LOG_ERROR(m_logger, "The polytropic index n must be less than 5.0 and greater than 0.0. Currently it is {}", m_polytropicIndex);
throw std::runtime_error("The polytropic index n must be less than 5.0 and greater than 0.0. Currently it is " + std::to_string(m_polytropicIndex));
}
diffusionCoeff = std::make_unique<mfem::VectorFunctionCoefficient>(mesh.SpaceDimension(), polycoeff::diffusionCoeff);
nonlinearSourceCoeff = std::make_unique<mfem::FunctionCoefficient>(polycoeff::nonlinearSourceCoeff);
m_polytropicIndex = n;
m_feOrder = order;
assembleNonlinearForm();
ResourceManager& rm = ResourceManager::getInstance();
const Resource& resource = rm.getResource("mesh:polySphere");
const auto &meshIO = std::get<std::unique_ptr<MeshIO>>(resource);
meshIO->LinearRescale(polycoeff::x1(m_polytropicIndex));
m_mesh = std::make_unique<mfem::Mesh>(meshIO->GetMesh());
// Use feOrder - 1 for the RT space to satisfy Brezzi-Babuska condition
// for the H1 and RT spaces
m_fecH1 = std::make_unique<mfem::H1_FECollection>(m_feOrder, m_mesh->SpaceDimension());
m_fecRT = std::make_unique<mfem::RT_FECollection>(m_feOrder - 1, m_mesh->SpaceDimension());
m_feTheta = std::make_unique<mfem::FiniteElementSpace>(m_mesh.get(), m_fecH1.get());
m_fePhi = std::make_unique<mfem::FiniteElementSpace>(m_mesh.get(), m_fecRT.get());
m_theta = std::make_unique<mfem::GridFunction>(m_feTheta.get());
m_phi = std::make_unique<mfem::GridFunction>(m_fePhi.get());
assembleBlockSystem();
}
PolySolver::~PolySolver() {}
void PolySolver::assembleNonlinearForm() {
// Add the \int_{\Omega}\nabla v\cdot\nabla\theta d\Omegaterm
auto wrappedDiffusionIntegrator = std::make_unique<polyMFEMUtils::BilinearIntegratorWrapper>(
new mfem::DiffusionIntegrator(*diffusionCoeff)
);
compositeIntegrator->add_integrator(wrappedDiffusionIntegrator.release());
void PolySolver::assembleBlockSystem() {
// Add the \int_{\Omega}v\theta^{n} d\Omega term
auto nonlinearIntegrator = std::make_unique<polyMFEMUtils::NonlinearPowerIntegrator>(*nonlinearSourceCoeff, n);
compositeIntegrator->add_integrator(nonlinearIntegrator.release());
// Start by defining the block structure of the system
// Block 0: Theta (scalar space, uses m_feTheta)
// Block 1: Phi (vector space, uses m_fePhi)
mfem::Array<mfem::FiniteElementSpace*> feSpaces;
feSpaces.Append(m_feTheta.get());
feSpaces.Append(m_fePhi.get());
// Add the contraint term \gamma(\nabla \theta(0)\cdot\nabla v(0))^{2}
double gamma = config.get<double>("Poly:Solver:Constraint:Gamma", 1e4);
auto constraintIntegrator = std::make_unique<polyMFEMUtils::BilinearIntegratorWrapper>(
new polyMFEMUtils::ConstraintIntegrator(gamma, &mesh)
);
compositeIntegrator->add_integrator(constraintIntegrator.release());
// Create the block offsets. These define the start of each block in the combined vector.
// Block offsets will be [0, thetaDofs, thetaDofs + phiDofs]
mfem::Array<int> blockOffsets;
blockOffsets.SetSize(3);
blockOffsets[0] = 0;
blockOffsets[1] = feSpaces[0]->GetVSize();
blockOffsets[2] = feSpaces[1]->GetVSize();
blockOffsets.PartialSum();
// Coefficients
mfem::ConstantCoefficient negOneCoeff(-1.0);
mfem::ConstantCoefficient oneCoeff(1.0);
mfem::Vector negOneVec(mfem::Vector(3));
mfem::Vector oneVec(mfem::Vector(3));
negOneVec = -1.0;
oneVec = 1.0;
mfem::VectorConstantCoefficient negOneVCoeff(negOneVec);
mfem::VectorConstantCoefficient oneVCoeff(oneVec);
// Add integrators to block form one by one
// We add integrators cooresponding to each term in the weak form
// The block form of the residual matrix
// ⎡ 0 -M ⎤ ⎡ θ ⎤ + ⎡f(θ)⎤ = ⎡ 0 ⎤ = R(X)
// ⎣ -Q D ⎦ ⎣ Φ ⎦ ⎣ 0 ⎦ ⎣ 0 ⎦
// This then simplifies to
// ⎡f(θ) - MΘ ⎤ = ⎡ 0 ⎤ = R
// ⎣ -Qɸ Dθ ⎦ ⎣ 0 ⎦
// Here M, Q, and D are
// M = ∫∇ψᶿ·Nᵠ dV (MixedVectorWeakDivergenceIntegrator)
// D = ∫ψᵠ·Nᵠ dV (VectorFEMassIntegrator)
// Q = ∫ψᵠ·∇Nᶿ dV (MixedVectorGradientIntegrator)
// f(θ) = ∫ψᶿ·θⁿ dV (NonlinearPowerIntegrator)
// Here ψᶿ and ψᵠ are the test functions for the theta and phi spaces, respectively
// Nᵠ and Nᶿ are the basis functions for the theta and phi spaces, respectively
// A full derivation of the weak form can be found in the 4DSSE documentation
// --- Assemble the MixedBilinear and Bilinear forms (M, D, and Q) ---
auto Mform = std::make_unique<mfem::MixedBilinearForm>(m_feTheta.get(), m_fePhi.get());
auto Qform = std::make_unique<mfem::MixedBilinearForm>(m_fePhi.get(), m_feTheta.get());
auto Dform = std::make_unique<mfem::BilinearForm>(m_fePhi.get());
// TODO: Check the sign on all of the integrators
Mform->AddDomainIntegrator(new mfem::MixedVectorWeakDivergenceIntegrator(negOneCoeff));
Qform->AddDomainIntegrator(new mfem::MixedVectorGradientIntegrator(negOneVCoeff));
Dform->AddDomainIntegrator(new mfem::VectorFEMassIntegrator(oneCoeff));
Mform->Assemble();
Mform->Finalize();
Qform->Assemble();
Qform->Finalize();
Dform->Assemble();
Dform->Finalize();
// --- Assemble the NonlinearForm (f) ---
auto fform = std::make_unique<mfem::NonlinearForm>(m_feTheta.get());
fform->AddDomainIntegrator(new polyMFEMUtils::NonlinearPowerIntegrator(oneCoeff, m_polytropicIndex));
// TODO: Add essential boundary conditions to the nonlinear form
// -- Build the BlockOperator --
m_polytropOperator = std::make_unique<PolytropeOperator>(
std::move(Mform),
std::move(Qform),
std::move(Dform),
std::move(fform),
blockOffsets
);
nonlinearForm->AddDomainIntegrator(compositeIntegrator.release());
}
void PolySolver::solve(){
// --- Set the initial guess for the solution ---
mfem::FunctionCoefficient initCoeff (
[this](const mfem::Vector &x) {
double r = x.Norml2();
// double theta = laneEmden::thetaSerieseExpansion(r, n, 10);
// return theta;
double radius = Probe::getMeshRadius(mesh);
double u = 1/radius;
setInitialGuess();
return -std::pow((u*r), 2)+1.0;
}
);
u->ProjectCoefficient(initCoeff);
if (config.get<bool>("Poly:Solver:ViewInitialGuess", false)) {
Probe::glVisView(*u, mesh, "initialGuess");
}
// mfem::DenseMatrix centerPoint(mesh.SpaceDimension(), 7);
mfem::DenseMatrix centerPoint(mesh.SpaceDimension(), 1);
centerPoint(0, 0) = 0.0;
centerPoint(1, 0) = 0.0;
centerPoint(2, 0) = 0.0;
// --- Set the essential true dofs for the operator ---
mfem::Array<int> theta_ess_tdof_list, phi_ess_tdof_list;
std::tie(theta_ess_tdof_list, phi_ess_tdof_list) = getEssentialTrueDof();
m_polytropOperator->SetEssentialTrueDofs(theta_ess_tdof_list, phi_ess_tdof_list);
mfem::Array<int> elementIDs;
mfem::Array<mfem::IntegrationPoint> ips;
mesh.FindPoints(centerPoint, elementIDs, ips);
mfem::Array<int> centerDofs;
mfem::Array<int> tempDofs;
for (int i = 0; i < elementIDs.Size(); i++) {
feSpace->GetElementDofs(elementIDs[i], tempDofs);
centerDofs.Append(tempDofs);
}
mfem::Array<int> ess_tdof_list;
mfem::Array<int> ess_brd(mesh.bdr_attributes.Max());
ess_brd = 1;
feSpace->GetEssentialTrueDofs(ess_brd, ess_tdof_list);
// combine the essential dofs with the center dofs
ess_tdof_list.Append(centerDofs);
nonlinearForm->SetEssentialTrueDofs(ess_tdof_list);
// Set the center elemID to be the Dirichlet boundary
// --- Load configuration parameters ---
double newtonRelTol = m_config.get<double>("Poly:Solver:Newton:RelTol", 1e-7);
double newtonAbsTol = m_config.get<double>("Poly:Solver:Newton:AbsTol", 1e-7);
int newtonMaxIter = m_config.get<int>("Poly:Solver:Newton:MaxIter", 200);
int newtonPrintLevel = m_config.get<int>("Poly:Solver:Newton:PrintLevel", 1);
// double alpha = config.get<double>("Poly:Solver:Newton:Alpha", 1e2);
double newtonRelTol = config.get<double>("Poly:Solver:Newton:RelTol", 1e-7);
double newtonAbsTol = config.get<double>("Poly:Solver:Newton:AbsTol", 1e-7);
int newtonMaxIter = config.get<int>("Poly:Solver:Newton:MaxIter", 200);
int newtonPrintLevel = config.get<int>("Poly:Solver:Newton:PrintLevel", 1);
double gmresRelTol = m_config.get<double>("Poly:Solver:GMRES:RelTol", 1e-10);
double gmresAbsTol = m_config.get<double>("Poly:Solver:GMRES:AbsTol", 1e-12);
int gmresMaxIter = m_config.get<int>("Poly:Solver:GMRES:MaxIter", 2000);
int gmresPrintLevel = m_config.get<int>("Poly:Solver:GMRES:PrintLevel", 0);
double gmresRelTol = config.get<double>("Poly:Solver:GMRES:RelTol", 1e-10);
double gmresAbsTol = config.get<double>("Poly:Solver:GMRES:AbsTol", 1e-12);
int gmresMaxIter = config.get<int>("Poly:Solver:GMRES:MaxIter", 2000);
int gmresPrintLevel = config.get<int>("Poly:Solver:GMRES:PrintLevel", 0);
LOG_DEBUG(m_logger, "Newton Solver (relTol: {:0.2E}, absTol: {:0.2E}, maxIter: {}, printLevel: {})", newtonRelTol, newtonAbsTol, newtonMaxIter, newtonPrintLevel);
LOG_DEBUG(m_logger, "GMRES Solver (relTol: {:0.2E}, absTol: {:0.2E}, maxIter: {}, printLevel: {})", gmresRelTol, gmresAbsTol, gmresMaxIter, gmresPrintLevel);
LOG_INFO(logger, "Newton Solver (relTol: {:0.2E}, absTol: {:0.2E}, maxIter: {}, printLevel: {})", newtonRelTol, newtonAbsTol, newtonMaxIter, newtonPrintLevel);
LOG_INFO(logger, "GMRES Solver (relTol: {:0.2E}, absTol: {:0.2E}, maxIter: {}, printLevel: {})", gmresRelTol, gmresAbsTol, gmresMaxIter, gmresPrintLevel);
std::vector<double> zeroSlopeCoordinate = {0.0, 0.0, 0.0};
// polyMFEMUtils::ZeroSlopeNewtonSolver newtonSolver(alpha, zeroSlopeCoordinate);
// --- Set up the Newton solver ---
mfem::NewtonSolver newtonSolver;
newtonSolver.SetRelTol(newtonRelTol);
newtonSolver.SetAbsTol(newtonAbsTol);
newtonSolver.SetMaxIter(newtonMaxIter);
newtonSolver.SetPrintLevel(newtonPrintLevel);
newtonSolver.SetOperator(*nonlinearForm);
newtonSolver.SetOperator(*m_polytropOperator);
mfem::GMRESSolver gmresSolver;
gmresSolver.SetRelTol(gmresRelTol);
gmresSolver.SetAbsTol(gmresAbsTol);
@@ -181,24 +221,96 @@ void PolySolver::solve(){
newtonSolver.SetSolver(gmresSolver);
// newtonSolver.SetAdaptiveLinRtol();
mfem::Vector B(feSpace->GetTrueVSize());
mfem::Vector B(m_feTheta->GetTrueVSize());
B = 0.0;
newtonSolver.Mult(B, *u);
newtonSolver.Mult(B, *m_theta);
Probe::glVisView(*u, mesh, "solution");
// --- Save and view the solution ---
saveAndViewSolution();
}
std::pair<mfem::Array<int>, mfem::Array<int>> PolySolver::getEssentialTrueDof() {
mfem::Array<int> theta_ess_tdof_list;
mfem::Array<int> phi_ess_tdof_list;
mfem::Array<int> centerDofs = findCenterElement();
phi_ess_tdof_list.Append(centerDofs);
mfem::Array<int> ess_brd(m_mesh->bdr_attributes.Max());
ess_brd = 1;
m_feTheta->GetEssentialTrueDofs(ess_brd, theta_ess_tdof_list);
// combine the essential dofs with the center dofs
theta_ess_tdof_list.Append(centerDofs);
return std::make_pair(theta_ess_tdof_list, phi_ess_tdof_list);
}
mfem::Array<int> PolySolver::findCenterElement() {
mfem::Array<int> centerDofs;
mfem::DenseMatrix centerPoint(m_mesh->SpaceDimension(), 1);
centerPoint(0, 0) = 0.0;
centerPoint(1, 0) = 0.0;
centerPoint(2, 0) = 0.0;
mfem::Array<int> elementIDs;
mfem::Array<mfem::IntegrationPoint> ips;
m_mesh->FindPoints(centerPoint, elementIDs, ips);
mfem::Array<int> tempDofs;
for (int i = 0; i < elementIDs.Size(); i++) {
m_feTheta->GetElementDofs(elementIDs[i], tempDofs);
centerDofs.Append(tempDofs);
}
return centerDofs;
}
void PolySolver::setInitialGuess() {
// --- Set the initial guess for the solution ---
mfem::FunctionCoefficient thetaInitGuess (
[this](const mfem::Vector &x) {
double r = x.Norml2();
double radius = Probe::getMeshRadius(*m_mesh);
double u = 1/radius;
return -std::pow((u*r), 2)+1.0;
}
);
mfem::VectorFunctionCoefficient phiInitGuess (m_mesh->SpaceDimension(),
[this](const mfem::Vector &x, mfem::Vector &v) {
double radius = Probe::getMeshRadius(*m_mesh);
double u = -1/std::pow(radius,2);
v(0) = 2*std::abs(x(0))*u;
v(1) = 2*std::abs(x(1))*u;
v(2) = 2*std::abs(x(2))*u;
}
);
m_theta->ProjectCoefficient(thetaInitGuess);
m_phi->ProjectCoefficient(phiInitGuess);
if (m_config.get<bool>("Poly:Solver:ViewInitialGuess", false)) {
Probe::glVisView(*m_theta, *m_mesh, "initialGuess");
}
}
void PolySolver::saveAndViewSolution() {
bool doView = m_config.get<bool>("Poly:Output:View", false);
if (doView) {
Probe::glVisView(*m_theta, *m_mesh, "solution");
}
// --- Extract the Solution ---
bool write11DSolution = config.get<bool>("Poly:Output:1D:Save", true);
bool write11DSolution = m_config.get<bool>("Poly:Output:1D:Save", true);
if (write11DSolution) {
std::string solutionPath = config.get<std::string>("Poly:Output:1D:Path", "polytropeSolution_1D.csv");
double rayCoLatitude = config.get<double>("Poly:Output:1D:RayCoLatitude", 0.0);
double rayLongitude = config.get<double>("Poly:Output:1D:RayLongitude", 0.0);
int raySamples = config.get<int>("Poly:Output:1D:RaySamples", 100);
std::string solutionPath = m_config.get<std::string>("Poly:Output:1D:Path", "polytropeSolution_1D.csv");
double rayCoLatitude = m_config.get<double>("Poly:Output:1D:RayCoLatitude", 0.0);
double rayLongitude = m_config.get<double>("Poly:Output:1D:RayLongitude", 0.0);
int raySamples = m_config.get<int>("Poly:Output:1D:RaySamples", 100);
std::vector rayDirection = {rayCoLatitude, rayLongitude};
Probe::getRaySolution(*u, *feSpace, rayDirection, raySamples, solutionPath);
Probe::getRaySolution(*m_theta, *m_feTheta, rayDirection, raySamples, solutionPath);
}
}