The correct way to calculate and plot the residual?

variational formulation
1

Hi all, I am using dolfin 2019.1.0 installed via pip. In this particular example I’m solving the Poisson equation (as described in the classic example https://fenicsproject.org/olddocs/dolfin/1.5.0/python/demo/documented/poisson/python/documentation.html). I am trying to figure out the correct way to evaluate the residual (specifically the difference between the left hand side and the ride hand side of the equation that I am trying to solve). What I find is that if I don’t enforce the boundary conditions correctly, I get extremely high values for the residual.

Here is the minimal code:

import dolfin as d
import matplotlib.pyplot as plt

# Creating the mesh and function space:
mesh = d.UnitSquareMesh(32,32)
V = d.FunctionSpace(mesh, "Lagrange", 1)

# Defining the boundary conditions:
def boundary(x):
    return x[0] < d.DOLFIN_EPS or x[0] > 1.0 - d.DOLFIN_EPS

u0 = d.Constant(0.0)
bc = d.DirichletBC(V, u0, boundary)

# Defining the Variational problem:
u_sol = d.TrialFunction(V)
v = d.TestFunction(V)

f = d.Expression("10*exp(-(pow(x[0] - 0.5, 2) + pow(x[1] - 0.5, 2)) / 0.02)", degree = 1)
g = d.Expression("sin(5*x[0])", degree = 1)
a = d.inner(d.grad(u_sol), d.grad(v))*d.dx
L = f*v*d.dx + g*v*d.ds
u_sol = d.Function(V)

# Compute the solution:
d.solve(a == L, u_sol, bc)

# Evaluate the residual (i.e. the difference between the LHS and the RHS):
RES = d.inner(d.grad(u_sol), d.grad(v))*d.dx - f*v*d.dx - g*v*d.ds
res1 = d.assemble(RES)
res2 = d.assemble(RES)

# Enforcing the boundary conditions only for res1 and not for res2:
bc.apply(res1, u_sol.vector())

# The max and min values of the residual for the two cases differ greatly:
print("Residual max with BC applied:", res1.max())
print("Residual min with BC applied:", res1.min())
print("Residual max without BC applied:", res2.max())
print("Residual min without BC applied:", res2.min())

# Plotting the value of the residual of the entire domain:
resF1 = d.Function(V,res1)

plt.figure()
col = d.plot(resF1)
plt.colorbar(col)
plt.title('Residual with BC applied')
plt.show()

# The value of the residual with no boundary conditions applied: 
resF2 = d.Function(V,res2)

plt.figure()
col = d.plot(resF2)
plt.colorbar(col)
plt.title('Residual without BC applied')
plt.show()

So my questions are as follows:

  1. What does bc.apply() do exactly?
  2. Am I calculating the residual (i.e. res1) correctly? And why is res2 so large?
  3. Is there an easier way to evaluate the residual at any given point of the domain?

Thanks!

2

I would consider having a look at this very similar post: Residual for Stokes equation not equal to zero - #6 by kamensky