### Neumann Boundary Conditions not at Correct Value

96
views
2
3 months ago by

I am attempting to apply Neumann boundary conditions to a PDE but I am having some difficulties. See the code below,

from fenics import *
import numpy as np

def atanh(x):
return ln((1 + x) / (1 - x)) / 2

def run(sys):
ds = globals()['ds']

n = 200

mesh = UnitIntervalMesh(n)
V = FunctionSpace(mesh, 'CG', 1)
v = TestFunction(V)
rho = Function(V)

class LeftEdge(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and near(x[0], 0)

class RightEdge(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and near(x[0], 1)

boundaryLeft = LeftEdge()
boundaryRight = RightEdge()

boundaries = FacetFunction("size_t", mesh)
boundaries.set_all(0)
boundaryLeft.mark(boundaries,1)
boundaryRight.mark(boundaries,2)

ds = ds(subdomain_data = boundaries)

G = Constant(1)

F = eval(sys)

solve(F == 0, rho)

rho = rho.vector().get_local()

print('Resulting BC', np.gradient(rho, 1./n)[[0,-1]], 'Applied BC', G.values()[0])

run('''
- v * G * ds(1) \
+ v * G * ds(2)
''')

run('''
- v * G * ds(1) \
+ v * G * ds(2)
''')

run('''
- 1000 * v * G * ds(1) \
+ 1000 * v * G * ds(2)
''')

run('''
(1000 * grad(rho)[0] * grad(v)[0] + 20000000 * v * atanh(rho)) * dx \
- 1000 * v * G * ds(1) \
+ 1000 * v * G * ds(2)
''')

run('''
(1000 * grad(rho)[0] * grad(v)[0] + 20000000 * v * atanh(rho)) * dx \
- 1.4 * 1000 * v * G * ds(1) \
+ 1.4 * 1000 * v * G * ds(2)
''')

This code outputs the following,

*** -------------------------------------------------------------------------
*** Warning: FacetFunction<T>(mesh, value) has been deprecated in FEniCS version 2017.2.0.
*** It will (likely) be removed in the next FEniCS release.
*** Use MeshFunction<T>(mesh, mesh->topology().dim() - 1, value)
*** -------------------------------------------------------------------------

No Jacobian form specified for nonlinear variational problem.
Differentiating residual form F to obtain Jacobian J = F'.
Solving nonlinear variational problem.
Newton iteration 0: r (abs) = 1.414e+00 (tol = 1.000e-10) r (rel) = 1.000e+00 (tol = 1.000e-09)
Newton iteration 1: r (abs) = 4.121e-13 (tol = 1.000e-10) r (rel) = 2.914e-13 (tol = 1.000e-09)
Newton solver finished in 1 iterations and 1 linear solver iterations.
Resulting BC [ 1.  1.] Applied BC 1.0
*** -------------------------------------------------------------------------
*** Warning: FacetFunction<T>(mesh, value) has been deprecated in FEniCS version 2017.2.0.
*** It will (likely) be removed in the next FEniCS release.
*** Use MeshFunction<T>(mesh, mesh->topology().dim() - 1, value)
*** -------------------------------------------------------------------------

No Jacobian form specified for nonlinear variational problem.
Differentiating residual form F to obtain Jacobian J = F'.
Solving nonlinear variational problem.
Newton iteration 0: r (abs) = 1.414e+00 (tol = 1.000e-10) r (rel) = 1.000e+00 (tol = 1.000e-09)
Newton iteration 1: r (abs) = 1.173e-13 (tol = 1.000e-10) r (rel) = 8.296e-14 (tol = 1.000e-09)
Newton solver finished in 1 iterations and 1 linear solver iterations.
Resulting BC [ 0.001  0.001] Applied BC 1.0
*** -------------------------------------------------------------------------
*** Warning: FacetFunction<T>(mesh, value) has been deprecated in FEniCS version 2017.2.0.
*** It will (likely) be removed in the next FEniCS release.
*** Use MeshFunction<T>(mesh, mesh->topology().dim() - 1, value)
*** -------------------------------------------------------------------------

No Jacobian form specified for nonlinear variational problem.
Differentiating residual form F to obtain Jacobian J = F'.
Solving nonlinear variational problem.
Newton iteration 0: r (abs) = 1.414e+03 (tol = 1.000e-10) r (rel) = 1.000e+00 (tol = 1.000e-09)
Newton iteration 1: r (abs) = 8.799e-11 (tol = 1.000e-10) r (rel) = 6.222e-14 (tol = 1.000e-09)
Newton solver finished in 1 iterations and 1 linear solver iterations.
Resulting BC [ 1.  1.] Applied BC 1.0
*** -------------------------------------------------------------------------
*** Warning: FacetFunction<T>(mesh, value) has been deprecated in FEniCS version 2017.2.0.
*** It will (likely) be removed in the next FEniCS release.
*** Use MeshFunction<T>(mesh, mesh->topology().dim() - 1, value)
*** -------------------------------------------------------------------------

No Jacobian form specified for nonlinear variational problem.
Differentiating residual form F to obtain Jacobian J = F'.
Solving nonlinear variational problem.
Newton iteration 0: r (abs) = 1.414e+03 (tol = 1.000e-10) r (rel) = 1.000e+00 (tol = 1.000e-09)
Newton iteration 1: r (abs) = 5.617e-03 (tol = 1.000e-10) r (rel) = 3.971e-06 (tol = 1.000e-09)
Newton iteration 2: r (abs) = 2.430e-11 (tol = 1.000e-10) r (rel) = 1.718e-14 (tol = 1.000e-09)
Newton solver finished in 2 iterations and 2 linear solver iterations.
Resulting BC [ 0.71300512  0.71300512] Applied BC 1.0
*** -------------------------------------------------------------------------
*** Warning: FacetFunction<T>(mesh, value) has been deprecated in FEniCS version 2017.2.0.
*** It will (likely) be removed in the next FEniCS release.
*** Use MeshFunction<T>(mesh, mesh->topology().dim() - 1, value)
*** -------------------------------------------------------------------------

No Jacobian form specified for nonlinear variational problem.
Differentiating residual form F to obtain Jacobian J = F'.
Solving nonlinear variational problem.
Newton iteration 0: r (abs) = 1.980e+03 (tol = 1.000e-10) r (rel) = 1.000e+00 (tol = 1.000e-09)
Newton iteration 1: r (abs) = 1.541e-02 (tol = 1.000e-10) r (rel) = 7.784e-06 (tol = 1.000e-09)
Newton iteration 2: r (abs) = 2.037e-11 (tol = 1.000e-10) r (rel) = 1.029e-14 (tol = 1.000e-09)
Newton solver finished in 2 iterations and 2 linear solver iterations.
Resulting BC [ 0.99820451  0.99820451] Applied BC 1.0

In the first example the result are what I expect, however in the second example I would still expect to have the correct boundary conditions applied as one and one. I was able to figure out how to get this working again in the third example, but is there a better way to do this? Then in the fourth example the boundary conditions are again not correct. I can add a constant to fix them as in the last example, but I would like to have this calculated automatically. What can I change to make sure the boundary conditions are applied correctly?

Community: FEniCS Project

4
3 months ago by
In the first three examples, you get the correct answer. If I'm not mistaken, you want to solve
$\frac{\mathrm{d}^2 \rho}{\mathrm{d} x^2} = 0, \qquad \text{with} \qquad \left. \frac{\mathrm{d} \rho}{\mathrm{d} x} \right|_{x = \{0, 1\}} = G.$
$0 = \int_0^1 \frac{\mathrm{d}^2 \rho}{\mathrm{d} x^2} v \,\mathrm{d} x = \left[ \frac{\mathrm{d} \rho}{\mathrm{d} x} v \right]_0^1 - \int_0^1 \frac{\mathrm{d} \rho}{\mathrm{d} x} \frac{\mathrm{d} v}{\mathrm{d} x} \,\mathrm{d}x,$
or
$\int_0^1 \frac{\mathrm{d} \rho}{\mathrm{d} x} \frac{\mathrm{d} v}{\mathrm{d} x} \,\mathrm{d}x = G \left( v(1) - v(0) \right).$
That's exactly what you have in case 1 and 3 (both sides being multiplied by 1000 in case 3). In the second case, you multiplied by 1000 only the left-hand side, i.e., you lowered the constant G by 1000. That's why the result is 0.001.

In the last two examples, you've got a different equation and the multiplication by a constant is more of a trial & error approach.
The OP is applying Neumann boundary conditions. Does your suggestion of
DirichletBC​
still apply?
written 3 months ago by ricopicone
1
My bad, I missed that there's a gradient of rho. I've edited the answer accordingly.
written 3 months ago by Adam Janecka
The nonlinear case is more problematic, there's a problem with boundary layer. In other words, $\rho = 0$ 'is' the solution in majority of the domain, except for the boundary where it needs to rapidly satisfy the boundary conditions -- and it fails to do so. The remedy is either using a finer mesh or lowering the magnitude of the nonlinearity.
written 3 months ago by Adam Janecka