I have encountered an issue when computing L2 norms for solution fields in different FEniCS Legacy installations. My flow problem employs mixed finite element spaces: a Discontinuous Galerkin (DG) element for variable D, and a Brezzi-Douglas-Marini (BDM) element for variable sig. The L2 norms are computed using the standard formulations assemble(inner(D, D) * dx) and assemble(inner(sig, sig) * dx).
In my primary development environment (Linux with apt-installed FEniCS, versions DOLFIN: 2019.2.0.64.dev0, UFL: 2024.2.0), these computations work correctly for both variables. However, in a secondary Conda/Mamba environment (versions DOLFIN: 2019.1.0, UFL: 2019.1.0), I encounter a compilation error specifically when assembling the quadratic form for the BDM field sig.
Traceback (most recent call last):
File "/.../Stokes-2D.py", line 94, in <module>
sig_L2_norm = assemble(inner(sig,sig)*dx)
^^^^^^^^^^^^^^^^^^^^^^^^^^^
............. a lot of errors of the same kind..................
File "/.../mambaforge/envs/FEniCS/lib/python3.11/site-packages/ufl/algorithms/apply_function_pullbacks.py", line 151, in apply_single_function_pullbacks
assert len(gsh) == 1
^^^^^^^^^^^^^
AssertionError
Can anyone help me, please ?
Best regards,
"""
This code solves a mixed FEM for a stationary Stokes for an incompressible flow in
symmetric-pseudostress-velocity-vorticity formulation
"""
# <----- preamble ----->
from fenics import *
import numpy as np
import sympy as sp
# create classes for defining boundary
def gamma(x,on_boundary):
return on_boundary
# <----- some useful classes ----->
def skew(scalar):
return as_tensor( [[0,scalar], [-scalar,0]] )
def trnull(vector):
return as_tensor( [[vector[0],vector[1]], [vector[2],-vector[0]]] )
def conv(v1,v2):
return as_tensor( [[v1[0]*v2[0],v1[0]*v2[1]], [v1[1]*v2[0],v1[1]*v2[1]]] )
def Div(tensor):
return as_vector( [tensor[0,0].dx(0) + tensor[0,1].dx(1), tensor[1,0].dx(0) + tensor[1,1].dx(1)] )
def fsig(tensor,const):
return as_tensor( [[tensor[0,0]+const, tensor[0,1]], [tensor[1,0], tensor[1,1]+const]] )
def Tran(tensor):
return as_tensor([ [tensor[0,0],tensor[1,0]],[tensor[0,1],tensor[1,1]] ])
# <----- model parameters and auxiliary symbolic expressions ----->
Id = Identity(2)
x,y = sp.symbols('x[0] x[1]')
nu = 1.0
# <----- Data and boundary conditions ----->
lbd = (-8.0*sp.pi)/(nu**(-1) + sp.sqrt(nu**(-2) + 16.*sp.pi**2))
ue_c1 = 1. - sp.exp(lbd*x)*sp.cos(2.*sp.pi*y)
ue_c2 = (lbd/(2*sp.pi))*sp.exp(lbd*x)*sp.sin(2.*sp.pi*y)
ff1 = 8.0*sp.pi*sp.exp(-16.0*sp.pi*x/(1.0 + sp.sqrt(1.0 + 16.0*sp.pi**2)))/(1.0 + sp.sqrt(1.0 + 16.0*sp.pi**2)) - 2.0*sp.pi**2*sp.exp(-8.0*sp.pi*x/(1.0 + sp.sqrt(1.0 + 16.0*sp.pi**2)))*sp.cos(2*sp.pi*y) + 32.0*sp.pi**2*sp.exp(-8.0*sp.pi*x/(1.0 + sp.sqrt(1.0 + 16.0*sp.pi**2)))*sp.cos(2*sp.pi*y)/(1.0 + sp.sqrt(1.0 + 16.0*sp.pi**2))**2
ff2 = -8.0*sp.pi**2*sp.exp(-8.0*sp.pi*x/(1.0 + sp.sqrt(1.0 + 16.0*sp.pi**2)))*sp.sin(2*sp.pi*y)/(1.0 + sp.sqrt(1.0 + 16.0*sp.pi**2)) + 128.0*sp.pi**2*sp.exp(-8.0*sp.pi*x/(1.0 + sp.sqrt(1.0 + 16.0*sp.pi**2)))*sp.sin(2*sp.pi*y)/(1.0 + sp.sqrt(1.0 + 16.0*sp.pi**2))**3
ue = Expression(( sp.printing.ccode(ue_c1), sp.printing.ccode(ue_c2) ), degree = 5)
ff = Expression(( sp.printing.ccode(ff1), sp.printing.ccode(ff2) ), degree = 5)
# <---- Stokes solver ------>
mesh = RectangleMesh(Point(-1./2., 0.), Point(3./2., 2.), 32, 32)
n = FacetNormal(mesh)
# <----- Polynomial order degree ----->
order = 0
# <----- define function spaces ----->
HD = VectorElement("DG" ,mesh.ufl_cell(), order+1, dim=3) # vector P_k discontinuous
Hsig = VectorElement("BDM",mesh.ufl_cell(), order+1) # tensor BDM_k+1
Hu = VectorElement("DG" ,mesh.ufl_cell(), order) # vector P_k discontinuous
Hgamm = FiniteElement("DG" ,mesh.ufl_cell(), order) # scalar P_k discontinuous
L = FiniteElement("R" ,mesh.ufl_cell(), 0) # Lagrange multiplier space
Xh = FunctionSpace(mesh,MixedElement([HD,Hsig,Hu,Hgamm,L]))
# <----- defining trial and test functions ----->
(D, sig, u, gamm, l1) = TrialFunctions(Xh)
(E, tau, v, delt, t1) = TestFunctions(Xh)
# <---- involved forms ----->
A = nu*inner(trnull(D),trnull(E))*dx - inner(sig,trnull(E))*dx - inner(tau,trnull(D))*dx
B2 = -dot(u,Div(tau))*dx - inner(skew(gamm),tau)*dx
B2t = -dot(v,Div(sig))*dx - inner(skew(delt),sig)*dx
M1 = l1*tr(tau)*dx
M2 = t1*tr(sig)*dx
AA = A + B2 + B2t + M1 + M2
FD1 = -dot(tau*n,ue)*ds
Ff2 = dot(ff,v)*dx
FF = FD1 + Ff2
# <---- solving the problem ---->
sol = Function(Xh)
solve(AA == FF, sol, solver_parameters = {"linear_solver": "umfpack"})
(D, sig, u, gamm, l1) = sol.split()
# <----- L2 norms ----->
u_L2_norm = sp.sqrt( assemble(inner(u,u)*dx) )
D_L2_norm = sp.sqrt( assemble(inner(D,D)*dx) )
gamm_L2_norm = sp.sqrt( assemble(inner(gamm,gamm)*dx) )
l1_L2_norm = sp.sqrt( assemble(inner(l1,l1)*dx) )
#sig_L2_norm = assemble(inner(sig,sig)*dx)
# <----- functions type ----->
print("# <----- functions type ----->")
print("u type", type(u))
print("D type", type(D))
print("gamm type", type(gamm))
print("l1 type", type(l1))
print("sig type", type(sig))
# <----- functions value shape ----->
print("# <----- functions value shape ----->")
print("u shape",u.function_space().ufl_element().value_shape())
print("D shape",D.function_space().ufl_element().value_shape())
print("gamm shape",gamm.function_space().ufl_element().value_shape())
print("l1 shape",l1.function_space().ufl_element().value_shape())
print("sig shape",sig.function_space().ufl_element() .value_shape())
# <----- FEniCS version ----->
import dolfin
import ufl
print("# <----- FEniCS version ----->")
print(f'DOLFIN: {dolfin.__version__}')
print(f'UFL: {ufl.__version__}')