import underworld3 as uw
import numpy as np
import sympyNotebook 4: Poisson Equation
In this example, we solve a steady, non-linear, scalar diffusion equation and show how to recover the fluxes.
mesh = uw.meshing.UnstructuredSimplexBox(
minCoords=(0.0, 0.0), maxCoords=(1.0, 1.0), cellSize=1 / 12, qdegree=3
)
x, y = mesh.X
mesh.view()
Mesh # 0: .meshes/uw_simplexbox_minC(0.0, 0.0)_maxC(1.0, 1.0)_csize0.08333333333333333_regFalse.msh
No variables are defined on the mesh
| Boundary Name | ID | Min Size | Max Size |
| ------------------------------------------------------ |
| Bottom | 11 | 23 | 23 |
| Top | 12 | 23 | 23 |
| Right | 13 | 23 | 23 |
| Left | 14 | 23 | 23 |
| All_Boundaries | 1001 | 48 | 48 |
| UW_Boundaries | -- | 92 | 92 |
| ------------------------------------------------------ |
DM Object: uw_.meshes/uw_simplexbox_minC(0.0, 0.0)_maxC(1.0, 1.0)_csize0.08333333333333333_regFalse.msh 1 MPI process
type: plex
uw_.meshes/uw_simplexbox_minC(0.0, 0.0)_maxC(1.0, 1.0)_csize0.08333333333333333_regFalse.msh in 2 dimensions:
Number of 0-cells per rank: 198
Number of 1-cells per rank: 543
Number of 2-cells per rank: 346
Labels:
depth: 3 strata with value/size (0 (198), 1 (543), 2 (346))
All_Boundaries: 1 strata with value/size (1001 (48))
Bottom: 1 strata with value/size (11 (23))
Elements: 1 strata with value/size (99999 (496))
Left: 1 strata with value/size (14 (23))
Right: 1 strata with value/size (13 (23))
Top: 1 strata with value/size (12 (23))
celltype: 3 strata with value/size (0 (198), 1 (543), 3 (346))
UW_Boundaries: 4 strata with value/size (11 (23), 12 (23), 13 (23), 14 (23))# mesh variables
t_soln = uw.discretisation.MeshVariable("T", mesh, 1, degree=3)
dTdY = uw.discretisation.MeshVariable(
r"dTdY", mesh, 1, degree=2,
varsymbol=r"\partial T/\partial y")The Poisson Solver
There are a number of pre-defined solver systems defined in underworld3 which are templates for orchestrating the underlying PETSc objects. A solver requires us to define the unknown in the form of meshVariables, provide boundary conditions, a constitutive model, and provide uw.functions to define constitutive properties, and driving terms.
We will use the Poisson solver for the diffusion equation, and we will use a Projection solver to compute the vertical gradient term.
The solver classes themselves are documented, so we can figure out what is needed before we define anything:
uw.systems.Poisson.view()This class provides functionality for a discrete representation of the Poisson equation
\[ \nabla \cdot \color{Blue}{\underbrace{\Bigl[ \boldsymbol\kappa \nabla u \Bigr]}_{\mathbf{F}}} = \color{Maroon}{\underbrace{\Bigl[ f \Bigl] }_{\mathbf{f}}} \]
The term \(\mathbf{F}\) relates the flux to gradients in the unknown \(u\)
Properties
The unknown is \(u\)
The diffusivity tensor, \(\kappa\) is provided by setting the
constitutive_modelproperty to one of the scalaruw.constitutive_modelsclasses and populating the parameters. It is usually a constant or a function of position / time and may also be non-linear or anisotropic.\(f\) is a volumetric source term
uw.systems.Projection.view()Solves \(u = \tilde{f}\) where \(\tilde{f}\) is a function that can be evaluated within an element and \(u\) is a meshVariable with associated shape functions. Typically, the projection is used to obtain a continuous representation of a function that is not well defined at the mesh nodes. For example, functions of the spatial derivatives of one or more meshVariable (e.g. components of fluxes) can be mapped to continuous variables with a projection. More broadly it is a projection from one basis to another and its limitations should be evaluated within that context.
The projection implemented by creating a solver for this problem
\[ \nabla \cdot \color{Blue}{\underbrace{\Bigl[ \boldsymbol\alpha \nabla u \Bigr]}_{\mathbf{F}}} - \color{Maroon}{\underbrace{\Bigl[ u - \tilde{f} \Bigl] }_{\mathbf{f}}} = 0 \]
Where the term \(\mathbf{F}\) provides a smoothing regularization. \(\alpha\) can be zero.
Constitutive Models
Most of the solvers require a constitutive model to be provided and its parameters populated. This is to allow flexibility in defining / redefining solvers during a model calculation.
We need a diffusion model. We can look at the documentation first.
uw.constitutive_models.DiffusionModel.view()class DiffusionModel(Constitutive_Model)
...diffusion_model = DiffusionModel(dim)
diffusion_model.material_properties = diffusion_model.Parameters(diffusivity=diffusivity_fn)
scalar_solver.constititutive_model = diffusion_model\[q_{i} = \kappa_{ij} \cdot \frac{\partial \phi}{\partial x_j}\]
where \(\kappa\) is a diffusivity, a scalar constant, sympy function, underworld mesh variable or any valid combination of those types. Access the constitutive model using:
flux = diffusion_model.flux(gradient_matrix)Problem setup
We will make the problem “interesting” by defining a non-linear diffusivity
\[ \kappa = \kappa_0 + \left| \nabla T \right |^2 \]
and apply a zero source term, non-zero boundary conditions on the top and bottom surfaces. The names of the boundaries of the mesh can be found by looking at mesh.view() above.
# The steady-state diffusion
poisson = uw.systems.Poisson(mesh, u_Field=t_soln)
poisson.constitutive_model = uw.constitutive_models.DiffusionModel
# Non-linear diffusivity
delT = mesh.vector.gradient(t_soln.sym)
k = 1 - (delT.dot(delT))
abs_r2 = x**2 + y**2
poisson.f = -abs_r2
# Essential boundary conditions
poisson.add_essential_bc([0], "Bottom")
poisson.add_essential_bc([1], "Top" )
# Create Projection solver for gradient terms
gradient = uw.systems.Projection(mesh, dTdY)
gradient.uw_function = t_soln.sym.diff(y)
gradient.smoothing = 1.0e-6
# Linear model - starting guess
poisson.constitutive_model.Parameters.diffusivity = 1
poisson.solve(zero_init_guess=True)
# Solve time
poisson.constitutive_model.Parameters.diffusivity = k
poisson.solve(zero_init_guess=False)
# Solve the gradient
gradient.solve()Non-linear problems, Jacobian
This is a non-linear problem and has a non-zero Jacobian (\(\partial F / \partial \nabla T )\). It’s a little bit hidden, but is available if you feel like checking:
display(poisson._G3)\(\displaystyle \left[\begin{matrix}- 3 {T}_{,0}^{ 2 }(\mathbf{x}) - {T}_{,1}^{ 2 }(\mathbf{x}) + 1 & - 2 {T}_{,0}(\mathbf{x}) {T}_{,1}(\mathbf{x})\\- 2 {T}_{,0}(\mathbf{x}) {T}_{,1}(\mathbf{x}) & - {T}_{,0}^{ 2 }(\mathbf{x}) - 3 {T}_{,1}^{ 2 }(\mathbf{x}) + 1\end{matrix}\right]\)
from mpi4py import MPI
if MPI.COMM_WORLD.size == 1:
import pyvista as pv
import underworld3.visualisation as vis
pvmesh = vis.mesh_to_pv_mesh(mesh)
pvmesh.point_data["T"] = vis.scalar_fn_to_pv_points(pvmesh, t_soln.sym)
pvmesh.point_data["dTdY"] = vis.scalar_fn_to_pv_points(pvmesh, dTdY.sym)
pl = pv.Plotter(window_size=(1000, 500), shape=(1, 2))
pl.subplot(0, 0)
pl.add_mesh(
pvmesh,
cmap="coolwarm",
edge_color="Black",
show_edges=True,
scalars="T",
use_transparency=False,
opacity=1,
show_scalar_bar=True,
)
pl.subplot(0, 1)
pl.add_mesh(
pvmesh,
cmap="coolwarm",
edge_color="Black",
show_edges=True,
scalars="dTdY",
use_transparency=False,
opacity=1,
scalar_bar_args=dict(vertical=False)
)
pl.export_html("html5/temperature_plot.html")from IPython.display import IFrame
IFrame(src="html5/temperature_plot.html", width=600, height=300)Interactive Image: Square mesh of triangular elements on which we evaluated a non-linear diffusion problem and the vertical gradient of the solution
Exercise 4.1
Since diffusivity, \(\kappa\), depends on the gradient of temperature, it also needs to be calculated using a projection to the mesh nodes. You can use the computation for \(\partial T / \partial y\) as a template to calculate and plot \(\kappa\), or you can compute the vertical flux term and plot that.