import underworld3 as uw
import numpy as np
import sympyNotebook 4: Poisson Equation
Introduces the various solver templates that are available in Underworld, starting with a steady-state diffusion problem. The template requires you to set some constitutive properties and define the unknowns. These are handled through subsitution into symbolic forms and the template equation can be inspected before you need to supply concrete expressions.
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
Number of cells: 346
No variables are defined on the mesh
| Boundary Name | ID |
| -------------------------------- |
| Bottom | 11 |
| Top | 12 |
| Right | 13 |
| Left | 14 |
| Null_Boundary | 666 |
| All_Boundaries | 1001 |
| All_Boundaries | 1001 |
| UW_Boundaries | -- |
| -------------------------------- |
Use view(1) to view detailed mesh information.
# 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 (Diffusion) 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 Projection solvers compute the meshVariable representation of the right hand side of the equation (potentially requiring some smoothing to regularise the problem). This is required when computing numerical values of gradients, for example, where the values are often not formally defined on the mesh points of the meshVariable. This can also be used to transfer information from discontinous variables to continuous representations when you need to visualise a field.
The solver classes themselves (the templates for specific solvers) are separately documented, so we can figure out what is needed before we define anything:
uw.systems.Poisson.view()Poisson Equation Solver
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{h}}} \]
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()Projection Solver
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{h}}} = 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 { \hspace{ 0.02pt } {T} }_{,0}^{ 2 }(\mathbf{x}) - { \hspace{ 0.02pt } {T} }_{,1}^{ 2 }(\mathbf{x}) + 1 & - 2 { \hspace{ 0.02pt } {T} }_{,0}(\mathbf{x}) { \hspace{ 0.02pt } {T} }_{,1}(\mathbf{x})\\- 2 { \hspace{ 0.02pt } {T} }_{,0}(\mathbf{x}) { \hspace{ 0.02pt } {T} }_{,1}(\mathbf{x}) & - { \hspace{ 0.02pt } {T} }_{,0}^{ 2 }(\mathbf{x}) - 3 { \hspace{ 0.02pt } {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.