Computing norm of solution obtained by LUSolver

General
1

Hello everyone,
I am trying implement a condition (in Fenics-Dolfin 2019.1.0 version) to see if and when the solution reaches a steady state. In order to achieve this I am computing the norm of the difference in the solutions (which are obtained by LUSolver) of the current and previous time steps. Following is my code :

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 = TrialFunction(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)
L = rhs(Res)

def calculateNorm(old_u,u):
    eps_u = norm((u.vector() - old_u.vector()),'l2')
    return eps_u

res = 100
eps = 1e-5

u = old_u
# main loop and plot solution each step
i = 0
while res > eps:
  t = t + tau
  i = i + 1
  b = assemble(L)
  solver.solve(u.vector(),b)
  res = calculateNorm(old_u,u)
  print(res)
  old_u = u

In the first iteration itself I obtain the value of res to be zero and the iterations stop. However without the LUSolver the value of res gradually reduces every iteration until it reaches below eps . What exactly am I missing out on?

Actually, this post is an extension of this post , wherein I have implemented the LUSolver. However since the focus was more on computing norms and understanding the nature of solution obtained when TrialFunction and solvers are used, I therefore started a new post.

2

You need to create an additional Function object. The statement

u = old_u

assigns the Python name u to point to the same underlying object as the Python name old_u. Therefore, in your code, there is only one Function object, and of course the residual between a function and itself will always be zero. What you need to do is create a new function and copy the values from the Function object old_u to the Function object u. The following modification (untested) should work for you:

# Create a new function `u` and assign initial values from `old_u`
u = Function(U)
old_u.assign(u)

# main loop and plot solution each step
i = 0
while res > eps:
  t = t + tau
  i = i + 1
  b = assemble(L)
  solver.solve(u.vector(),b)
  res = calculateNorm(old_u,u)
  print(res)

  # Update the `old_u` solution by copying the solution
  # vector from `u` to `old_u`
  u.assign(old_u)
3 Likes
3

I am extremely thankful to you for promptly and appropriately responding.

However, I observed that in Fenics 2019.1.0, the functionality of assign is the other way around. The vector which gets copied is written in the parentheses of the assign function and the vector on which the copy is imposed, is written before the assign function by a dot separator. For example:

<solution_that_gets_written>.assign(<solution_that_writes>)

I would also like to ask that if the calculateNorm function, receives a vector having various fields (for example, temperature, pressure, entropy and velocity), then why does the split function needs to contain a boolean argument True? For example:

v,p,T,H = u.split(True)

I am asking this because, without True as an argument the norm doesn’t get calculated properly and yields fluctuating values. I had not used the True argument in the case without the LUSolver, and norm was observed to steadily decrease with each iteration.

4

You are correct. My mistake.

Can you provide a code example that demonstrates this behavior? The code you have provided in this thread does not use split, and the code you have provided in your previous post does not contain an implementation of calculateNorm. Could you also provide your implementation of the “case without the LUSolver”? It’s difficult for me to determine why you are not getting the expected results without having a code that can produce the expected results.

1 Like
5

Below I have pasted the implementations of the cases with and without the LUSolver. Both the cases work fine, and I get the desired results. All I am trying to understand is the reason behind the different syntax in the function calculateNorm for the cases with and without the use LUSolver.

The case where I have not used the LUSolver: Here I have to project back the solution on the function spaces, in the function calculateNorm (otherwise I get an error " ‘Sum’ object has no attribute ‘vector’ " even if I give True as an argument to split function )

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

mesh = Mesh("../meshes/annulus.xml")
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(1e-2) # Reynolds number 
m       = 0
debye   = 0.1
kappa   = 1.0

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)

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_u_walls,bc_u_farfield,bc_c_farfield,bc_q_farfield,bc_phi_walls,bc_phi_farfield]

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

def evolve(old_var,dt):
    var,dvar = Function(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 - (old_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.assign(old_var)
    solve(Res==0,var,bcs)
    return var

beta=1.0
var_init = Expression(("0","0","0","2","0","-beta*x[0]"),beta=beta,degree=2)
old_var = interpolate(var_init,W)

def calculateNorm(old_var,var):
    old_u,old_p,old_c,old_q,old_phi = old_var.split()
    u,p,c,q,phi = var.split()
    U = project((u-old_u),VectorFunctionSpace(mesh, "CG", 2))
    C = project((c-old_c),FunctionSpace(mesh, "CG", 1))
    charge = project((q-old_q),FunctionSpace(mesh, "CG", 1))
    eps_u = norm(U.vector())
    eps_c = norm(C.vector())
    eps_q = norm(charge.vector())
    return max(set([eps_u,eps_c,eps_q]))

t = 0
dt = 0.01
res = 1
eps = 1e-8
i = 0
while res > eps:
    i = i+1
    t = t+dt
    st = time.time()
    var = evolve(old_var,dt)
    et = time.time()
    res = calculateNorm(old_var,var)
    print(res)
    old_var.assign(var)

The case where I have implemented the LUSolver. Unless I don’t put True as an argument to the split function (in the calculateNorm function), norm isn’t computed correctly.

import matplotlib.pyplot as plt
from dolfin import *
from mshr import *
import logging
import time
# Importing mesh
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 
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
dt = 0.01
res = 1
eps = 1e-8
i = 0

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(old_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(old_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(old_phi),du)*dx + dt*dp*div(u)*dx  

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

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

varTilde = Function(W)
varTilde.assign(old_var)

def calculateNorm(old_var,var):
    (old_u,old_p,old_c,old_q,old_phi) = old_var.split(True) 
    (u,p,c,q,phi) = var.split(True)

    eps_u = norm((u.vector()-old_u.vector()),'l2')
    eps_c = norm((c.vector()-old_c.vector()),'l2')
    eps_q = norm((q.vector()-old_q.vector()),'l2')  
        
    return max(set([eps_u,eps_c,eps_q]))

while res > eps:
    i = i+1
    t = t+dt
    st = time.time()
    b = assemble(ll)
    [bc.apply(b) for bc in bcs]
    solver.solve(varTilde.vector(),b)
    et = time.time()
    res = calculateNorm(old_var,varTilde)
    print(res)
    old_var.assign(varTilde)

Pardon me if I have missed out on anything, and I am thankful to you for your time and patience.

6

Regarding the need to use split(deepcopy=True): when you call split without creating a “deep copy”, the .vector() method of the sub-functions will be a shallow copy, i.e. it is just a pointer to the vector() of its parent function. The vector of the parent function will include all dofs of the mixed space, so when you do (u,p,c,q,phi) = var.split(False) and then take norm(u.vector()), you are actually doing the same as norm(var.vector()). Using deepcopy=True causes each of the sub-functions to have its own PETScVector, so the norm is correctly computed only over the dofs of the subspace. Check out the following example which illustrates this:

from dolfin import *

mesh = UnitSquareMesh(4, 4)
elem_s = FiniteElement("P", mesh.ufl_cell(), 1)
elem_v = VectorElement("P", mesh.ufl_cell(), 1)
elem_m = MixedElement([elem_s, elem_v])

V = FunctionSpace(mesh, elem_m)
m = Function(V)

print("m.vector()")
print(m.vector().id(), m.vector().size())

print("\ndeepcopy=False")
s0, v0 = m.split()
print(s0.vector().id(), s0.vector().size())
print(v0.vector().id(), v0.vector().size())

print("\ndeepcopy=True")
s1, v1 = m.split(deepcopy=True)
print(s1.vector().id(), s1.vector().size())
print(v1.vector().id(), v1.vector().size())

producing

m.vector()
7 75

deepcopy=False
7 75
7 75

deepcopy=True
20 25
25 50

You can see that with deepcopy=False, the vectors returned by s.vector() and v.vector() have the same id, i.e. they are the same object.

Regarding the need to use project: you are calling norm with two different syntaxes. For the non-LU case, you are calling it with the difference of two Function objects, while in the LU case, you are calling it with the difference of two PETScVector objects. When you subtract one function from another, this creates a UFL symbolic algebra object, rather than numerically calculating the difference of the two solutions. norm() does not accept UFL objects as input. If you copy the definition of calculateNorm from the LU case to the non-LU case, it should work (and would be preferred, since you would no longer have to calculate the solution of three extra linear problems at every iteration).

1 Like
7

I am extremely thankful to you for explaining it to me with such clarity and detail. I have updated the code in the case without the LUSolver. It worked!

I wish to explore and know more about Fenics hence, is there any channel or documentation you would recommend?

Thanks once again!!

8

My top resources are:

  1. This forum - there are lots of great answers to search from when you run into an issue.

  2. The legacy FEniCS source code (Bitbucket). A brief introduction to the source:

    1. The fast C++ “backend” code is located in the “dolfin” directory (Bitbucket)
    2. The Python “frontend” code is located in the “python/dolfin” directory (Bitbucket)
    3. The wrapper layer that allows C++ code to be called from Python is found in “python/src” (Bitbucket). Any imports in the Python layer that include .cpp are defined in this layer.

    I often use the source to find the correct call signature for functions.

  3. The Python dir() function - to find out what methods are accessible from the objects in my code.

  4. The legacy FEniCS demos: https://fenicsproject.org/olddocs/dolfin/2019.1.0/python/demos.html (although your knowledge of FEniCS appears to exceed the point where the demos would be useful)

  5. The UFL manual: User manual — Unified Form Language (UFL) 2021.1.0 documentation

I will also include here an advertisement for DOLFINx with a gentle encouragement to make the switch :slight_smile:

2 Likes
9

Thank you very much for sharing such a comprehensive set of documentation!

Is DOLFINx/ FEniCSx an update of FEniCS? If not then what are the advantages?

10

See for instance this post The new DOLFINx solver is now recommended over DOLFIN
and
https://fenicsproject.org/roadmap/