Regarding implementation of the LUSolver

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.

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

1 Like

Thank you so much for the prompt response!!!

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.

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

2 Likes

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?

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

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

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

1 Like

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.

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

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:

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):

3 Likes

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!!