From 1305ce3ba9a8254e31b5d849e8529e7baf91cb93 Mon Sep 17 00:00:00 2001 From: Tyagi Date: Thu, 19 Mar 2026 13:04:10 +1100 Subject: [PATCH 1/4] Fix Stokes velocity FE dual-space setup --- src/underworld3/cython/petsc_generic_snes_solvers.pyx | 2 ++ 1 file changed, 2 insertions(+) diff --git a/src/underworld3/cython/petsc_generic_snes_solvers.pyx b/src/underworld3/cython/petsc_generic_snes_solvers.pyx index 0c052a22d..f7df4b77d 100644 --- a/src/underworld3/cython/petsc_generic_snes_solvers.pyx +++ b/src/underworld3/cython/petsc_generic_snes_solvers.pyx @@ -3405,6 +3405,8 @@ class SNES_Stokes_SaddlePt(SolverBaseClass): options = PETSc.Options() options.setValue("private_{}_u_petscspace_degree".format(self.petsc_options_prefix), u_degree) # for private variables + options.setValue("private_{}_u_petscdualspace_lagrange_continuity".format(self.petsc_options_prefix), self.u.continuous) + options.setValue("private_{}_u_petscdualspace_lagrange_node_endpoints".format(self.petsc_options_prefix), False) self.petsc_fe_u = PETSc.FE().createDefault(mesh.dim, mesh.dim, mesh.isSimplex, mesh.qdegree, "private_{}_u_".format(self.petsc_options_prefix), PETSc.COMM_SELF) self.petsc_fe_u.setName("velocity") self.petsc_fe_u_id = self.dm.getNumFields() From b5e7858afef7976c702506f122b75cdf31dd0c92 Mon Sep 17 00:00:00 2001 From: Tyagi Date: Thu, 19 Mar 2026 13:07:04 +1100 Subject: [PATCH 2/4] Add box MMS reproducer for high-order Stokes --- examples/stokes_box_mms_simplex.py | 156 +++++++++++++++++++++++++++++ 1 file changed, 156 insertions(+) create mode 100644 examples/stokes_box_mms_simplex.py diff --git a/examples/stokes_box_mms_simplex.py b/examples/stokes_box_mms_simplex.py new file mode 100644 index 000000000..e18d3ea55 --- /dev/null +++ b/examples/stokes_box_mms_simplex.py @@ -0,0 +1,156 @@ +""" +Simple manufactured-solution Stokes example on a triangular box. + +This script is useful for checking higher-order simplex velocity / pressure +pairs on a straight-edged domain without curved-geometry effects. + +Default setup: + python examples/stokes_box_mms_simplex.py + +Useful comparison: + python examples/stokes_box_mms_simplex.py --vdegree 2 --pdegree 1 + python examples/stokes_box_mms_simplex.py --vdegree 3 --pdegree 2 +""" + +import argparse +import os + +os.environ.setdefault("OPENBLAS_NUM_THREADS", "1") +os.environ.setdefault("OMP_NUM_THREADS", "1") +os.environ.setdefault("MKL_NUM_THREADS", "1") +os.environ.setdefault("VECLIB_MAXIMUM_THREADS", "1") +os.environ["SYMPY_USE_CACHE"] = "no" + +import numpy as np +import sympy as sp + +import underworld3 as uw +from underworld3.systems import Stokes + + +def parse_args(): + parser = argparse.ArgumentParser() + parser.add_argument("--cellsize", type=float, default=0.125) + parser.add_argument("--vdegree", type=int, default=3) + parser.add_argument("--pdegree", type=int, default=2) + parser.add_argument("--pcont", action="store_true", default=True) + parser.add_argument("--stokes-tol", type=float, default=1.0e-10) + return parser.parse_args() + + +def subtract_pressure_mean(mesh, pressure_var): + p_int = uw.maths.Integral(mesh, pressure_var.sym[0]).evaluate() + volume = uw.maths.Integral(mesh, 1.0).evaluate() + pressure_var.data[:, 0] -= p_int / volume + + +def relative_l2_error(mesh, err_expr, ana_expr): + if isinstance(err_expr, sp.MatrixBase): + err_expr = err_expr.dot(err_expr) + ana_expr = ana_expr.dot(ana_expr) + else: + err_expr = err_expr * err_expr + ana_expr = ana_expr * ana_expr + + err_I = uw.maths.Integral(mesh, err_expr) + ana_I = uw.maths.Integral(mesh, ana_expr) + + return np.sqrt(err_I.evaluate()) / np.sqrt(ana_I.evaluate()) + + +def main(): + args = parse_args() + qdegree = max(2 * args.vdegree, args.vdegree + args.pdegree) + + mesh = uw.meshing.UnstructuredSimplexBox( + minCoords=(0.0, 0.0), + maxCoords=(1.0, 1.0), + cellSize=args.cellsize, + regular=True, + qdegree=qdegree, + filename=f"/tmp/stokes_box_mms_simplex_{args.vdegree}_{args.pdegree}_{args.cellsize}.msh", + ) + + x, y = mesh.X + + psi = x**2 * y**2 + v_ana_expr = sp.Matrix([sp.diff(psi, y), -sp.diff(psi, x)]) + p_ana_expr = x**2 - sp.Rational(1, 3) + + bodyforce = sp.Matrix( + [ + -(sp.diff(v_ana_expr[0], x, 2) + sp.diff(v_ana_expr[0], y, 2)) + + sp.diff(p_ana_expr, x), + -(sp.diff(v_ana_expr[1], x, 2) + sp.diff(v_ana_expr[1], y, 2)) + + sp.diff(p_ana_expr, y), + ] + ) + + v_soln = uw.discretisation.MeshVariable( + varname="Velocity", + mesh=mesh, + degree=args.vdegree, + vtype=uw.VarType.VECTOR, + ) + + p_soln = uw.discretisation.MeshVariable( + varname="Pressure", + mesh=mesh, + degree=args.pdegree, + vtype=uw.VarType.SCALAR, + continuous=args.pcont, + ) + + stokes = Stokes(mesh, velocityField=v_soln, pressureField=p_soln) + stokes.constitutive_model = uw.constitutive_models.ViscousFlowModel + stokes.constitutive_model.Parameters.viscosity = 1.0 + stokes.saddle_preconditioner = 1.0 + stokes.bodyforce = bodyforce + + for boundary_name in ( + mesh.boundaries.Bottom.name, + mesh.boundaries.Top.name, + mesh.boundaries.Left.name, + mesh.boundaries.Right.name, + ): + stokes.add_essential_bc(v_ana_expr, boundary_name) + + stokes.tolerance = args.stokes_tol + stokes.petsc_options["snes_type"] = "ksponly" + stokes.petsc_options["ksp_type"] = "fgmres" + stokes.petsc_options["ksp_rtol"] = args.stokes_tol + stokes.petsc_options["ksp_atol"] = 0.0 + stokes.petsc_options["ksp_monitor"] = None + stokes.petsc_options["ksp_monitor_true_residual"] = None + stokes.petsc_options["ksp_converged_reason"] = None + + stokes.petsc_options.setValue("fieldsplit_velocity_pc_mg_type", "kaskade") + stokes.petsc_options.setValue("fieldsplit_velocity_pc_mg_cycle_type", "w") + stokes.petsc_options["fieldsplit_velocity_mg_coarse_pc_type"] = "svd" + stokes.petsc_options["fieldsplit_velocity_ksp_type"] = "fcg" + stokes.petsc_options["fieldsplit_velocity_mg_levels_ksp_type"] = "chebyshev" + stokes.petsc_options["fieldsplit_velocity_mg_levels_ksp_max_it"] = 5 + stokes.petsc_options["fieldsplit_velocity_mg_levels_ksp_converged_maxits"] = None + stokes.petsc_options.setValue("fieldsplit_pressure_pc_type", "mg") + stokes.petsc_options.setValue("fieldsplit_pressure_pc_mg_type", "multiplicative") + stokes.petsc_options.setValue("fieldsplit_pressure_pc_mg_cycle_type", "v") + + stokes.solve(verbose=False, debug=False) + + subtract_pressure_mean(mesh, p_soln) + + v_err_expr = sp.Matrix(v_soln.sym).T - v_ana_expr + p_err_expr = p_soln.sym[0] - p_ana_expr + + v_err_l2 = relative_l2_error(mesh, v_err_expr, v_ana_expr) + p_err_l2 = relative_l2_error(mesh, p_err_expr, p_ana_expr) + + uw.pprint("cellsize:", args.cellsize) + uw.pprint("vdegree:", args.vdegree) + uw.pprint("pdegree:", args.pdegree) + uw.pprint("Relative velocity L2 error:", v_err_l2) + uw.pprint("Relative pressure L2 error:", p_err_l2) + + +if __name__ == "__main__": + main() From 684e89983aebbe714fbebbdd2b9bf7150aef95bb Mon Sep 17 00:00:00 2001 From: Tyagi Date: Thu, 19 Mar 2026 13:14:00 +1100 Subject: [PATCH 3/4] Convert box MMS reproducer to notebook style --- examples/stokes_box_mms_simplex.py | 294 ++++++++++++++++------------- 1 file changed, 165 insertions(+), 129 deletions(-) diff --git a/examples/stokes_box_mms_simplex.py b/examples/stokes_box_mms_simplex.py index e18d3ea55..3fcb0c306 100644 --- a/examples/stokes_box_mms_simplex.py +++ b/examples/stokes_box_mms_simplex.py @@ -1,7 +1,7 @@ """ -Simple manufactured-solution Stokes example on a triangular box. +Notebook-style manufactured-solution Stokes example on a triangular box. -This script is useful for checking higher-order simplex velocity / pressure +This example is useful for checking higher-order simplex velocity / pressure pairs on a straight-edged domain without curved-geometry effects. Default setup: @@ -27,130 +27,166 @@ import underworld3 as uw from underworld3.systems import Stokes - -def parse_args(): - parser = argparse.ArgumentParser() - parser.add_argument("--cellsize", type=float, default=0.125) - parser.add_argument("--vdegree", type=int, default=3) - parser.add_argument("--pdegree", type=int, default=2) - parser.add_argument("--pcont", action="store_true", default=True) - parser.add_argument("--stokes-tol", type=float, default=1.0e-10) - return parser.parse_args() - - -def subtract_pressure_mean(mesh, pressure_var): - p_int = uw.maths.Integral(mesh, pressure_var.sym[0]).evaluate() - volume = uw.maths.Integral(mesh, 1.0).evaluate() - pressure_var.data[:, 0] -= p_int / volume - - -def relative_l2_error(mesh, err_expr, ana_expr): - if isinstance(err_expr, sp.MatrixBase): - err_expr = err_expr.dot(err_expr) - ana_expr = ana_expr.dot(ana_expr) - else: - err_expr = err_expr * err_expr - ana_expr = ana_expr * ana_expr - - err_I = uw.maths.Integral(mesh, err_expr) - ana_I = uw.maths.Integral(mesh, ana_expr) - - return np.sqrt(err_I.evaluate()) / np.sqrt(ana_I.evaluate()) - - -def main(): - args = parse_args() - qdegree = max(2 * args.vdegree, args.vdegree + args.pdegree) - - mesh = uw.meshing.UnstructuredSimplexBox( - minCoords=(0.0, 0.0), - maxCoords=(1.0, 1.0), - cellSize=args.cellsize, - regular=True, - qdegree=qdegree, - filename=f"/tmp/stokes_box_mms_simplex_{args.vdegree}_{args.pdegree}_{args.cellsize}.msh", - ) - - x, y = mesh.X - - psi = x**2 * y**2 - v_ana_expr = sp.Matrix([sp.diff(psi, y), -sp.diff(psi, x)]) - p_ana_expr = x**2 - sp.Rational(1, 3) - - bodyforce = sp.Matrix( - [ - -(sp.diff(v_ana_expr[0], x, 2) + sp.diff(v_ana_expr[0], y, 2)) - + sp.diff(p_ana_expr, x), - -(sp.diff(v_ana_expr[1], x, 2) + sp.diff(v_ana_expr[1], y, 2)) - + sp.diff(p_ana_expr, y), - ] - ) - - v_soln = uw.discretisation.MeshVariable( - varname="Velocity", - mesh=mesh, - degree=args.vdegree, - vtype=uw.VarType.VECTOR, - ) - - p_soln = uw.discretisation.MeshVariable( - varname="Pressure", - mesh=mesh, - degree=args.pdegree, - vtype=uw.VarType.SCALAR, - continuous=args.pcont, - ) - - stokes = Stokes(mesh, velocityField=v_soln, pressureField=p_soln) - stokes.constitutive_model = uw.constitutive_models.ViscousFlowModel - stokes.constitutive_model.Parameters.viscosity = 1.0 - stokes.saddle_preconditioner = 1.0 - stokes.bodyforce = bodyforce - - for boundary_name in ( - mesh.boundaries.Bottom.name, - mesh.boundaries.Top.name, - mesh.boundaries.Left.name, - mesh.boundaries.Right.name, - ): - stokes.add_essential_bc(v_ana_expr, boundary_name) - - stokes.tolerance = args.stokes_tol - stokes.petsc_options["snes_type"] = "ksponly" - stokes.petsc_options["ksp_type"] = "fgmres" - stokes.petsc_options["ksp_rtol"] = args.stokes_tol - stokes.petsc_options["ksp_atol"] = 0.0 - stokes.petsc_options["ksp_monitor"] = None - stokes.petsc_options["ksp_monitor_true_residual"] = None - stokes.petsc_options["ksp_converged_reason"] = None - - stokes.petsc_options.setValue("fieldsplit_velocity_pc_mg_type", "kaskade") - stokes.petsc_options.setValue("fieldsplit_velocity_pc_mg_cycle_type", "w") - stokes.petsc_options["fieldsplit_velocity_mg_coarse_pc_type"] = "svd" - stokes.petsc_options["fieldsplit_velocity_ksp_type"] = "fcg" - stokes.petsc_options["fieldsplit_velocity_mg_levels_ksp_type"] = "chebyshev" - stokes.petsc_options["fieldsplit_velocity_mg_levels_ksp_max_it"] = 5 - stokes.petsc_options["fieldsplit_velocity_mg_levels_ksp_converged_maxits"] = None - stokes.petsc_options.setValue("fieldsplit_pressure_pc_type", "mg") - stokes.petsc_options.setValue("fieldsplit_pressure_pc_mg_type", "multiplicative") - stokes.petsc_options.setValue("fieldsplit_pressure_pc_mg_cycle_type", "v") - - stokes.solve(verbose=False, debug=False) - - subtract_pressure_mean(mesh, p_soln) - - v_err_expr = sp.Matrix(v_soln.sym).T - v_ana_expr - p_err_expr = p_soln.sym[0] - p_ana_expr - - v_err_l2 = relative_l2_error(mesh, v_err_expr, v_ana_expr) - p_err_l2 = relative_l2_error(mesh, p_err_expr, p_ana_expr) - - uw.pprint("cellsize:", args.cellsize) - uw.pprint("vdegree:", args.vdegree) - uw.pprint("pdegree:", args.pdegree) - uw.pprint("Relative velocity L2 error:", v_err_l2) - uw.pprint("Relative pressure L2 error:", p_err_l2) - - -if __name__ == "__main__": - main() +# %% [markdown] +# # Stokes Box MMS on a Simplex Mesh +# +# This notebook-style example solves a manufactured Stokes problem on the unit +# square using triangular elements. It is intentionally simple and useful for +# checking higher-order simplex pairs such as `P2/P1` and `P3/P2`. + +# %% [markdown] +# ## Runtime parameters +# +# The script still accepts command-line arguments so it can be run directly. + +# %% +parser = argparse.ArgumentParser() +parser.add_argument("--cellsize", type=float, default=0.125) +parser.add_argument("--vdegree", type=int, default=3) +parser.add_argument("--pdegree", type=int, default=2) +parser.add_argument("--pcont", action="store_true", default=True) +parser.add_argument("--stokes-tol", type=float, default=1.0e-10) +args = parser.parse_args() + +qdegree = max(2 * args.vdegree, args.vdegree + args.pdegree) + +# %% [markdown] +# ## Mesh + +# %% +mesh = uw.meshing.UnstructuredSimplexBox( + minCoords=(0.0, 0.0), + maxCoords=(1.0, 1.0), + cellSize=args.cellsize, + regular=True, + qdegree=qdegree, + filename=f"/tmp/stokes_box_mms_simplex_{args.vdegree}_{args.pdegree}_{args.cellsize}.msh", +) + +x, y = mesh.X + +# %% [markdown] +# ## Manufactured solution +# +# Stream function: +# `psi = x^2 y^2` +# +# Exact velocity: +# `u = (dpsi/dy, -dpsi/dx)` +# +# Exact pressure: +# `p = x^2 - 1/3` + +# %% +psi = x**2 * y**2 +v_ana_expr = sp.Matrix([sp.diff(psi, y), -sp.diff(psi, x)]) +p_ana_expr = x**2 - sp.Rational(1, 3) + +bodyforce = sp.Matrix( + [ + -(sp.diff(v_ana_expr[0], x, 2) + sp.diff(v_ana_expr[0], y, 2)) + + sp.diff(p_ana_expr, x), + -(sp.diff(v_ana_expr[1], x, 2) + sp.diff(v_ana_expr[1], y, 2)) + + sp.diff(p_ana_expr, y), + ] +) + +# %% [markdown] +# ## Discretisation + +# %% +v_soln = uw.discretisation.MeshVariable( + varname="Velocity", + mesh=mesh, + degree=args.vdegree, + vtype=uw.VarType.VECTOR, +) + +p_soln = uw.discretisation.MeshVariable( + varname="Pressure", + mesh=mesh, + degree=args.pdegree, + vtype=uw.VarType.SCALAR, + continuous=args.pcont, +) + +# %% [markdown] +# ## Stokes system + +# %% +stokes = Stokes(mesh, velocityField=v_soln, pressureField=p_soln) +stokes.constitutive_model = uw.constitutive_models.ViscousFlowModel +stokes.constitutive_model.Parameters.viscosity = 1.0 +stokes.saddle_preconditioner = 1.0 +stokes.bodyforce = bodyforce + +for boundary_name in ( + mesh.boundaries.Bottom.name, + mesh.boundaries.Top.name, + mesh.boundaries.Left.name, + mesh.boundaries.Right.name, +): + stokes.add_essential_bc(v_ana_expr, boundary_name) + +stokes.tolerance = args.stokes_tol +stokes.petsc_options["snes_type"] = "ksponly" +stokes.petsc_options["ksp_type"] = "fgmres" +stokes.petsc_options["ksp_rtol"] = args.stokes_tol +stokes.petsc_options["ksp_atol"] = 0.0 +stokes.petsc_options["ksp_monitor"] = None +stokes.petsc_options["ksp_monitor_true_residual"] = None +stokes.petsc_options["ksp_converged_reason"] = None + +stokes.petsc_options.setValue("fieldsplit_velocity_pc_mg_type", "kaskade") +stokes.petsc_options.setValue("fieldsplit_velocity_pc_mg_cycle_type", "w") +stokes.petsc_options["fieldsplit_velocity_mg_coarse_pc_type"] = "svd" +stokes.petsc_options["fieldsplit_velocity_ksp_type"] = "fcg" +stokes.petsc_options["fieldsplit_velocity_mg_levels_ksp_type"] = "chebyshev" +stokes.petsc_options["fieldsplit_velocity_mg_levels_ksp_max_it"] = 5 +stokes.petsc_options["fieldsplit_velocity_mg_levels_ksp_converged_maxits"] = None +stokes.petsc_options.setValue("fieldsplit_pressure_pc_type", "mg") +stokes.petsc_options.setValue("fieldsplit_pressure_pc_mg_type", "multiplicative") +stokes.petsc_options.setValue("fieldsplit_pressure_pc_mg_cycle_type", "v") + +# %% [markdown] +# ## Solve + +# %% +stokes.solve(verbose=False, debug=False) + +# %% [markdown] +# ## Pressure gauge + +# %% +p_int = uw.maths.Integral(mesh, p_soln.sym[0]).evaluate() +volume = uw.maths.Integral(mesh, 1.0).evaluate() +p_soln.data[:, 0] -= p_int / volume + +# %% [markdown] +# ## Relative `L2` errors + +# %% +v_err_expr = sp.Matrix(v_soln.sym).T - v_ana_expr +p_err_expr = p_soln.sym[0] - p_ana_expr + +v_err_sq_expr = v_err_expr.dot(v_err_expr) +v_ana_sq_expr = v_ana_expr.dot(v_ana_expr) +p_err_sq_expr = p_err_expr * p_err_expr +p_ana_sq_expr = p_ana_expr * p_ana_expr + +v_err_l2 = np.sqrt(uw.maths.Integral(mesh, v_err_sq_expr).evaluate()) / np.sqrt( + uw.maths.Integral(mesh, v_ana_sq_expr).evaluate() +) +p_err_l2 = np.sqrt(uw.maths.Integral(mesh, p_err_sq_expr).evaluate()) / np.sqrt( + uw.maths.Integral(mesh, p_ana_sq_expr).evaluate() +) + +# %% [markdown] +# ## Report + +# %% +uw.pprint("cellsize:", args.cellsize) +uw.pprint("vdegree:", args.vdegree) +uw.pprint("pdegree:", args.pdegree) +uw.pprint("Relative velocity L2 error:", v_err_l2) +uw.pprint("Relative pressure L2 error:", p_err_l2) From 2f67421f4f12c0412be7d66de12909238683e03c Mon Sep 17 00:00:00 2001 From: Louis Moresi Date: Thu, 19 Mar 2026 21:05:14 +1100 Subject: [PATCH 4/4] Add dual-space options to mesh coordinate FE and comment Stokes fix MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit The Stokes velocity FE fix (previous commit) is mirrored here for the mesh coordinate projection FE in discretisation_mesh.py, which had the same omission. This only affects higher-order (curved) coordinate meshes on simplices — P1 meshes work by accident with PETSc defaults. Added comments to both sites explaining why these options are required. Underworld development team with AI support from Claude Code --- src/underworld3/cython/petsc_generic_snes_solvers.pyx | 7 ++++++- src/underworld3/discretisation/discretisation_mesh.py | 5 +++++ 2 files changed, 11 insertions(+), 1 deletion(-) diff --git a/src/underworld3/cython/petsc_generic_snes_solvers.pyx b/src/underworld3/cython/petsc_generic_snes_solvers.pyx index f7df4b77d..20d5c7e2a 100644 --- a/src/underworld3/cython/petsc_generic_snes_solvers.pyx +++ b/src/underworld3/cython/petsc_generic_snes_solvers.pyx @@ -3403,8 +3403,13 @@ class SNES_Stokes_SaddlePt(SolverBaseClass): print(f"{uw.mpi.rank}: Building FE / quadrature for {self.name}", flush=True) + # Both degree AND dual-space options must be set before createDefault(). + # The dual-space options control node placement on simplices: without them, + # PETSc may choose a different nodal basis than the solver's weak form expects. + # For P2 the defaults happen to be correct, but P3+ gets wrong node placement. + options = PETSc.Options() - options.setValue("private_{}_u_petscspace_degree".format(self.petsc_options_prefix), u_degree) # for private variables + options.setValue("private_{}_u_petscspace_degree".format(self.petsc_options_prefix), u_degree) options.setValue("private_{}_u_petscdualspace_lagrange_continuity".format(self.petsc_options_prefix), self.u.continuous) options.setValue("private_{}_u_petscdualspace_lagrange_node_endpoints".format(self.petsc_options_prefix), False) self.petsc_fe_u = PETSc.FE().createDefault(mesh.dim, mesh.dim, mesh.isSimplex, mesh.qdegree, "private_{}_u_".format(self.petsc_options_prefix), PETSc.COMM_SELF) diff --git a/src/underworld3/discretisation/discretisation_mesh.py b/src/underworld3/discretisation/discretisation_mesh.py index cb339c636..eb0784916 100644 --- a/src/underworld3/discretisation/discretisation_mesh.py +++ b/src/underworld3/discretisation/discretisation_mesh.py @@ -1096,8 +1096,13 @@ def nuke_coords_and_rebuild( # later where we call the interpolation routines to project from the linear # mesh coordinates to other mesh coordinates. + # Dual-space options control node placement on simplices and must be set + # before createDefault(). Currently only P1 coordinate meshes are used, + # but these are needed for higher-order (curved) coordinate meshes. options = PETSc.Options() options.setValue(f"meshproj_{self.mesh_instances}_petscspace_degree", self.degree) + options.setValue(f"meshproj_{self.mesh_instances}_petscdualspace_lagrange_continuity", True) + options.setValue(f"meshproj_{self.mesh_instances}_petscdualspace_lagrange_node_endpoints", False) self.petsc_fe = PETSc.FE().createDefault( self.dim,