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 < d.DOLFIN_EPS or x > 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.5, 2) + pow(x - 0.5, 2)) / 0.02)", degree = 1) g = d.Expression("sin(5*x)", 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:
- What does bc.apply() do exactly?
- Am I calculating the residual (i.e. res1) correctly? And why is res2 so large?
- Is there an easier way to evaluate the residual at any given point of the domain?