How to perform modulo operation with ufl

General
1

Hi all!

I’m trying to perform a modulo operation with ufl:


flux = ufl.conditional(
     t % 10 < 5, 1, 0
)

Where t is a ufl.Constant.

Obviously this doesn’t work, but I couldn’t find the correct ufl function. Any pointers?

2

As far as I can tell UFL does not support this operator (ufl/ChangeLog.md at 2125dcec5df6dd005cf54175687d43fc64564bed · FEniCS/ufl · GitHub,
ufl/ufl/core/expr.py at 2125dcec5df6dd005cf54175687d43fc64564bed · FEniCS/ufl · GitHub)
I think what you would have to do instead is to

t = dolfinx.fem.Constant(mesh, some_value
flux = dolfinx.fem.Constant(mesh, some_other_value)
# do something
flux.value = 1 if float(t) % 10 < 5 else 0
# Doe something to update
t.value += float(dt)
flux.value = 1 if float(t) % 10 < 5 else 0
2 Likes
3

Without having tested it myself, could an option be to use numpy.mod? For example:

flux = ufl.conditional(
     ufl.lt(numpy.mod(t.value, 10), 5), 1, 0
)

Let me know if this worked/was of help!

4

Thanks @dokken I see how this will achieve the correct behaviour. Is it worth opening an issue somewhere?

5

This one is a bit scary, as you will loose the control over t, i.e.

from mpi4py import MPI
import dolfinx
import ufl
import numpy as np
mesh =dolfinx.mesh.create_unit_interval(MPI.COMM_WORLD, 1)

def condition(t):
    return dolfinx.default_scalar_type(1 if float(t) % 10 < 5 else 0)

t = dolfinx.fem.Constant(mesh, 16.0)
flux = dolfinx.fem.Constant(mesh, condition(t))
flux_form = dolfinx.fem.form(flux*ufl.dx(domain=mesh))
print(dolfinx.fem.assemble_scalar(flux_form))

# Alternative flux
flux2 = ufl.conditional(
     ufl.lt(np.mod(t.value, 10), 5), 1, 0
)
flux2_form = dolfinx.fem.form(flux2*ufl.dx(domain=mesh))

print(dolfinx.fem.assemble_scalar(flux2_form)) 

t.value = 12.0
# First one updates
flux.value = condition(t)
print(dolfinx.fem.assemble_scalar(flux_form))

# Second one doesnt
print(dolfinx.fem.assemble_scalar(flux2_form)) 


0.0
0.0
1.0
0.0

You would then have to recompile the form each time.

1 Like
6

@dokken could you explain what is dolfinx.default_scalar_type here?

What would happen if we don’t use it and just return a float?

7

Ah, that is a subtle yet very important point, thanks for showcasing that!

8

dolfinx.default_scalar_type is numpy float64 if PETSc is not installed.
If PETSc is installed, it is whatever is the default type in PETSc (it could be np.float32, np.float64, np.complex64, np.complex128.

DOLFINx can use any floating time at any point, while PETSc is limited to a specific one (one needs one petsc installation per floating type).

So long story short:
If you use PETSc, then this is a wrapper around PETSc.ScalarType.

Python floating types have arbitrary precision, and must be cast to an appropriate type for working in C++.

1 Like