Skip to content

Commit 684e899

Browse files
gthyagigthyagi
committed
Convert box MMS reproducer to notebook style
1 parent b5e7858 commit 684e899

1 file changed

Lines changed: 165 additions & 129 deletions

File tree

‎examples/stokes_box_mms_simplex.py‎

Lines changed: 165 additions & 129 deletions
Original file line numberDiff line numberDiff line change
@@ -1,7 +1,7 @@
11
"""
2-
Simple manufactured-solution Stokes example on a triangular box.
2+
Notebook-style manufactured-solution Stokes example on a triangular box.
33
4-
This script is useful for checking higher-order simplex velocity / pressure
4+
This example is useful for checking higher-order simplex velocity / pressure
55
pairs on a straight-edged domain without curved-geometry effects.
66
77
Default setup:
@@ -27,130 +27,166 @@
2727
import underworld3 as uw
2828
from underworld3.systems import Stokes
2929

30-
31-
def parse_args():
32-
parser = argparse.ArgumentParser()
33-
parser.add_argument("--cellsize", type=float, default=0.125)
34-
parser.add_argument("--vdegree", type=int, default=3)
35-
parser.add_argument("--pdegree", type=int, default=2)
36-
parser.add_argument("--pcont", action="store_true", default=True)
37-
parser.add_argument("--stokes-tol", type=float, default=1.0e-10)
38-
return parser.parse_args()
39-
40-
41-
def subtract_pressure_mean(mesh, pressure_var):
42-
p_int = uw.maths.Integral(mesh, pressure_var.sym[0]).evaluate()
43-
volume = uw.maths.Integral(mesh, 1.0).evaluate()
44-
pressure_var.data[:, 0] -= p_int / volume
45-
46-
47-
def relative_l2_error(mesh, err_expr, ana_expr):
48-
if isinstance(err_expr, sp.MatrixBase):
49-
err_expr = err_expr.dot(err_expr)
50-
ana_expr = ana_expr.dot(ana_expr)
51-
else:
52-
err_expr = err_expr * err_expr
53-
ana_expr = ana_expr * ana_expr
54-
55-
err_I = uw.maths.Integral(mesh, err_expr)
56-
ana_I = uw.maths.Integral(mesh, ana_expr)
57-
58-
return np.sqrt(err_I.evaluate()) / np.sqrt(ana_I.evaluate())
59-
60-
61-
def main():
62-
args = parse_args()
63-
qdegree = max(2 * args.vdegree, args.vdegree + args.pdegree)
64-
65-
mesh = uw.meshing.UnstructuredSimplexBox(
66-
minCoords=(0.0, 0.0),
67-
maxCoords=(1.0, 1.0),
68-
cellSize=args.cellsize,
69-
regular=True,
70-
qdegree=qdegree,
71-
filename=f"/tmp/stokes_box_mms_simplex_{args.vdegree}_{args.pdegree}_{args.cellsize}.msh",
72-
)
73-
74-
x, y = mesh.X
75-
76-
psi = x**2 * y**2
77-
v_ana_expr = sp.Matrix([sp.diff(psi, y), -sp.diff(psi, x)])
78-
p_ana_expr = x**2 - sp.Rational(1, 3)
79-
80-
bodyforce = sp.Matrix(
81-
[
82-
-(sp.diff(v_ana_expr[0], x, 2) + sp.diff(v_ana_expr[0], y, 2))
83-
+ sp.diff(p_ana_expr, x),
84-
-(sp.diff(v_ana_expr[1], x, 2) + sp.diff(v_ana_expr[1], y, 2))
85-
+ sp.diff(p_ana_expr, y),
86-
]
87-
)
88-
89-
v_soln = uw.discretisation.MeshVariable(
90-
varname="Velocity",
91-
mesh=mesh,
92-
degree=args.vdegree,
93-
vtype=uw.VarType.VECTOR,
94-
)
95-
96-
p_soln = uw.discretisation.MeshVariable(
97-
varname="Pressure",
98-
mesh=mesh,
99-
degree=args.pdegree,
100-
vtype=uw.VarType.SCALAR,
101-
continuous=args.pcont,
102-
)
103-
104-
stokes = Stokes(mesh, velocityField=v_soln, pressureField=p_soln)
105-
stokes.constitutive_model = uw.constitutive_models.ViscousFlowModel
106-
stokes.constitutive_model.Parameters.viscosity = 1.0
107-
stokes.saddle_preconditioner = 1.0
108-
stokes.bodyforce = bodyforce
109-
110-
for boundary_name in (
111-
mesh.boundaries.Bottom.name,
112-
mesh.boundaries.Top.name,
113-
mesh.boundaries.Left.name,
114-
mesh.boundaries.Right.name,
115-
):
116-
stokes.add_essential_bc(v_ana_expr, boundary_name)
117-
118-
stokes.tolerance = args.stokes_tol
119-
stokes.petsc_options["snes_type"] = "ksponly"
120-
stokes.petsc_options["ksp_type"] = "fgmres"
121-
stokes.petsc_options["ksp_rtol"] = args.stokes_tol
122-
stokes.petsc_options["ksp_atol"] = 0.0
123-
stokes.petsc_options["ksp_monitor"] = None
124-
stokes.petsc_options["ksp_monitor_true_residual"] = None
125-
stokes.petsc_options["ksp_converged_reason"] = None
126-
127-
stokes.petsc_options.setValue("fieldsplit_velocity_pc_mg_type", "kaskade")
128-
stokes.petsc_options.setValue("fieldsplit_velocity_pc_mg_cycle_type", "w")
129-
stokes.petsc_options["fieldsplit_velocity_mg_coarse_pc_type"] = "svd"
130-
stokes.petsc_options["fieldsplit_velocity_ksp_type"] = "fcg"
131-
stokes.petsc_options["fieldsplit_velocity_mg_levels_ksp_type"] = "chebyshev"
132-
stokes.petsc_options["fieldsplit_velocity_mg_levels_ksp_max_it"] = 5
133-
stokes.petsc_options["fieldsplit_velocity_mg_levels_ksp_converged_maxits"] = None
134-
stokes.petsc_options.setValue("fieldsplit_pressure_pc_type", "mg")
135-
stokes.petsc_options.setValue("fieldsplit_pressure_pc_mg_type", "multiplicative")
136-
stokes.petsc_options.setValue("fieldsplit_pressure_pc_mg_cycle_type", "v")
137-
138-
stokes.solve(verbose=False, debug=False)
139-
140-
subtract_pressure_mean(mesh, p_soln)
141-
142-
v_err_expr = sp.Matrix(v_soln.sym).T - v_ana_expr
143-
p_err_expr = p_soln.sym[0] - p_ana_expr
144-
145-
v_err_l2 = relative_l2_error(mesh, v_err_expr, v_ana_expr)
146-
p_err_l2 = relative_l2_error(mesh, p_err_expr, p_ana_expr)
147-
148-
uw.pprint("cellsize:", args.cellsize)
149-
uw.pprint("vdegree:", args.vdegree)
150-
uw.pprint("pdegree:", args.pdegree)
151-
uw.pprint("Relative velocity L2 error:", v_err_l2)
152-
uw.pprint("Relative pressure L2 error:", p_err_l2)
153-
154-
155-
if __name__ == "__main__":
156-
main()
30+
# %% [markdown]
31+
# # Stokes Box MMS on a Simplex Mesh
32+
#
33+
# This notebook-style example solves a manufactured Stokes problem on the unit
34+
# square using triangular elements. It is intentionally simple and useful for
35+
# checking higher-order simplex pairs such as `P2/P1` and `P3/P2`.
36+
37+
# %% [markdown]
38+
# ## Runtime parameters
39+
#
40+
# The script still accepts command-line arguments so it can be run directly.
41+
42+
# %%
43+
parser = argparse.ArgumentParser()
44+
parser.add_argument("--cellsize", type=float, default=0.125)
45+
parser.add_argument("--vdegree", type=int, default=3)
46+
parser.add_argument("--pdegree", type=int, default=2)
47+
parser.add_argument("--pcont", action="store_true", default=True)
48+
parser.add_argument("--stokes-tol", type=float, default=1.0e-10)
49+
args = parser.parse_args()
50+
51+
qdegree = max(2 * args.vdegree, args.vdegree + args.pdegree)
52+
53+
# %% [markdown]
54+
# ## Mesh
55+
56+
# %%
57+
mesh = uw.meshing.UnstructuredSimplexBox(
58+
minCoords=(0.0, 0.0),
59+
maxCoords=(1.0, 1.0),
60+
cellSize=args.cellsize,
61+
regular=True,
62+
qdegree=qdegree,
63+
filename=f"/tmp/stokes_box_mms_simplex_{args.vdegree}_{args.pdegree}_{args.cellsize}.msh",
64+
)
65+
66+
x, y = mesh.X
67+
68+
# %% [markdown]
69+
# ## Manufactured solution
70+
#
71+
# Stream function:
72+
# `psi = x^2 y^2`
73+
#
74+
# Exact velocity:
75+
# `u = (dpsi/dy, -dpsi/dx)`
76+
#
77+
# Exact pressure:
78+
# `p = x^2 - 1/3`
79+
80+
# %%
81+
psi = x**2 * y**2
82+
v_ana_expr = sp.Matrix([sp.diff(psi, y), -sp.diff(psi, x)])
83+
p_ana_expr = x**2 - sp.Rational(1, 3)
84+
85+
bodyforce = sp.Matrix(
86+
[
87+
-(sp.diff(v_ana_expr[0], x, 2) + sp.diff(v_ana_expr[0], y, 2))
88+
+ sp.diff(p_ana_expr, x),
89+
-(sp.diff(v_ana_expr[1], x, 2) + sp.diff(v_ana_expr[1], y, 2))
90+
+ sp.diff(p_ana_expr, y),
91+
]
92+
)
93+
94+
# %% [markdown]
95+
# ## Discretisation
96+
97+
# %%
98+
v_soln = uw.discretisation.MeshVariable(
99+
varname="Velocity",
100+
mesh=mesh,
101+
degree=args.vdegree,
102+
vtype=uw.VarType.VECTOR,
103+
)
104+
105+
p_soln = uw.discretisation.MeshVariable(
106+
varname="Pressure",
107+
mesh=mesh,
108+
degree=args.pdegree,
109+
vtype=uw.VarType.SCALAR,
110+
continuous=args.pcont,
111+
)
112+
113+
# %% [markdown]
114+
# ## Stokes system
115+
116+
# %%
117+
stokes = Stokes(mesh, velocityField=v_soln, pressureField=p_soln)
118+
stokes.constitutive_model = uw.constitutive_models.ViscousFlowModel
119+
stokes.constitutive_model.Parameters.viscosity = 1.0
120+
stokes.saddle_preconditioner = 1.0
121+
stokes.bodyforce = bodyforce
122+
123+
for boundary_name in (
124+
mesh.boundaries.Bottom.name,
125+
mesh.boundaries.Top.name,
126+
mesh.boundaries.Left.name,
127+
mesh.boundaries.Right.name,
128+
):
129+
stokes.add_essential_bc(v_ana_expr, boundary_name)
130+
131+
stokes.tolerance = args.stokes_tol
132+
stokes.petsc_options["snes_type"] = "ksponly"
133+
stokes.petsc_options["ksp_type"] = "fgmres"
134+
stokes.petsc_options["ksp_rtol"] = args.stokes_tol
135+
stokes.petsc_options["ksp_atol"] = 0.0
136+
stokes.petsc_options["ksp_monitor"] = None
137+
stokes.petsc_options["ksp_monitor_true_residual"] = None
138+
stokes.petsc_options["ksp_converged_reason"] = None
139+
140+
stokes.petsc_options.setValue("fieldsplit_velocity_pc_mg_type", "kaskade")
141+
stokes.petsc_options.setValue("fieldsplit_velocity_pc_mg_cycle_type", "w")
142+
stokes.petsc_options["fieldsplit_velocity_mg_coarse_pc_type"] = "svd"
143+
stokes.petsc_options["fieldsplit_velocity_ksp_type"] = "fcg"
144+
stokes.petsc_options["fieldsplit_velocity_mg_levels_ksp_type"] = "chebyshev"
145+
stokes.petsc_options["fieldsplit_velocity_mg_levels_ksp_max_it"] = 5
146+
stokes.petsc_options["fieldsplit_velocity_mg_levels_ksp_converged_maxits"] = None
147+
stokes.petsc_options.setValue("fieldsplit_pressure_pc_type", "mg")
148+
stokes.petsc_options.setValue("fieldsplit_pressure_pc_mg_type", "multiplicative")
149+
stokes.petsc_options.setValue("fieldsplit_pressure_pc_mg_cycle_type", "v")
150+
151+
# %% [markdown]
152+
# ## Solve
153+
154+
# %%
155+
stokes.solve(verbose=False, debug=False)
156+
157+
# %% [markdown]
158+
# ## Pressure gauge
159+
160+
# %%
161+
p_int = uw.maths.Integral(mesh, p_soln.sym[0]).evaluate()
162+
volume = uw.maths.Integral(mesh, 1.0).evaluate()
163+
p_soln.data[:, 0] -= p_int / volume
164+
165+
# %% [markdown]
166+
# ## Relative `L2` errors
167+
168+
# %%
169+
v_err_expr = sp.Matrix(v_soln.sym).T - v_ana_expr
170+
p_err_expr = p_soln.sym[0] - p_ana_expr
171+
172+
v_err_sq_expr = v_err_expr.dot(v_err_expr)
173+
v_ana_sq_expr = v_ana_expr.dot(v_ana_expr)
174+
p_err_sq_expr = p_err_expr * p_err_expr
175+
p_ana_sq_expr = p_ana_expr * p_ana_expr
176+
177+
v_err_l2 = np.sqrt(uw.maths.Integral(mesh, v_err_sq_expr).evaluate()) / np.sqrt(
178+
uw.maths.Integral(mesh, v_ana_sq_expr).evaluate()
179+
)
180+
p_err_l2 = np.sqrt(uw.maths.Integral(mesh, p_err_sq_expr).evaluate()) / np.sqrt(
181+
uw.maths.Integral(mesh, p_ana_sq_expr).evaluate()
182+
)
183+
184+
# %% [markdown]
185+
# ## Report
186+
187+
# %%
188+
uw.pprint("cellsize:", args.cellsize)
189+
uw.pprint("vdegree:", args.vdegree)
190+
uw.pprint("pdegree:", args.pdegree)
191+
uw.pprint("Relative velocity L2 error:", v_err_l2)
192+
uw.pprint("Relative pressure L2 error:", p_err_l2)

0 commit comments

Comments
 (0)