# Nonlinear diffusion equation: Periodic boundary conditions not satisfied (BUG?)

it's me again.

Strange stuff to report today. I suspect I've found a bug! Here is the nonlinear diffusion equation direct from the Mathematica documentation for FEM.

c = 1/Sqrt[(1 + Grad[u[x, y], {x, y}].Grad[u[x, y], {x, y}])];
Cu = {{{{c, 0}, {0, c}}}};
eqn = {Inactive[Div][
Cu[[1, 1]].Inactive[Grad][u[x, y], {x, y}], {x, y}] == 0};


And a simple mesh to solve it with:

Needs["NDSolveFEM"];
mesh = ToElementMesh[FullRegion[2], {{-1, 1}, {-1, 1}}];
Show[mesh["Wireframe"], Frame -> True]


Note the exact solution (the diffusion tensor is constant for this case):

uA[x_, y_] = y;


Our boundary conditions will "target" this solution, using mixed Dirichlet and periodic boundary conditions (can do it with pure Dirichlet but that misses the point of this post):

bcs = {DirichletCondition[u[x, y] == uA[x, y], -1 < x < 1],
PeriodicBoundaryCondition[u[x, y], x == 1, # - {2, 0} &]};


We provide the solver an initial guess (seed) that agrees on the boundary with the exact solution, but deviates inside. (This is not important, but we want the solver to work a bit for the solution.)

uSeed[x_, y_] = (1 - 0.3 (1 - x^2) (1 - y^2)) uA[x, y];


Now we solve this problem with NDSolveValue:

{ufA} = NDSolveValue[Join[eqn, bcs], {u}, Element[{x, y}, mesh],
InitialSeeding -> {u[x, y] == uSeed[x, y]}];
Plot3D[ufA[x, y], Element[{x, y}, mesh]]


Oh dear! This doesn't look good! We wanted it to look this of course:

Plot3D[uA[x, y], Element[{x, y}, mesh]]


But, to the point now. It does not even satisfy the periodic boundary condition, on the target boundary x == 1! That is the problem, simply stated. What is going on here?

I will scratch a little deeper to gather some clues, using FEM programming. Just mostly copying code from the documentation here:

iSeeding = {uSeed[x, y]};
vd = NDSolveVariableData[{"DependentVariables", "Space"} -> {{u}, {x, y}}];
sd = NDSolveSolutionData[{"DependentVariables",
"Space"} -> {iSeeding, ToNumericalRegion[mesh]}];
coefficients = {"DiffusionCoefficients" -> Cu};
initCoeffs = InitializePDECoefficients[vd, sd, coefficients];
initBCs = InitializeBoundaryConditions[vd, sd, bcs] ;
methodData =
InitializePDEMethodData[vd, sd, Method -> {"FiniteElement"}];
linearizedPDECoeffs = LinearizePDECoefficients[initCoeffs, vd, sd];
SplitPDECoefficients[linearizedPDECoeffs, vd, sd];
sdU = EvaluateInitialSeeding[methodData, vd, sd];
linear = DiscretizePDE[linearizedPDECoeffs, methodData,
sdU]; {linearLoad, linearStiffness, linearDamping, linearMass} =
linear["SystemMatrices"];
linearBCs = DiscretizeBoundaryConditions[initBCs, methodData, sdU];
seed = NDSolveSolutionDataComponent[sdU, "DependentVariables"];


All standard stuff. Now we come to something interesting. We call DeployDirichletConditions on the seed data that we just created. The way we set up the boundary conditions, this should do nothing because the seed already satisfies the boundary conditions. It requires no modification. However it is indeed modified quite signifiantly:

{DeployDirichletConditions[seed, linearBCs],
Norm@(seed - seedOLD)/Norm[seedOLD]}


{Null, 0.175549}

Now let's visualize the modified seed data:

uSeedf = ElementMeshInterpolation[mesh, seed];
Plot3D[uSeedf[x, y], {x, -1, 1}, {y, -1, 1}, AxesLabel -> Automatic]


This seems like an important clue. The seed has been modified so that the values at x==1 (the target of PeriodicBoundaryCondition) are now all zero (rather than periodic as they should be)! It seems something has gone wrong. To finish the solution, we need two functions femJacobian and femRHS, copied from the documentation, and I'll give their definitions at the end of this post for reference; you'll have to execute that them first. Then we run FindRoot to get the solution:

root = U /.
FindRoot[femRHS[U], {U, seedOLD}, Jacobian -> femJacobian[U],
Method -> {"AffineCovariantNewton"}];
NDSolveSetSolutionDataComponent[sdU, "DependentVariables", root];
{uf} = ProcessPDESolutions[methodData, sdU];
Plot3D[uf[x, y], Element[{x, y}, mesh]]


The solution agrees with the one that came from NDSolveValue, as expected -- i.e. it's wrong. However, let's try FindRoot again, but this time circumvent the effect of DeployDirichletConditions by using the unmodified seed, seedOLD:

root = U /.
FindRoot[femRHS[U], {U, seedOLD}, Jacobian -> femJacobian[U],
Method -> {"AffineCovariantNewton"}];
NDSolveSetSolutionDataComponent[sdU, "DependentVariables", root];
{uf} = ProcessPDESolutions[methodData, sdU];
Plot3D[uf[x, y], Element[{x, y}, mesh]]


This looks great! Time to celebrate? Sorry, not so fast. There are more problems. The solver seems to do OK if the initial seed agrees with the final solution on the target boundary (x==1). This is rather artificial. For many problems, we won't know what the solution will be on the boundary. For example, if we try the following seed function, things go very sour again:

uSeed[x_, y_] = (1 - 0.8 (1 - y^2)) uA[x, y];


This seed function is similar to the original, but it deviates from the exact solution when Abs[y] < 1, i.e. on the boundaries x==-1 and x==1. If we solve again (we have to go back to the definition of iSeed above), the standard way, with DeployDirichletConditions, we obtain the solution that violates periodicity (agres with the original output of NDSolveValue, uA). If we try our new "trick" and skip DeployDirichletConditions, things get interesting again:

Plot3D[uf[x, y], Element[{x, y}, mesh]]


If we look at the solution near the x == 1 boundary it seems there is a remnant of the seed function. Indeed if we subtract the seed we find

Plot3D[uf[x, y] - uSeed[x, y], Element[{x, y}, mesh], PlotRange -> All]

Instead of enforcing periodic BC, the solver is effectively forcing the solution to be equal to the seed function uSeed at the target boundary x == 1. This is very curious behavior! I really hope someone has an idea about this. @user21?

Below are the functions you need copied from the Mathematica documentation. Thanks for reading.

femRHS[u_?VectorQ] :=
NDSolveSetSolutionDataComponent[sdU, "DependentVariables", u];
nonlinear =
nonlinear = Null;
(*subtract the linear Robin boundary value*)
nonlinearBCs =
DiscretizeBoundaryConditions[initBCs, methodData, sdU,
"Nonlinear"];

femJacobian[u_?VectorQ] :=
Block[{stiffness, nonlinear, nonlinearStiffness, nonlinearBCs},
NDSolveSetSolutionDataComponent[sdU, "DependentVariables", u];
nonlinear =
DiscretizePDE[linStiffnessPDEC, methodData, sdU, "Nonlinear"];
nonlinearStiffness = nonlinear["StiffnessMatrix"];
nonlinear = Null;
stiffness = linearStiffness + nonlinearStiffness;
nonlinearStiffness = Null;
nonlinearBCs =
DiscretizeBoundaryConditions[initBCs, methodData, sdU,
"Nonlinear"];
DeployPartialBoundaryConditions[{Null, stiffness}, nonlinearBCs];
DeployPartialBoundaryConditions[{Null, stiffness}, linearBCs];
stiffness];


I'm in contact with Mathematica support about this. Meanwhile, I can offer a workaround. The code looks long below, but it is mostly just copied from above, with very few changes.

We need to define new functions PfemJacobian and PfemRHS to provide to FindRoot at the solution stage. These are alternatives to femJacobian and femRHS, provided in the documentation.

Needs["NDSolveFEM"];
PfemRHS[uV_?VectorQ] :=
dof}, NDSolveSetSolutionDataComponent[sdU, "DependentVariables",
uV];
nonlinear =
nonlinear = Null;
(*subtract the linear Robin boundary value*)
nonlinearBCs =
DiscretizeBoundaryConditions[initBCs, methodData, sdU,
"Nonlinear"];
stiffnessDummy = SparseArray[{}, {dof, dof}];
linearBCsPartial];
PfemJacobian[uV_?VectorQ] :=
Block[{stiffness, nonlinear, nonlinearStiffness, nonlinearBCs,
NDSolveSetSolutionDataComponent[sdU, "DependentVariables", uV];
nonlinear =
DiscretizePDE[linStiffnessPDEC, methodData, sdU, "Nonlinear"];
nonlinearStiffness = nonlinear["StiffnessMatrix"];
nonlinear = Null;
stiffness = linearStiffness + nonlinearStiffness;
nonlinearStiffness = Null;
nonlinearBCs =
DiscretizeBoundaryConditions[initBCs, methodData, sdU,
"Nonlinear"];
dof = Length[stiffness];
DeployPartialBoundaryConditions[{Null, stiffness}, nonlinearBCs];
linearBCsPartial];
stiffness];


Here is an explanation. The only difference with femRHS and femJacobian is that a second call to DeployPartialBoundaryConditions is replaced with a call to DeployBoundaryConditions (traditional way to deploy BCs when solving Linear PDEs), with globally defined discretized BC data named linearBCsPartial.

By inspecting the behavior of DeployPartialBoundaryConditions I concluded that it was not implementing the expected DirichletCondition because it had already been enforced on the seed data. Each iteration of the solver produces a change to the previous solution, and this change should have a zero Dirichlet condition on the target boundary, if the new solution is going to satisfy the desired Dirichlet condition of the full problem.

With these definitions, we continue mostly as before. I repeat the code from above so it is self-contained in this post. Defining the problem as before:

c = 1/Sqrt[(1 + Grad[u[x, y], {x, y}].Grad[u[x, y], {x, y}])];
Cu = {{{{c, 0}, {0, c}}}};
mesh = ToElementMesh[FullRegion[2], {{-1, 1}, {-1, 1}}];
uA[x_, y_] = y; (* Target solution *)


Now we define several separated boundary conditions

bcs = {DirichletCondition[u[x, y] == uA[x, y], -1 < x < 1],
PeriodicBoundaryCondition[u[x, y], x == 1, # - {2, 0} &]};
bcsDirichlet = {DirichletCondition[u[x, y] == uA[x, y], -1 < x < 1]};
bcsPartial = {DirichletCondition[u[x, y] == 0, -1 < x < 1],
PeriodicBoundaryCondition[u[x, y], x == 1, # - {2, 0} &]};


Note the zero Dirichlet condition for bcsPartial. Continuing as before:

uSeed[x_, y_] = (1 - 0.8 (1 - y^2)) uA[x, y];
iSeeding = {uSeed[x, y]};
vd = NDSolveVariableData[{"DependentVariables",
"Space"} -> {{u}, {x, y}}];
sd = NDSolveSolutionData[{"DependentVariables",
"Space"} -> {iSeeding, ToNumericalRegion[mesh]}];
coefficients = {"DiffusionCoefficients" -> Cu};
initCoeffs = InitializePDECoefficients[vd, sd, coefficients];


Here are the new statements to initialize the separated boundary conditions.

initBCs = InitializeBoundaryConditions[vd, sd, bcs] ;
initBCsDirichlet =
InitializeBoundaryConditions[vd, sd, bcsDirichlet] ;
initBCsPartial = InitializeBoundaryConditions[vd, sd, bcsPartial] ;


Continuing...

methodData =
InitializePDEMethodData[vd, sd, Method -> {"FiniteElement"}];
linearizedPDECoeffs = LinearizePDECoefficients[initCoeffs, vd, sd];
SplitPDECoefficients[linearizedPDECoeffs, vd, sd];
sdU = EvaluateInitialSeeding[methodData, vd, sd];
linear = DiscretizePDE[linearizedPDECoeffs, methodData, sdU];
linear["SystemMatrices"];


Here are the new statements to discretize the separated boundary conditions

linearBCs = DiscretizeBoundaryConditions[initBCs, methodData, sdU];
linearBCsDirichlet = DiscretizeBoundaryConditions[initBCsDirichlet, methodData, sdU];
linearBCsPartial = DiscretizeBoundaryConditions[initBCsPartial, methodData, sdU];


Because linearBCsDirichlet contains only the Dirichlet conditions, we can deploy this part using DeployDirichletConditions without worrying about ill effects due to PeriodicBoundaryCondition. (Although in this case it is not needed because the seed already satisfies the Dirichlet conditions.)

seed = NDSolveSolutionDataComponent[sdU, "DependentVariables"];
DeployDirichletConditions[seed, linearBCsDirichlet];


Finally, to solve, we call FindRoot with the new functions defined above PfemRHS and PfemJacobian.

root = U /.
FindRoot[PfemRHS[U], {U, seed}, Jacobian -> PfemJacobian[U],
Method -> {"AffineCovariantNewton"}];
NDSolveSetSolutionDataComponent[sdU, "DependentVariables", root];
{uf} = ProcessPDESolutions[methodData, sdU];
Plot3D[uf[x, y], Element[{x, y}, mesh]]


I'm not sure how general this workaround is, but it may be helpful for some.

• We have discussed several times some issue connected with compatibility between PeriodicBoundaryCondition and DirichletCondition in Mathematica FEM. Also your solution is helpful for some reason, and I gave you +1. Commented Nov 16, 2020 at 11:37

As alternative method we can use linear FEM to solve this problem

uSeed[x_, y_] = (1 - 0.3 (1 - x^2) (1 - y^2)) uA[x, y];
U[0][x_, y_] := uSeed[x, y]; n = 4;
Do[c1 = 1/
Sqrt[(1 +
Cu1 = {{{{c1, 0}, {0, c1}}}};
eqn1 = {Inactive[Div][
Cu1[[1, 1]].Inactive[Grad][u[x, y], {x, y}], {x, y}] == 0};
U[i] = NDSolveValue[{eqn1, {DirichletCondition[
u[x, y] == uA[x, y], -1 < x < 1],
PeriodicBoundaryCondition[u[x, y], x == 1, # - {2, 0} &]}}, u,
Element[{x, y}, mesh]];, {i, 1, n}]


Visualization of numerical solution and error on every step

Table[{Plot3D[U[i][x, y], Element[{x, y}, mesh],
AxesLabel -> Automatic, PlotRange -> All],
Plot3D[U[i][x, y] - uA[x, y], Element[{x, y}, mesh],
AxesLabel -> Automatic, PlotRange -> All, PlotLabel -> i]}, {i, n}]


As Figure 1 shows the error not decreasing with number of iteration increasing for i>2. Unfortunately this is the problem of compatibility of DirichletCondition[] and PeriodicBoundaryCondition[]. For instance, if we plot error=uf[x,y]-y for numerical solution from Will.Mo answer, then we got this picture with the same big error in the corner points: From the other side, if we exclude PeriodicBoundaryCondition[] from the code above, then we got higher precision numerical solution for n=30

Do[c1 = 1/
Sqrt[(1 +
Cu1 = {{{{c1, 0}, {0, c1}}}};
eqn1 = {Inactive[Div][
Cu1[[1, 1]].Inactive[Grad][u[x, y], {x, y}], {x, y}] == 0};
U[i] = NDSolveValue[{eqn1,
DirichletCondition[
u[x, y] == uA[x, y], (y == -1 || y == 1) && -1 <= x <= 1]}, u,
Element[{x, y}, mesh]];, {i, 1, 30}]

Table[Plot3D[U[i][x, y] - uA[x, y], Element[{x, y}, mesh],
AxesLabel -> Automatic, PlotRange -> All, PlotLabel -> i], {i, 25,
30}]



• Hi and thanks. I believe your code accomplishes something similar to my workaround in that NDSolveValue will use DeployBoundaryConditions for linear systems. However, your linearization is a bit different from what Mathematica is doing with nonlinear systems. Still, I agree that this also seems to work. You bring up a more general issue of compatibility between PeriodicBoundaryCondition and DirichletCondition. Note that this issue seems to go away with somewhat different BCs: Commented Nov 16, 2020 at 7:30
• Namely: {DirichletCondition[u[x, y] == uA[x, y], Or[y == 1, y == -1]], PeriodicBoundaryCondition[u[x, y], And[x == 1, -1 < y < 1], # - {2, 0} &]}}`. Also, the stray points where the BC fails may be related to a recent post of mine: mathematica.stackexchange.com/questions/230157/… Commented Nov 16, 2020 at 7:33
• @Will.Mo Thank you very much for the post linked. Commented Nov 16, 2020 at 11:36