Regarding implementation of the LUSolver

[atalasik]

Hello everyone,
I am using the Fenics-Dolfin 2019.1.0 version and trying to implement the LUSolver for solving the parabolic PDE. Following in my code and the error.

from fenics import *
from matplotlib import pyplot as plt
from mshr import *

# parameters, function spaces, initial data
n_steps = 20
tau     = 5e-4
t       = 0.0


#mesh    = RectangleMesh(Point(0,0),Point(1,1),128,128)
domain  = Circle(Point(0, 0),  1)
mesh    = generate_mesh(domain, 64)

U       = FunctionSpace(mesh,'CG',1)
old_u   = interpolate(Expression("2*(float)(rand())-1 + tanh(100*(x[0]-0.5))",degree=2),U)

u = Function(U)
v = TestFunction(U)
Res = (u-old_u)*v*dx + tau*inner(grad(u),grad(v))*dx
a = lhs(Res)
A = assemble(a)
solver = LUSolver(A)
solver.parameters['reuse_factorization'] = True
L = rhs(Res)

u.assign(old_u)
# main loop and plot solution each step
for i in range(1,n_steps+1):
  plt.figure()
  plot(old_u)
  t = t + tau
  b = assemble(L)
  solver.solve(u.vector(),b)
  print(et-st)
  old_u.assign(u)

The error I receive:

Form has no parts with arity 2.

---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
<ipython-input-21-b018cbdfd804> in <module>
     21 a = lhs(Res)
     22 A = assemble(a)
---> 23 solver = LUSolver(A)
     24 solver.parameters['reuse_factorization'] = True
     25 L = rhs(Res)

TypeError: __init__(): incompatible constructor arguments. The following argument types are supported:
    1. dolfin.cpp.la.LUSolver()
    2. dolfin.cpp.la.LUSolver(arg0: str)
    3. dolfin.cpp.la.LUSolver(A: dolfin.cpp.la.GenericLinearOperator, method: str = 'default')
    4. dolfin.cpp.la.LUSolver(comm: MPICommWrapper, A: dolfin.cpp.la.GenericLinearOperator, method: str = 'default')

Invoked with: 0.0

However, when I run the code on for advection-diffusion on the Fenics Bitbucket repsoitory(link to the respository) it runs without any error (contains a similar syntax for invoking the LUSolver).

I am very new to Fenics and I could be completely missing out on something, so any help would be greatly appreciated.

[dokken]

@atalasik

u = Function(U)

When you are solving a linear problem, you need to use a trial function for u (u = TrialFunction(V)),
as the form a has to be assembled into a matrix

[atalasik]

Thank you so much for the prompt response!!!

[atalasik]

This worked well also for linearised Nerst-Planck-Poisson system of equations.

However when I try to extend this to the Navier Stokes equations (using Taylor Hood finite elements as well as linearised convective term : using velocity at the previous time step) for the case of obtaining Von-Karman vortex street, I get an error while implementing boundary conditions (which otherwise works without LUSolver).

Following is my code:

import matplotlib.pyplot as plt
from dolfin import *
from mshr import *
import logging
import time

logging.getLogger('FFC').setLevel(logging.WARNING)
logging.getLogger('UFL').setLevel(logging.WARNING)

# numerical parameters
N_T  = 10000           # number of time steps
dt   = 0.0025          # time step size
Re   = Constant(250.0) # Reynolds number

# Create geometry and mesh
L  = 10.0 # channel length
H  =  4.0 # channel height
Xr =  2.0 # position of hole (assume centered)
Rr =  0.5 # radius of hole

channel  = Rectangle(Point(0, 0), Point(L, H))
p        = Point(Xr, H/2)
cylinder = Circle(p, Rr)
domain   = channel - cylinder
mesh     = generate_mesh(domain, 40)
R0 = H*0.5-Rr

for i in range(2):
    cell_markers = MeshFunction("bool", mesh, mesh.topology().dim())
    for c in cells(mesh):
        if c.midpoint().distance(p) < Rr + R0:
            cell_markers[c] = True
        else:
            cell_markers[c] = False
    mesh = refine(mesh, cell_markers)
    R0 = R0/2
#plt.figure(figsize=(30,30))
#plot(mesh)

# Define Taylor Hood element & function spaces
V = VectorElement("CG", mesh.ufl_cell(), 2)
Q = FiniteElement("CG", mesh.ufl_cell(), 1)
W = FunctionSpace(mesh, V*Q)

# Define boundaries
tol = 0.05
def B1(x, on_boundary):
    return on_boundary and (near(x[0],0,tol) or near(x[1],0,tol) or near(x[0],L,tol) or near(x[1],H,tol))

def B2(x, on_boundary):
    return on_boundary and near((x[0]-Xr)**2+(x[1]-H/2)**2,Rr*Rr,tol)

bc1 = DirichletBC( W.sub(0), Constant((1.0,0.0)), B1 )
bc2 = DirichletBC( W.sub(0), Constant((0.0,0.0)), B2 )
boundary = [bc1,bc2]

parameters["form_compiler"]["cpp_optimize"] = True
parameters["form_compiler"]["quadrature_degree"] = 3

qinit = Expression(("1.0","0.0","0.0"),degree=2)
old_q = interpolate(qinit,W)

file1 = File("output/vel.pvd")
file2 = File("output/pre.pvd")

q,dq = TrialFunction(W),TestFunction(W)
u,p  = split(q)
v,r  = split(dq)
old_u,old_p = split(old_q) 
Res = inner(u-old_u,v)*dx + dt*inner(v,grad(old_u)*old_u)*dx - dt*div(v)*p*dx
Res += (dt/Re)*inner(grad(u),grad(v))*dx + dt*r*div(u)*dx

q = old_q
a = lhs(Res)
L = rhs(Res)
A = assemble(a)

solver = LUSolver(A)
[bc.apply(A) for bc in boundary]
# bc2.apply(A)

t=0
for i in range(N_T):
    t = t+dt
    st = time.time()
    b = assemble(L)
    [bc.apply(b,q.vector()) for bc in boundary]
    solver.solve(q.vector(),b)
    et = time.time()
    print(i,(et-st))
    old_q = q

and the error I receive is:

TypeError                                 Traceback (most recent call last)
<ipython-input-30-d992e8ee9125> in <module>
     82 
     83 solver = LUSolver(A)
---> 84 [bc.apply(A) for bc in boundary]
     85 # bc2.apply(A)
     86 

<ipython-input-30-d992e8ee9125> in <listcomp>(.0)
     82 
     83 solver = LUSolver(A)
---> 84 [bc.apply(A) for bc in boundary]
     85 # bc2.apply(A)
     86 

/usr/local/lib/python3.6/dist-packages/dolfin/fem/dirichletbc.py in inside(self, x, on_boundary)
     45             return self.inside_function(x)
     46         else:
---> 47             return self.inside_function(x, on_boundary)
     48 
     49 

<ipython-input-30-d992e8ee9125> in B1(x, on_boundary)
     51 #     e4 = on_boundary and near(x[1],4.0,tol)
     52 #     e = e1 or e2 or e3 or e4
---> 53     return on_boundary and (near(x[0],0,tol) or near(x[1],0,tol) or near(x[0],L,tol) or near(x[1],H,tol))
     54 
     55 def B2(x, on_boundary):

TypeError: near(): incompatible function arguments. The following argument types are supported:
    1. (x0: float, x1: float, eps: float = 3e-16) -> bool

Is there any conflict between the boundary condition implementing for-loop and the function near or is it due to the formulation (Res)?

I would be obliged if you could provide any leads. Thanking you for your time and patience.

[conpierce8]

@atalasik

L = rhs(Res)

Before you apply the bcs, you have redefined L. Changing this to l = rhs(Res) and b = assemble(l) should fix the issue.

[atalasik]

I am thankful to you for the prompt response. It worked and the code ran.

However the boundary and initial conditions have been reversed. I am getting a velocity field around the obstacle where I have implemented the no slip boundary condition (as bc2 in the code)

This is what I get without the LUSolver (which further develops into a Von Karmann vortex street)

What else should be changed?

[conpierce8]

@atalasik

[bc.apply(b,q.vector()) for bc in boundary]

This form of bc.apply is for nonlinear problems (see DirichletBC.h). You should instead use

[bc.apply(b) for bc in boundary]

[atalasik]

I am extremely thankful to you for the prompt and appropriate response!

[atalasik]

When I couple the linearised Navier Stokes (NS) equations with linearised Nerst-Planck-Poisson (NPP) equations, the LUSolver runs but gives a wrong results, especially the velocity field.

The LUSolver yeilds correct results when the NPP and NS equations are solved independently, as well as the entire problem (NS+NPP) runs correctly without the LUSolver. Following is the code when I couple the NPP and NS equations and try to solve it with the LUSolver:

import matplotlib.pyplot as plt
from dolfin import *
from mshr import *
import logging
import time
# Importing mesh
# !dolfin-convert annulus.msh annulus.xml
mesh = Mesh("../meshes/annulus.xml")

# Defining function space for fluid dynamics
V = VectorElement("Lagrange", mesh.ufl_cell(), 2)
Q = FiniteElement("Lagrange", mesh.ufl_cell(), 1)
W = FunctionSpace(mesh,MixedElement([V,Q,Q,Q,Q]))

Re   = Constant(1.0) # Reynolds number 
t       = 0.0
m       = 0
debye   = 0.1
kappa   = 1.0

# Boundaries
tol=0.1
tolout = 100.0
domain_end = 1000.0
def walls(x,on_boundary):
    return on_boundary and near(x[0]**2 + x[1]**2,1.0,tol)

def farfield(x,on_boundary):
    return on_boundary and near(sqrt(x[0]**2 + x[1]**2),domain_end,tolout)

# Boundary conditions
beta = 1
bc_u_walls = DirichletBC(W.sub(0),Constant((0.0,0.0)),walls)
bc_u_farfield = DirichletBC(W.sub(0),Constant((0.0,0.0)),farfield)

bc_phi_walls = DirichletBC(W.sub(4),Constant(0.0),walls)
bc_phi_farfield = DirichletBC(W.sub(4),Expression("-beta*x[0]",beta=1,degree=2),farfield) 
bc_c_farfield = DirichletBC(W.sub(2),Constant(2.0),farfield)
bc_q_farfield = DirichletBC(W.sub(3),Constant(0.0),farfield)

bcs = [bc_c_farfield,bc_q_farfield,bc_phi_walls,bc_phi_farfield,bc_u_walls,bc_u_farfield]

parameters["form_compiler"]["cpp_optimize"] = True
parameters["form_compiler"]["quadrature_degree"] = 3
    
beta=1.0
var_init = Expression(("0","0","0","2","0","-beta*x[0]"),beta=beta,degree=2)
old_var = interpolate(var_init,W)

t=0
N_T=10000
dt = 0.01
res = 1
old_res = 100
eps = 1e-8
i = 0

file1 = File("coupled_output_LU/vel.pvd")
file2 = File("coupled_output_LU/pre.pvd")
file3 = File("coupled_output_LU/c.pvd")
file4 = File("coupled_output_LU/q.pvd")
file5 = File("coupled_output_LU/phi.pvd")
file6 = File("coupled_output_LU/Ephi.pvd")

var,dvar = TrialFunction(W),TestFunction(W)
u,p,c,q,phi = split(var)
du,dp,dc,dq,dphi = split(dvar)
old_u,old_p,old_c,old_q,old_phi = split(old_var)
    
Res = (c-old_c)*dc*dx + dt*inner(old_u,grad(c))*dc*dx + kappa*dt*(inner(grad(c),grad(dc))*dx +\
                                                                      old_q*inner(grad(phi),grad(dc))*dx)
Res += (q-old_q)*dq*dx + dt*inner(old_u,grad(q))*dq*dx + kappa*dt*(inner(grad(q),grad(dq))*dx +\
                                                                      old_c*inner(grad(phi),grad(dq))*dx)
Res += inner(grad(phi),grad(dphi))*dx - (q/(2*debye**2))*dphi*dx
    
Res += Re*(inner(u-old_u,du)*dx + inner(du,grad(old_u)*old_u)*dt*dx) - p*div(du)*dt*dx +\
            dt*inner(grad(u),grad(du))*dx - old_q*dt*inner(grad(phi),du)*dx + dt*dp*div(u)*dx  

var = old_var
a = lhs(Res)
ll = rhs(Res)

A = assemble(a)
[bc.apply(A) for bc in bcs]
solver = LUSolver(A)

for i in range(N_T):
    i = i+1
    t = t+dt
    st = time.time()
    b = assemble(ll)
    [bc.apply(b) for bc in bcs]
    solver.solve(var.vector(),b)
    et = time.time()
    print(i,t,(et-st))
    if (i % 20 ==0):
        u,p,c,q,phi  = var.split()
        u.rename("v","velocity")
        p.rename("p","pressure")
        c.rename("c","concentration")
        q.rename("q","charge")
        phi.rename("phi","potential")
        E = project(-grad(phi),VectorFunctionSpace(mesh, "CG", 1))
        E.rename("E","electric field")
        
        file1 << (u,t)
        file2 << (p,t)
        file3 << (c,t)
        file4 << (q,t)
        file5 << (phi,t)
        file6 << (E,t)
        
    old_var = var

I am unable to figure out what exactly I have implemented wrongly, or if I have missed out on anything. Any help would be greatly appreciated.

[conpierce8]

Can you provide the mathematical statement of the problem? Can you also show the results you are getting and the results you expect?

[atalasik]

The mathematical statements:

In classical formulation :

Variational formulation:

The velocity field I expect (the one I obtain without LUSolver) :

The velocity field I obtain with LUSolver:

[conpierce8]

I changed a few instances of phi to old_phi in your code to match the variational form you have supplied (hopefully my extensive formatting of the code doesn’t obscure the changes):

Res = (
    (c - old_c) * dc * dx
    + dt * inner(old_u, grad(c)) * dc * dx
    + kappa * dt * (
        inner(grad(c), grad(dc)) * dx +
        old_q * inner(grad(old_phi), grad(dc)) * dx  # changed `phi` to `old_phi`
    )
)
Res += (
    (q - old_q) * dq * dx
    + dt * inner(old_u, grad(q)) * dq * dx
    + kappa * dt * (
        inner(grad(q), grad(dq)) * dx
        + old_c * inner(grad(old_phi), grad(dq)) * dx  # changed `phi` to `old_phi`
    )
)
Res += inner(grad(phi), grad(dphi))*dx - (q / (2 * debye**2)) * dphi * dx
Res += (
    Re * (
        inner(u - old_u, du) * dx
        + dt * inner(du, grad(old_u) * old_u) * dx
    )
    - dt * p * div(du) * dx
    + dt * inner(grad(u), grad(du)) * dx
    - dt * old_q * inner(grad(old_phi), du) * dx  # changed `phi` to `old_phi`
    + dt * dp * div(u) * dx
)

With these changes, I obtain (at timestep i=500 and with a smaller domain having an outer radius of 100 instead of 1000):

[atalasik]

Thank you so much for pointing it out.

I would have completely missed out on this as, I was under the impression that old_c, old_q were enough to linearise the weak formulation, and even though I had documented the weak formulation with phi^n, I simply plugged in phi (instead of old_phi), thinking that my formulation has already been linearised.

Thanks again!!