Cross product of normal and test function in the variational form using ufl's FacetNormal?

General
,
1

For the variational form of the (inhomogeneous) shallow water equations, I need to find

u^\bot = \textbf{k}\times \textbf{u}

where \textbf{k} is the unit vector pointing normal to the surface. I tried to use ufl.FacetNormal taking inspiration from this example to create the unit vector normal to the mesh k = ufl.FacetNormal(mesh), and from this, I thought simply running u_n = ufl.cross(k,u) would suffice. Running the following compiled

from mpi4py import MPI
import meshio
import numpy as np
import ufl
import basix.ufl
import dolfinx as df
import scifem
from petsc4py import PETSc
from dolfinx.fem.petsc import apply_lifting, assemble_matrix, assemble_vector, set_bc
from Projector import Projector

order = 1
exact_case = 2
tmp_data = meshio.dolfin.read("supplement/examples/mixed-poisson/hdiv-l2/meshes/sphere_ico4.xml")

mesh = df.mesh.create_mesh(
    comm=MPI.COMM_WORLD,
    cells=tmp_data.cells_dict["triangle"],
    x=tmp_data.points,
    e=ufl.Mesh(basix.ufl.element("Lagrange", "triangle", 1, shape=(3,))),
)

def create_function_spaces(mesh, family, order):
    S = df.fem.functionspace(mesh, (family, order))
    V = df.fem.functionspace(mesh, ("DG", order-1))
    PlottingMesh = df.fem.functionspace(mesh, ("CG", order))
    return (S, V, PlottingMesh)

def initial_conditions(S, V):
    u0 = df.fem.Function(S)
    u0.interpolate(lambda x: np.zeros(x[:S.mesh.topology.dim].shape))
    f = ufl.exp(-((-x[1]-1)/0.1)**2)
    petsc_options = {"ksp_type": "preonly", "pc_type": "lu"}
    V_projector = Projector(V, petsc_options=petsc_options)
    D0 = V_projector.project(f)
    return (u0, D0)

# Define global normal
x = ufl.SpatialCoordinate(mesh)
global_normal = x # Not sure exactly what this is used for
family = "RT"
order = 1

# Define mixed mesh
(S,V,PlottingMesh) = create_function_spaces(mesh,family,order)
W = ufl.MixedFunctionSpace(*(S, V))

# Constants
dt = df.default_scalar_type(0.05)
H = 1.0
g = 1.0

# forcing term
f = x[2]

# Initial values of solutions u and D
u_, D_ = initial_conditions(S,V)

# normal unit vector
k = ufl.FacetNormal(mesh)

# u_n = k x u
u_n = ufl.cross(k,u)

F = (ufl.inner(u - u_, w) - dt*ufl.div(w)*g*D_mid + dt*f*ufl.inner(u_n, w)
    + ufl.inner(D - D_, phi) + dt*H*ufl.div(u_mid)*phi)*ufl.dx

(a, L) = ufl.system(F)

where the mesh is a 2D spherical manifold of radius 1 and the specific mesh used is found in the following file (make sure to extract the .gz file containing the mesh), and the code for the projector can be found here.

However, the code breaks when attempting to setup the problem

petsc_options = {
    "ksp_type": "preonly",
    "pc_type": "lu",
    "pc_factor_mat_solver_type": "mumps",
    "ksp_error_if_not_converged": True,
}
problem = df.fem.petsc.LinearProblem(
    ufl.extract_blocks(a),
    L,
    bcs=[],
    petsc_options=petsc_options,
    petsc_options_prefix="mixed_poisson_",
    kind="mpi",
)

with error

ValueError: Integral of type cell cannot contain a ReferenceNormal.

which I presume to mean that the cross product with k being in the variational form is causing the error. Are there things I can do to get around this, or should ufl.FacetNormal generally be avoided for what I’m trying to do?

2

You need the ufl.CellNormal, as you are integrating over cells, not facets.

1 Like
3

Okay that resolves that issue. Now I’m running into another being a timeout error which tracebacked to attempting create_form for L.

Running df.fem.form(a) also causes a timeout with the following traceback

---------------------------------------------------------------------------
FileExistsError                           Traceback (most recent call last)
File /dolfinx-env/lib/python3.12/site-packages/ffcx/codegeneration/jit.py:98, in get_cached_module(module_name, object_names, cache_dir, timeout)
     96 try:
     97     # Create C file with exclusive access
---> 98     with open(c_filename, "x"):
     99         pass

FileExistsError: [Errno 17] File exists: '/root/.cache/fenics/libffcx_forms_ef21e46178e68bfbcd19443e673c35419b88f2a3.c'

During handling of the above exception, another exception occurred:

TimeoutError                              Traceback (most recent call last)
Cell In[16], line 1
----> 1 df.fem.form(a)

File /usr/local/dolfinx-real/lib/python3.12/dist-packages/dolfinx/fem/forms.py:438, in form(form, dtype, form_compiler_options, jit_options, entity_maps)
    435     else:
    436         return form
--> 438 return _create_form(form)

File /usr/local/dolfinx-real/lib/python3.12/dist-packages/dolfinx/fem/forms.py:430, in form.<locals>._create_form(form)
    420 """Recursively convert ufl.Forms to dolfinx.fem.Form.
    421 
    422 Args:
   (...)    427     A ``dolfinx.fem.Form`` or a list of ``dolfinx.fem.Form``.
    428 """
    429 if isinstance(form, ufl.Form):
--> 430     return _form(form)
    431 elif isinstance(form, ufl.ZeroBaseForm):
    432     return _zero_form(form)

File /usr/local/dolfinx-real/lib/python3.12/dist-packages/dolfinx/fem/forms.py:339, in form.<locals>._form(form)
    337 if msh is None:
    338     raise RuntimeError("Expecting to find a Mesh in the form.")
--> 339 ufcx_form, module, code = jit.ffcx_jit(
    340     msh.comm, form, form_compiler_options=form_compiler_options, jit_options=jit_options
    341 )
    343 # For each argument in form extract its function space
    344 V = [arg.ufl_function_space()._cpp_object for arg in form.arguments()]

File /usr/local/dolfinx-real/lib/python3.12/dist-packages/dolfinx/jit.py:61, in mpi_jit_decorator.<locals>.mpi_jit(comm, *args, **kwargs)
     57 @functools.wraps(local_jit)
     58 def mpi_jit(comm, *args, **kwargs):
     59     # Just call JIT compiler when running in serial
     60     if comm.size == 1:
---> 61         return local_jit(*args, **kwargs)
     63     # Default status (0 == ok, 1 == fail)
     64     status = 0

File /usr/local/dolfinx-real/lib/python3.12/dist-packages/dolfinx/jit.py:216, in ffcx_jit(ufl_object, form_compiler_options, jit_options)
    214 # Switch on type and compile, returning cffi object
    215 if isinstance(ufl_object, ufl.Form):
--> 216     r = ffcx.codegeneration.jit.compile_forms([ufl_object], options=p_ffcx, **p_jit)
    217 elif isinstance(ufl_object, tuple) and isinstance(ufl_object[0], ufl.core.expr.Expr):
    218     r = ffcx.codegeneration.jit.compile_expressions([ufl_object], options=p_ffcx, **p_jit)

File /dolfinx-env/lib/python3.12/site-packages/ffcx/codegeneration/jit.py:189, in compile_forms(forms, options, cache_dir, timeout, cffi_extra_compile_args, cffi_verbose, cffi_debug, cffi_libraries, visualise)
    187 if cache_dir is not None:
    188     cache_dir = Path(cache_dir)
--> 189     obj, mod = get_cached_module(module_name, form_names, cache_dir, timeout)
    190     if obj is not None:
    191         return obj, mod, (None, None)

File /dolfinx-env/lib/python3.12/site-packages/ffcx/codegeneration/jit.py:123, in get_cached_module(module_name, object_names, cache_dir, timeout)
    121     logger.info(f"Waiting for {ready_name} to appear.")
    122     time.sleep(1)
--> 123 raise TimeoutError(
    124     "JIT compilation timed out, probably due to a failed previous compile. "
    125     f"Try cleaning cache (e.g. remove {c_filename}) or increase timeout option."
    126 )

TimeoutError: JIT compilation timed out, probably due to a failed previous compile. Try cleaning cache (e.g. remove /root/.cache/fenics/libffcx_forms_ef21e46178e68bfbcd19443e673c35419b88f2a3.c) or increase timeout option.

Do you know what might be causing this?

4

There is so much code missing from your example, that anything could be wrong at this point.

For instance, you never define most of the variables.

You also haven’t specified which version of DOLFINx, ufl etc that you are using.

Here is a code that runs in the nightly docker image:

from mpi4py import MPI
import meshio
import numpy as np
import ufl
import basix.ufl
import dolfinx as df
from petsc4py import PETSc

from typing import Optional

import numpy as np
import dolfinx.fem.petsc


class Projector:
    """
    Projector for a given function.
    Solves Ax=b, where

    .. highlight:: python
    .. code-block:: python

        u, v = ufl.TrialFunction(Space), ufl.TestFunction(space)
        dx = ufl.Measure("dx", metadata=metadata)
        A = inner(u, v) * dx
        b = inner(function, v) * dx(metadata=metadata)

    Args:
        function: UFL expression of function to project
        space: Space to project function into
        petsc_options: Options to pass to PETSc
        jit_options: Options to pass to just in time compiler
        form_compiler_options: Options to pass to the form compiler
        metadata: Data to pass to the integration measure
    """

    _A: PETSc.Mat  # The mass matrix
    _b: PETSc.Vec  # The rhs vector
    _lhs: dolfinx.fem.Form  # The compiled form for the mass matrix
    _ksp: PETSc.KSP  # The PETSc solver
    _x: dolfinx.fem.Function  # The solution vector
    _dx: ufl.Measure  # Integration measure

    def __init__(
        self,
        space: dolfinx.fem.FunctionSpace,
        petsc_options: Optional[dict] = None,
        jit_options: Optional[dict] = None,
        form_compiler_options: Optional[dict] = None,
        metadata: Optional[dict] = None,
    ):
        petsc_options = {} if petsc_options is None else petsc_options
        jit_options = {} if jit_options is None else jit_options
        form_compiler_options = (
            {} if form_compiler_options is None else form_compiler_options
        )

        # Assemble projection matrix once
        u = ufl.TrialFunction(space)
        v = ufl.TestFunction(space)
        self._dx = ufl.Measure("dx", domain=space.mesh, metadata=metadata)
        a = ufl.inner(u, v) * self._dx(metadata=metadata)
        self._lhs = dolfinx.fem.form(
            a, jit_options=jit_options, form_compiler_options=form_compiler_options
        )
        self._A = dolfinx.fem.petsc.assemble_matrix(self._lhs)
        self._A.assemble()

        # Create vectors to store right hand side and the solution
        self._x = dolfinx.fem.Function(space)
        self._b = dolfinx.fem.Function(space)

        # Create Krylov Subspace solver
        self._ksp = PETSc.KSP().create(space.mesh.comm)
        self._ksp.setOperators(self._A)

        # Set PETSc options
        prefix = f"projector_{id(self)}"
        opts = PETSc.Options()
        opts.prefixPush(prefix)
        for k, v in petsc_options.items():
            opts[k] = v
        opts.prefixPop()
        self._ksp.setFromOptions()
        for opt in opts.getAll().keys():
            del opts[opt]

        # Set matrix and vector PETSc options
        self._A.setOptionsPrefix(prefix)
        self._A.setFromOptions()
        self._b.x.petsc_vec.setOptionsPrefix(prefix)
        self._b.x.petsc_vec.setFromOptions()

    def reassemble_lhs(self):
        dolfinx.fem.petsc.assemble_matrix(self._A, self._lhs)
        self._A.assemble()

    def assemble_rhs(self, h: ufl.core.expr.Expr):
        """
        Assemble the right hand side of the problem
        """
        v = ufl.TestFunction(self._b.function_space)
        rhs = ufl.inner(h, v) * self._dx
        rhs_compiled = dolfinx.fem.form(rhs)
        self._b.x.array[:] = 0.0
        dolfinx.fem.petsc.assemble_vector(self._b.x.petsc_vec, rhs_compiled)
        self._b.x.petsc_vec.ghostUpdate(
            addv=PETSc.InsertMode.ADD_VALUES, mode=PETSc.ScatterMode.REVERSE
        )
        self._b.x.petsc_vec.ghostUpdate(
            addv=PETSc.InsertMode.INSERT_VALUES, mode=PETSc.ScatterMode.FORWARD
        )

    def project(self, h: ufl.core.expr.Expr) -> dolfinx.fem.Function:
        """
        Compute projection using a PETSc KSP solver

        Args:
            assemble_rhs: Re-assemble RHS and re-apply boundary conditions if true
        """
        self.assemble_rhs(h)
        self._ksp.solve(self._b.x.petsc_vec, self._x.x.petsc_vec)
        return self._x

    def __del__(self):
        self._A.destroy()
        self._ksp.destroy()


order = 1
exact_case = 2
tmp_data = meshio.dolfin.read(
    "supplement/examples/mixed-poisson/hdiv-l2/meshes/sphere_ico4.xml"
)

mesh = df.mesh.create_mesh(
    comm=MPI.COMM_WORLD,
    cells=tmp_data.cells_dict["triangle"],
    x=tmp_data.points,
    e=ufl.Mesh(basix.ufl.element("Lagrange", "triangle", 1, shape=(3,))),
)


def create_function_spaces(mesh, family, order):
    S = df.fem.functionspace(mesh, (family, order))
    V = df.fem.functionspace(mesh, ("DG", order - 1))
    PlottingMesh = df.fem.functionspace(mesh, ("CG", order))
    return (S, V, PlottingMesh)


def initial_conditions(S, V):
    u0 = df.fem.Function(S)
    u0.x.array[:] = 0.0
    f = ufl.exp(-(((-x[1] - 1) / 0.1) ** 2))
    petsc_options = {"ksp_type": "preonly", "pc_type": "lu"}
    V_projector = Projector(V, petsc_options=petsc_options)
    D0 = V_projector.project(f)
    return (u0, D0)


# Define global normal
x = ufl.SpatialCoordinate(mesh)
global_normal = x  # Not sure exactly what this is used for
family = "RT"
order = 1

# Define mixed mesh
(S, V, PlottingMesh) = create_function_spaces(mesh, family, order)
W = ufl.MixedFunctionSpace(*(S, V))

# Constants
dt = df.default_scalar_type(0.05)
H = 1.0
g = 1.0

# forcing term
f = x[2]

# Initial values of solutions u and D
u_, D_ = initial_conditions(S, V)
w, phi = ufl.TestFunctions(W)
D = dolfinx.fem.Function(V)
u = dolfinx.fem.Function(S)

D_mid = 0.5 * (D + D_)
u_mid = 0.5 * (u_ + u)
# normal unit vector
k = ufl.CellNormal(mesh)


# u_n = k x u
u_n = ufl.cross(k, u)

F = (
    ufl.inner(u - u_, w)
    - dt * ufl.div(w) * g * D_mid
    + dt * f * ufl.inner(u_n, w)
    + ufl.inner(D - D_, phi)
    + dt * H * ufl.div(u_mid) * phi
) * ufl.dx

(a, L) = ufl.system(F)

For further questions, please make an effort in making your example reproducible.

5

You’re right I realize I forgot to copy paste my code segment and run to make sure it worked before posting, otherwise I would have seen that I forgot to define many variables. It’s roughly 4 AM where I am, I should have saved my questions for the next day and I would have less likely made these mistakes.

And yes, not everyone here knows that I’m on the nightly version, so specifying my version should be something that I’ll make known in the future. I will add that to my question in an edit.

Yes the code you’ve created fills in several of the gaps. The only difference is the definitions of u and D

(u, D) = ufl.TrialFunctions(W)
(w, phi) = ufl.TestFunctions(W)

So the full code from start to finish for the code which successfully compiles is

from mpi4py import MPI
import meshio
import numpy as np
import ufl
import basix.ufl
import dolfinx as df
from petsc4py import PETSc

from typing import Optional

import numpy as np
import dolfinx.fem.petsc


class Projector:
    """
    Projector for a given function.
    Solves Ax=b, where

    .. highlight:: python
    .. code-block:: python

        u, v = ufl.TrialFunction(Space), ufl.TestFunction(space)
        dx = ufl.Measure("dx", metadata=metadata)
        A = inner(u, v) * dx
        b = inner(function, v) * dx(metadata=metadata)

    Args:
        function: UFL expression of function to project
        space: Space to project function into
        petsc_options: Options to pass to PETSc
        jit_options: Options to pass to just in time compiler
        form_compiler_options: Options to pass to the form compiler
        metadata: Data to pass to the integration measure
    """

    _A: PETSc.Mat  # The mass matrix
    _b: PETSc.Vec  # The rhs vector
    _lhs: dolfinx.fem.Form  # The compiled form for the mass matrix
    _ksp: PETSc.KSP  # The PETSc solver
    _x: dolfinx.fem.Function  # The solution vector
    _dx: ufl.Measure  # Integration measure

    def __init__(
        self,
        space: dolfinx.fem.FunctionSpace,
        petsc_options: Optional[dict] = None,
        jit_options: Optional[dict] = None,
        form_compiler_options: Optional[dict] = None,
        metadata: Optional[dict] = None,
    ):
        petsc_options = {} if petsc_options is None else petsc_options
        jit_options = {} if jit_options is None else jit_options
        form_compiler_options = (
            {} if form_compiler_options is None else form_compiler_options
        )

        # Assemble projection matrix once
        u = ufl.TrialFunction(space)
        v = ufl.TestFunction(space)
        self._dx = ufl.Measure("dx", domain=space.mesh, metadata=metadata)
        a = ufl.inner(u, v) * self._dx(metadata=metadata)
        self._lhs = dolfinx.fem.form(
            a, jit_options=jit_options, form_compiler_options=form_compiler_options
        )
        self._A = dolfinx.fem.petsc.assemble_matrix(self._lhs)
        self._A.assemble()

        # Create vectors to store right hand side and the solution
        self._x = dolfinx.fem.Function(space)
        self._b = dolfinx.fem.Function(space)

        # Create Krylov Subspace solver
        self._ksp = PETSc.KSP().create(space.mesh.comm)
        self._ksp.setOperators(self._A)

        # Set PETSc options
        prefix = f"projector_{id(self)}"
        opts = PETSc.Options()
        opts.prefixPush(prefix)
        for k, v in petsc_options.items():
            opts[k] = v
        opts.prefixPop()
        self._ksp.setFromOptions()
        for opt in opts.getAll().keys():
            del opts[opt]

        # Set matrix and vector PETSc options
        self._A.setOptionsPrefix(prefix)
        self._A.setFromOptions()
        self._b.x.petsc_vec.setOptionsPrefix(prefix)
        self._b.x.petsc_vec.setFromOptions()

    def reassemble_lhs(self):
        dolfinx.fem.petsc.assemble_matrix(self._A, self._lhs)
        self._A.assemble()

    def assemble_rhs(self, h: ufl.core.expr.Expr):
        """
        Assemble the right hand side of the problem
        """
        v = ufl.TestFunction(self._b.function_space)
        rhs = ufl.inner(h, v) * self._dx
        rhs_compiled = dolfinx.fem.form(rhs)
        self._b.x.array[:] = 0.0
        dolfinx.fem.petsc.assemble_vector(self._b.x.petsc_vec, rhs_compiled)
        self._b.x.petsc_vec.ghostUpdate(
            addv=PETSc.InsertMode.ADD_VALUES, mode=PETSc.ScatterMode.REVERSE
        )
        self._b.x.petsc_vec.ghostUpdate(
            addv=PETSc.InsertMode.INSERT_VALUES, mode=PETSc.ScatterMode.FORWARD
        )

    def project(self, h: ufl.core.expr.Expr) -> dolfinx.fem.Function:
        """
        Compute projection using a PETSc KSP solver

        Args:
            assemble_rhs: Re-assemble RHS and re-apply boundary conditions if true
        """
        self.assemble_rhs(h)
        self._ksp.solve(self._b.x.petsc_vec, self._x.x.petsc_vec)
        return self._x

    def __del__(self):
        self._A.destroy()
        self._ksp.destroy()


order = 1
exact_case = 2
tmp_data = meshio.dolfin.read(
    "supplement/examples/mixed-poisson/hdiv-l2/meshes/sphere_ico4.xml"
)

mesh = df.mesh.create_mesh(
    comm=MPI.COMM_WORLD,
    cells=tmp_data.cells_dict["triangle"],
    x=tmp_data.points,
    e=ufl.Mesh(basix.ufl.element("Lagrange", "triangle", 1, shape=(3,))),
)


def create_function_spaces(mesh, family, order):
    S = df.fem.functionspace(mesh, (family, order))
    V = df.fem.functionspace(mesh, ("DG", order - 1))
    PlottingMesh = df.fem.functionspace(mesh, ("CG", order))
    return (S, V, PlottingMesh)


def initial_conditions(S, V):
    u0 = df.fem.Function(S)
    u0.x.array[:] = 0.0
    f = ufl.exp(-(((-x[1] - 1) / 0.1) ** 2))
    petsc_options = {"ksp_type": "preonly", "pc_type": "lu"}
    V_projector = Projector(V, petsc_options=petsc_options)
    D0 = V_projector.project(f)
    return (u0, D0)


# Define global normal
x = ufl.SpatialCoordinate(mesh)
global_normal = x  # Not sure exactly what this is used for
family = "RT"
order = 1

# Define mixed mesh
(S, V, PlottingMesh) = create_function_spaces(mesh, family, order)
W = ufl.MixedFunctionSpace(*(S, V))

# Constants
dt = df.default_scalar_type(0.05)
H = 1.0
g = 1.0

# forcing term
f = x[2]

# Initial values of solutions u and D
u_, D_ = initial_conditions(S, V)
(u, D) = ufl.TrialFunctions(W)
(w, phi) = ufl.TestFunctions(W)

D_mid = 0.5 * (D + D_)
u_mid = 0.5 * (u_ + u)
# normal unit vector
k = ufl.CellNormal(mesh)


# u_n = k x u
u_n = ufl.cross(k, u)

F = (
    ufl.inner(u - u_, w)
    - dt * ufl.div(w) * g * D_mid
    + dt * f * ufl.inner(u_n, w)
    + ufl.inner(D - D_, phi)
    + dt * H * ufl.div(u_mid) * phi
) * ufl.dx

(a, L) = ufl.system(F)

Sorry again for flubbing that so badly.

1 Like
6

My current versions are

fenics-dolfinx             0.10.0.dev0
fenics-ufl                 2025.2.0.dev0

which is accessible on the Dolfinx nightly version

7

To avoid having to grab the mesh from another file, I’ve decided to write the code to create a generic spherical 2D manifold of radius 1 to work on the problem with using gmsh, and provide some needed additional context for the variational form.

import gmsh
import dolfinx as df
from mpi4py import MPI

gmsh.initialize()
gmsh.model.occ.add_sphere(0,0,0,radius = 1, tag = 1)
gmsh.model.occ.synchronize()
## Add Physical Group (important step for Dolfinx meshing)
areas = gmsh.model.getEntities(dim=2) 
fluid_marker = 11 # What do you want the fluid marker to be named?
gmsh.model.addPhysicalGroup(areas[0][0], [areas[0][1]], fluid_marker) # Create the physical group
gmsh.model.setPhysicalName(areas[0][0], fluid_marker, "Fluid area")

## Mesh refinement
res = 1/10
gmsh.option.setNumber("Mesh.CharacteristicLengthMin", res)
gmsh.option.setNumber("Mesh.CharacteristicLengthMax", res)

## Generate mesh
gmsh.model.occ.synchronize()
gmsh.model.mesh.generate(dim = 2)

# Create Dolfinx mesh
model_rank = 0
mesh = df.io.gmshio.model_to_mesh(gmsh.model, MPI.COMM_WORLD, model_rank)[0]

On this spherical mesh, everything still compiles up to creation of the variational form, however I still get the same timeout issue when attempting to create

bilinear_form = df.fem.form(a)
linear_form = df.fem.form(L)

necessary for using dolfinx.fem.petsc.assemble_matrix() on the bilinear form for latter computations.

Now, to provide additional context, the system of shallow water equations which need to be solved are

(1)\frac{\partial \textbf{u}}{\partial t} + f \textbf{u}^\perp = -g \nabla D

(2)\frac{\partial D}{\partial t} + H \nabla \cdot \textbf{u} = 0

where g and H are constants, and \textbf{u} and D are vector and scalar solutions, respectively, and I am following the following code written in legacy Dolfin which solves these equations

"""
Linear shallow water code for mimetic element pairs.

The equations are the linear shallow water equations.  Let u be
velocity, D be depth _perturbation_, g = gravity, H = base layer depth

du/dt + fu^\perp = -g grad(D)
dD/dt + H div(u) = 0

In the continuous case, if we pick a velocity field u = curl(psi),
then div(u) = 0, so dD/dt is initially 0.  If you then choose D to
satisfy fu^\perp = -g grad(D), du/dt is initially 0 also.  Then both u
and D remain fixed for all time.  This is known as "geostrophic
balance".  One of the special properties of the spatial discretisation
is that it represents this balance.

With implicit midpoint timestepping, the scheme conserves energy -
this is verified in.
"""

# Copyright (C) 2013 Andrew McRae and Imperial College London
# Written by Andrew McRae and Colin J. Cotter
# Modified by Marie E. Rognes (meg@simula.no), 2013

from dolfin import *

parameters["form_compiler"]["cpp_optimize"] = True
parameters["form_compiler"]["optimize"] = True

outward_normals_code = """
class OutwardNormal : public Expression
{
public:

boost::shared_ptr<Mesh> mesh;

OutwardNormal() : Expression(3)
{
}

void eval(Array<double>& values, const Array<double>& x,
          const ufc::cell& c) const
  {
   assert(mesh);
   const Cell cell(*mesh, c.index);

   // Extract local normal and orientation (already computed)
   const Point& local_normal(cell.cell_normal());
   const int orientation = mesh->cell_orientations()[c.index];

   // Set values (1 means flipped, so then we flip back)
   if (orientation == 0)
   {
     values[0] = local_normal.x();
     values[1] = local_normal.y();
     values[2] = local_normal.z();
   } else
   {
     values[0] = - local_normal.x();
     values[1] = - local_normal.y();
     values[2] = - local_normal.z();
   }
  }
};"""

def initialize_parameters(verbose=True):
    # Create application parameters
    params = Parameters("Discretization")
    params.add("family", "RT")
    params.add("order", 1)
    params.add("meshid", "sphere_ico4")
    params.add("store_pvds", True)
    params.add("T", 10.0)
    params.add("dt", 0.05)
    params.add("directory", "default")

    # Read parameters from the command-line
    params.parse()
    if params["directory"] == "default":
        default = "%s%d_%s" % (params["family"], params["order"],
                               params["meshid"])
        params["directory"] = default

    info(params, verbose)

    # Return parameters
    return params

def create_function_spaces(mesh, family, order):
    S = FunctionSpace(mesh, family, order)
    V = FunctionSpace(mesh, "DG", order-1)
    PlottingMesh = VectorFunctionSpace(mesh, "CG", order)
    return (S, V, PlottingMesh)

def myperp(u):
    return cross(en_exp, u)

def initial_conditions(S, V):
    u0 = interpolate(Constant((0.0,0.0,0.0)), S)
    D0 = project(Expression("exp(-pow((-x[1]-1.0)/0.1, 2))"), V)
    return (u0, D0)

def energy(u, D, H, g):
    kinetic = 0.5*H*dot(u, u)*dx
    potential = 0.5*g*D*D*dx
    return (kinetic, potential)

def output_energy(energyfile, t, kinetic, potential):
    energyfile.write("times.append(%g);\n" % t)
    energyfile.write("kinetics.append(%.18g);\n" % kinetic)
    energyfile.write("potentials.append(%.18g);\n" % potential)

def main():

    # Setup parameters
    params = initialize_parameters()

    # Store parameters
    storage = "results/%s" % params["directory"]
    paramsfile = File("%s/application_parameters.xml" % storage)
    paramsfile << params

    # Set-up storage for energies
    energyfilename = "results/%s/energies.py" % params["directory"]
    info_blue("Storing energies to %s" % energyfilename)
    energyfile = open(energyfilename, "w")
    energyfile.write("# Results from linear shallow water test case.\n")
    energyfile.write("kinetics = [];\n")
    energyfile.write("potentials = [];\n")
    energyfile.write("times = [];\n")

    # Initialize mesh
    mesh = Mesh("../mixed-poisson/hdiv-l2/meshes/%s.xml.gz" % params["meshid"])
    global_normal = Expression(("x[0]", "x[1]", "x[2]"))
    mesh.init_cell_orientations(global_normal)

    # Set-up function spaces
    (S, V, PlottingMesh) = create_function_spaces(mesh, params["family"], params["order"])
    W = S * V

    # Setting up perp operator
    outward_normals = Expression(outward_normals_code)
    outward_normals.mesh = mesh
    perp = lambda u: cross(outward_normals, u)

    # Extract initial conditions
    (u_, D_) = initial_conditions(S, V)

    # Extract some parameters for discretization
    dt = Constant(params["dt"])
    f = Expression("x[2]")
    H = 1.0
    g = 1.0

    # Implicit midoint scheme discretization in time
    (u, D) = TrialFunctions(W)
    (w, phi) = TestFunctions(W)
    u_mid = 0.5*(u + u_)
    D_mid = 0.5*(D + D_)
    F = (inner(u - u_, w) - dt*div(w)*g*D_mid + dt*f*inner(perp(u_mid), w)
         + inner(D - D_, phi) + dt*H*div(u_mid)*phi)*dx
    (a, L) = system(F)

    # Preassemble matrix (because we can)
    A = assemble(a)

    # Define energy functional
    (kinetic_func, potential_func) = energy(u_, D_, H, g)

    # Setup solution function for current time
    uD = Function(W)

    # Predefine b (for the sake of reuse of memory)
    b = Vector(W.dim())

    # Set-up linear solver (so that we can reuse LU)
    solver = LUSolver(A)
    solver.parameters["same_nonzero_pattern"] = True
    solver.parameters["reuse_factorization"] = True

    # Set-up some pvd storage
    if params["store_pvds"]:
        ufile = File("%s/pvds/u.pvd" % storage)
        Dfile = File("%s/pvds/D.pvd" % storage)

    # Set-up some time related variables
    k = 0
    t = 0.0
    T = params["T"]
    dt = float(dt)

    # Output initial energy
    E_k = assemble(kinetic_func)
    E_p = assemble(potential_func)
    print t, E_k, E_p, E_k + E_p, D_.vector().min(), D_.vector().max()
    output_energy(energyfile, t, E_k, E_p)

    # Time loop
    while(t < (T - 0.5*dt)):

        # Assemble right-hand side, reuse b
        assemble(L, tensor=b)

        # Solve system
        solver.solve(uD.vector(), b)

        # Update previous solution
        u, D = uD.split()
        u_.assign(u)
        D_.assign(D)

        # Update time and counter
        t += dt
        k += 1

        # Output current energy and max/min depth
        E_k = assemble(kinetic_func)
        E_p = assemble(potential_func)
        print t, E_k, E_p, E_k + E_p, D_.vector().min(), D_.vector().max()
        output_energy(energyfile, t, E_k, E_p)

        # Store solutions to xml and pvd
        solutionfile = File("%s/xmls/uD_%d.xml.gz" % (storage, k))
        solutionfile << uD

        if params["store_pvds"]:
            uplot = project(u,PlottingMesh)
            ufile << uplot
            Dfile << D

if __name__ == "__main__":

    tic()
    main()
    print toc()

To derive the variational form as it’s defined here, we first take the inner product between (1) and the test function \boldsymbol w and the inner product between (2) and \phi (and leveraging integration by parts so that \int _\Omega \boldsymbol w \cdot g \nabla D dx = -\int _\Omega g D \nabla \cdot\boldsymbol w dx), and adding both together, we get the variational form which needs to be solved

(3) <\boldsymbol w, \frac{\partial \textbf{u}}{\partial t}> + <\boldsymbol w, f \textbf{u}^\perp> + <\phi,\frac{\partial D}{\partial t}> + <\phi,H \text{div}\textbf{u}>= <\text{div} \boldsymbol w, gD>

where <a,b> is the inner product between a and b. We then break the time derivatives down into a discretized form using something called the Implicit midoint scheme discretization in time (I’ve not heard of or used a midpoint discretization scheme before) where the definition of this scheme, from observation, appears to be first approximating the temporal derivatives via backwards difference scheme

\frac{\partial \textbf{u}}{\partial t} = \frac{\textbf u^{(m+1)}-\textbf u^{(m)}}{\Delta t}

\frac{\partial D}{\partial t} = \frac{D^{(m+1)}-D^{(m)}}{\Delta t}

and then assigning

\textbf u:= \textbf u_{mid} = \frac{1}{2}(\textbf u^{(m+1)}+\textbf u^{(m)})

D:= D_{mid} = \frac{1}{2}(D^{(m+1)}+D^{(m)}).

Our “current time” variational form at time step m=0 is now

<\boldsymbol w, (\textbf{u}-\textbf{u}_0)> + \Delta t <\boldsymbol w, f(\hat{\textbf k}\times \textbf u_{mid})> + <\phi,(D-D_0)> + \Delta t <\phi,H\text{div} \textbf u_{mid}> = \Delta t <\text{div}\textbf{w}, g D_{mid}>

which reflects how we’ve initially defined F.

Hopefully this post helps anyone else reading the initial post!

8

I think I figured it out now. It wasn’t just the wrong normal, FacetNormal, being used, but creating a MixedFunctionSpace was also getting in the way of creating a cross product. Opting to create a mixed element to create a functionspace instead did not raise any errors.

from mpi4py import MPI
import meshio
import numpy as np
import ufl
import basix.ufl
import dolfinx as df
from petsc4py import PETSc

from typing import Optional

import numpy as np
import dolfinx.fem.petsc


class Projector:
    """
    Projector for a given function.
    Solves Ax=b, where

    .. highlight:: python
    .. code-block:: python

        u, v = ufl.TrialFunction(Space), ufl.TestFunction(space)
        dx = ufl.Measure("dx", metadata=metadata)
        A = inner(u, v) * dx
        b = inner(function, v) * dx(metadata=metadata)

    Args:
        function: UFL expression of function to project
        space: Space to project function into
        petsc_options: Options to pass to PETSc
        jit_options: Options to pass to just in time compiler
        form_compiler_options: Options to pass to the form compiler
        metadata: Data to pass to the integration measure
    """

    _A: PETSc.Mat  # The mass matrix
    _b: PETSc.Vec  # The rhs vector
    _lhs: dolfinx.fem.Form  # The compiled form for the mass matrix
    _ksp: PETSc.KSP  # The PETSc solver
    _x: dolfinx.fem.Function  # The solution vector
    _dx: ufl.Measure  # Integration measure

    def __init__(
        self,
        space: dolfinx.fem.FunctionSpace,
        petsc_options: Optional[dict] = None,
        jit_options: Optional[dict] = None,
        form_compiler_options: Optional[dict] = None,
        metadata: Optional[dict] = None,
    ):
        petsc_options = {} if petsc_options is None else petsc_options
        jit_options = {} if jit_options is None else jit_options
        form_compiler_options = (
            {} if form_compiler_options is None else form_compiler_options
        )

        # Assemble projection matrix once
        u = ufl.TrialFunction(space)
        v = ufl.TestFunction(space)
        self._dx = ufl.Measure("dx", domain=space.mesh, metadata=metadata)
        a = ufl.inner(u, v) * self._dx(metadata=metadata)
        self._lhs = dolfinx.fem.form(
            a, jit_options=jit_options, form_compiler_options=form_compiler_options
        )
        self._A = dolfinx.fem.petsc.assemble_matrix(self._lhs)
        self._A.assemble()

        # Create vectors to store right hand side and the solution
        self._x = dolfinx.fem.Function(space)
        self._b = dolfinx.fem.Function(space)

        # Create Krylov Subspace solver
        self._ksp = PETSc.KSP().create(space.mesh.comm)
        self._ksp.setOperators(self._A)

        # Set PETSc options
        prefix = f"projector_{id(self)}"
        opts = PETSc.Options()
        opts.prefixPush(prefix)
        for k, v in petsc_options.items():
            opts[k] = v
        opts.prefixPop()
        self._ksp.setFromOptions()
        for opt in opts.getAll().keys():
            del opts[opt]

        # Set matrix and vector PETSc options
        self._A.setOptionsPrefix(prefix)
        self._A.setFromOptions()
        self._b.x.petsc_vec.setOptionsPrefix(prefix)
        self._b.x.petsc_vec.setFromOptions()

    def reassemble_lhs(self):
        dolfinx.fem.petsc.assemble_matrix(self._A, self._lhs)
        self._A.assemble()

    def assemble_rhs(self, h: ufl.core.expr.Expr):
        """
        Assemble the right hand side of the problem
        """
        v = ufl.TestFunction(self._b.function_space)
        rhs = ufl.inner(h, v) * self._dx
        rhs_compiled = dolfinx.fem.form(rhs)
        self._b.x.array[:] = 0.0
        dolfinx.fem.petsc.assemble_vector(self._b.x.petsc_vec, rhs_compiled)
        self._b.x.petsc_vec.ghostUpdate(
            addv=PETSc.InsertMode.ADD_VALUES, mode=PETSc.ScatterMode.REVERSE
        )
        self._b.x.petsc_vec.ghostUpdate(
            addv=PETSc.InsertMode.INSERT_VALUES, mode=PETSc.ScatterMode.FORWARD
        )

    def project(self, h: ufl.core.expr.Expr) -> dolfinx.fem.Function:
        """
        Compute projection using a PETSc KSP solver

        Args:
            assemble_rhs: Re-assemble RHS and re-apply boundary conditions if true
        """
        self.assemble_rhs(h)
        self._ksp.solve(self._b.x.petsc_vec, self._x.x.petsc_vec)
        return self._x

    def __del__(self):
        self._A.destroy()
        self._ksp.destroy()


order = 1
exact_case = 2
tmp_data = meshio.dolfin.read("supplement/examples/mixed-poisson/hdiv-l2/meshes/sphere_ico4.xml")

mesh = df.mesh.create_mesh(
    comm=MPI.COMM_WORLD,
    cells=tmp_data.cells_dict["triangle"],
    x=tmp_data.points,
    e=ufl.Mesh(basix.ufl.element("Lagrange", "triangle", 1, shape=(3,))),
)


def create_function_spaces(mesh, family, order):
    S_el = basix.ufl.element(family, mesh.basix_cell(), order)
    S = df.fem.functionspace(mesh, S_el)
    V_el = basix.ufl.element("DG", mesh.basix_cell(), order - 1)
    V = df.fem.functionspace(mesh, V_el)
    W_el = basix.ufl.mixed_element([S_el, V_el])
    W = df.fem.functionspace(mesh, W_el)
    PlottingMesh = df.fem.functionspace(mesh, ("CG", order))
    return (S, V, W, PlottingMesh)


def initial_conditions(S, V):
    u0 = df.fem.Function(S)
    u0.x.array[:] = 0.0
    f = ufl.exp(-(((-x[1] - 1) / 0.1) ** 2))
    petsc_options = {"ksp_type": "preonly", "pc_type": "lu"}
    V_projector = Projector(V, petsc_options=petsc_options)
    D0 = V_projector.project(f)
    return (u0, D0)


# Define global normal
x = ufl.SpatialCoordinate(mesh)
global_normal = x  # Not sure exactly what this is used for
family = "RT"
order = 1

# Define mixed mesh
S, V, W, PlottingMesh = create_function_spaces(mesh, family, order)

# Constants
dt = df.default_scalar_type(0.05)
H = 1.0
g = 1.0

# forcing term
f = x[2]

# Initial values of solutions u and D
u_, D_ = initial_conditions(S, V)
(u, D) = ufl.TrialFunctions(W)
(w, phi) = ufl.TestFunctions(W)

D_mid = 0.5 * (D + D_)
u_mid = 0.5 * (u_ + u)
# normal unit vector
k = ufl.CellNormal(mesh)


# u_n = k x u
u_n = ufl.cross(k, u)

F = (
    ufl.inner(u - u_, w)
    - dt * ufl.div(w) * g * D_mid
    + dt * f * ufl.inner(u_n, w)
    + ufl.inner(D - D_, phi)
    + dt * H * ufl.div(u_mid) * phi
) * ufl.dx

a, L = ufl.system(F)

(where the changes were made to def create_function_spaces()) still compiles, and now

bilinear_form = df.fem.form(a)
linear_form = df.fem.form(L)

successfully compiles, as well. Unless you see any potential issues for doing it this way, I think I can now close this problem and move on to tackling the next one.

9

This seems rather strange to me, could you elaborate with a minimal example?

10

Actually, testing my variational form when the functionspace is a MixedFunctionSpace type, this assertion I’ve made about the cross product causing the error is incorrect. In case I haven’t said this, the error comes at the step of reading the variational form into some PETSc’s KSP solver or dolfinx.fem.petsc.LinearProblem when working with a mixed problem, however it appears as thought the cross product being involved (or ufl.perp(u) in the 2D case) has nothing to do with it. Here is a shortened version (no projection included, and instead opting to interpolate despite not yielding correct results) on a square where the variational form is that of the Shallow Water Equations with all constants set to 1

\frac{\partial \textbf{u}}{\partial t} + f (\hat{\textbf{k}}\times\textbf{u}) + \nabla D = 0

\frac{\partial D}{\partial t} + \nabla \cdot \textbf{u} = 0

<\boldsymbol{w},\textbf{u}_t> + <\boldsymbol{w} ,f (\hat{\textbf{k}}\times\textbf{u})>-<\nabla \cdot \boldsymbol w, D> + <\phi,D_t> + <\nabla \cdot \textbf{u}>=0

from mpi4py import MPI
import numpy as np
from dolfinx import fem, mesh
from Projector import Projector
import ufl

def initial_conditions(S, V):
    u0 = fem.Function(S)
    u0.interpolate(lambda x: np.zeros(x[:S.mesh.topology.dim].shape))
    D0 = fem.Function(V)
    D0.interpolate(lambda x: np.exp(-((-x[1]-1)/0.1)**2))
    return (u0, D0)

domain = mesh.create_unit_square(MPI.COMM_WORLD, 32, 32, mesh.CellType.quadrilateral)
x = ufl.SpatialCoordinate(domain)

S = fem.functionspace(domain, ("RT",1))
V = fem.functionspace(domain, ("DG",1))
W = ufl.MixedFunctionSpace(*(S, V))

# Initial values of solutions u and D
u_, D_ = initial_conditions(S, V)
(u, D) = ufl.TrialFunctions(W)
(w, phi) = ufl.TestFunctions(W)

# u_n = k x u
u_n = ufl.perp(u)

dt = 0.05 # delta t
f = x[1]# coriolis parameter

F = (
    ufl.inner(u - u_, w)
    - ufl.div(w) * D
    + dt * f * ufl.inner(u_n, w)
    + ufl.inner(D - D_, phi)
    + dt * ufl.div(u) * phi
) * ufl.dx

(a, L) = ufl.system(F)

petsc_options = {
    "ksp_type": "preonly",
    "pc_type": "lu",
    "pc_factor_mat_solver_type": "mumps",
    "ksp_error_if_not_converged": True,
}
problem = fem.petsc.LinearProblem(
    a,
    L,
    bcs=[],
    petsc_options=petsc_options,
    petsc_options_prefix="mixed_poisson_",
    kind="mpi",
)

returns the error:

---------------------------------------------------------------------------
ArityMismatch                             Traceback (most recent call last)
Cell In[13], line 7
     1 petsc_options = {
     2     "ksp_type": "preonly",
     3     "pc_type": "lu",
     4     "pc_factor_mat_solver_type": "mumps",
     5     "ksp_error_if_not_converged": True,
     6 }
----> 7 problem = fem.petsc.LinearProblem(
     8     ufl.extract_blocks(a),
     9     L,
    10     bcs=[],
    11     petsc_options=petsc_options,
    12     petsc_options_prefix="mixed_poisson_",
    13     kind="mpi",
    14 )

File /usr/local/dolfinx-real/lib/python3.12/dist-packages/dolfinx/fem/petsc.py:829, in LinearProblem.__init__(self, a, L, petsc_options_prefix, bcs, u, P, kind, petsc_options, form_compiler_options, jit_options, entity_maps)
   743 """Initialize solver for a linear variational problem.
   744 
   745 By default, the underlying KSP solver uses PETSc's default
  (...)    820         entity ``i`` in the integration domain mesh.
   821 """
   822 self._a = _create_form(
   823     a,
   824     dtype=PETSc.ScalarType,  # type: ignore[attr-defined]
  (...)    827     entity_maps=entity_maps,
   828 )
--> 829 self._L = _create_form(
   830     L,
   831     dtype=PETSc.ScalarType,  # type: ignore[attr-defined]
   832     form_compiler_options=form_compiler_options,
   833     jit_options=jit_options,
   834     entity_maps=entity_maps,
   835 )
   836 self._A = create_matrix(self._a, kind=kind)
   837 self._preconditioner = _create_form(
   838     P,  # type: ignore[arg-type]
   839     dtype=PETSc.ScalarType,  # type: ignore[attr-defined]
  (...)    842     entity_maps=entity_maps,
   843 )

File /usr/local/dolfinx-real/lib/python3.12/dist-packages/dolfinx/fem/forms.py:438, in form(form, dtype, form_compiler_options, jit_options, entity_maps)
   435     else:
   436         return form
--> 438 return _create_form(form)

File /usr/local/dolfinx-real/lib/python3.12/dist-packages/dolfinx/fem/forms.py:430, in form.<locals>._create_form(form)
   420 """Recursively convert ufl.Forms to dolfinx.fem.Form.
   421 
   422 Args:
  (...)    427     A ``dolfinx.fem.Form`` or a list of ``dolfinx.fem.Form``.
   428 """
   429 if isinstance(form, ufl.Form):
--> 430     return _form(form)
   431 elif isinstance(form, ufl.ZeroBaseForm):
   432     return _zero_form(form)

File /usr/local/dolfinx-real/lib/python3.12/dist-packages/dolfinx/fem/forms.py:339, in form.<locals>._form(form)
   337 if msh is None:
   338     raise RuntimeError("Expecting to find a Mesh in the form.")
--> 339 ufcx_form, module, code = jit.ffcx_jit(
   340     msh.comm, form, form_compiler_options=form_compiler_options, jit_options=jit_options
   341 )
   343 # For each argument in form extract its function space
   344 V = [arg.ufl_function_space()._cpp_object for arg in form.arguments()]

File /usr/local/dolfinx-real/lib/python3.12/dist-packages/dolfinx/jit.py:61, in mpi_jit_decorator.<locals>.mpi_jit(comm, *args, **kwargs)
    57 @functools.wraps(local_jit)
    58 def mpi_jit(comm, *args, **kwargs):
    59     # Just call JIT compiler when running in serial
    60     if comm.size == 1:
---> 61         return local_jit(*args, **kwargs)
    63     # Default status (0 == ok, 1 == fail)
    64     status = 0

File /usr/local/dolfinx-real/lib/python3.12/dist-packages/dolfinx/jit.py:216, in ffcx_jit(ufl_object, form_compiler_options, jit_options)
   214 # Switch on type and compile, returning cffi object
   215 if isinstance(ufl_object, ufl.Form):
--> 216     r = ffcx.codegeneration.jit.compile_forms([ufl_object], options=p_ffcx, **p_jit)
   217 elif isinstance(ufl_object, tuple) and isinstance(ufl_object[0], ufl.core.expr.Expr):
   218     r = ffcx.codegeneration.jit.compile_expressions([ufl_object], options=p_ffcx, **p_jit)

File /dolfinx-env/lib/python3.12/site-packages/ffcx/codegeneration/jit.py:206, in compile_forms(forms, options, cache_dir, timeout, cffi_extra_compile_args, cffi_verbose, cffi_debug, cffi_libraries, visualise)
   203     for name in form_names:
   204         decl += form_template.format(name=name)
--> 206     impl = _compile_objects(
   207         decl,
   208         forms,
   209         form_names,
   210         module_name,
   211         p,
   212         cache_dir,
   213         cffi_extra_compile_args,
   214         cffi_verbose,
   215         cffi_debug,
   216         cffi_libraries,
   217         visualise=visualise,
   218     )
   219 except Exception as e:
   220     try:
   221         # remove c file so that it will not timeout next time

File /dolfinx-env/lib/python3.12/site-packages/ffcx/codegeneration/jit.py:331, in _compile_objects(decl, ufl_objects, object_names, module_name, options, cache_dir, cffi_extra_compile_args, cffi_verbose, cffi_debug, cffi_libraries, visualise)
   327 libraries = _libraries + cffi_libraries if cffi_libraries is not None else _libraries
   329 # JIT uses module_name as prefix, which is needed to make names of all struct/function
   330 # unique across modules
--> 331 _, code_body = ffcx.compiler.compile_ufl_objects(
   332     ufl_objects, prefix=module_name, options=options, visualise=visualise
   333 )
   335 # Raise error immediately prior to compilation if no support for C99
   336 # _Complex. Doing this here allows FFCx to be used for complex codegen on
   337 # Windows.
   338 if sys.platform.startswith("win32"):

File /dolfinx-env/lib/python3.12/site-packages/ffcx/compiler.py:108, in compile_ufl_objects(ufl_objects, options, object_names, prefix, visualise)
   106 # Stage 1: analysis
   107 cpu_time = time()
--> 108 analysis = analyze_ufl_objects(ufl_objects, options["scalar_type"])  # type: ignore
   109 _print_timing(1, time() - cpu_time)
   111 # Stage 2: intermediate representation

File /dolfinx-env/lib/python3.12/site-packages/ffcx/analysis.py:96, in analyze_ufl_objects(ufl_objects, scalar_type)
    93     else:
    94         raise TypeError("UFL objects not recognised.")
---> 96 form_data = tuple(_analyze_form(form, scalar_type) for form in forms)
    97 for data in form_data:
    98     elements += data.unique_sub_elements

File /dolfinx-env/lib/python3.12/site-packages/ffcx/analysis.py:96, in <genexpr>(.0)
    93     else:
    94         raise TypeError("UFL objects not recognised.")
---> 96 form_data = tuple(_analyze_form(form, scalar_type) for form in forms)
    97 for data in form_data:
    98     elements += data.unique_sub_elements

File /dolfinx-env/lib/python3.12/site-packages/ffcx/analysis.py:182, in _analyze_form(form, scalar_type)
   179 complex_mode = np.issubdtype(scalar_type, np.complexfloating)
   181 # Compute form metadata
--> 182 form_data: ufl.algorithms.formdata.FormData = ufl.algorithms.compute_form_data(
   183     form,
   184     do_apply_function_pullbacks=True,
   185     do_apply_integral_scaling=True,
   186     do_apply_geometry_lowering=True,
   187     preserve_geometry_types=(ufl.classes.Jacobian,),
   188     do_apply_restrictions=True,
   189     do_append_everywhere_integrals=False,  # do not add dx integrals to dx(i) in UFL
   190     complex_mode=complex_mode,
   191 )
   193 # Determine unique quadrature degree and quadrature scheme
   194 # per each integral data
   195 for id, integral_data in enumerate(form_data.integral_data):
   196     # Iterate through groups of integral data. There is one integral
   197     # data for all integrals with same domain, itype, subdomain_id
  (...)    201     # all integrals in this integral data group, i.e. must be the
   202     # same for for the same (domain, itype, subdomain_id)

File /dolfinx-env/lib/python3.12/site-packages/ufl/algorithms/compute_form_data.py:482, in compute_form_data(form, do_apply_function_pullbacks, do_apply_integral_scaling, do_apply_geometry_lowering, preserve_geometry_types, do_apply_default_restrictions, do_apply_restrictions, do_estimate_degrees, do_append_everywhere_integrals, do_replace_functions, coefficients_to_split, complex_mode, do_remove_component_tensors)
   479 preprocessed_form = reconstruct_form_from_integral_data(self.integral_data)
   481 # TODO: Test how fast this is
--> 482 check_form_arity(preprocessed_form, self.original_form.arguments(), complex_mode)
   484 # TODO: This member is used by unit tests, change the tests to
   485 # remove this!
   486 self.preprocessed_form = preprocessed_form

File /dolfinx-env/lib/python3.12/site-packages/ufl/algorithms/check_arities.py:213, in check_form_arity(form, arguments, complex_mode)
   211 """Check the arity of a form."""
   212 for itg in form.integrals():
--> 213     check_integrand_arity(itg.integrand(), arguments, complex_mode)

File /dolfinx-env/lib/python3.12/site-packages/ufl/algorithms/check_arities.py:194, in check_integrand_arity(expr, arguments, complex_mode)
   192 arguments = tuple(sorted(set(arguments), key=lambda x: (x.number(), x.part())))
   193 rules = ArityChecker(arguments)
--> 194 arg_tuples = map_expr_dag(rules, expr, compress=False)
   195 args = tuple(a[0] for a in arg_tuples)
   196 if args != arguments:

File /dolfinx-env/lib/python3.12/site-packages/ufl/corealg/map_dag.py:35, in map_expr_dag(function, expression, compress, vcache, rcache)
    15 def map_expr_dag(function, expression, compress=True, vcache=None, rcache=None):
    16     """Apply a function to each subexpression node in an expression DAG.
    17 
    18     If the same function is called multiple times in a transformation
  (...)     33         The result of the final function call
    34     """
---> 35     (result,) = map_expr_dags(
    36         function, [expression], compress=compress, vcache=vcache, rcache=rcache
    37     )
    38     return result

File /dolfinx-env/lib/python3.12/site-packages/ufl/corealg/map_dag.py:103, in map_expr_dags(function, expressions, compress, vcache, rcache)
   101     r = handlers[v._ufl_typecode_](v)
   102 else:
--> 103     r = handlers[v._ufl_typecode_](v, *(vcache[u] for u in v.ufl_operands))
   105 # Optionally check if r is in rcache, a memory optimization
   106 # to be able to keep representation of result compact
   107 if compress:

File /dolfinx-env/lib/python3.12/site-packages/ufl/algorithms/check_arities.py:58, in ArityChecker.sum(self, o, a, b)
    56 """Apply to sum."""
    57 if a != b:
---> 58     raise ArityMismatch(
    59         f"Adding expressions with non-matching form arguments "
    60         f"{tuple(map(_afmt, a))} vs {tuple(map(_afmt, b))}."
    61     )
    62 return a

ArityMismatch: Adding expressions with non-matching form arguments ('v_0^1',) vs ('v_0^0',).

If we use the special case where f=0, and therefore no cross product is present, by replacing the last lines with

F = (
    ufl.inner(u - u_, w)
    - ufl.div(w) * D
    + ufl.inner(D - D_, phi)
    + dt * ufl.div(u) * phi
) * ufl.dx

(a, L) = ufl.system(F)

petsc_options = {
    "ksp_type": "preonly",
    "pc_type": "lu",
    "pc_factor_mat_solver_type": "mumps",
    "ksp_error_if_not_converged": True,
}
problem = fem.petsc.LinearProblem(
    a,
    L,
    bcs=[],
    petsc_options=petsc_options,
    petsc_options_prefix="mixed_poisson_",
    kind="mpi",
)

I get similar error as before

---------------------------------------------------------------------------
ArityMismatch                             Traceback (most recent call last)
Cell In[1], line 47
     39 (a, L) = ufl.system(F)
     41 petsc_options = {
     42     "ksp_type": "preonly",
     43     "pc_type": "lu",
     44     "pc_factor_mat_solver_type": "mumps",
     45     "ksp_error_if_not_converged": True,
     46 }
---> 47 problem = fem.petsc.LinearProblem(
     48     a,
     49     L,
     50     bcs=[],
     51     petsc_options=petsc_options,
     52     petsc_options_prefix="mixed_poisson_",
     53     kind="mpi",
     54 )

File /usr/local/dolfinx-real/lib/python3.12/dist-packages/dolfinx/fem/petsc.py:822, in LinearProblem.__init__(self, a, L, petsc_options_prefix, bcs, u, P, kind, petsc_options, form_compiler_options, jit_options, entity_maps)
    728 def __init__(
    729     self,
    730     a: typing.Union[ufl.Form, Sequence[Sequence[ufl.Form]]],
   (...)    741     entity_maps: typing.Optional[Sequence[_EntityMap]] = None,
    742 ) -> None:
    743     """Initialize solver for a linear variational problem.
    744 
    745     By default, the underlying KSP solver uses PETSc's default
   (...)    820             entity ``i`` in the integration domain mesh.
    821     """
--> 822     self._a = _create_form(
    823         a,
    824         dtype=PETSc.ScalarType,  # type: ignore[attr-defined]
    825         form_compiler_options=form_compiler_options,
    826         jit_options=jit_options,
    827         entity_maps=entity_maps,
    828     )
    829     self._L = _create_form(
    830         L,
    831         dtype=PETSc.ScalarType,  # type: ignore[attr-defined]
   (...)    834         entity_maps=entity_maps,
    835     )
    836     self._A = create_matrix(self._a, kind=kind)

File /usr/local/dolfinx-real/lib/python3.12/dist-packages/dolfinx/fem/forms.py:438, in form(form, dtype, form_compiler_options, jit_options, entity_maps)
    435     else:
    436         return form
--> 438 return _create_form(form)

File /usr/local/dolfinx-real/lib/python3.12/dist-packages/dolfinx/fem/forms.py:430, in form.<locals>._create_form(form)
    420 """Recursively convert ufl.Forms to dolfinx.fem.Form.
    421 
    422 Args:
   (...)    427     A ``dolfinx.fem.Form`` or a list of ``dolfinx.fem.Form``.
    428 """
    429 if isinstance(form, ufl.Form):
--> 430     return _form(form)
    431 elif isinstance(form, ufl.ZeroBaseForm):
    432     return _zero_form(form)

File /usr/local/dolfinx-real/lib/python3.12/dist-packages/dolfinx/fem/forms.py:339, in form.<locals>._form(form)
    337 if msh is None:
    338     raise RuntimeError("Expecting to find a Mesh in the form.")
--> 339 ufcx_form, module, code = jit.ffcx_jit(
    340     msh.comm, form, form_compiler_options=form_compiler_options, jit_options=jit_options
    341 )
    343 # For each argument in form extract its function space
    344 V = [arg.ufl_function_space()._cpp_object for arg in form.arguments()]

File /usr/local/dolfinx-real/lib/python3.12/dist-packages/dolfinx/jit.py:61, in mpi_jit_decorator.<locals>.mpi_jit(comm, *args, **kwargs)
     57 @functools.wraps(local_jit)
     58 def mpi_jit(comm, *args, **kwargs):
     59     # Just call JIT compiler when running in serial
     60     if comm.size == 1:
---> 61         return local_jit(*args, **kwargs)
     63     # Default status (0 == ok, 1 == fail)
     64     status = 0

File /usr/local/dolfinx-real/lib/python3.12/dist-packages/dolfinx/jit.py:216, in ffcx_jit(ufl_object, form_compiler_options, jit_options)
    214 # Switch on type and compile, returning cffi object
    215 if isinstance(ufl_object, ufl.Form):
--> 216     r = ffcx.codegeneration.jit.compile_forms([ufl_object], options=p_ffcx, **p_jit)
    217 elif isinstance(ufl_object, tuple) and isinstance(ufl_object[0], ufl.core.expr.Expr):
    218     r = ffcx.codegeneration.jit.compile_expressions([ufl_object], options=p_ffcx, **p_jit)

File /dolfinx-env/lib/python3.12/site-packages/ffcx/codegeneration/jit.py:206, in compile_forms(forms, options, cache_dir, timeout, cffi_extra_compile_args, cffi_verbose, cffi_debug, cffi_libraries, visualise)
    203     for name in form_names:
    204         decl += form_template.format(name=name)
--> 206     impl = _compile_objects(
    207         decl,
    208         forms,
    209         form_names,
    210         module_name,
    211         p,
    212         cache_dir,
    213         cffi_extra_compile_args,
    214         cffi_verbose,
    215         cffi_debug,
    216         cffi_libraries,
    217         visualise=visualise,
    218     )
    219 except Exception as e:
    220     try:
    221         # remove c file so that it will not timeout next time

File /dolfinx-env/lib/python3.12/site-packages/ffcx/codegeneration/jit.py:331, in _compile_objects(decl, ufl_objects, object_names, module_name, options, cache_dir, cffi_extra_compile_args, cffi_verbose, cffi_debug, cffi_libraries, visualise)
    327 libraries = _libraries + cffi_libraries if cffi_libraries is not None else _libraries
    329 # JIT uses module_name as prefix, which is needed to make names of all struct/function
    330 # unique across modules
--> 331 _, code_body = ffcx.compiler.compile_ufl_objects(
    332     ufl_objects, prefix=module_name, options=options, visualise=visualise
    333 )
    335 # Raise error immediately prior to compilation if no support for C99
    336 # _Complex. Doing this here allows FFCx to be used for complex codegen on
    337 # Windows.
    338 if sys.platform.startswith("win32"):

File /dolfinx-env/lib/python3.12/site-packages/ffcx/compiler.py:108, in compile_ufl_objects(ufl_objects, options, object_names, prefix, visualise)
    106 # Stage 1: analysis
    107 cpu_time = time()
--> 108 analysis = analyze_ufl_objects(ufl_objects, options["scalar_type"])  # type: ignore
    109 _print_timing(1, time() - cpu_time)
    111 # Stage 2: intermediate representation

File /dolfinx-env/lib/python3.12/site-packages/ffcx/analysis.py:96, in analyze_ufl_objects(ufl_objects, scalar_type)
     93     else:
     94         raise TypeError("UFL objects not recognised.")
---> 96 form_data = tuple(_analyze_form(form, scalar_type) for form in forms)
     97 for data in form_data:
     98     elements += data.unique_sub_elements

File /dolfinx-env/lib/python3.12/site-packages/ffcx/analysis.py:96, in <genexpr>(.0)
     93     else:
     94         raise TypeError("UFL objects not recognised.")
---> 96 form_data = tuple(_analyze_form(form, scalar_type) for form in forms)
     97 for data in form_data:
     98     elements += data.unique_sub_elements

File /dolfinx-env/lib/python3.12/site-packages/ffcx/analysis.py:182, in _analyze_form(form, scalar_type)
    179 complex_mode = np.issubdtype(scalar_type, np.complexfloating)
    181 # Compute form metadata
--> 182 form_data: ufl.algorithms.formdata.FormData = ufl.algorithms.compute_form_data(
    183     form,
    184     do_apply_function_pullbacks=True,
    185     do_apply_integral_scaling=True,
    186     do_apply_geometry_lowering=True,
    187     preserve_geometry_types=(ufl.classes.Jacobian,),
    188     do_apply_restrictions=True,
    189     do_append_everywhere_integrals=False,  # do not add dx integrals to dx(i) in UFL
    190     complex_mode=complex_mode,
    191 )
    193 # Determine unique quadrature degree and quadrature scheme
    194 # per each integral data
    195 for id, integral_data in enumerate(form_data.integral_data):
    196     # Iterate through groups of integral data. There is one integral
    197     # data for all integrals with same domain, itype, subdomain_id
   (...)    201     # all integrals in this integral data group, i.e. must be the
    202     # same for for the same (domain, itype, subdomain_id)

File /dolfinx-env/lib/python3.12/site-packages/ufl/algorithms/compute_form_data.py:482, in compute_form_data(form, do_apply_function_pullbacks, do_apply_integral_scaling, do_apply_geometry_lowering, preserve_geometry_types, do_apply_default_restrictions, do_apply_restrictions, do_estimate_degrees, do_append_everywhere_integrals, do_replace_functions, coefficients_to_split, complex_mode, do_remove_component_tensors)
    479 preprocessed_form = reconstruct_form_from_integral_data(self.integral_data)
    481 # TODO: Test how fast this is
--> 482 check_form_arity(preprocessed_form, self.original_form.arguments(), complex_mode)
    484 # TODO: This member is used by unit tests, change the tests to
    485 # remove this!
    486 self.preprocessed_form = preprocessed_form

File /dolfinx-env/lib/python3.12/site-packages/ufl/algorithms/check_arities.py:213, in check_form_arity(form, arguments, complex_mode)
    211 """Check the arity of a form."""
    212 for itg in form.integrals():
--> 213     check_integrand_arity(itg.integrand(), arguments, complex_mode)

File /dolfinx-env/lib/python3.12/site-packages/ufl/algorithms/check_arities.py:194, in check_integrand_arity(expr, arguments, complex_mode)
    192 arguments = tuple(sorted(set(arguments), key=lambda x: (x.number(), x.part())))
    193 rules = ArityChecker(arguments)
--> 194 arg_tuples = map_expr_dag(rules, expr, compress=False)
    195 args = tuple(a[0] for a in arg_tuples)
    196 if args != arguments:

File /dolfinx-env/lib/python3.12/site-packages/ufl/corealg/map_dag.py:35, in map_expr_dag(function, expression, compress, vcache, rcache)
     15 def map_expr_dag(function, expression, compress=True, vcache=None, rcache=None):
     16     """Apply a function to each subexpression node in an expression DAG.
     17 
     18     If the same function is called multiple times in a transformation
   (...)     33         The result of the final function call
     34     """
---> 35     (result,) = map_expr_dags(
     36         function, [expression], compress=compress, vcache=vcache, rcache=rcache
     37     )
     38     return result

File /dolfinx-env/lib/python3.12/site-packages/ufl/corealg/map_dag.py:103, in map_expr_dags(function, expressions, compress, vcache, rcache)
    101     r = handlers[v._ufl_typecode_](v)
    102 else:
--> 103     r = handlers[v._ufl_typecode_](v, *(vcache[u] for u in v.ufl_operands))
    105 # Optionally check if r is in rcache, a memory optimization
    106 # to be able to keep representation of result compact
    107 if compress:

File /dolfinx-env/lib/python3.12/site-packages/ufl/algorithms/check_arities.py:58, in ArityChecker.sum(self, o, a, b)
     56 """Apply to sum."""
     57 if a != b:
---> 58     raise ArityMismatch(
     59         f"Adding expressions with non-matching form arguments "
     60         f"{tuple(map(_afmt, a))} vs {tuple(map(_afmt, b))}."
     61     )
     62 return a

ArityMismatch: Adding expressions with non-matching form arguments ('v_0^0', 'v_1^1') vs ('v_0^1', 'v_1^1').

. The errors go away when replacing the lines creating the functionspaces with

S_el = basix.ufl.element("RT", domain.basix_cell(), 1)
S = fem.functionspace(domain, S_el)
V_el = basix.ufl.element("DG", domain.basix_cell(), 0)
V = fem.functionspace(domain, V_el)
W_el = basix.ufl.mixed_element([S_el, V_el])
W = fem.functionspace(domain, W_el)

This might help clear some things up:

I originally tried to follow your example here as a template for handling this time varying problem as they both use a similar mixed function space, however examining it again, I just am now noticing some differences in your variational form from the one I tried to set up (which might hold the answer for this “non-matching form” error).

a = (
    ufl.inner(sigma_, tau_) + ufl.div(sigma_) * v + ufl.div(tau_) * u + r * v + t * u
) * ufl.dx
L = [ufl.ZeroBaseForm((tau,)), g * v * ufl.dx, ufl.ZeroBaseForm((t,))]

Your L is a list of 3 items rather than simply a convolution of the RHS, g, and the test function, v. Your list instead, contains that integral as well as the other 2 test functions by themselves, being the tau vector and t scalar; I don’t know what ZeroBaseForm means or what it does. Is my Linear term on the RHS supposed to be like it is here; does it have anything to do with the error about mismatching forms?

11

That is because you need to call ufl.extract_blocks on a and L before passing them into LinearProblem. See for instance dolfinx/python/demo/demo_navier-stokes.py at 213bf7cbc09dff556f5d69f3ad50adc9314e2f6d · FEniCS/dolfinx · GitHub

If the rhs (L in your case) do not have test functions from all sub spaces, one has to provide the sequence of forms manually as input to Linear problem. This is so the correct rhs vector can be created

1 Like
12

Oh, I didn’t know I needed to do that for both. It compiles now. Pressing shift+tab to pull up the docs for extract_blocks,

Docstring:
Extract blocks.

Given a linear or bilinear form on a mixed space,
extract the block corresponding to the indices ix, iy.

Example:
   a = inner(grad(u), grad(v))*dx + div(u)*q*dx + div(v)*p*dx
   extract_blocks(a, 0, 0) -> inner(grad(u), grad(v))*dx
   extract_blocks(a) -> [inner(grad(u), grad(v))*dx, div(v)*p*dx, div(u)*q*dx, 0]

The line
extract_blocks(a) -> [inner(grad(u), grad(v))*dx, div(v)*p*dx, div(u)*q*dx, 0]
tells me that each of these “blocks” extracted from the linear or bilinear form over the mixed space returns a “list” of the integrals inside (and there’s a 0 for whatever reason?). If the linear or bilinear form does not contain a convolution integral for each and every test function, then we need to build the “block” by creating and filling this list of integrals ourself, and ZeroBaseForm accomplishes this. I’m guessing I could also fill in the list with convolution integrals with 0 and the test function as this would amount to adding 0 to the RHS and also include the missing test functions?

I saw the NS problem in the tutorials, but I never actually worked on it. Looking at it now, I notice several similarities; time dependent, 2 variables of interest (one a scalar and the other a vector), dual form/system of 2 equations, mixed functionspace of 2 functionspaces. The variational form looks quite different with things like “jump” and “outer”. There is hopefully some additional value to be gained from this as well if not to only illustrate just this.

Anyways, I’m at the point where the shallow water code is all compiling now, but the values are off. I think this covers the different points of confusion I was having in this post.

13

The zero is there to illustrate that there is no (p, q) block in that particular bi-linear form.

You could add a dolfinx.fem.Constant(mesh,0.)*test_function*ufl.dx for each missing test function in the linear form (or a similar one for missing entry in the bilinear form, here with both test and trial functions). However, using ufl.ZeroBaseForm(tuple_of_test_and_trial_functions) is a more efficient way of instructing dolfinx about this pattern.

Please make a new topic with a more tailored question if you need further guidance. This has strayed far away from the original topic of the post.

1 Like