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.

meshVariables can be composed with other symbolic objects and evaluated numerically when required. They can also be differentiated. Most importantly, sympy can manipulate expressions, simplify them and cancel terms.

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

#|  echo: false 
# This is required to fix pyvista 
# (visualisation) crashes in interactive notebooks (including on binder)

import nest_asyncio
nest_asyncio.apply()

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). (If you take out the [0] in this case, sympy will complain, but you are allowed to do this: s + t + v[:,0] - try it !)

s[0] + t[0] + v[0]

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.diff(x) + s.diff(y)

s.diff(x)

# grad
s.gradient()

v.curl()

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

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

Underworld Expressions

We often want to define symbols that represent complicated expressions which 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.

Note: like MeshVariables expressions have unique names and are persistent objects. They are containers for expressions (in the way MeshVariables are containers for numerical information. This means that we can write an expression and change what it represents while it remains part of another expression.

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

curl_sv.view()
display(curl_sv.sym) 

curl_sv + 1

# If you run the cell above multiple times, the object id should not change
id(curl_sv)

Underworld Sub-expressions

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

\[ \eta_\textrm{eff} = \displaystyle \min\left({ {\eta_0} }, \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

    uw.unwrap(expression, 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.

viscosity = uw.function.expression(
                    r"\eta_0",
                    sympy.sympify(1),
                    rf"viscosity"
                ) 

tau_yield = uw.function.expression(
                    r"\tau_y",
                    sympy.sympify(10),
                    rf"yield stress"
                ) 

edot_II = uw.function.expression(
                    r"\dot{\varepsilon}_{II}",
                    sympy.Symbol('epsilon'),
                    rf"strain rate invariant"
                ) 


viscosity_eff = uw.function.expression(
                    r"\eta_\mathrm{eff}",
                    sympy.Min(viscosity, tau_yield / (2 * edot_II)),
                    rf"Curl of {v.symbol}"
                )

# Try this with True / False in the keep_constants argument !

uw.unwrap(viscosity_eff, keep_constants=True)

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

s.array[...] = uw.function.evaluate(sympy.cos(3 * sympy.pi * x)**2 , s.coords)
t.array[...] = uw.function.evaluate(sympy.sin(3 * sympy.pi * y)**2 , t.coords)

fn  = sympy.cos(3 * sympy.pi * x)**2 + s.sym[0] + t.sym[0]
a = uw.function.evaluate(fn , np.atleast_2d(mesh._get_domain_centroids()[uw.mpi.rank]))
shape = np.atleast_2d(a).shape

viscosity_eff.sym

viscosity_eff.sym.diff(viscosity)

uw.function.derivative(viscosity_eff, viscosity)

uw.function.derivative(viscosity_eff, viscosity, evaluate=False)

b = uw.function.derivative(viscosity_eff, viscosity, evaluate=False)

b.sym

b.doit()

# 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")

#| fig-cap: "Interactive Image: Square mesh of triangular elements on which we evaluated a simple `sympy` function of position"

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

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
               

R2 = R.sym + 1

S1

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

uw.unwrap(R1) # Lazy evaluation

uw.unwrap(S1)

uw.unwrap(R2)

Derivatives of expressions & Lazy evaluation

sympy can differentiate expressions

viscosity_eff.sym

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 uw.unwrap(stokes_solver.constitutive_model.viscosity). 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 !