SOLVE method error with TensorElement using Mixed Formulation

Errors
1

Hello, I’m soling a set of coupled partial differential equations, similar to the Navier-Stokes, but with tensor and vector functions. I’ve been trying to implement a modified version of the IPCS for my particular problem, and right now I’m stuck in the first step, in which I try to solve for a tensor equation. The error I get is this: AttributeError: 'ListTensor' object has no attribute 'vector'.
Could you please help me to understand my issue?
Here I provide the code in which I’m currently working.
(Sorry if it’s too long)

Thanks!

## Define vector elements and function spaces
QE = TensorElement("P", mesh.ufl_cell(), 2)
UE = VectorElement("P", mesh.ufl_cell(), 2)
ME = MixedElement([QE, UE])
QU = FunctionSpace(mesh, ME)
Q = QU.sub(0)
U = QU.sub(1)
## Define boundaries
inflow   = 'near(x[0], 0)'                     
outflow  = 'near(x[0], 1.0)'
walls    = 'near(x[1], 0) || near(x[1], 0.2)'

## Define inflow profile for the velocity
inflow_profile = ('4.0*1.5*x[1]*(0.1 - x[1]) / pow(0.1, 2)', 
                  '0')

## Define boundary conditions for the velocity
## For mixed systems, we should specify which subsystem we set
## the boundary condition. This is done via the sub(0) and sub(1)
bcu_inflow = DirichletBC(U, 
                         Expression(inflow_profile, 
                                    degree=2), 
                         inflow)
bcu_walls = DirichletBC(U, 
                        Constant((0, 0)), 
                        walls)
bcu = [bcu_inflow, bcu_walls]         ## Colect all boundary conditions for the velocity
## Define boundary conditions for the order parameter
## Define initial value for the order parameter
initial_order = Expression((('0.8', '0.8'), ('0.8', '-0.8')), degree=2)

bcq_inflow = DirichletBC(Q, 
                          initial_order, 
                          outflow)
bcq_walls = DirichletBC(Q, 
                        initial_order, 
                        walls)
bcq = [bcq_inflow, bcq_walls]        ## Colect all boundary conditions for the order parameter
## Define trial and test functions
qu = TrialFunction(QU)
q, u = split(qu)
rv   = TestFunction(QU)
r, v = split(rv)
# Define functions for solutions at previous and current time steps
qu_n     = Function(QU)
q_n, u_n = split(qu_n)
qu_      = Function(QU)
q_, u_   = split(qu_)
## Define epsilon tensor
def epsilon(u):
    return sym(nabla_grad(u))

## Define vorticity tensor
def omega(u):
    return skew(nabla_grad(u))

## Define L tensor
def L(u):
    L1 = (1/3.) * dot((lamb*epsilon(u) + omega(u)), Identity(len(u)))
    L2 = (1/3.) * dot((lamb*epsilon(u) - omega(u)), Identity(len(u)))
    return L1 + L2

## Define N tensor
def N(q, u):
    N1 = dot(lamb*epsilon(u) + omega(u), q)
    N2 = dot(q, lamb*epsilon(u) - omega(u))
    N3 = -2*lamb*(q + Identity(mesh.geometry().dim())/3.) * inner(q, nabla_grad(u))
    return N1 + N2 + N3
## Define expressions used in variational forms
k1 = Constant(dt)
k2 = Constant(Pe)
## Define variational problem for step 1
F1 = inner((q - q_n) / k1, r)*dx                      \
    + inner(dot(u_n, nabla_grad(q_n)),r)*dx           \
    - (1/k2) * inner(nabla_grad(q), nabla_grad(r))*dx \
    - inner(L(u_n), r)*dx                             \
    - inner(N(q_n, u_n), r)*dx
a1 = lhs(F1)
L1 = rhs(F1)

# Assemble matrices
A1 = assemble(a1)

# Apply boundary conditions to matrices
[bc.apply(A1) for bc in bcq]
# Time-stepping
t = 0
for n in range(num_steps):
    
    # Update current time
    t += dt

    # Step 1: Tentative velocity step
    b1 = assemble(L1)
    [bc.apply(b1) for bc in bcq]
    solve(A1, q_.vector(), b1, 'bicgstab', 'hypre_amg')
2

If you want to use an IPCS like method, you shouldn’t create a mixed element and a function space based on that element. The mixed element means that any matrix based on the trial or test functions from such a space should generate a matrix with number of dofs = dofs in QE + dofs in UE, while what you want is a single matrix for dofs in QE and then another one for UE. See:
https://bitbucket.org/fenics-project/dolfin/src/master/python/demo/documented/navier-stokes/demo_navier-stokes.py

1 Like
3

Thanks for the answer!
I’ve changed that and set the code into two different spaces, just like you said, and now, my code works, but I have another question, I need that the Tensor Space only admits tensors symmetric and traceless. I found that you can add symmetry = True in the declaration of the TensorFunctionSpace but I’m not sure if that is correct, could you help me with that?
Sorry, I know that it’s a silly question, but I’m a rookie using FEniCS.

4

What do you mean by traceless?
Do you mean Trace(A)=0 or A_ii = 0 for all i?

5

I mean, Trace(A) = 0.
The problem that I’ve been trying to solve is a modified version of the dynamic equations for active nematics and in this case, the order parameter, which in 2D is a second order symmetric and with Trace(Q) = 0, is coupled with a velocity field.
That’s why I want to found a way to implement that in FEniCS.
Thank you!

6

That would be an additional constraint to add to your variational form.

You could try to do something like:
https://fenicsproject.org/qa/1281/symmetric-zero-trace-tensor-space/
The equivalent in 2D has been shown in https://www.mn.uio.no/math/english/research/projects/feec-a/events/conferences-and-workshops/complex-materials/dassbach.pdf

7

Awesome, thank you for the recommendations!