import underworld3 as uw
import numpy as np
import sympyNotebook 5: Stokes Equation
Mesh with embedded internal surface where we apply harmonic thermal forcing. This allows us to compute topography and gravity kernels. We’ll use a cylindrical annulus this timeand we’ll use free-slip boundary conditions throughout.
res = 0.075
r_o = 1.0
r_int = 0.825
r_i = 0.55meshball = uw.meshing.AnnulusInternalBoundary(radiusOuter=r_o,
radiusInternal=r_int,
radiusInner=r_i,
cellSize_Inner=res,
cellSize_Internal=res*0.5,
cellSize_Outer=res,
centre=False,)
meshball.view()
x, y = meshball.CoordinateSystem.X
r, th = meshball.CoordinateSystem.xR
unit_rvec = meshball.CoordinateSystem.unit_e_0
# Orientation of surface normals
Gamma_N = unit_rvec
Mesh # 0: .meshes/uw_annulus_internalBoundary_rO1.0rInt0.825_rI0.55_csize0.1_csizefs0.075.msh
No variables are defined on the mesh
| Boundary Name | ID | Min Size | Max Size |
| ------------------------------------------------------ |
| Lower | 1 | 94 | 94 |
| Internal | 2 | 278 | 278 |
| Upper | 3 | 166 | 166 |
| Centre | 10 | 0 | 0 |
| All_Edges | 1000 | 0 | 0 |
| All_Boundaries | 1001 | 132 | 132 |
| UW_Boundaries | -- | 538 | 538 |
| ------------------------------------------------------ |
DM Object: uw_.meshes/uw_annulus_internalBoundary_rO1.0rInt0.825_rI0.55_csize0.1_csizefs0.075.msh 1 MPI process
type: plex
uw_.meshes/uw_annulus_internalBoundary_rO1.0rInt0.825_rI0.55_csize0.1_csizefs0.075.msh in 2 dimensions:
Number of 0-cells per rank: 1142
Number of 1-cells per rank: 3294
Number of 2-cells per rank: 2152
Labels:
depth: 3 strata with value/size (0 (1142), 1 (3294), 2 (2152))
All_Boundaries: 1 strata with value/size (1001 (132))
Elements: 1 strata with value/size (666666 (3022))
Internal: 1 strata with value/size (2 (278))
Lower: 1 strata with value/size (1 (94))
Upper: 1 strata with value/size (3 (166))
celltype: 3 strata with value/size (0 (1142), 1 (3294), 3 (2152))
UW_Boundaries: 3 strata with value/size (1 (94), 2 (278), 3 (166))Solver setup
We can obtain unit vectors in the natural coordinate system (here \(r\), \(\theta\)) as mesh.CoordinateSystem.unit_e_0, mesh.CoordinateSystem.unit_e_1. There is a null space if we apply the boundary conditions exactly, and so we define a function to represent the null space.
We can set solver options via the petsc_options interface on the solver.
# Create a density structure / buoyancy force
# gravity will vary linearly from zero at the centre
# of the sphere to (say) 1 at the surface
# Null space in velocity (constant v_theta) expressed in x,y coordinates
v_theta_fn_xy = r * meshball.CoordinateSystem.rRotN.T * sympy.Matrix((0,1))
# Mesh variables for the unknowns
v_soln = uw.discretisation.MeshVariable("V0", meshball, 2, degree=2, varsymbol=r"{v_0}")
p_soln = uw.discretisation.MeshVariable("p", meshball, 1, degree=1, continuous=True)
stokes = uw.systems.Stokes(
meshball,
velocityField=v_soln,
pressureField=p_soln,
)
stokes.constitutive_model = uw.constitutive_models.ViscousFlowModel
stokes.constitutive_model.Parameters.shear_viscosity_0 = 1
stokes.tolerance = 1.0e-6
stokes.petsc_options.setValue("ksp_monitor", None)
stokes.petsc_options.setValue("snes_monitor", None)
stokes.petsc_options["fieldsplit_velocity_mg_coarse_pc_type"] = "svd"Natural boundary conditions
In FEM, natural boundary conditions are fluxes specified at a boundary (through surface integrals). We can also apply integrals to internal surfaces. We need to compute the vector components of the boundary condition (normal / tangential) and supply them in the Cartesian frame.
The bouyancy force on the internal surface is already radial, so this translates to
stokes.add_natural_bc(-t_init * unit_rvec, "Internal")To set a no-normal-flow boundary condition, we need to penalise the radial velocity at the boundaries. Symbolically this is
\[ \mathbf{f} = \lambda \left( \mathbf{v} \cdot \Gamma_N \right) \,\, \Gamma_N \]
where \(\lambda\) is a large penalty value. This translates into sympy code as
stokes.add_natural_bc(
penalty * Gamma_N.dot(v_soln.sym) * Gamma_N, "Lower"
)
stokes.bodyforce = sympy.Matrix([0,0])
# Note, the thermal bouyancy field is localised in the radius using a
# gaussian solely for the purposes of plotting. It is automatically
# a delta function when applied through a surface integral
t_init = sympy.sin(5*th) * sympy.exp(-1000.0 * ((r - r_int) ** 2))
stokes.add_natural_bc(-t_init * unit_rvec, "Internal")
stokes.add_natural_bc(10000 * Gamma_N.dot(v_soln.sym) * Gamma_N, "Upper")
if r_i != 0.0:
stokes.add_natural_bc(10000 * Gamma_N.dot(v_soln.sym) * Gamma_N, "Lower")
stokes.solve() 0 SNES Function norm 2.307507502179e-01
Residual norms for stokes_ solve.
0 KSP Residual norm 8.742204758135e+00
1 KSP Residual norm 7.603223176805e-07
1 SNES Function norm 4.751452850114e-08Removal of the null space
We can use the uw.maths.Integral to compute the inner product of the null space and the velocity solution. It is not zero, so we remove it !
# Null space evaluation
I0 = uw.maths.Integral(meshball, v_theta_fn_xy.dot(v_theta_fn_xy))
vnorm = I0.evaluate()
I0.fn = v_theta_fn_xy.dot(v_soln.sym)
ns = I0.evaluate()
print(ns/vnorm, vnorm)
with meshball.access(v_soln):
dv = uw.function.evaluate(ns * v_theta_fn_xy, v_soln.coords) / vnorm
v_soln.data[...] -= dv
ns = I0.evaluate()
print(ns/vnorm, vnorm)0.0005624988949826505 1.4249504223200373
1.9675625281382883e-19 1.4249504223200373
if uw.mpi.size == 1:
import pyvista as pv
import underworld3.visualisation as vis
pvmesh = vis.mesh_to_pv_mesh(meshball)
pvmesh.point_data["P"] = vis.scalar_fn_to_pv_points(pvmesh, p_soln.sym)
pvmesh.point_data["V"] = vis.vector_fn_to_pv_points(pvmesh, v_soln.sym)
pvmesh.point_data["T"] = vis.scalar_fn_to_pv_points(pvmesh, t_init)
skip = 1
points = np.zeros((meshball._centroids[::skip].shape[0], 3))
points[:, 0] = meshball._centroids[::skip, 0]
points[:, 1] = meshball._centroids[::skip, 1]
point_cloud = pv.PolyData(points)
pvstream = pvmesh.streamlines_from_source(
point_cloud, vectors="V",
integration_direction="both",
integrator_type=45,
surface_streamlines=True,
initial_step_length=0.01,
max_time=0.5,
max_steps=500,
)
pl = pv.Plotter(window_size=(750, 750))
pl.add_mesh(
pvmesh,
cmap="coolwarm",
edge_color="Grey",
edge_opacity=0.33,
scalars="T",
show_edges=True,
use_transparency=False,
opacity=1.0,
show_scalar_bar=False
)
pl.add_mesh(pvstream, opacity=0.3, show_scalar_bar=False, cmap="Greens", render_lines_as_tubes=False)
pl.export_html("html5/stokes_annulus_plot.html")from IPython.display import IFrame
IFrame(src="html5/stokes_annulus_plot.html", width=500, height=400)Interactive Image: Annulus mesh of triangular elements on which we evaluated Stokes flow driven by an internal delta function buoyancy. Boundary conditions are free slip, imposed using a penalty on the radial velocity at the boundary