Skip to content

Add in riesz map and remove pc update #713

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

Merged
jshipton merged 13 commits into from
Mar 24, 2026
Merged

Add in riesz map and remove pc update #713

Show file tree
Hide file tree
Changes from 5 commits
Commits
Show all changes
13 commits
Select commit Hold shift + click to select a range
335d6bf
Add in riesz map and remove pc update
atb1995 Mar 12, 2026
fbd42ea
Intermediate commit - need to pass solver stuff down
atb1995 Mar 18, 2026
baa395a
Passing through of solver settings
atb1995 Mar 19, 2026
8cdfab1
Merge branch 'main' into compressible_fix
atb1995 Mar 19, 2026
30d8284
Fix of settings unpack
atb1995 Mar 19, 2026
e29a439
Update gusto/solvers/preconditioners.py
atb1995 Mar 19, 2026
0caafbf
Tidy up of options ready in using JHC suggestions
atb1995 Mar 19, 2026
3ee0a5e
Update kgos
atb1995 Mar 20, 2026
7f2bf74
Merge branch 'main' into compressible_fix
atb1995 Mar 20, 2026
aa884e8
Updating kgos
atb1995 Mar 20, 2026
85d18c6
Remove a tol
atb1995 Mar 24, 2026
2f9349a
Updating kgos
atb1995 Mar 24, 2026
765f869
Merge branch 'main' into compressible_fix
atb1995 Mar 24, 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
66 changes: 28 additions & 38 deletions gusto/solvers/preconditioners.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -2,7 +2,7 @@

from firedrake import (dot, jump, dS_h, ds_b, ds_t, ds,
FacetNormal, Tensor, AssembledVector,
AuxiliaryOperatorPC, PETSc)
AuxiliaryOperatorPC, PETSc, RieszMap)

from firedrake.preconditioners import PCBase
from firedrake.matrix_free.operators import ImplicitMatrixContext
Expand Down Expand Up @@ -460,35 +460,6 @@ class CompressibleHybridisedSCPC(PCBase):

_prefix = "compressible_hybrid_scpc"

scpc_parameters = {
'mat_type': 'matfree',
'ksp_type': 'preonly',
'pc_type': 'python',
'pc_python_type': 'firedrake.SCPC',
'pc_sc_eliminate_fields': '0, 1',
# The reduced operator is not symmetric
'condensed_field': {
'ksp_type': 'fgmres',
'ksp_rtol': 1.0e-8,
'ksp_atol': 1.0e-8,
'ksp_max_it': 100,
'pc_type': 'gamg',
'pc_gamg_sym_graph': None,
'mg_levels': {
'ksp_type': 'gmres',
'ksp_max_it': 5,
'pc_type': 'bjacobi',
'sub_pc_type': 'ilu'
}
}
}

cg_ilu_parameters = {
'ksp_type': 'cg',
'pc_type': 'bjacobi',
'sub_pc_type': 'ilu'
}

def initialize(self, pc):
"""
Set up problem and solver
Expand All @@ -507,15 +478,19 @@ def initialize(self, pc):
raise ValueError("Expecting PC type python")

self._process_context(pc)
prefix = pc.getOptionsPrefix()
opts = PETSc.Options()

self.hybridised_scpc_parameters = self.scpc_parameters
self.rho_avg_solver_parameters = self.cg_ilu_parameters
self.theta_solver_parameters = self.cg_ilu_parameters
self.exner_avg_solver_parameters = self.cg_ilu_parameters
# Unpack sub-solver parameters from the options prefix
self.theta_solver_parameters = self._get_sub_params(opts, prefix + "theta_backsub")
self.exner_avg_solver_parameters = self._get_sub_params(opts, prefix + "exner_ave")
self.rho_avg_solver_parameters = self._get_sub_params(opts, prefix + "rho_ave")
self.riesz_map_parameters = self._get_sub_params(opts, prefix + "riesz_map")
self.scpc_solve_parameters = self._get_sub_params(opts, prefix + "scpc_solve")

@JHopeCollins JHopeCollins Mar 19, 2026

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This isn't necessary, just give each solver the right prefix and they will grab these options themselves.

@atb1995 atb1995 Mar 19, 2026

Copy link
Copy Markdown
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

So something like:

      self.exner_avg_solver = LinearVariationalSolver(
            options_prefix=pc.getOptionsPrefix()+'exner_ave'
        )

What about RieszMap which has no options_prefix as an argument?

@JHopeCollins JHopeCollins Mar 19, 2026 •
edited
Loading

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Arguably that is an oversight. Nevertheless, you can ask PETSc to get all the options with a particular prefix:

prefixed_options = PETSc.Options("prefix_")  # This is another Options object
options_dict = PETSc.Options("prefix_").getAll()  # this is a python dict for solver_parameters

If PETSc.Options() contains "prefix_obj_option" then PETSc.Options("prefix_") will contain "obj_option"`.

@JHopeCollins JHopeCollins Mar 19, 2026

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

But yes for the other solvers that is exactly what I meant.


if logger.isEnabledFor(DEBUG):
self.hybridised_scpc_parameters['ksp_monitor_true_residual'] = None
self.hybridised_scpc_parameters['ksp_converged_reason'] = None
self.scpc_solve_parameters['ksp_monitor_true_residual'] = None
self.scpc_solve_parameters['ksp_converged_reason'] = None

# Equations and parameters
equations = self.equations
Expand Down Expand Up @@ -548,6 +523,11 @@ def initialize(self, pc):
self.xrhs = Function(self.W)
self.y = Function(self.W)

# Riesz map for the dual of the original function space
self.riesz_map = RieszMap(self.W.dual(),
solver_parameters=self.riesz_map_parameters,
constant_jacobian=True)

self.y_hybrid = Function(self.W_hyb)

# Split input functions
Expand Down Expand Up @@ -691,7 +671,7 @@ def L_tr(f):
aeqn, Leqn, self.y_hybrid, constant_jacobian=True
)
self.hybridized_solver = LinearVariationalSolver(
hybridized_prb, solver_parameters=self.hybridised_scpc_parameters,
hybridized_prb, solver_parameters=self.scpc_solve_parameters,
options_prefix=pc.getOptionsPrefix()+self._prefix, appctx=appctx
)

Expand Down Expand Up @@ -753,7 +733,7 @@ def apply(self, pc, x, y):
with self.xstar.dat.vec_wo as xv:
x.copy(xv)

self.xrhs.assign(self.xstar.riesz_representation())
self.xrhs.assign(self.riesz_map.__call__(self.xstar))
Comment thread
atb1995 marked this conversation as resolved.
Outdated
# Solve hybridized system
self.hybridized_solver.solve()

Expand Down Expand Up @@ -813,3 +793,13 @@ def applyTranspose(self, pc, x, y):
"""

raise NotImplementedError("The transpose application of the PC is not implemented.")

def _get_sub_params(self, opts, prefix):
result = {}
for key, val in opts.getAll().items():
if key.startswith(prefix):
stripped = key[len(prefix):]
# Remove leading underscore
stripped = stripped.lstrip('_')
result[stripped] = val if val != "" else None
return result
63 changes: 63 additions & 0 deletions gusto/solvers/solver_presets.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -67,12 +67,75 @@ def trace_nullsp(T):
# Compressible Euler equations - (u, rho, theta) system. We use a
# bespoke hybridization preconditioner for this system owing to the
# nonlinear pressure gradient term.
r_tol = 1e-8
a_tol = 1e-8

@JHopeCollins JHopeCollins Mar 19, 2026

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The right atol depends on the dimension of the problem, so hardcoding it may lead to some unexpected behaviour if someone tries a problem with particularly different scales.

@atb1995 atb1995 Mar 19, 2026

Copy link
Copy Markdown
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

You would suggest just not setting it at all?

@JHopeCollins JHopeCollins Mar 19, 2026 •
edited
Loading

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Personally no, but that's it is opinion, I was just pointing it out. If atol was part of the options before then it isn't up to me to decide gusto policy!
If the majority of use cases have scales O(earth) then it isn't unreasonable, it just needs bearing in mind.

@jshipton jshipton Mar 23, 2026

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

If both rtol and atol are set, does it stop as soon as either are met? I wonder which one we're hitting...

@JHopeCollins JHopeCollins Mar 23, 2026 •
edited
Loading

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yes, it will stop as soon as any of atol, rtol, divtol, stol, or max_it are hit. If you are hitting divtol or stol you are definitely in trouble but it tends to be obvious.

The message from either ksp_converged_reason, ksp_converged_rate, or snes_converged_reason will say what convergence criteria was met.

atol might mean you stop prematurely if you have something with small physical scales that hasn't been nondimensionalised.

@jshipton jshipton Mar 23, 2026

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@atb1995 @tommbendall do either of you have any logs that contain the convergence reasons? In general, I agree with Josh that we should not be setting atol but I am also wary of changing too much in one go!

@JHopeCollins JHopeCollins Mar 23, 2026

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I imagine it is unlikely for you to be hitting atol for anything on the globe. But maybe for something local or nondimensionalised if the initial residual is O(1) or smaller.

@atb1995 atb1995 Mar 23, 2026

Copy link
Copy Markdown
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

It is always hitting r_tol as far as I am aware, we are getting nowhere near 1e-8 atol for a lot of the runs I am doing at least

@jshipton jshipton Mar 23, 2026

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

let's get rid then


theta_backsub_settings = {
'ksp_type': 'cg',
'pc_type': 'bjacobi',
'sub_pc_type': 'ilu',
'ksp_rtol': r_tol,
'ksp_atol': a_tol,
}

exner_ave_settings = {
'ksp_type': 'cg',
'pc_type': 'bjacobi',
'sub_pc_type': 'ilu',
'ksp_rtol': r_tol,
'ksp_atol': a_tol,
}

rho_ave_settings = {
'ksp_type': 'cg',
'pc_type': 'bjacobi',
'sub_pc_type': 'ilu',
'ksp_rtol': r_tol,
'ksp_atol': a_tol,
}

riesz_map_settings = {
'ksp_type': 'cg',
'pc_type': 'bjacobi',
'sub_pc_type': 'ilu',
'ksp_rtol': r_tol,
'ksp_atol': a_tol,
}

scpc_solve_settings = {
'mat_type': 'matfree',
'ksp_type': 'preonly',
'pc_type': 'python',
'pc_python_type': 'firedrake.SCPC',
'pc_sc_eliminate_fields': '0,1',
# The reduced operator is not symmetric
'condensed_field': {
'ksp_type': 'fgmres',
'ksp_rtol': r_tol,
'ksp_atol': a_tol,
'ksp_max_it': 100,
'pc_type': 'gamg',
'pc_gamg_sym_graph': None,
'mg_levels': {
'ksp_type': 'gmres',
'ksp_max_it': 5,
'pc_type': 'bjacobi',
'sub_pc_type': 'ilu'
}
}
}

settings = {
'ksp_monitor': None,
'ksp_type': 'preonly',
'mat_type': 'matfree',
'pc_type': 'python',
'pc_python_type': 'gusto.CompressibleHybridisedSCPC',
'theta_backsub': theta_backsub_settings,
'exner_ave': exner_ave_settings,
'rho_ave': rho_ave_settings,
'riesz_map': riesz_map_settings,
'scpc_solve': scpc_solve_settings,
}

# We pass the implicit weighting parameter (alpha) and tau_values to the
Expand Down
7 changes: 0 additions & 7 deletions gusto/timestepping/semi_implicit_quasi_newton.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -20,7 +20,6 @@
from gusto.time_discretisation.time_discretisation import ExplicitTimeDiscretisation
from gusto.timestepping.timestepper import BaseTimestepper
from gusto.solvers.solver_presets import hybridised_solver_parameters
from gusto.equations.compressible_euler_equations import CompressibleEulerEquations


__all__ = ["SemiImplicitQuasiNewton", "Forcing", "QuasiNewtonLinearSolver"]
Expand Down Expand Up @@ -937,12 +936,6 @@ def update_reference_profiles(self):
self.equation.update_reference_profiles()
self.solver.invalidate_jacobian()

# TODO: Issue #686 is to address this reference profile update bug (pythonPC update not called)
# this line forces it to update for now
pc = self.solver.snes.getKSP().getPC()
if (isinstance(self.equation, CompressibleEulerEquations) and pc.getType() == "python"):
pc.getPythonContext().update(pc)

def zero_non_prognostics(self, equation, xrhs, field_names, prognostic_names):
"""
Zero rhs term F(x) for non-prognostics.
Expand Down
Loading