Equivalent to dolfin.fem.norms.norm with dolfinx

Hi,
Is there an equivalent to dolfin.fem.norms.norm with dolfinx?
Im trying to get the L2 norm of a difference of two solutions: ur-ul. Is there a simple way to do that?
Thanks

See for instance: Error control: Computing convergence rates — FEniCSx tutorial

Thanks, but the error_L2 gives the following error:

    degree = uh.function_space.ufl_element().degree()
AttributeError: 'ListTensor' object has no attribute 'function_space'

Without a minimal reproducible example, I cannot help you

Yes, but it is difficult to get a minimal one. Sorry.
Maybe the complete code is not to complicated:

import numpy as np
import numpy

import ufl
from dolfinx import fem, io, mesh, plot,  nls, log
from ufl import div, dx, grad, inner

from mpi4py import MPI
from petsc4py.PETSc import ScalarType
from petsc4py import PETSc


from scipy.optimize import minimize
import scipy.linalg
import math
import logging
 
from logging.handlers import RotatingFileHandler

msh = mesh.create_rectangle(comm=MPI.COMM_WORLD,
                            points=((0.0, 0.0), (1.0, 1.0)), n=(40, 40),
                            cell_type=mesh.CellType.triangle,)

nu = 0.001

# Function to mark x = 0, x = 1, y = 0 
def noslip_boundary(x):
    return np.logical_or(np.logical_or(np.isclose(x[0], 0.0), np.isclose(x[0], 1.0)), np.isclose(x[1], 0.0))

# Function to mark the lid (y = 1)
def lid(x):
    return np.isclose(x[1], 1.0)

# Lid velocity#
def lid_velocity_expression(x):
    return np.stack((np.ones(x.shape[1]), np.zeros(x.shape[1])))

P2 = ufl.VectorElement("Lagrange", msh.ufl_cell(), 2)
P1 = ufl.FiniteElement("Lagrange", msh.ufl_cell(), 1)
V, Q = fem.FunctionSpace(msh, P2), fem.FunctionSpace(msh, P1)

# Create the function space
TH = P2 * P1
W = fem.FunctionSpace(msh, TH)
W0, _ = W.sub(0).collapse()

# No slip boundary condition
noslip = fem.Function(V)
facets = mesh.locate_entities_boundary(msh, 1, noslip_boundary)
dofs = fem.locate_dofs_topological((W.sub(0), V), 1, facets)
bc0 = fem.dirichletbc(noslip, dofs, W.sub(0))

# Driving velocity condition u = (1, 0) on top boundary (y = 1)
lid_velocity = fem.Function(W0)
lid_velocity.interpolate(lid_velocity_expression)
facets = mesh.locate_entities_boundary(msh, 1, lid)
dofs = fem.locate_dofs_topological((W.sub(0), V), 1, facets)
bc1 = fem.dirichletbc(lid_velocity, dofs, W.sub(0))

# Since for this problem the pressure is only determined up to a constant, we pin the pressure at the point (0, 0)
zero = fem.Function(Q)
zero.x.set(0.0)
dofs = fem.locate_dofs_geometrical(
    (W.sub(1), Q), lambda x: np.isclose(x.T, [0, 0, 0]).all(axis=1))
bc2 = fem.dirichletbc(zero, dofs, W.sub(1))

# Collect Dirichlet boundary conditions
bcs = [bc0, bc1, bc2]

# Define variational problem
w = ufl.TrialFunction(W)

(u, p) = ufl.split(w)
(v, q) = ufl.TestFunctions(W)

log.set_log_level(log.LogLevel.INFO)
u_k = fem.Function(V)
def initU(x):
    k = [1.0, 0.0]
    return np.vstack((
        x[0] * k[0],
        x[1] * k[1]
    ))
u_k.interpolate(initU)  # previous (known) u

# Tentative velocity step
abilin = inner(grad(u)*u_k, v)*dx + nu*inner(grad(u), grad(v))*dx - div(v)*p*dx - q*div(u)*dx

Fext = fem.Function(V)
Fext.interpolate(initU)  # 

ZeroC = fem.Constant(msh,0.0)

Lrhs = inner(Fext,v) * dx + ZeroC*q*dx


def bilin(uk_loc):
    return inner(grad(u)*uk_loc, v)*dx + nu*inner(grad(u), grad(v)) * \
    dx - div(v)*p*dx - q*div(u)*dx

def g(uk):
    ak = bilin(uk)
    problem_loc = fem.petsc.LinearProblem(ak, Lrhs, bcs=bcs, petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
    w_loc = problem_loc.solve()
    return w_loc


# Picard iterations
eps = 1.0           # error measure ||u-u_k||
tol = 1.0E-6     # tolerance
maxiter = 60        # max no of iterations allowed

xit0 = u_k

############### log  ###########################
# création de l'objet logger qui va nous servir à écrire dans les logs
logger = logging.getLogger()
# on met le niveau du logger à DEBUG, comme ça il écrit tout
logger.setLevel("INFO")


def error_L2(uh, u_ex, degree_raise=3):
    # Create higher order function space
    degree = uh.function_space.ufl_element().degree()
    family = uh.function_space.ufl_element().family()
    mesh = uh.function_space.mesh
    W = fem.FunctionSpace(mesh, (family, degree+ degree_raise))
    # Interpolate approximate solution
    u_W = fem.Function(W)
    u_W.interpolate(uh)

    # Interpolate exact solution, special handling if exact solution
    # is a ufl expression or a python lambda function
    u_ex_W = fem.Function(W)
    if isinstance(u_ex, ufl.core.expr.Expr):
        u_expr = fem.Expression(u_ex, W.element.interpolation_points)
        u_ex_W.interpolate(u_expr)
    else:
        u_ex_W.interpolate(u_ex)
    
    # Compute the error in the higher order function space
    e_W = fem.Function(W)
    e_W.x.array[:] = u_W.x.array - u_ex_W.x.array
    
    # Integrate the error
    error = fem.form(ufl.inner(e_W, e_W) * ufl.dx)
    error_local = fem.assemble_scalar(error)
    error_global = mesh.comm.allreduce(error_local, op=MPI.SUM)
    return np.sqrt(error_global)

pointfixe = True
anderson = True
if(pointfixe):
    ## point fixe
    iter=0
    while eps > tol and iter < maxiter:
        iter += 1
        w_kp = g(u_k)
        (u_kp, p_kp) = ufl.split(w_kp)
        eps = error_L2(u_kp,u_k)
        u_k = u_kp.copy()   # update for next iteration
        logger.info("======================")
        logger.info("Point fixe")
        logger.info("iteration: "+ str(iter))
        logger.info("||r(k)|| = "+str(eps))

use u_kp, p_kp = w_kp.split(deepcopy=True) instead, as the ufl.split does not separate the function spaces of the underlying object.

Thank you. First I got:

TypeError: Function.split() got an unexpected keyword argument 'deepcopy'

then I removed the keyword, and I got:

 File "/home/mc/anaconda3/envs/fenicsx-env/lib/python3.10/site-packages/dolfinx/fem/function.py", line 324, in _
    self._cpp_object.interpolate(u._cpp_object, cells)
RuntimeError: Interpolation: elements have different value dimensions

Here is an improved and generalized version of the error-norm:

def error_L2(uh, u_ex, degree_raise=3):
    # Create higher order function space
    ufl_el = uh.function_space.ufl_element()
    new_el = ufl_el.reconstruct(degree=degree_raise+ufl_el.degree())
    mesh = uh.function_space.mesh
    W = fem.FunctionSpace(mesh, new_el)
    # Interpolate approximate solution
    u_W = fem.Function(W)
    u_W.interpolate(uh)

    # Interpolate exact solution, special handling if exact solution
    # is a ufl expression or a python lambda function
    u_ex_W = fem.Function(W)
    if isinstance(u_ex, ufl.core.expr.Expr):
        u_expr = fem.Expression(u_ex, W.element.interpolation_points())
        u_ex_W.interpolate(u_expr)
    else:
        u_ex_W.interpolate(u_ex)

    # Compute the error in the higher order function space
    e_W = fem.Function(W)
    e_W.x.array[:] = u_W.x.array - u_ex_W.x.array

    # Integrate the error
    error = fem.form(ufl.inner(e_W, e_W) * ufl.dx)
    error_local = fem.assemble_scalar(error)
    error_global = mesh.comm.allreduce(error_local, op=MPI.SUM)
    return np.sqrt(error_global)

Thank you. Now, I have a “nan” in the above code inside the loop…

That has nothing to do with the norm, but the fact that

u_kp is Nan/inf.
Please revise your variational form/g function.

1 Like