Evaluate ufl expression, element CFL and Peclet

dolfinx
1

Hi everyone,

it took me ages to figure out the many steps needed to compute a local element-wise ufl expression, in this case CFL and Peclet. can you confirm this is the correct (and easiest way) of doing it?

    # Compute Peclet and CFL number locally
    vnorm = ufl.sqrt(ufl.dot(vel, vel))  # local norm of the velocity
    h = ufl.CellDiameter(mesh)           # local mesh size
    Peclet = h*vnorm/ufl.sqrt(ufl.inner(D, D)/3) # Peclet number
    cfl = 2*dt*vnorm/h                # local time step
    
    DG = dolfinx.fem.FunctionSpace(mesh, ("DG", 0))
    cfl_expr = dolfinx.fem.Expression(cfl,DG.element.interpolation_points())
    cfl_fun = dolfinx.fem.Function(DG)
    cfl_fun.interpolate(cfl_expr)
    maxcfl = np.max(cfl_fun.x.array)
    print("maxcfl = ", maxcfl)
2

You could do it this way, or using

v = ufl.TrialFunction(DG)
cfl_vec = dolfinx.fem.petsc.assemble_vector(dolfinx.fem.form(cfl_fun*v/ufl.CellVolume(mesh)*ufl.dx))

as you are using a DG-0 function space for the cfl function.
Either is fine.
Consider for instance the comparison:

import dolfinx
import ufl
from mpi4py import MPI
import numpy as np

mesh = dolfinx.mesh.create_unit_square(MPI.COMM_WORLD, 10, 10)
V = dolfinx.fem.VectorFunctionSpace(mesh, ("Lagrange", 2))
Q = dolfinx.fem.FunctionSpace(mesh, ("DG", 0))
vel = dolfinx.fem.Function(V)
D = dolfinx.fem.Function(Q)
vel.interpolate(lambda x: (np.sin(x[0]), np.cos(x[1])))
D.interpolate(lambda x: 1 + x[0]+x[1])
dt = dolfinx.fem.Constant(mesh, 0.01)

# Compute Peclet and CFL number locally
vnorm = ufl.sqrt(ufl.dot(vel, vel))  # local norm of the velocity
h = ufl.CellDiameter(mesh)           # local mesh size
Peclet = h*vnorm/ufl.sqrt(ufl.inner(D, D)/3) # Peclet number
cfl = 2*dt*vnorm/h                # local time step

DG = dolfinx.fem.FunctionSpace(mesh, ("DG", 0))
cfl_expr = dolfinx.fem.Expression(cfl,DG.element.interpolation_points())
cfl_fun = dolfinx.fem.Function(DG)
cfl_fun.interpolate(cfl_expr)
maxcfl = np.max(cfl_fun.x.array)
print("maxcfl = ", maxcfl)


v = ufl.TrialFunction(DG)
cfl_vec = dolfinx.fem.petsc.assemble_vector(dolfinx.fem.form(cfl_fun*v/ufl.CellVolume(mesh)*ufl.dx))
print(cfl_vec.array-cfl_fun.x.array)
1 Like