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
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