# Neumann boundary condition ignored

Im trying to solve a system of two PDE (im[t,x] and ia[t,x]) that are dependent on time and distance. I’m using my own anisotropic mesh (thanks to @TimLaska) that is aimed to represent an interface between a membrane and a liquid. The interface is located at L/2 where L = thickness of the system. “im” takes positive values in the interval 0 <= x <= L/2 and is equal to 0 at x > L/2. “ia” behaves similarly.

The problem that I’m encountering is that when I specify a Neumann boundary value at L/2, Mathematica doesn’t find this value in the mesh and the Neumann value is effectively ignored.

These are the functions for creating anisotropic meshes (kindly provided by @TimLaska in one of my previous questions)

(*Import required FEM package*)
Needs["NDSolveFEM"];
(*Define Some Helper Functions For Structured Meshes*)
pointsToMesh[data_] := MeshRegion[Transpose[{data}], Line@Table[{i, i + 1}, {i,Length[data]- 1}]];
unitMeshGrowth[n_, r_] :=  Table[(r^(j/(-1 + n)) - 1.)/(r - 1.), {j, 0, n - 1}]
meshGrowth[x0_, xf_, n_, r_] := (xf - x0) unitMeshGrowth[n, r] + x0
firstElmHeight[x0_, xf_, n_, r_] := Abs@First@Differences@meshGrowth[x0, xf, n, r]
lastElmHeight[x0_, xf_, n_, r_] := Abs@Last@Differences@meshGrowth[x0, xf, n, r]
findGrowthRate[x0_, xf_, n_, fElm_] :=  Quiet@Abs@ FindRoot[firstElmHeight[x0, xf, n, r] -    fElm, {r, 1.0001, 100000},Method -> "Brent"][[1, 2]]
meshGrowthByElm[x0_, xf_, n_, fElm_] := N@Sort@Chop@meshGrowth[x0, xf, n, findGrowthRate[x0, xf, n, fElm]]
meshGrowthByElm0[len_, n_, fElm_] := meshGrowthByElm[0, len, n, fElm]
flipSegment[l_] := (#1 - #2) & @@ {First[#], #} &@Reverse[l];
leftSegmentGrowth[len_, n_, fElm_] := meshGrowthByElm0[len, n, fElm]
rightSegmentGrowth[len_, n_, fElm_] := Module[{seg}, seg = leftSegmentGrowth[len, n, fElm];
flipSegment[seg]]
extendMesh[mesh_, newmesh_] := Union[mesh, Max@mesh + newmesh]


Here I define a couple of constants and create my own mesh

(* Constants *)
L = 0.1; dim = dia = 1.*^-6; kf = kr = 1; kr =
1./kf;
imo[i_] := Piecewise[{{0.1, 0 <= i <= L/2.}, {0, i > L/2.}}]
iao[i_] := Piecewise[{{0, 0 <= i <= L/2.}, {0.1, i > L/2.}}]
(* Mesh generation *)
seg1 = rightSegmentGrowth[L/2., 100, L/1000];
seg2 = leftSegmentGrowth[L/2., 100, L/1000];
totalseg = extendMesh[seg1, seg2];
rh = pointsToMesh@totalseg ;
crd = MeshCoordinates[rh];
mesh = ToElementMesh[crd];


And here is the simple code I use to find the solution alongside the error I get

(*PDE system*)

eqsim = {D[im[t, x], t] - dim D[im[t, x], x, x] == NeumannValue[kf*ia[t, x], x == L/2.] + NeumannValue[0, x == L/2.],im[0, x] == imo[x]};
eqsia = {D[ia[t, x], t] - dia D[ia[t, x], x, x] ==  NeumannValue[kr*im[t, x], x == L/2.] + NeumannValue[0, x ==   L/2.], ia[0, x] == iao[x]};
sol = NDSolve[{eqsim, eqsia}, {im, ia},x \[Element] mesh, {t, 0, 60}, Method ->{"MethodOfLines", "SpatialDiscretization" -> {"FiniteElement"}}]
(*NDSolve::bcnop: No places were found on the boundary where x==0.05 was True, so NeumannValue[ia,x==0.05] will effectively be ignored.*)
(*NDSolve::bcnop: No places were found on the boundary where x==0.05 was True, so NeumannValue[1. im,x==0.05] will effectively be ignored.*)
(* NDSolve::bcnop: No places were found on the boundary where x==0.05 was True, so NeumannValue[0,x==0.05] will effectively be ignored.*)
(*General::stop: Further output of NDSolve::bcnop will be suppressed during this calculation.*)


I first thought that the discretisation in the mesh didn’t include the point L/2 but after double checking I found out that it’s there

Position[crd, L/2.]
(*{{100, 1}}*)


Similar problems to this have been reported when using Dirichlet conditions here DiscretizeRegion does not include the boundary specified in ImplicitRegion (10.1) , here PeriodicBoundaryConditions: missing points (a simpler example) and here Solving Laplace's equation in 2D using region primitives , but after reading the solutions provided I haven’t been able to find a solution to my issue.

Any help, feedback or advice is more than welcome (I’m using v 12.2.0.0).

• Neumann Values need to be applied to boundaries (0 or 0.1) only. They cannot be applied to internal nodes. The Neumann value takes into account the surface normal which does not exist for an internal node. Commented Jan 18, 2021 at 18:08
• You also may want to take a look at my answer to "how to model diffusion through a membrane?" here 219140. It is a 2D problem but the same principle can be used for a 1D problem. Commented Jan 18, 2021 at 22:41
• Thanks a lot Tim.
– Afmo
Commented Jan 19, 2021 at 7:01
• There is an example of a phase change in the documentation. See PDEModels/tutorial/HeatTransfer/HeatTransfer#135129751 or online here Commented Jan 19, 2021 at 9:51
• Thanks @user21. I'm working now with the interphase mass transfer model provided in the documentation which I find very useful (reference.wolfram.com/language/PDEModels/tutorial/MassTransport/…). I just thought in my ignorance that by including a NeumannValue at the interface, I could coupled fluxes and avoid the use of the interfacial region used in the documented example. I now realise that there was a reason for that. I wonder if there is another mathematical way to include a flux coupling at certain point in the mesh apart from the boundaries, maybe there isn't.
– Afmo
Commented Jan 19, 2021 at 10:00

Here's one way to solve the problem by explicitly modeling the interface region (note that this is explained in more detail here). We start 1st by defining some helper functions.

# Mesh helper functions (edited)

These helper functions are similar to those described in the OP but have a couple of additions.

(*Import required FEM package*)
Needs["NDSolveFEM"];
(* Define Some Helper Functions For Structured Meshes*)
pointsToMesh[data_] :=
MeshRegion[Transpose[{data}],
Line@Table[{i, i + 1}, {i, Length[data] - 1}]];
unitMeshGrowth[n_, r_] :=
Table[(r^(j/(-1 + n)) - 1.)/(r - 1.), {j, 0, n - 1}]
meshGrowth[x0_, xf_, n_, r_] := (xf - x0) unitMeshGrowth[n, r] + x0
firstElmHeight[x0_, xf_, n_, r_] :=
Abs@First@Differences@meshGrowth[x0, xf, n, r]
lastElmHeight[x0_, xf_, n_, r_] :=
Abs@Last@Differences@meshGrowth[x0, xf, n, r]
findGrowthRate[x0_, xf_, n_, fElm_] :=
Quiet@Abs@
FindRoot[firstElmHeight[x0, xf, n, r] - fElm, {r, 1.0001, 1/fElm},
Method -> "Brent"][[1, 2]]
meshGrowthByElm[x0_, xf_, n_, fElm_] :=
N@Sort@Chop@meshGrowth[x0, xf, n, findGrowthRate[x0, xf, n, fElm]]
meshGrowthByElm0[len_, n_, fElm_] := meshGrowthByElm[0, len, n, fElm]
flipSegment[l_] := (#1 - #2) & @@ {First[#], #} &@Reverse[l];
leftSegmentGrowth[len_, n_, fElm_] := meshGrowthByElm0[len, n, fElm]
rightSegmentGrowth[len_, n_, fElm_] := Module[{seg},
seg = leftSegmentGrowth[len, n, fElm];
flipSegment[seg]
]
reflectRight[pts_] := With[{rt = ReflectionTransform[{1}, {Last@pts}]},
Union[pts, Flatten[rt /@ Partition[pts, 1]]]]
reflectLeft[pts_] :=
With[{rt = ReflectionTransform[{-1}, {First@pts}]},
Union[pts, Flatten[rt /@ Partition[pts, 1]]]]
extendMesh[mesh_, newmesh_] := Union[mesh, Max@mesh + newmesh]


# Build 3 region mesh

The following workflow will build a 3 material region mesh representing the left side, a thin interface region, and the right side.

(*Constants*)
l = 0.1;
hl = l/2;
int = l/1000;
hint = int/2;
hlrem = hl - hint;
nhint = 5;
nhl = 100;
intelm = hint/nhint;
(*Build right-hand side of the mesh*)
hsegint = Subdivide[0, hint, nhint];
rint = pointsToMesh@hsegint
hseg = leftSegmentGrowth[hlrem, 100, intelm/10];
rhseg = pointsToMesh@hseg
(*Extend and reflect mesh about origin*)
r = pointsToMesh@reflectLeft@extendMesh[hsegint, hseg]
(* Extract Coords from RegionProduct *)
crd = MeshCoordinates[r];
(* grab line element incidents RegionProduct mesh *)
inc = Delete[0] /@ MeshCells[r, 1];
mesh = ToElementMesh["Coordinates" -> crd,
"MeshElements" -> {LineElement[inc]}];
(* Association for Clearer Region Assignment *)
reg = <|"left" -> 1, "int" -> 2, "right" -> 3|>;
(*Define regions for region marker function*)
Ωleft = ImplicitRegion[-hl <= x <= -hint, {x}];
Ωint = ImplicitRegion[-hint <= x <= hint, {x}];
Ωright = ImplicitRegion[hint <= x <= hl, {x}];
(* RegionMember Function *)
rmfl = RegionMember[Ωleft];
rmfi = RegionMember[Ωint];
rmfr = RegionMember[Ωright];
regmarkerfn =
Piecewise[{{reg["left"], rmfl[#1]}, {reg["int"], rmfi[#1]}},
reg["right"]] &;
(* Get mean coordinate of each hexa for region marker assignment *)
mean = Mean /@ GetElementCoordinates[mesh["Coordinates"], #] & /@
ElementIncidents[mesh["MeshElements"]] // First;
regmarkers = regmarkerfn /@ mean;
mesh = ToElementMesh["Coordinates" -> mesh["Coordinates"],
"MeshElements" -> {LineElement[inc, regmarkers]}];


# Set up and solve PDE system

We will use the MassTransportPDEComponentto create a PDE operator that depends on region dependent coefficients as shown below:

(*Mass transport parameters*)
dl = 10^-6;
dr = 10^-6;
kf = 0.1; kr = 1/8 kf;
K = kf/kr;
k0 = 10^4;
(*Initial condition function*)
imo[i_] := Piecewise[{{0.1, i <= -hint}, {0, i > hint}}]
iao[i_] := Piecewise[{{0, i < hint}, {0.1, i >= hint}}]
(*Define switch to activate reactions in membrane region*)
sigma = Evaluate[
Piecewise[{{1, ElementMarker == reg["int"]}, {0, True}}]];
(*Defined region dependent diffusion coefficients*)
diffleft =
If[ElementMarker == reg["left"] || ElementMarker == reg["int"], dl,
0];
diffright =
If[ElementMarker == reg["right"] || ElementMarker == reg["int"],
dr, 0];;
(*PDE system*)
op = MassTransportPDEComponent[{{im[t, x], ia[t, x]}, t, {x}}, <|
"DiffusionCoefficient" -> {{diffleft, 0}, {0, diffright}},
"MassReactionRate" -> {{-sigma k0 (K ia[t, x] - im[t, x])/
im[t, x], 0}, {0,
sigma k0 (K ia[t, x] - im[t, x])/ia[t, x]}}|>];
eqns = op == {0, 0};
{imFun, iaFun} =
NDSolveValue[{eqns, im[0, x] == imo[x], ia[0, x] == iao[x]}, {im,
ia}, x ∈ mesh, {t, 0, 60}];
imgs = Plot[{imFun[#, x], iaFun[#, x]}, x ∈ mesh,
PlotPoints -> 500,
PlotRange -> {0, Max@imFun["ValuesOnGrid"]}] & /@
Subdivide[0, 60, 60];
ListAnimate@imgs


It is straightforward to change parameters and initial conditions and estimate their effect on the solution. Here, we will initialize the left-hand side at 0 concentration and give the right-hand side a much higher diffusion coefficient similar to a gas and watch the system evolve.

(*Mass transport parameters*)
dl = 10^-6;
dr = 10^-3;
kf = 0.1; kr = 1/8 kf;
K = kf/kr;
k0 = 10^4;
(*Defined region dependent diffusion coefficients*)
diffleft =
If[ElementMarker == reg["left"] || ElementMarker == reg["int"], dl,
0];
diffright =
If[ElementMarker == reg["right"] || ElementMarker == reg["int"], dr,
0];
(*Initial condition function*)
imo[i_] := Piecewise[{{0, i <= -hint}, {0, i > hint}}]
iao[i_] := Piecewise[{{0, i < hint}, {0.1, i >= hint}}]
op = MassTransportPDEComponent[{{im[t, x], ia[t, x]}, t, {x}}, <|
"DiffusionCoefficient" -> {{diffleft, 0}, {0, diffright}},
"MassReactionRate" -> {{-sigma k0 (K ia[t, x] - im[t, x])/
im[t, x], 0}, {0,
sigma k0 (K ia[t, x] - im[t, x])/ia[t, x]}}|>];
eqns = op == {0, 0};
{imFun, iaFun} =
NDSolveValue[{eqns, im[0, x] == imo[x], ia[0, x] == iao[x]}, {im,
ia}, x ∈ mesh, {t, 0, 60}];
imgs = Plot[{imFun[#, x], iaFun[#, x]}, x ∈ mesh,
PlotPoints -> 500,
PlotRange -> {0, Max@imFun["ValuesOnGrid"]}] & /@
Subdivide[0, 60, 60];
ListAnimate@imgs


# Modeling very thin membranes

@Afmo had a question in the comments about the numerical stability of modeling a 1 nm membrane on a 1 cm domain (a 7 order of magnitude difference). Now, if you're not interested in what happens within the membrane material and just want to use it as a means to model the jump condition, the membrane does not need to be 1 nm. One can simply merge the 2 fluid domains as a postprocessing step and ignore the membrane itself. However, in this case, the anisotropic meshing procedure still appears to be numerically stable over a 7 order magnitude ratio of domain size to membrane size as shown in the following workflow.

## Very thin membrane mesh

To maintain numerical stability, it is generally a good practice to at least match the element thickness at the interface. An example is shown with the following workflow (please note that I adjusted the mesh helper functions to accommodate these huge aspect ratio differences):

(*Constants*)
l = 0.1;
hl = l/2;
(*7 order magnitude difference between membrane thickness and domain \
size*)
int = l/10000000;
hint = int/2;
hlrem = hl - hint;
nhint = 1;
nhl = 100;
intelm = hint/nhint;
(*Build right-hand side of the mesh*)
hsegint = Subdivide[0, hint, nhint];
rint = pointsToMesh@hsegint
hseg = leftSegmentGrowth[hlrem, nhl, intelm];
rhseg = pointsToMesh@hseg
(*Extend and reflect mesh about origin*)
r = pointsToMesh@reflectLeft@extendMesh[hsegint, hseg]
(* Extract Coords from RegionProduct *)
crd = MeshCoordinates[r];
(* grab line element incidents RegionProduct mesh *)
inc = Delete[0] /@ MeshCells[r, 1];
mesh = ToElementMesh["Coordinates" -> crd,
"MeshElements" -> {LineElement[inc]}];
(* Association for Clearer Region Assignment *)
reg = <|"left" -> 1, "int" -> 2, "right" -> 3|>;
(*Define regions for region marker function*)
Ωleft = ImplicitRegion[-hl <= x <= -hint, {x}];
Ωint = ImplicitRegion[-hint <= x <= hint, {x}];
Ωright = ImplicitRegion[hint <= x <= hl, {x}];
(* RegionMember Function *)
rmfl = RegionMember[Ωleft];
rmfi = RegionMember[Ωint];
rmfr = RegionMember[Ωright];
regmarkerfn =
Piecewise[{{reg["left"], rmfl[#1]}, {reg["int"], rmfi[#1]}},
reg["right"]] &;
(* Get mean coordinate of each hexa for region marker assignment *)
mean = Mean /@ GetElementCoordinates[mesh["Coordinates"], #] & /@
ElementIncidents[mesh["MeshElements"]] // First;
regmarkers = regmarkerfn /@ mean;
mesh = ToElementMesh["Coordinates" -> mesh["Coordinates"],
"MeshElements" -> {LineElement[inc, regmarkers]}];


## Set up and solve PDE system

Set up and solve the PDE system as before:

(*Mass transport parameters*)
dl = 10^-6;
dr = 10^-6;
kf = 0.1; kr = 1/8 kf;
K = kf/kr;
k0 = 10^4;
(*Initial condition function*)
imo[i_] := Piecewise[{{0.1, i <= -hint}, {0, i > hint}}]
iao[i_] := Piecewise[{{0, i < hint}, {0.1, i >= hint}}]
(*Define switch to activate reactions in membrane region*)
sigma = Evaluate[
Piecewise[{{1, ElementMarker == reg["int"]}, {0, True}}]];
(*Defined region dependent diffusion coefficients*)
diffleft =
If[ElementMarker == reg["left"] || ElementMarker == reg["int"], dl,
0];
diffright =
If[ElementMarker == reg["right"] || ElementMarker == reg["int"],
dr, 0];;
(*PDE system*)
op = MassTransportPDEComponent[{{im[t, x], ia[t, x]}, t, {x}}, <|
"DiffusionCoefficient" -> {{diffleft, 0}, {0, diffright}},
"MassReactionRate" -> {{-sigma k0 (K ia[t, x] - im[t, x])/
im[t, x], 0}, {0,
sigma k0 (K ia[t, x] - im[t, x])/ia[t, x]}}|>];
eqns = op == {0, 0};
{imFun, iaFun} =
NDSolveValue[{eqns, im[0, x] == imo[x], ia[0, x] == iao[x]}, {im,
ia}, x ∈ mesh, {t, 0, 60},
StartingStepSize -> 0.000000001];
imgs = Plot[{imFun[#, x], iaFun[#, x]}, x ∈ mesh,
PlotPoints -> 500,
PlotRange -> {0, Max@imFun["ValuesOnGrid"]}] & /@
Subdivide[0, 60, 60];
ListAnimate@imgs


Obviously, one will ultimately reach a limit, but we been able to extend the anisotropic meshing to an impressive 7 order magnitude difference without evidence of great numerical instability.

• Thanks a lot Tim, It's excellent. I was just worried about running into numerical problems when using an interfacial region in cases where the simulation requires a very long domain and a very thin interface. For instance, if the interface is 1 nm and the domain is 1 cm. Your approach works very well though.
– Afmo
Commented Jan 20, 2021 at 7:24
• @Afmo I made a slight modification to the mesh helper functions and is able to extend the anisotropic meshing to handle a 1 nm interface without much evidence of numerical instability. Commented Jan 21, 2021 at 1:53
• Excellent, thanks again Tim.
– Afmo
Commented Jan 21, 2021 at 7:46
• @Afmo, Tim made a proposal to add anisotropic meshing capabilities to the FEM mesh generation here. If you think this is worthwhile, consider giving this proposal an upvote Commented Feb 22, 2021 at 10:29