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('''
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
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. \]
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 
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  
Please login to add an answer/comment or follow this question.

Similar posts:
Search »