Poisson equation: Impose full gradient as boundary condition via Lagrange multipliers

Question

I have a physical problem governed by the Poisson equation in two dimensions $$ -\nabla^2 u = f(x,y), \; in \; \Omega $$ I have measurements of the two gradient components $\partial{u}/\partial{x}$ and $\partial{u}/\partial{y}$ along some part of the boundary, $\Gamma_m$, so would like to impose $$ \frac{\partial{u}}{\partial{x_i}}_0 = g_m, \; on \; \Gamma_m $$ and propagate into the far field.

The tangential gradient component, $\frac{\partial{u}}{\partial{x}}_{(t,0)}$, I can just integrate and then enforce through a Dirichlet condition, such that $$ \int_{\Gamma_m}\frac{\partial{u}}{\partial{x}}_{(t,0)} \, ds = u_0 $$ In order to simultaneously impose the normal component, $\frac{\partial{u}}{\partial{x}}_{(n,0)}$, I gathered I would have to go via Lagrange multipliers.

So I think the variational form is then $$ \int_\Omega \nabla{u} \cdot\nabla{v}\, dx - \lambda \int_{\Gamma_m} ( \frac{\partial{u}}{\partial{x}}_{(n,0)}-g_m ) v \,ds = \int_\Omega f\, v\, dx $$ I spent a long time trying to piece it together from the information on related problems such as https://answers.launchpad.net/fenics/+question/212434 https://answers.launchpad.net/fenics/+question/216323

but still cannot see where I am going wrong. My solution attempt so far is:

from dolfin import *

# Create mesh and define function space
mesh = UnitSquareMesh(64, 64)
V = FunctionSpace(mesh, "Lagrange", 1)
R = FunctionSpace(mesh, "R", 0)
W = V * R

# Create mesh function over cell facets
boundary_parts = MeshFunction("uint", mesh, mesh.topology().dim()-1)

# Mark left boundary facets as subdomain 0
class LeftBoundary(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and x[0] < DOLFIN_EPS

Gamma_Left = LeftBoundary()
Gamma_Left.mark(boundary_parts, 0)

class FarField(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and ( (x[0] > 1.0-DOLFIN_EPS) \
               or (x[1]<DOLFIN_EPS) or (x[1]> 1.0-DOLFIN_EPS) )

Gamma_FF = FarField()
Gamma_FF.mark(boundary_parts, 1)

# Define boundary condition
u0 = Expression("sin(x[1]*pi)")
bcs = [DirichletBC(V, u0, Gamma_Left)]

# Define variational problem
(u, lmbd) = TrialFunctions(W)
(v, d) = TestFunctions(W)

f = Expression("10*exp(-(pow(x[0] - 0.5, 2) + pow(x[1] - 0.5, 2)) / 0.02)")
g = Constant(0.0)
h = Constant(-4.0)
n = FacetNormal(mesh)

F = inner(grad(u), grad(v))*dx + d*dot(grad(u),n)*ds(0) + lmbd*dot(grad(v),n)*ds(0)-\
   (f*v*dx + g*v*ds(1) + h*d*ds(0) + lmbd*h*ds(0))

a = lhs(F)
L = rhs(F)

# Compute solution
A = assemble(a, exterior_facet_domains=boundary_parts)
b = assemble(L, exterior_facet_domains=boundary_parts)
for bc in bcs: bc.apply(A, b)

w = Function(W)
solve(A, w.vector(), b, 'lu')
(u,lmbd) = w.split()

# Plot solution
plot(u, interactive=True)

which runs but gives a noisy result not at all resembling a solution to a Poisson equation. It seems to have something to do with the combined function spaces, but I cannot find the mistake.
I would appreciate any help or pointers in the right direction - many thanks already!
Cheers
Markus

Answer 1

First, a general point: you cannot prescribe arbitrary boundary conditions for a partial differential operator and expect that the partial differential equation (which always includes both operator and boundary conditions) is well-posed, i.e., admits a unique solution that depends continuously on the data -- all of which is a necessary condition for actually trying to compute something.

Depending on the operator, there are often quite a number of valid conditions you can impose (to get a taste, see the three-volume monograph by Lions and Magenes). However, what you are trying to do (specify the full gradient, which is equivalent to both Dirichlet and Neumann conditions on the same (part of the) boundary for a second-order elliptic PDE) is not among them - this is known as a lateral Cauchy problem, and is ill-posed: there is no guarantee that a given pair of boundary data admits a solution, and even if one exists, there is no stability with respect to small perturbations in the data. (In fact, this is the original ill-posed problem in the sense of Hadamard, and the classical example why you cannot ignore boundary conditions when discussing well-posedness. You can find an explicit example in his Lectures on Cauchy's problem in linear partial differential equations from the 1920s.)

Now to your specific problem (which might be a textbook example of an XY problem). If I understand correctly, you wish to solve an exterior problem for the Poisson equation, and you take a two-dimensional computional domain $(r,R)\times(a,b)$. You have full data (Dirichlet - using the tangential derivative trick - and Neumann) on (say) $x=r$. If $R$ is sufficiently large, you can justify a radiation condition (which prescribes the asymptotic behavior of your solution a t $x\to \infty$); see, e.g., the book Numerical methods for exterior problems by Long'an Ying. The question is what you know about the remaining boundary parts, $y=a$ and $y=b$.

1. If you can impose boundary conditions (Neumann, Robin, Dirichlet - which you would need to fix the constant in the integration of the tangential derivative, by the way), then it is enough to use either the normal components of your gradient as Neumann condition (if you can fix the constant mode) or integrate the tangential components as a Dirichlet condition. Since both conditions presumably correspond to the same function, you get the same solution either way.

2. If you don't know the behavior at $y=a$ and $y=b$, you really have a lateral Cauchy problem, and you cannot compute a solution using standard finite element methods. The standard way of dealing with this is the method of quasi-reversibility (introduced by Lattés and Lions in the 1960s; see their book), which consists in approximating the second-order problem by a fourth-order problem (where you can - and need to - prescribe two boundary conditions). In your case, this would amount to replacing $-\Delta u=f$ by $-\Delta u + \varepsilon \Delta^2 u=f$ for some small $\varepsilon>0$. (This can also be interpreted as minimizing the residual in a suitable norm and adding an $H^2$ regularization term; related to your comment on treating the problem as an inverse problem.) You can show that for compatible data (i.e., a pair of boundary conditions that actually correspond to a solution $u$ of the Poisson equation), the quasi-reversibility solutions $u_\varepsilon$ converge to $u$ as $\varepsilon\to0$.

   Since this is now a fourth-order problem with solutions in $H^2$, the best way of solving it numerically is using a mixed finite element formulation such as the one described in the paper by Dardé, Hannukainen and Hyvönen. (There are also some slides online.) It should not be too difficult to implement this approach using FEniCS.

Answer 2

You can't expect that solution to your altered problem would be a solution to Poisson problem because you need to change the problem somehow to make it well-posed.

One could guess that possible problem formulation is to minimize $$F(u, \lambda) = \int_\Omega \frac{1}{2}|\nabla u|^2 \;\mathrm{d}x - \int_\Omega f u \;\mathrm{d}x - \int_{\Gamma_\mathrm{N}} g u \;\mathrm{d}S + \int_{\Gamma_\mathrm{N}} \lambda (u-u_\mathrm{D})\;\mathrm{d}S$$ over $(u,\lambda) \in V\times L^2(\Gamma_\mathrm{N})$ where $V = \{v\in H^1;v|_{\Gamma_\mathrm{D}}=0\}$ is space with zero trace on Dirichlet boundary $\Gamma_\mathrm{D}$ and $u_\mathrm{D}$ is Dirichlet value you would like to enforce additionaly on Neumann boundary $\Gamma_\mathrm{N}$. Without the last term minimization of $F(u)$ yields $$ 0 = \mathrm{D}F(u)[v] = \int_\Omega \nabla u\cdot\nabla v\;\mathrm{d}x - \int_\Omega fv\;\mathrm{d}x - \int_{\Gamma_\mathrm{N}} gv\;\mathrm{d}S \quad \forall v\in V,$$ which is exactly Poisson problem with Neumann condition on ${\Gamma_\mathrm{N}}$ and zero (it is simply adjusted to non-zero) Dirichlet condition on ${\Gamma_\mathrm{D}}$. With additional Lagrange multiplier term necessary condition of minimality reads $$ 0 = \mathrm{D}F(u,\lambda)[v,\mu] = \int_\Omega \nabla u\cdot\nabla v\;\mathrm{d}x - \int_\Omega fv\;\mathrm{d}x - \int_{\Gamma_\mathrm{N}} gv\;\mathrm{d}S \quad + \int_{\Gamma_\mathrm{N}} \lambda v\;\mathrm{d}S + \int_{\Gamma_\mathrm{N}} \mu (u-u_\mathrm{D})\;\mathrm{d}S \quad \forall (v,\mu)\in V\times L^2(\Gamma_\mathrm{N}),$$ so we see that this modified problem admits Poisson problem $-\Delta u=f$ but with modified Neumann condition $\frac{\partial u}{\partial\mathbf{n}} = g-\lambda$ on $\Gamma_\mathrm{N}$. Simply said - requiring additional Dirichlet condition on $\Gamma_\mathrm{N}$ changed original Neumann condition required there.

If you solve this and $\lambda\ll|g|$ holds then maybe your over-specified conditions are consistent with Poisson problem. I say 'maybe' because of stability issue @ChristianClason mentioned in his answer.

Conversely we could try to enforce additional Neumann condition as a constraint (via Lagrange multiplier) on Dirichlet boundary $\Gamma_\mathrm{D}$. But that would not work either as test function $v\in V$ is zero on $\Gamma_\mathrm{D}$ which rules out constraint.

Conclusion is you that you can't expect that second order PDE will admit two independent boundary conditions.

---

As this question relates to FEniCS, note that it easy in FEniCS to specify $F(u,\lambda)$ and let UFL calculate its Gateaux derivative $\mathrm{D}F(u,\lambda)[v,\mu]$ by derivative() function.

Unfortunately $F(u,\lambda)$ can't be easily defined in FEniCS yet, as library does not support function $\lambda\in L^2(\Gamma_\mathrm{N})$ restricted to boundary. But another minimizations problems with Lagrange multipliers $\lambda\in L^2(\Omega)$ or $\lambda\in\mathbb{R}$ can be easily implemented this way.

Answer 3

Your approach cannot work, definitely because of the implementation and probably because of your formulation.

Imposing Dirichlet conditions in dolfin, eventually sets the corresponding DOFs of your testspace to zero.

This is an excerpt from the fenics-manual:

Chapter 3.3.9 (end): The application of a Dirichlet boundary condition to a linear system will identify all degrees of freedom that should be set to the given value and modify the linear system such that its solution respects the boundary condition. This is accomplished by zeroing and inserting 1 on the diagonal of the rows of the matrix corresponding to Dirichlet values, and inserting the Dirichlet value in the corresponding entry of the right-hand side vector [...]

Thus, your $v$ will be set to zero on $\Gamma_m$, what makes also your term with the multiplier disappear.

In summary, using the default routine in dolfin you cannot apply both Dirichlet and other conditions on the same boundary.

However, before you try to fix this in your implementation, go find the right test spaces for your mathematical formulation (as @Jan Blechta just has mentioned.)