# Solving pde with boundary integral constraint

I want to solve a PDE, but I'm struggling a lot getting it into Mathematica.

This is what I currently have (thanks to @user21). However, there are still a couple of issues which I explain below.

Needs["NDSolveFEM"]
{dpde, dbc, vd, sd, md} =
ProcessPDEEquations[{D[C1[x, y]*u[x, y], x] + D[C2[x, y]*u[x, y], y] +
D[C3[x, y]*u[x, y], x, x] + D[C4[x, y]*u[x, y], x, y] +
D[C5[x, y]*u[x, y], y, y] ==
NeumannValue[0, x <= 0] + NeumannValue[0, y <= 0] +
NeumannValue[0, x + y >= 1],
DirichletCondition[u[x, y] == 0, x == 1 && y == 0]},
u, {x, y} \[Element] region,
Method -> {"PDEDiscretization" -> {"FiniteElement",
"MeshOptions" -> {"MaxCellMeasure" -> 0.000125}}}]
s = dpde["StiffnessMatrix"];
DeployBoundaryConditions[{l, s}, dbc];
mesh = md["ElementMesh"];
dof = md["DegreesOfFreedom"];

mkIntegralConstraint[level_?NumberQ, mesh_, dof_Integer] :=
Module[{s1, s2, d1, diriRows, diriValMat},
s1 = SparseArray[{}, {dof, 1}];
s2 = SparseArray[{}, {dof, dof}];
d1 = SparseArray[{ConstantArray[1., {dof}]}];
diriRows = {1};
(*level the pure neumann to a specific value*)
diriValMat =
SparseArray[{{N[
level*dof/Total[mesh["MeshElementMeasure"], 2]]}}];
DiscretizedBoundaryConditionData[{s1, s2, d1, diriRows,
diriValMat, {dof, 0, Length[diriRows]}}, 1]]

dc = mkIntegralConstraint[1, mesh, dof]
DeployBoundaryConditions[{l, s}, dc, "ConstraintMethod" -> "Append"];

sol = LinearSolve[s, l];
if2 = ElementMeshInterpolation[{mesh}, Take[sol, dof]];
Plot3D[if2[x, y], {x, y} \[Element] region, PlotRange -> All]


However, I believe for these parameter values, the peak of the pdf should not be occurring at the edge of the domain like this.

• You've said the derivative is zero on the boundary and there is no flux, but for this particular equation isn't $u(x,y)=0$ a valid solution as your source condition u[x,1] == 1 is a single point on the top left corner of the triangle? There are no source terms in your equation, and no non-zero boundary conditions. Jun 24 '20 at 13:43
• @flinty Ah yeah, thanks for reminding me. I left out another condition from the problem. I also need the condition that u integrates to 1 over the domain! Jun 24 '20 at 14:00
• @flinty Also, I made a mistake above, the absorbing state is when $x=1$, not $y=1$. (So $u[1,y]$, which on my domain is only the point $u[1,0]$). Jun 24 '20 at 14:40

Here I am sketching out a solution path. My feeling is that the equations are not quite right, or something with the setup. More on this later.

This setup solves the equation given. Make sure it is correct as the examples shown in the question differ slightly.

uval = NDSolveValue[{D[C1[x, y]*u[x, y], x] +
D[C2[x, y]*u[x, y], y] + (1/(2 NN))*
D[C3[x, y]*u[x, y], x, x] - (1/NN)*
D[C4[x, y]*u[x, y], x, y] + (1/(2*NN))*
D[C5[x, y]*u[x, y], y, y] == 0,
DirichletCondition[u[x, y] == 1, x == 1 && y == 0]},
u, {x, y} \[Element] region]


This gives a message:

If the equation and setup are indeed correct then there are two ways to keep the convection dominance of this at bay. One is to add artificial diffusion which will change the PDE and the other is to mesh it better. I have gone with the second approach; and because I did not want to think about the PDE too much I have chosen the mesh-it-do-death approach.

uval = NDSolveValue[{D[C1[x, y]*u[x, y], x] +
D[C2[x, y]*u[x, y], y] + (1/(2 NN))*
D[C3[x, y]*u[x, y], x, x] - (1/NN)*
D[C4[x, y]*u[x, y], x, y] + (1/(2*NN))*
D[C5[x, y]*u[x, y], y, y] == 0,
DirichletCondition[u[x, y] == 1, x == 1 && y == 0]},
u, {x, y} \[Element] region,
Method -> {"PDEDiscretization" -> {"FiniteElement",
"MeshOptions" -> {"MaxCellMeasure" -> 0.0001}}}];

Plot3D[uval[x, y], {x, y} \[Element] region, PlotRange -> All]


When we check the integral constraint we see that:

NIntegrate[uval[x, y], {x, y} \[Element] region]
NIntegrate::slwcon: Numerical integration converging too slowly; suspect one of the following: singularity, value of the integration is 0, highly oscillatory integrand, or WorkingPrecision too small.
(* 0.000345862 *)


which is not the requested value. How could it be, we did not improve a constraint. Now, currently NDSolve can not handle integral constraints but many moons back I gave a presentation where I showed how to do it. I am going to reproduce the relevant part here.

Parse and discretize the PDE:

Needs["NDSolveFEM"]
{dpde, dbc, vd, sd, md} =
ProcessPDEEquations[{D[C1[x, y]*u[x, y], x] +
D[C2[x, y]*u[x, y], y] + (1/(2 NN))*
D[C3[x, y]*u[x, y], x, x] - (1/NN)*
D[C4[x, y]*u[x, y], x, y] + (1/(2*NN))*
D[C5[x, y]*u[x, y], y, y] == 0,
DirichletCondition[u[x, y] == 1, x == 1 && y == 0]},
u, {x, y} \[Element] region,
Method -> {"PDEDiscretization" -> {"FiniteElement",
"MeshOptions" -> {"MaxCellMeasure" -> 0.0000125}}}];


Extract the system matrices and deploy the given boundary conditions:

l = dpde["LoadVector"];
s = dpde["StiffnessMatrix"];
DeployBoundaryConditions[{l, s}, dbc];


Extract the mesh and the degrees of freedom:

mesh = md["ElementMesh"];
dof = md["DegreesOfFreedom"];


Next, we write a little helper function to create a boundary condition to satisfy boundary integral constrains:

mkIntegralConstraint[level_?NumberQ, mesh_, dof_Integer] :=
Module[{s1, s2, d1, diriRows, diriValMat},
s1 = SparseArray[{}, {dof, 1}];
s2 = SparseArray[{}, {dof, dof}];
d1 = SparseArray[{ConstantArray[1., {dof}]}];
diriRows = {1};
(* level the pure neumann to a specific value *)

diriValMat =
SparseArray[{{N[
level*dof/Total[mesh["MeshElementMeasure"], 2]]}}];
DiscretizedBoundaryConditionData[{s1, s2, d1, diriRows,
diriValMat, {dof, 0, Length[diriRows]}}, 1]
]


This function works by creating a weighted value for every node in the mesh in such a way that once these are summed up they satisfy the boundary integral value.

Make use of that:

dc = mkIntegralConstraint[1, mesh, dof]


Now, we deploy that boundary condition in such that the values are appended to the system matrices as Lagrangian multipliers.

DeployBoundaryConditions[{l, s}, dc,
"ConstraintMethod" -> "Append" ];


Solve the equations and post process the equations:

sol = LinearSolve[s, l];
if2 = ElementMeshInterpolation[{mesh}, Take[sol, dof]];
Plot3D[if2[x, y], {x, y} \[Element] region, PlotRange -> All]


Check the integral:

NIntegrate[if2[x, y], {x, y} \[Element] region,
Method -> "FiniteElement"]


(* 0.990778 *)


This is not perfect but that's do to the nature of the PDE used that is convection dominated and unstable. In that light, I think, the result is OK.

• Hi, Thanks you very much for your response. I'll have a proper read through your answer, and compare it to an approximation I have. Thanks! Edit: Actually, in the example you gave, you've used Dirichlet boundary conditions. The majority of my boundary conditions I want to be neumann, as in the second exmaple in my post (at least, that's what I've tried to implement). Would you be able to show me how to put in the neumann conditions I require? (Is this what you meant about the equation not being quite right?) Jun 26 '20 at 10:02
• I've tried running your code as is. It produces the first plot, as it did for you. However, when I get to the point of running the last section, I get the following error messages: 1) Take::seqs: Sequence specification (+n, -n, {+n}, {-n}, {m, n}, or {m, n, s}) expected at position 2 in Take[LinearSolve[dpde[StiffnessMatrix],dpde[LoadVector]],md[DegreesOfFreedom]]. 2) Plot3D::idomdim: {x,y}[Element]region does not have a valid dimension as a plotting domain. Jun 26 '20 at 11:08
• Also, thinking about it a bit more, I might want the absorbing boundary condition to be u[1,0]=0,rather than u[1,0]=1. Jun 26 '20 at 12:36
• Do you have any thoughts on the issues I'm currently having? Jul 3 '20 at 13:16
• @user112495, make sure that you re-add definitions for NN and region. Jul 16 '20 at 10:18

Infered from the NeumannValue settings the region is

RegionPlot[
x + y >= 1 \[Or] x <= 0 \[Or] y <= 0, {x, -.1, 1.1}, {y, -.1, 1.1}]


This can be enhanced using PlotPoints option.

As Region it looks:

Region[MeshRegion[{{0, 0}, {1, 0}, {0, 1}}, Line[{{1, 2, 3, 1}}]]]


Error message of the original question is than:

Better work with

polygon = Region[Polygon[{{0, 0}, {1, 0}, {0, 1}}]]


From the fundamentals of mathemtics follows that without restriction all C_i functions can be set equal one.

Then the questions looks like this:

{dpde, dbc, vd, sd, md} =
ProcessPDEEquations[{D[u[x, y], x] + D[u[x, y], y] +
D[u[x, y], x, x] + D[u[x, y], x, y] + D[u[x, y], y, y] ==
NeumannValue[0, x <= 0] + NeumannValue[0, y <= 0] +
NeumannValue[0, x + y >= 1],
DirichletCondition[u[x, y] == 0, x == 1 && y == 0]},
u, {x, y} \[Element] polygon,
Method -> {"PDEDiscretization" -> {"FiniteElement",
"MeshOptions" -> {"MaxCellMeasure" -> 0.000125}}}]
s = dpde["StiffnessMatrix"];
DeployBoundaryConditions[{l, s}, dbc];
mesh = md["ElementMesh"];
dof = md["DegreesOfFreedom"];

mkIntegralConstraint[level_?NumberQ, mesh_, dof_Integer] :=
Module[{s1, s2, d1, diriRows, diriValMat},
s1 = SparseArray[{}, {dof, 1}];
s2 = SparseArray[{}, {dof, dof}];
d1 = SparseArray[{ConstantArray[1., {dof}]}];
diriRows = {1};
(*level the pure neumann to a specific value*)
diriValMat =
SparseArray[{{N[
level*dof/Total[mesh["MeshElementMeasure"], 2]]}}];
DiscretizedBoundaryConditionData[{s1, s2, d1, diriRows,
diriValMat, {dof, 0, Length[diriRows]}}, 1]]

dc = mkIntegralConstraint[1, mesh, dof]
DeployBoundaryConditions[{l, s}, dc, "ConstraintMethod" -> "Append"];

sol = LinearSolve[s, l];
if2 = ElementMeshInterpolation[{mesh}, Take[sol, dof]];
Plot3D[if2[x, y], {x, y} \[Element] region, PlotRange -> All]


This confirms the expection of the question in all manners:

MaxValue[if2[x, y], {x, y} \[Element] polygon]
if2[0, 1]


2.14897

2.14897

The minimum is in a corner too:

if2[1, 0]

-2.74605*10^-7


These confirm that the given boundary conditions hold only for the Dirichlet condition!

From the documentation of DirichletCondition

Typically, at least one Dirichlet-type boundary condition needs to be specified to make the differential equation uniquely solvable. Dirichlet conditions are also called essential boundary conditions.

Something different from the coordinate function is definitely wrong after this changes.