Solve Method with TensorFunctionSpace

General
1

Hi, everyone!
I’ve been working with FEniCS for some time, and I’m trying to implement a modified version of the IPCS scheme, related to the Navier-Stokes equations, to a continuum model for active nematics.
Having said that, my problem is the following, I’ve a dynamical equation for a symmetric and traceless 2x2 tensor. But when I use the method SOLVE I get the following error:
solve(): incompatible function arguments. The following argument types are supported:
1. (A: dolfin.cpp.la.GenericLinearOperator, x: dolfin.cpp.la.GenericVector, b: dolfin.cpp.la.GenericVector, method: str = ‘lu’, preconditioner: str = ‘none’) → int

Invoked with: <dolfin.cpp.la.Matrix object at 0x7f487ba38f90>, Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement(‘Lagrange’, triangle, 1), dim=2), 7), TensorElement(FiniteElement(‘Lagrange’, triangle, 2), shape=(2, 2), symmetry={(1, 0): (0, 1)})), 17), <dolfin.cpp.la.Vector object at 0x7f487ad19a90>, ‘bicgstab’, ‘hypre_amg’

Here I provide some parts of my code:

Vector = VectorElement(‘Lagrange’, mesh.ufl_cell(), 2);
T = TensorFunctionSpace(mesh, ‘Lagrange’, 2, symmetry = True);
V = FunctionSpace(mesh, Vector)

#—Define strain rate tensor symmetric gradient
def epsilon(u):
return sym(nabla_grad(u))

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

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

#—Define N_ij tensor
def N(epsilon, omega, Q):
N1 = dot(lamb * epsilon + omega, Q)
N2 = 1/3. *dot(lamb * epsilon - omega, Q)
N3 = -2/3. * inner(Q, nabla_grad(u)) * lamb * (Q + Identity(len(u)))
return N1 + N2 #+ N3

#—Define variational problem for step 1
F1 = inner((Q - Q_n) / k, psi)*dx + inner(dot(u_n, nabla_grad(Q_n)), psi)*dx \

  • inner(nabla_grad(Q_n), nabla_grad(psi))*dx \
  • inner(L(epsilon(u_n), omega(u_n)), psi)*dx \
  • inner(N(epsilon(u_n), omega(u_n), Q_n), psi)*dx

a1 = lhs(F1)
L1 = rhs(F1)
A1 = assemble(a1)

Time-stepping

t = 0
for n in range(ns):
# Update current time
t += dt

# Step 1: Tentative velocity step
b1 = assemble(L1)
solve(A1, Q_, b1, 'bicgstab', 'hypre_amg')

Hope you can help me.
Thanks!

2

Make sure that the code is formatted with markdown syntax, ie.

```python
# add code here

```

Also make sure that the code can reproduce the issue,

I cant spot what Q_ is.
If it is a function, could you try to send in Q_.vector() instead?

3

Yep, sorry with the syntax, I’ll fix it in the following discussions.
Here’s some context:

#---Define mesh and functions spaces
mesh = UnitSquareMesh(8,8)
Vector = VectorElement('Lagrange', mesh.ufl_cell(), 2);
T = TensorFunctionSpace(mesh, 'Lagrange', 2, symmetry = True);
V = FunctionSpace(mesh, Vector)

Then I define the trial and test functions over the functions spaces

#---Define trial and test functions
Q   = TrialFunction(T)
psi = TestFunction(T)
u   = TrialFunction(V)
v   = TestFunction(V)

#---Define functions for sulutions at previous and current time steps
Q_n = Function(T)
Q_   = Function(T)
u_n  = Function(V)
u_    = Function(V)

#---Define expressions used in variational forms
lamb = Constant(lamb)
xi   = Constant(xi)
eta  = Constant(eta)
k  = Constant(dt)

#---Define strain rate tensor symmetric gradient
def epsilon(u):
    return sym(nabla_grad(u))

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

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

#---Define N_ij tensor
def N(epsilon, omega, Q):
    N1 = dot(lamb * epsilon + omega, Q)
    N2 = 1/3. *dot(lamb * epsilon - omega, Q)  
    N3 = -2/3. * inner(Q, nabla_grad(u)) * lamb * (Q + Identity(len(u)))  
    return N1 + N2 #+ N3

#---Define variational problem for step 1
F1 = inner((Q - Q_n) / k, psi)*dx + inner(dot(u_n, nabla_grad(Q_n)), psi)*dx   \
   + inner(nabla_grad(Q_n), nabla_grad(psi))*dx  \
   - inner(L(epsilon(u_n), omega(u_n)), psi)*dx   \
   + inner(N(epsilon(u_n), omega(u_n), Q_n), psi)*dx   

a1 = lhs(F1)
L1 = rhs(F1)
A1 = assemble(a1)

# Time-stepping
t = 0
for n in range(ns):    
    # Update current time
    t += dt

    # Step 1: First variational problem
    b1 = assemble(G1)
    solve(A1, Q_, b1, 'bicgstab', 'hypre_amg')

And in relation with your question, I’m not sure of what the method vector() does to the equation, but I could try with that.

4

.vector() accesses the GenericVector that contains the degrees of freedom (unknowns) that you solve for.
This would make it match the input

5

That fix the output, but now I’ve a different problem; the output it’s just the null vector, for every step in time. Here’s the code fixed (with the new problem), sorry for the basic questions

# Sub domain for Periodic boundary condition
class PeriodicBoundary(SubDomain):

    # Left boundary is "target domain" G
    def inside(self, x, on_boundary):
        return bool(x[0] < DOLFIN_EPS and x[0] > -DOLFIN_EPS and on_boundary)

    # Map right boundary (H) to left boundary (G)
    def map(self, x, y):
        y[0] = x[0] - 1.0
        y[1] = x[1]
pbc = PeriodicBoundary()
#---Tensor-Vector elemet
#---Define the mesh
mesh = UnitSquareMesh(8,8)
#Tensor = TensorElement('Lagrange', mesh.ufl_cell(), 2, symmetry = True)
Vector = VectorElement('Lagrange', mesh.ufl_cell(), 2);
T = TensorFunctionSpace(mesh, 'Lagrange', 2, symmetry = True, constrained_domain=pbc);
V = FunctionSpace(mesh, Vector, constrained_domain=pbc)
#---Declare constants of the model
#alpha = ; alpha = 1 / Pe
#beta = ; beta = 1 / Chi
eta  = 0.1      #---viscosity
xi   = 0.1      #---activity parameter
lamb = 0.1      #---tumbling parameter
T    = 1.0      #---final time
ns   = 1000     #---number of time steps
dt   = T / ns   #---time step size
#---Define trial and test functions
Q   = TrialFunction(T)
psi = TestFunction(T)
u   = TrialFunction(V)
v   = TestFunction(V)

#---Define functions for sulutions at previous and current time steps

Q_n = Function(T)
Q_  = Function(T)
u_n = Function(V)
u_  = Function(V)
#---Define expressions used in variational forms
lamb = Constant(lamb)
xi   = Constant(xi)
eta  = Constant(eta)
k  = Constant(dt)

#---Define strain rate tensor symmetric gradient
def epsilon(u):
    return sym(nabla_grad(u))

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

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

#---Define N_ij tensor
def N(epsilon, omega, Q):
    N1 = dot(lamb * epsilon + omega, Q)
    N2 = 1/3. *dot(lamb * epsilon - omega, Q)  
    N3 = -2/3. * inner(Q, nabla_grad(u)) * lamb * (Q + Identity(len(u)))  
    return N1 + N2 #+ N3
#---Define variational problem for step 1
F1 = inner((Q - Q_n) / k, psi)*dx + inner(dot(u_n, nabla_grad(Q_n)), psi)*dx   \
   + inner(nabla_grad(Q_n), nabla_grad(psi))*dx  \
   - inner(L(epsilon(u_n), omega(u_n)), psi)*dx   \
   + inner(N(epsilon(u_n), omega(u_n), Q_n), psi)*dx   

a1 = lhs(F1)
L1 = rhs(F1)
A1 = assemble(a1)
# Time-stepping
t = 0
for n in range(ns):
    
    # Update current time
    t += dt

    # Step 1: First variational problem
    b1 = assemble(L1)
    #[bc.apply(b1) for bc in bcu]
    solve(A1, Q_.vector(), b1, 'bicgstab', 'hypre_amg')

    # Update previous solution
    Q_n.assign(Q_)
    u_n.assign(u_)

    # Update progress bar
    set_log_level(LogLevel.ERROR)
    print('Q max:', Q_.vector().get_local().max())

And the output is:

....
Q max: 0.0
Q max: 0.0
Q max: 0.0
Q max: 0.0
Q max: 0.0
Q max: 0.0
Q max: 0.0
Q max: 0.0
....
6

Im not really sure what you expect. The initial u_n and Q_n are Zero, so all terms with those are Zero in F1

Thus you are only left with inner(Q/k,psi)*dx, with no Dirichlet conditions, Which has the solution 0.

7

Oh, I see, then I should work with the Dirichlet conditions.
And finally, is there a way in which I start with u_n and Q_n different from zero?
In Navier-Stokes you have an input flow which makes sure those values are not zero, but in my problem, I’m not sure how to implement that.
Thanks a lot!

8

You can interpolate or project any function into u_n and Q_n.