Notebook 8: Particle Swarms

We used a particle swarm to track the flow in Example 7. We called this a “passive” swarm because the points did not influence the flow in any way but were simply carried along with the fluid.

Particle swarms are unstructured data objects that live within the computational domain. Their points can be moved arbitrarily through the domain and points may migrate from one process to another when the coordinates are changed. By default they carry only the particle location, but we can add scalar, vector and tensor variables to the swarm and they will be transported with the particles.

Particle transport is usually through a velocity or displacement field that incrementally changes the locations. This is a common use, but particles can be used to represent any unstructured field. For example, during mesh adaptation, the nodal points from the previous mesh are equivalent to a disconnected swarm from the point of view of the new mesh. The same is true when reading data save from one mesh to the MeshVariables on another.

import numpy as np
import sympy
import underworld3 as uw

res = 12
width = 8

mesh = uw.meshing.UnstructuredSimplexBox(
    cellSize=1 / res,
    minCoords=(0.0, 0.0),
    maxCoords=(width, 1.0),
    qdegree=3,
)

# Coordinate directions etc
x, y = mesh.CoordinateSystem.X

# original y coordinate

y0 = uw.discretisation.MeshVariable(
    "Y0", mesh, vtype=uw.VarType.SCALAR, varsymbol=r"y_0"
)

with mesh.access(y0):
    y0.data[:, 0] = uw.function.evaluate(y, y0.coords)

Deform the mesh

Move all nodes down to meet an undulating lower surface. The displacement field is smooth and continuous, so there is no particular need to remesh in this case. However, it is generally better to consider either deforming the mesh with gmsh before triangulation, or remeshing (both are possible with underworld3, but not as simple to demonstrate).

new_coords = mesh.data
dy = (1 - y) * (sympy.sin(sympy.pi * x) / 10)
new_coords[:, 1] = uw.function.evaluate(y - dy, mesh.data)

display(dy)

mesh.deform_mesh(new_coords)

\(\displaystyle \frac{\left(1 - \mathrm{y}\right) \sin{\left(\mathrm{x} \pi \right)}}{10}\)

# y coordinate after 1st deformation
y1 = uw.discretisation.MeshVariable(
    "Y1", mesh, vtype=uw.VarType.SCALAR, varsymbol=r"y_1"
)

with mesh.access(y1):
    y1.data[:, 0] = uw.function.evaluate(y, y1.coords)

Variable with name Y1 already exists on the mesh - Skipping.

Deform the mesh again

Now we tilt the lower surface down. The variable y1 is negative in regions where the surface was low in the deformed mesh so it is a good first approximation to where ponds of dense fluid might form (e.g. water under ice).

new_coords = mesh.data
dy = (1 - y) * x / 16
new_coords[:, 1] = uw.function.evaluate(y - dy, mesh.data)

mesh.deform_mesh(new_coords)

if uw.mpi.size == 1:
    import pyvista as pv
    import underworld3.visualisation as vis

    pvmesh = vis.mesh_to_pv_mesh(mesh)
    pvmesh.point_data["y0"] = vis.scalar_fn_to_pv_points(pvmesh, y0.sym)
    pvmesh.point_data["y1"] = vis.scalar_fn_to_pv_points(pvmesh, y1.sym)

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

    pl.add_mesh(
        pvmesh,
        scalars="y1",
        cmap="RdBu_r",
        edge_color="Grey",
        edge_opacity=0.33,
        show_edges=True,
        use_transparency=False,
        opacity=1,
        show_scalar_bar=True,
    )

    pl.camera.zoom(3)

Add a swarm to this mesh

A swarm object can be constructed either by adding local points (Example 7) or by filling the mesh with a given density. The density value that we provide (fill_param) adds particles on the Gaussian integration points: 0 will populate the centroids of the mesh elements A value of 1 provides three points per triangle, four per quad, four in a tetrahedron, eight in a hexahedron (the points that support linear interpolation in standard FEM).

Note: in the current implementation of underworld swarm, we an only start adding points when the swarm variables are all finalised.

swarm = uw.swarm.Swarm(mesh)

s = uw.swarm.SwarmVariable(
    "S",
    swarm,
    vtype=uw.VarType.SCALAR,
    proxy_degree=1,
    proxy_continuous=True,
    varsymbol=r"\Sigma",
)

## It is a limitation of the swarm that we
## can only start adding points when the swarm variables
## are all finalised.

swarm.populate(0)

## Set some values on s

with swarm.access(s):
    s.data[:, 0] = uw.function.evaluate(y1.sym, swarm.particle_coordinates.data)

Variable with name proxy_S already exists on the mesh - Skipping.

Proxy Variables

Swarm variables are completely unstructured (they are not located on a grid and they also have no neighbour information). We want to be able to mix their values with mesh-based variables in sympy expressions. To do this, each swarm variable has (by default) a proxy mesh variable which is constructed through a projection to nodal point values. It is accessed as SwarmVariable.sym and, symbolically, it is bracketed \(\left< \cdot \right>\) as a reminder that there is a projection in between the values and the representation on the mesh.

Resolution of proxy variables

The default projection method from Lagrangian swarm variables onto nodal point variables is an inverse-distance weighted average computed using the underworld radial basis function module. This is robust and reproducible, but does not preserve moments of the particle field. You can explore the effect of different swarm densities on the fidelity of the mesh representation by adjusting the number of particles in the swarm (in the cell above) and checking the resulting fields.

s.sym + y1.sym

\(\displaystyle \left[\begin{matrix}{ \hspace{ 0.02pt } {\left<\Sigma\right>} }(\mathbf{x}) + { \hspace{ 0.02pt } {y_1} }(\mathbf{x})\end{matrix}\right]\)

if uw.mpi.size == 1:
    import pyvista as pv
    import underworld3.visualisation as vis

    pvmesh = vis.mesh_to_pv_mesh(mesh)
    pvmesh.point_data["y0"] = vis.scalar_fn_to_pv_points(pvmesh, y0.sym)
    pvmesh.point_data["y1"] = vis.scalar_fn_to_pv_points(pvmesh, y1.sym)
    pvmesh.point_data["S"] = vis.scalar_fn_to_pv_points(pvmesh, s.sym)
    pvmesh.point_data["delta"] = vis.scalar_fn_to_pv_points(pvmesh, s.sym - y1.sym)

    swarm_points = vis.swarm_to_pv_cloud(swarm)

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

    pl.add_mesh(
        pvmesh,
        cmap="RdBu_r",
        edge_color="Grey",
        scalars="S",
        edge_opacity=0.33,
        show_edges=True,
        use_transparency=False,
        opacity=0.75,
        show_scalar_bar=True,
    )

    pl.add_points(
        swarm_points.points,
        point_size=3,
        style="points",
        color="Black",
        opacity=1,
    )

    pl.export_html(f"html5/terrain_flow_mesh_only.html")

from IPython.display import IFrame

IFrame(src=f"html5/terrain_flow_mesh_only.html", width=750, height=300)

Interactive Image: Convection model output

# Mesh variables for the unknowns

v_soln = uw.discretisation.MeshVariable("V0", mesh, 2, degree=2, varsymbol=r"{v_0}")
p_soln = uw.discretisation.MeshVariable("p", mesh, 1, degree=1, continuous=True)

fixed_base = True


stokes = uw.systems.Stokes(
    mesh,
    velocityField=v_soln,
    pressureField=p_soln,
)

stokes.bodyforce = sympy.Matrix((0, 0))

stokes.constitutive_model = uw.constitutive_models.ViscousFlowModel
stokes.constitutive_model.Parameters.shear_viscosity_0 = sympy.Piecewise(
    (sympy.sympify(1) / 100, s.sym[0] < 0.05),
    (
        1,
        True,
    ),
)

stokes.add_essential_bc((1, 0), "Top")

if fixed_base == True:
    stokes.add_essential_bc((0, 0), "Bottom")
else:

    # Free slip on the (deformed) base below a certain Level:
    # 1) Penalise normal velocity below
    # 2) Penalise the entire velocity vector above

    Gamma = mesh.Gamma
    GammaNorm = uw.function.expression(
        r"|\Gamma|", sympy.sqrt(Gamma.dot(Gamma)), "Scaling for surface normals"
    )
    GammaN = Gamma / GammaNorm

    stokes.add_natural_bc(10000 * v_soln.sym.dot(GammaN) * GammaN, "Bottom")

stokes.tolerance = 1.0e-4

stokes.constitutive_model.Parameters.shear_viscosity_0.sym

\(\displaystyle \begin{cases} \frac{1}{100} & \text{for}\: { \hspace{ 0.02pt } {\left<\Sigma\right>} }(\mathbf{x}) < 0.05 \\1 & \text{otherwise} \end{cases}\)

stokes.solve()

if uw.mpi.size == 1:
    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)
    pvmesh.point_data["eta"] = vis.scalar_fn_to_pv_points(
        pvmesh, stokes.constitutive_model.Parameters.shear_viscosity_0.sym
    )
    pvmesh.point_data["V"] = vis.vector_fn_to_pv_points(pvmesh, v_soln.sym)

    pvswarm = vis.swarm_to_pv_cloud(swarm)
    pvswarm.point_data["eta"] = vis.scalar_fn_to_pv_points(
        pvswarm, stokes.constitutive_model.Parameters.shear_viscosity_0.sym
    )

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

    pl.add_mesh(
        pvmesh,
        cmap="RdBu_r",
        edge_color="Grey",
        scalars="eta",
        edge_opacity=0.33,
        show_edges=True,
        use_transparency=False,
        opacity=0.05,
        show_scalar_bar=False,
    )

    pl.add_points(
        pvswarm,
        scalars="eta",
        cmap="Blues_r",
        opacity=0.33,
        point_size=5,
    )

    pl.add_arrows(
        pvmesh.points,
        pvmesh.point_data["V"],
        cmap="Greens",
        mag=0.5,
        show_scalar_bar=False,
    )

    pl.export_html(f"html5/terrain_flow_plot.html")

from IPython.display import IFrame

IFrame(src=f"html5/terrain_flow_plot.html", width=750, height=300)

Interactive Image: Flow model output

Exercise - 8.1

Look at the Stokes example and try with a free slip base and low / high viscosity for the trapped material. You could also see what mixing the boundary condition looks like: free slip in the troughs, fixed in the highs. This is possible with a single boundary condition if you use sympy.Piecewise to switch between different penalty conditions.

    ## Free slip base (conditional)

    Gamma = mesh.Gamma
    GammaNorm = uw.function.expression(r"|\Gamma|", sympy.sqrt(Gamma.dot(Gamma)), "Scaling for surface normals")
    GammaN = Gamma / GammaNorm
    
    bc_mask = sympy.Piecewise((1.0, y1.sym[0] < -0.05), (0.0, True))
    nbc = 10000 * bc_mask * GammaN.dot(v_soln.sym) *  GammaN + (1-bc_mask) * v_soln.sym            
    stokes.add_natural_bc(nbc, "Lower")

Give that a try !