### Neumann Boundary Conditions not at Correct Value

96

views

2

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('''
grad(rho)[0] * grad(v)[0] * dx \
- v * G * ds(1) \
+ v * G * ds(2)
''')
run('''
1000 * grad(rho)[0] * grad(v)[0] * dx \
- v * G * ds(1) \
+ v * G * ds(2)
''')
run('''
1000 * grad(rho)[0] * grad(v)[0] * dx \
- 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

### 1 Answer

4

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. \]

Then the weak formulation reads

\[ 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.

\[ \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. \]

Then the weak formulation reads

\[ 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

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

Please login to add an answer/comment or follow this question.