Skip to content

Allow integrating hypre AMS/ADS in FSPs #4412

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Draft
nmnobre wants to merge 3 commits into
base: devel
Choose a base branch
from
+53 −23
Draft

Allow integrating hypre AMS/ADS in FSPs #4412

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
fa0c190
Rename PC->Pc to avoid hiding typename (now mimics jac -> Jac)
nmnobre Feb 25, 2026
89899fe
Allow integrating hypre AMS/ADS in FSPs
nmnobre Feb 7, 2025
d01f48c
Prevent nested FSP as I dunno how to get matrix blocks
nmnobre Feb 25, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
44 changes: 36 additions & 8 deletions src/numerics/petsc_preconditioner.C
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -173,6 +173,10 @@ void PetscPreconditioner<T>::set_petsc_aux_data(PC & pc, System & sys, const uns
PCType pc_type = nullptr;
LibmeshPetscCallA(comm, PCGetType(pc, &pc_type));

// Get the nesting level of this PC's KSP
PetscInt level;
LibmeshPetscCallA(comm, PCGetKSPNestLevel(pc, &level));

// Check if hypre ams/ads, otherwise we quit with nothing to do
if (pc_type && std::string(pc_type) == PCHYPRE)
{
Expand All @@ -183,10 +187,6 @@ void PetscPreconditioner<T>::set_petsc_aux_data(PC & pc, System & sys, const uns
// If not ams/ads, we quit with nothing to do
if (std::string(hypre_type) == "ams")
{
// If multiple variables, we error out as senseless
libmesh_error_msg_if(sys.n_vars() > 1,
"Error applying hypre AMS to a system with multiple "
"variables");
// If not a 1st order Nédélec or a 2d 1st order Raviart-Thomas system, we
// error out as we do not support anything else at the moment
libmesh_error_msg_if(sys.variable(v).type() != FEType(1, NEDELEC_ONE) &&
Expand All @@ -199,10 +199,6 @@ void PetscPreconditioner<T>::set_petsc_aux_data(PC & pc, System & sys, const uns
}
else if (std::string(hypre_type) == "ads")
{
// If multiple variables, we error out as senseless
libmesh_error_msg_if(sys.n_vars() > 1,
"Error applying hypre ADS to a system with multiple "
"variables");
// If not a 3d 1st order Raviart-Thomas system, we error out as we do not
// support anything else at the moment
libmesh_error_msg_if(sys.variable(v).type() != FEType(1, RAVIART_THOMAS) ||
Expand All @@ -213,6 +209,38 @@ void PetscPreconditioner<T>::set_petsc_aux_data(PC & pc, System & sys, const uns
set_hypre_ads_data(pc, sys, v);
}
}
else if (pc_type && std::string(pc_type) == PCFIELDSPLIT && !level)
{
// Annoyingly, if using Schur complement preconditioning, we need to call PCSetUp()
auto pmatrix = cast_ptr<PetscMatrixBase<T> *>(sys.request_matrix("System Matrix"));
pmatrix->close();
Mat mat = pmatrix->mat();
LibmeshPetscCallA(comm, PCSetOperators(pc, mat, mat));
LibmeshPetscCallA(comm, PCSetUp(pc));

// Get sub KSP contexts for all splits
PetscInt n_splits;
KSP * subksps;
LibmeshPetscCallA(comm, PCFieldSplitGetSubKSP(pc, &n_splits, &subksps));

// We assume a one-to-one map between splits and variables
if (sys.n_vars() != n_splits)
return;

// Get the sub PC context for each split and recursively call this function
for (auto s : make_range(n_splits))
{
PC subpc;
LibmeshPetscCallA(comm, KSPGetPC(subksps[s], &subpc));
#if PETSC_VERSION_LESS_THAN(3, 25, 0)
LibmeshPetscCallA(comm, PCSetKSPNestLevel(subpc, level + 1));
#endif
set_petsc_aux_data(subpc, sys, s);
}

// Free the array of sub KSP contexts
LibmeshPetscCallA(comm, PetscFree(subksps));
}
}


Expand Down
32 changes: 17 additions & 15 deletions src/solvers/petsc_nonlinear_solver.C
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -432,7 +432,7 @@ extern "C"

NonlinearImplicitSystem & sys = solver->system();

PetscMatrixBase<Number> * const PC = pc ? PetscMatrixBase<Number>::get_context(pc, sys.comm()) : nullptr;
PetscMatrixBase<Number> * const Pc = pc ? PetscMatrixBase<Number>::get_context(pc, sys.comm()) : nullptr;
PetscMatrixBase<Number> * Jac = jac ? PetscMatrixBase<Number>::get_context(jac, sys.comm()) : nullptr;
PetscVector<Number> & X_sys = *cast_ptr<PetscVector<Number> *>(sys.solution.get());
PetscVector<Number> X_global(x, sys.comm());
Expand Down Expand Up @@ -462,13 +462,13 @@ extern "C"
Jac->close();

if (pc && (p_is_shell == PETSC_TRUE))
PC->close();
Pc->close();

PetscFunctionReturn(LIBMESH_PETSC_SUCCESS);
}

// Set the dof maps
PC->attach_dof_map(sys.get_dof_map());
Pc->attach_dof_map(sys.get_dof_map());
Jac->attach_dof_map(sys.get_dof_map());

// Use the systems update() to get a good local version of the parallel solution
Expand All @@ -484,36 +484,42 @@ extern "C"
sys.get_dof_map().enforce_constraints_exactly(sys, sys.current_local_solution.get());

if (solver->_zero_out_jacobian)
PC->zero();
Pc->zero();


if (solver->jacobian != nullptr)
solver->jacobian(*sys.current_local_solution.get(), *PC, sys);
solver->jacobian(*sys.current_local_solution.get(), *Pc, sys);

else if (solver->jacobian_object != nullptr)
solver->jacobian_object->jacobian(*sys.current_local_solution.get(), *PC, sys);
solver->jacobian_object->jacobian(*sys.current_local_solution.get(), *Pc, sys);

else if (solver->matvec != nullptr)
solver->matvec(*sys.current_local_solution.get(), nullptr, PC, sys);
solver->matvec(*sys.current_local_solution.get(), nullptr, Pc, sys);

else
libmesh_error_msg("Error! Unable to compute residual and/or Jacobian!");

PC->close();
Pc->close();
if (solver->_exact_constraint_enforcement)
{
sys.get_dof_map().enforce_constraints_on_jacobian(sys, PC);
PC->close();
sys.get_dof_map().enforce_constraints_on_jacobian(sys, Pc);
Pc->close();
}

if (Jac != PC)
if (Jac != Pc)
{
// Assume that shells know what they're doing
libmesh_assert(!solver->_exact_constraint_enforcement || (j_is_mffd == PETSC_TRUE) ||
(j_is_shell == PETSC_TRUE));
Jac->close();
}

KSP kspc;
PC pcc;
LibmeshPetscCall2(sys.comm(), SNESGetKSP(snes, &kspc));
LibmeshPetscCall2(sys.comm(), KSPGetPC(kspc, &pcc));
PetscPreconditioner<Number>::set_petsc_aux_data(pcc, sys);

PetscFunctionReturn(LIBMESH_PETSC_SUCCESS);
}

Expand Down Expand Up @@ -985,10 +991,6 @@ PetscNonlinearSolver<T>::solve (SparseMatrix<T> & pre_in, // System Preconditi
#endif
LibmeshPetscCall(SNESSetFromOptions(_snes));

PC pc;
LibmeshPetscCall(KSPGetPC(ksp, &pc));
PetscPreconditioner<T>::set_petsc_aux_data(pc, this->system());

#if defined(LIBMESH_HAVE_PETSC_HYPRE) && PETSC_VERSION_LESS_THAN(3, 23, 0) && \
!PETSC_VERSION_LESS_THAN(3, 12, 0) && defined(PETSC_HAVE_HYPRE_DEVICE)
{
Expand Down