Notebook 3: Symbolic forms

Underworld is deeply integrated with sympy (www.sympy.org) so that any mesh variable can also be used in a sympy expression. We already saw sympy expressions for the coordinates and coordinate directions.

In the examples below, we use a simple 2D, Cartesian mesh because it is much simpler to see the various changes.

import underworld3 as uw
import numpy as np
import sympy

mesh = uw.meshing.UnstructuredSimplexBox(
    minCoords = (-1.0, -1.0),
    maxCoords = (+1.0, +1.0),
    cellSize = 0.05,
    regular=True,
    verbose=False,
)

x,y = mesh.CoordinateSystem.X

As before, we add discrete variables

# mesh variable example / test

s = uw.discretisation.MeshVariable(
    varname="S",
    mesh=mesh, 
    vtype = uw.VarType.SCALAR,
    varsymbol=r"\cal{S}"
)

t = uw.discretisation.MeshVariable(
    varname="T",
    mesh=mesh, 
    vtype = uw.VarType.SCALAR,
    varsymbol=r"\cal{T}"
)

v = uw.discretisation.MeshVariable(
    varname="V1",
    mesh=mesh, 
    degree=2,
    vtype = uw.VarType.VECTOR,
    varsymbol=r"\mathbf{v}",
)

Symbolic forms, derivatives

Variables can be part of complicated sympy expressions. It is important to note that all symbols are matrices and sympy can be fussy when it comes to operations with other matrices (scalars are not entirely equivalent to \(1 \times 1\) matrices).

s.sym[0]+t.sym[0] + v.sym[0]

\(\displaystyle {\cal{S}}(\mathbf{x}) + {\cal{T}}(\mathbf{x}) + {\mathbf{v}}_{ 0 }(\mathbf{x})\)

Derivatives can be handled explicitly, but the mesh also provides vector operators and these are generally better because they are automatically consistent with the underlying coordinate system for the mesh.

For compound expressions of variables, use mesh.vector.curl(expression) but for individual variables, variable.curl() is an equivalent shorthand.

# grad by hand
s.sym[0].diff(x) + s.sym[0].diff(y)

\(\displaystyle {\cal{S}}_{,0}(\mathbf{x}) + {\cal{S}}_{,1}(\mathbf{x})\)

# grad
s.gradient()

\(\displaystyle \left[\begin{matrix}{\cal{S}}_{,0}(\mathbf{x}) & {\cal{S}}_{,1}(\mathbf{x})\end{matrix}\right]\)

v.curl()

\(\displaystyle - {\mathbf{v}}_{ 0,1}(\mathbf{x}) + {\mathbf{v}}_{ 1,0}(\mathbf{x})\)

# curl
mesh.vector.curl(s.sym * v.sym)

\(\displaystyle - {\cal{S}}(\mathbf{x}) {\mathbf{v}}_{ 0,1}(\mathbf{x}) + {\cal{S}}(\mathbf{x}) {\mathbf{v}}_{ 1,0}(\mathbf{x}) + {\cal{S}}_{,0}(\mathbf{x}) {\mathbf{v}}_{ 1 }(\mathbf{x}) - {\cal{S}}_{,1}(\mathbf{x}) {\mathbf{v}}_{ 0 }(\mathbf{x})\)

# v dot grad (scalar)... 
v.sym.dot(mesh.vector.gradient(s.sym))

\(\displaystyle {\cal{S}}_{,0}(\mathbf{x}) {\mathbf{v}}_{ 0 }(\mathbf{x}) + {\cal{S}}_{,1}(\mathbf{x}) {\mathbf{v}}_{ 1 }(\mathbf{x})\)

Underworld Expressions

We often want to define symbols that represent complicated expressions which we want do not want to expand when we probe the mathematical formulation.

An example might be a constitutive model that has a number of conditional expressions, or a concept such as a timestep which we want to refer to as \(\delta t\) regardless of its current numerical value.

Underworld expressions are objects that have a sympy symbolic representation that is only expanded at the time numerical evaluations are required. How about the example above in expression form:

    curl_sv = uw.function.expression(
                    r"\nabla \times \left(\cal{S} \mathbf{v}\right)",
                    mesh.vector.curl(s.sym * v.sym),
                    f"Curl of {v.symbol}"
                ) 

creates an expression object which displays as \(\nabla \times \left(\cal{S} \mathbf{v}\right)\) but which also has curl_sv.sym which is the full expression.

curl_sv = uw.function.expression(
                    r"\nabla \times \left(\cal{S} \mathbf{v}\right)",
                    mesh.vector.curl(s.sym * v.sym),
                    rf"Curl of {v.symbol}"
                ) 

curl_sv.view()
display(curl_sv.sym) 

curl_sv + 1

Class: <class ‘underworld3.function.expressions.UWexpression’>

\(\quad\)\(\displaystyle { \nabla \times \left(\cal{S} \mathbf{v}\right) \hspace{ 0.0pt } }\)\(=\)\(\displaystyle - {\cal{S}}(\mathbf{x}) {\mathbf{v}}_{ 0,1}(\mathbf{x}) + {\cal{S}}(\mathbf{x}) {\mathbf{v}}_{ 1,0}(\mathbf{x}) + {\cal{S}}_{,0}(\mathbf{x}) {\mathbf{v}}_{ 1 }(\mathbf{x}) - {\cal{S}}_{,1}(\mathbf{x}) {\mathbf{v}}_{ 0 }(\mathbf{x})\)

\(\quad\)Description: Curl of {}

\(\displaystyle - {\cal{S}}(\mathbf{x}) {\mathbf{v}}_{ 0,1}(\mathbf{x}) + {\cal{S}}(\mathbf{x}) {\mathbf{v}}_{ 1,0}(\mathbf{x}) + {\cal{S}}_{,0}(\mathbf{x}) {\mathbf{v}}_{ 1 }(\mathbf{x}) - {\cal{S}}_{,1}(\mathbf{x}) {\mathbf{v}}_{ 0 }(\mathbf{x})\)

\(\displaystyle { \nabla \times \left(\cal{S} \mathbf{v}\right) \hspace{ 0.0pt } } + 1\)

Underworld Sub-expressions

The viscosity for a material with a yield stress might look like this:

\[ \eta_\textrm{eff} = \displaystyle \min\left({ {\eta} }, \frac{\max\left({ {\tau_{y}} }, { {\tau_{y, \mathrm{min}}} }\right)}{2 { \dot\varepsilon_{II} }}\right) \]

where \(\tau_y\) is a scalar yield stress, and \(\dot\varepsilon_{II}\) is the second invariant of the strain-rate tensor. Each of these would be a nested Underworld expression.

To expand the expression and see inside, we can use

    expression.unwrap(keep_constants=True)

which reaches down the hierarchy and recursively replaces any expression with its expanded sympy representation. If a symbol represents a constant value (float, int, or sympy number), then it is only replaced if keep_constants is set to False.

Symbolic forms can be evaluated at points in the (meshed domain) using uw.function.evaluate. Pure sympy functions can be used to set values in the data container of a meshVariable object

with mesh.access(s, t):
    s.data[:,0] = uw.function.evaluate(sympy.cos(3 * sympy.pi * x)**2 , s.coords)
    t.data[:,0] = uw.function.evaluate(sympy.sin(3 * sympy.pi * y)**2 , t.coords)

s.view()

Class: <class ‘underworld3.discretisation._MeshVariable’>

MeshVariable:

symbol: \({\cal{S}}\)

shape: \((1, 1)\)

degree: \(1\)

continuous: True

type: SCALAR

FE Data:

PETSc field id: \(0\)

PETSc field name: S

array([[1.        ],
       [1.        ],
       [1.        ],
       ...,
       [0.79389263],
       [0.79389263],
       [0.79389263]])

# Visualise it / them

import pyvista as pv
import underworld3.visualisation as vis

pvmesh = vis.mesh_to_pv_mesh(mesh)
pvmesh.point_data["s"] = vis.scalar_fn_to_pv_points(pvmesh, s.sym[0])
pvmesh.point_data["t"] = vis.scalar_fn_to_pv_points(pvmesh, t.sym[0])
pvmesh.point_data["sxt"] = vis.scalar_fn_to_pv_points(pvmesh, s.sym[0]*t.sym[0])

pvmesh.warp_by_scalar(scalars="sxt", factor=0.3, normal=(0,0,1), inplace=True)

# pvmesh.plot(show_edges=True, show_scalar_bar=False)

pl = pv.Plotter(window_size=(750, 750))

pl.add_mesh(pvmesh, 
            show_edges=True,
            edge_color="#4455FF",
            cmap="Greys",
            scalars="sxt", 
            show_scalar_bar=False)

# Save and show the mesh

pl.camera_position = 'yz'
pl.camera.azimuth = 45
pl.camera.elevation = 45

pl.export_html("html5/sine_squared.html")

from IPython.display import IFrame
IFrame(src="html5/sine_squared.html", width=600, height=400)

Interactive Image: Square mesh of triangular elements on which we evaluated a simple sympy function of position

Lazy evaluation of expressions

Generally speaking, we use expressions as placeholders for parameters or functions that we know will be needed when it comes to solve a problem, but we can’t be sure that we can specify them at the time we set up the solver.

In the code for our solvers, for example, we set up a template with expressions that describe where the parameters of the problem will be and we expect you to fill the details in when we create a new solver-object. These template expressions are protected so that assignment just changes the value that the expression holds, it does not change the symbol or the description (though you can do this if you want).

This is a rather contrived example:

R = uw.function.expression(r"r\left( \mathbf{x} \right)",
                           sympy.sqrt(x**2+y**2), 
                           "distance from origin")

R1 = R + 1

S = sympy.sqrt(x**2+y**2)
S1 = S + 1
               

R1 = R + 1
R2 = R.sym + 1
S1 = S + 1

S1

\(\displaystyle \sqrt{\mathrm{x}^{2} + \mathrm{y}^{2}} + 1\)

R.sym = S = sympy.sqrt((x-1)**2 + (y-1)**2)

uw.function.expressions.unwrap(R1) # Lazy evaluation 

\(\displaystyle \sqrt{\left(\mathrm{x} - 1\right)^{2} + \left(\mathrm{y} - 1\right)^{2}} + 1\)

uw.function.expressions.unwrap(S1)

\(\displaystyle \sqrt{\mathrm{x}^{2} + \mathrm{y}^{2}} + 1\)

uw.function.expressions.unwrap(R2)

\(\displaystyle \sqrt{\mathrm{x}^{2} + \mathrm{y}^{2}} + 1\)

Exercise 3.1

Have a look at the visco-plastic constitutive model (template) for Stokes equation

stokes_solver = uw.systems.solvers.SNES_Stokes(mesh)
stokes_solver.constitutive_model = uw.constitutive_models.ViscoPlasticFlowModel
stokes_solver.constitutive_model.Parameters.shear_viscosity_0 = 1
stokes_solver.constitutive_model.Parameters.yield_stress = 100

You can examine this expression in more detail using the view method of the stokes_solver.constitutive_model, and you can expand the expression to see how it reduces when sympy needs to evaluate this expression at one or more locations in the domain.

Exercise 3.2

Assignment to an expression object replaces the sympy value but does not change the rest of the object. This is the concept of lazy evaluation which we introduced earlier.

Validate this using the constititutive model above. Try changing the yield stress or the shear viscosity and see how the expression for the apparent viscosity changes.

e.g. 

    stokes_solver.constitutive_model.Parameters.shear_viscosity_0 = sympy.exp(-10 * t.sym[0])
    display(stokes_solver.constitutive_model.viscosity)

Actually, that won’t look different, will it ? You need to expand out the expressions a bit to see it. Try using stokes_solver.constitutive_model.viscosity.unwrap(). This function makes substitutions of all (underworld) sub-expressions in any sympy expression.

You can pass keep_constants=False if you want to expand all the numerical values as well. An expression is considered to be a constant if it contains no sub-expressions, MeshVariables, or mesh-coordinates. We usually do not want to have long floating point numbers all over the place when we check an expression except if it’s the values that we want to check.

Try it !