Solution of PDE system does no satisfy boundary conditions

Question

I have a problem with solution of the 2D PDE system. It appears that solution does not match boundary condition at all.

I've tried using separately generated mesh, and it had not helped, just making error oscillation frequency larger, but not lowering the error amplitude.

Adding "Method -> {"ExplicitRungeKutta", "SpatialDiscretization" -> {"FiniteElement", "MeshOptions" -> {"MaxCellMeasure" -> 0.00001}}}" does not make any difference.

The region is a rectangle (I tried oval too) with two round holes. There is no such problem in the case of just one hole.

Do you have any ideas making it better?

    system = {nxU[x, y] == D[HU[x, y], x] + (D[nyU[x, y], x, y] + D[nxU[x, y], x, x]),
    nyU[x, y] == D[HU[x, y], y] + (D[nyU[x, y], y, y] + D[nxU[x, y], x, y]),
    D[nyU[x, y], y] + D[nxU[x, y], x] == D[HU[x, y], y, y] + D[HU[x, y],x,x]};

    FuncLst = {nxU[x, y], nyU[x, y], HU[x, y]};

    n0=0.3;hIn0=-0.5;rOut=30;r0=1.;Centre1={-10,0};Centre2={10,0};

    RegionOut=Rectangle[{-rOut,-rOut/2},{rOut,rOut/2}];
    BoundaryOut=RegionBoundary[RegionOut];

    RegionIn1=Disk[Centre1,r0];
    RegionIn2=Disk[Centre2,r0];
    RegionIn=RegionUnion[RegionIn1,RegionIn2];
    BoundaryIn1=RegionBoundary[RegionIn1];
    BoundaryIn2=RegionBoundary[RegionIn2];
    BoundaryIn=RegionBoundary[RegionIn];

    RegionAll=RegionUnion[RegionDifference[RegionOut,RegionIn],BoundaryIn];

    n0func1x[x_,y_]=n0 (x-Centre1[[1]])/r0;
    n0func1y[x_,y_]=n0 (y-Centre1[[2]])/r0;
    n0func2x[x_,y_]=n0 (x-Centre2[[1]])/r0;
    n0func2y[x_,y_]=n0 (y-Centre2[[2]])/r0;

    BndCond={DirichletCondition[nxU[x,y]==0,{x,y}\[Element]BoundaryOut],
    DirichletCondition[nyU[x,y]==0,{x,y}\[Element]BoundaryOut],
    DirichletCondition[HU[x,y]==0,{x,y}\[Element]BoundaryOut],
    DirichletCondition[nxU[x,y]==n0func1x[x,y],{x,y}\[Element]BoundaryIn1],
    DirichletCondition[nyU[x,y]==n0func1y[x,y],{x,y}\[Element]BoundaryIn1],
    DirichletCondition[nxU[x,y]==n0func2x[x,y],{x,y}\[Element]BoundaryIn2],
    DirichletCondition[nyU[x,y]==n0func2y[x,y],{x,y}\[Element]BoundaryIn2],
    DirichletCondition[HU[x,y]==hIn,{x,y}\[Element]BoundaryIn1],
    DirichletCondition[HU[x,y]==hIn,{x,y}\[Element]BoundaryIn2]};

    solPDE = NDSolve[Join[system, (BndCond /. hIn -> hIn0)], Join[FuncLst],{x, y} \[Element] RegionAll][[1]];
    nxUs[x_, y_] = Evaluate[nxU[x, y] /. solPDE];
    xs[t_] = Centre1[[1]] + Cos[t];
    ys[t_] = Centre1[[2]] + Sin[t];
    Plot[{nxUs[xs[t], ys[t]], n0func1x[xs[t], ys[t]]}, {t, 0, 2 \[Pi]}]

Answer 1

This is not an answer (yet), just a bunch of observations and suspicions that will hopefully lead a savvier person than myself to find a solution.

The problem, as far as I can tell, stems from the InterpolationOrder of the results returned by NDSolve. By default, the solutions returned by NDSolve seem to be second-order interpolating functions in both directions:

Needs["DifferentialEquations`InterpolatingFunctionAnatomy`"];
InterpolatingFunctionInterpolationOrder[Head[nxU[x, y] /. solPDE]]

(* {2, 2} *)

What's odd is that NDSolve is supposed to allow you to specify the InterpolationOrder of the function it returns; but for this problem (maybe for all finite-element methods?), NDSolve just seems to ignore this option:

solPDE2 = NDSolve[Join[system, (BndCond /. hIn -> hIn0)], Join[FuncLst], {x, y} \[Element] mesh, InterpolationOrder -> 1][[1]];
InterpolatingFunctionInterpolationOrder[Head[nxU[x, y] /. solPDE2]]

(* {2, 2} *)

OK... let's force a first-order mesh instead, which should yield a first-order interpolating function:

Needs["NDSolve`FEM`"]
mesh = ToElementMesh[RegionAll, "MeshOrder" -> 1];
solPDE3 = NDSolve[Join[system, (BndCond /. hIn -> hIn0)],    
                  Join[FuncLst], {x, y} \[Element] mesh][[1]];

(* {1, 1} *)

Excellent. What do the boundaries look like?

nxUsnew[x_, y_] = Evaluate[nxU[x, y] /. solPDE3];
Plot[{nxUsnew[xs[t], ys[t]], n0func1x[xs[t], ys[t]]}, {t, 0, 2 \[Pi]}]

.... OK then.

We can see that at most of the actual points on the mesh, the boundary conditions are pretty well satisfied. This is probably true for the original function as well: the result of NDSolve is in contact with the true boundary function at several points, presumably the mesh points; and it's only the regions between the mesh points (where interpolation is necessary) that disagree. Also, running this version of the code with "MeshOrder" -> 2 yields something similar (though not identical) to your original result.

However, I have no idea why the points at $\theta = 0$ and $\theta = \pi$ (i.e., when $y = 0$) are off of the curve. On my system (Mac OS, MM 10.4) these points are aberrant both when using the ToElementMesh method and when using your original method to find the solution. However, your original graphs don't seem to have this problem. Perhaps you might have better luck using the ToElementMesh method on your system.

Still, it would be helpful to know the answers to two questions:

• Is there a way to get NDSolve to use a first-order mesh on a finite-element problem, without resorting to loading additional packages?
• What's up with the points at $y = 0$ being off the curve, regardless of the method used?