Adding expressions with non-matching form arguments ('v_0^1', 'v_1^0') vs ('v_0^0', 'v_1^1')

variational formulation
1

Hi All,

This bothers me. I’ve read similar posts (1) two unknowns in a form, (2) form a is not a bilinear, (3) no trial function for nonlinear problem. I tried multiple times but couldn’t figure out. A MWE is listed here:

from fenics import *

h = 3.0e6
lmbda = 6.0e6
eta0 = 5.0e22
a = ln(2.5)
dT = 3.0e3
b = ln(2.0)
d = 1.5e5/h
#unit vector downwards
k = Constant((0.0, -1.0))
#mass diffusivity
kch = 1.0e-10
#thermal diffusivity
kth = 1.0e-6
#mass diffusivity/thermal diffusivity
kct = kch/kth
#Large number for stabilization
ah = 50.0

#real numbers: less efficient, for those less updated 
aspect = lmbda/h
alpha = 2.0e-5
g = 10.0
rho0 = 3.1e3
drho = 1.85e2

nx = 160
ny = 80

mesh = RectangleMesh(Point(0.0, 0.0), Point(aspect, 1.0), nx, ny)
V = VectorFunctionSpace(mesh, 'Lagrange', 2)
Q = FunctionSpace(mesh, 'Lagrange', 1)
W = MixedFunctionSpace(V, Q)
(u, p) = TrialFunctions(W)
(u_n, p_n) = TrialFunctions(W)
(u_p, p_p) = TrialFunctions(W)
(v, q) = TestFunctions(W)
up = Function(W)

#DG for temperature and concentration in advection-dominated flow
TC = FunctionSpace(mesh, 'DG', 1)
T = TrialFunction(TC)
C = TrialFunction(TC)
T_n = TrialFunction(TC)
C_n = TrialFunction(TC)
T_p = TrialFunction(TC)
C_p = TrialFunction(TC)
T_t = TestFunction(TC)
C_t = TestFunction(TC)

hsize = CellDiameter(mesh)
x = SpatialCoordinate(mesh)
n = FacetNormal(mesh)

#Boundaries
sides = 'near(x[0], 0) || near(x[0], aspect)'
top = 'near(x[1], 1)'
bottom = 'near(x[1], 0)'

#Boundary conditions
bc_Ttop = DirichletBC(TC, Constant(0.0), top)
bc_Tbottom = DirichletBC(TC, Constant(1.0), bottom)
bc_Cbottom = DirichletBC(TC, Constant(1.0), bottom)
bc_vbottom = DirichletBC(W.sub_space(0), Constant((0.0, 0.0)), bottom)

Tbcs = [bc_Ttop, bc_Tbottom]
vbcs = [bc_vbottom]
Cbcs = [bc_Cbottom]

#Constant objects: fast evaluation, for those updated frequently
RaT = alpha*g*rho0*dT*h**3/(kth*eta0)
RaC = drho*g*h**3/(kth*eta0)


def eta(T):
    #T and pressure dependent viscosity
    return eta0*exp(-a*T/dT + b*(h-x[1])/h)

def epsilon(u):
    return (as_tensor([[u[0].dx(0), 0.5*u[0].dx(1) + 0.5*u[1].dx(0)], [0.5*u[1].dx(0) + 0.5*u[0].dx(1), u[1].dx(1)]]))

def sigma(T,u):
    return 2.0*eta(T)*epsilon(u)

def Tini(aspect, RaT):
    #Initial T    
    q0 = aspect**(7.0/3.0)*(RaT/(2.0*sqrt(pi)))**(2.0/3.0)/(1.0+aspect**4.0)**(2.0/3.0)
    Q0 = 2.0*sqrt(aspect/(pi*q0))
    Tu = Expression('0.5*erf(0.5*(1.0-x[1])*sqrt(q0/x[0]))', degree=0, q0=q0)
    Tl = Expression('1.0 - 0.5*erf(0.5*x[1]*sqrt(q0/(aspect-x[0])))', degree=0, q0=q0, aspect=aspect)
    Tr = Expression('0.5 + 0.5*Q0*sqrt(q0/(pi+pi*x[1]))*exp(-pow(x[0],2)*q0/(4.0*x[1]+4.0))', degree=0, Q0=Q0, q0=q0)
    Ts = Expression('0.5 - 0.5*Q0*sqrt(q0/(2.0*pi-x[1]*pi))*exp(-pow(aspect-x[0],2)*q0/(8.0-4.0*x[1]))', degree=0, Q0=Q0, q0=q0, aspect=aspect)
    T0 = Expression('Tu + Tl + Tr + Ts - 1.5', degree=1, Tu=Tu, Tl=Tl, Tr=Tr, Ts=Ts)
    #Mapping all T in [0, 1]
    T0 = Expression('T0 < 0.0 ? 0.0 : T0', degree=0, T0=T0)
    T0 = Expression('T0 > 1.0 ? 1.0 : T0', degree=0, T0=T0)
    
    return T0

def Cini(d):
    #Initial C
    C0 = Expression('x[1] <= d ? 1.0 : 0.0', degree=0, d=d)#degree=2 will give too much interpolation.
    
    return C0

def Stokes(T,C):
    #Form of Stokes system: symmetric but not positive-definite
    a = inner(grad(u), grad(v))*dx#inner(sigma(T,u), epsilon(v))*dx##
    a += p*div(v)*dx
    a += q*div(u)*dx 
    L = (RaT*T - RaC*C)*inner(k, v)*dx
    
    #preconditioner
    pc0 = inner(grad(u), grad(v))*dx + p*q*dx
    
    return a, L, pc0
    
#Initialize Temperature and Concentration
T0 = Tini(aspect, RaT)
T_n = interpolate(T0, TC)
C0 = Cini(d)
C_n = project(C0, TC)

#Compute initial velocity u0 from T0, C0
solver = KrylovSolver('tfqmr', 'petsc_amg')#or mines
a, L, pc0 = Stokes(T_n,C_n)
A, b = assemble_system(a, L, vbcs) 
pc, ptmp = assemble_system(pc0, L, vbcs)
solver.set_operators(A, pc)
solver.solve(up.vector(), b)
u_n, p_n = up.split()

This code is a part of Chapter 31 from The FEniCS book (Dynamic simulations of convection in the Earth’s mantle). The error is thrown at line A, b = assemble_system(a, L, vbcs): Adding expressions with non-matching form arguments (‘v_0^1’, ‘v_1^0’) vs (‘v_0^0’, ‘v_1^1’). I guess the error could be either of previous 3 posts but I cannot tell. Any ideas? Thanks!

2

It seems to me that you are trying to solve the stokes problem once T and C are already computed. If that is the case, then

must not use the trail functions T and C, but use instead the solutions T_n, C_n.

3

Thanks! Yes, T and C have been initialized as T0 and C0 which should be interpolated to Function objects T_n and C_n. I just realized T_n and C_n in my old code are TrialFunction objects, which is not right. After the correction, i.e. from

(u_n, p_n) = TrialFunctions(W)
(u_p, p_p) = TrialFunctions(W)

to

up_n = Function(W)
up_p = Function(W)
(u_n, p_n) = up_n.split()
(u_p, p_p) = up_n.split()

and from

T_n = TrialFunction(TC)
C_n = TrialFunction(TC)
T_p = TrialFunction(TC)
C_p = TrialFunction(TC)

to

T_n = Function(TC)
C_n = Function(TC)
T_p = Function(TC)
C_p = Function(TC)

Now up_n, up_p, T_n, C_n, T_p, C_p are all Function objects (numerical approximations). The line a, L, pc0 = Stokes(T_n,C_n) after initialization of T and C should pass Function objects T_n and C_n to function Stokes, giving form a, L and preconditioned pc0. It looks fine for now. When the next line is executed A, b = assemble_system(a, L, vbcs), it gives the same error:
Adding expressions with non-matching form arguments ('v_0^1', 'v_1^0') vs ('v_0^0', 'v_1^1'). What does v_0^1 mean exactly? I have no further clue to debug this, any suggestions?

4

Also I tested the a demo code of Stokes from FEniCS project, it gave me the same error:

raise ArityMismatch("Adding expressions with non-matching form arguments {0} vs {1}.".format(_afmt(a), _afmt(b)))

ufl_legacy.algorithms.check_arities.ArityMismatch: Adding expressions with non-matching form arguments ('v_0^1', 'v_1^0') vs ('v_0^0', 'v_1^1').
5

Please note that you are looking at a very old version of the demo (1.6.0). Please consider at least: https://fenicsproject.org/olddocs/dolfin/2019.1.0/python/demos/stokes-iterative/demo_stokes-iterative.py.html?highlight=stokes%20iterative
I cannot reproduce any error messages with the code above.

It seems like you are trying to replicate very old dolfin code. Changing your finite element function space definition to:

el_u = VectorElement("Lagrange", mesh.ufl_cell(), 2)
el_p = FiniteElement("Lagrange", mesh.ufl_cell(), 1)
W = FunctionSpace(mesh, MixedElement([el_u, el_p]))

and

bc_vbottom = DirichletBC(W.sub(0), Constant((0.0, 0.0)), bottom)

the code does not throw an error.

6

Thanks a lot but why! Just because it’s very old?

7

Yes, you are looking at the API from over 10 years ago. Alot changed since then.

1 Like