1D mesh generation for PDE solution

Question

I'm trying to solve a system of two PDE that are dependent on time and distance (H[x,t] and P[x,t]). I'm solving the problem using the method of lines but I want to generate myself the mesh and introduce it in NDsolve. The mesh that I want to generate is the following

I need a mesh like this one because the values for one of the functions (P[x,t]) changes with time only very close to x = 0, whereas H[x,t] changes over the whole region 0< x < xmax. Below is an example of the code I'm using

(* Constants *)
f = 38.94; logL = -2;  
Ls = 10^logL; a = 0.5;            
C1 = 1*^-5; dH = 1*^-6;       
Ea = 0.1; 
tmax = 40; (* Time in seconds *)
xmax = 10 Sqrt[dH] Sqrt[tmax]; (* Maximum distance to simulate. cm *)


(* PDE system *)
 eqsH = {D[H[x, t], t] - dH D[H[x, t], x, x] ==  NeumannValue[Ls Exp[a f Ea ] P[x, t] - Ls Exp[-a f Ea ] H[x, t], 
    x == 0], H[x, 0] == 1};
eqsP = {D[P[x, t], t] == NeumannValue[-Ls Exp[a f Ea ] P[x, t] + Ls Exp[-a f Ea ] H[x, t], 
     x == 0], P[x, 0] == 1};

(*Solution of the differential equations*)
prec = 7;
msf = 0.001; 

sol = NDSolve[{eqsH, eqsP}, {H, P}, {x, 0, xmax}, {t, 0, tmax}, 
     AccuracyGoal -> prec, PrecisionGoal -> prec, 
     Method -> {"MethodOfLines", 
       "SpatialDiscretization" -> {"FiniteElement"}}] // First // 
   Quiet;

Can I please get some help in how to create the mesh and introduce it in NDSolve? thanks in advance !

Answer 1

Here is an alternate approach using a graded mesh.

Define some helper functions for a graded mesh

Here are some functions that I've used to create 1d to 3D anisotropic meshes. Not all functions are used.

(*Import required FEM package*)
Needs["NDSolve`FEM`"];
(* 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]
  ]
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]

Create a graded horizontal mesh segment

The following will create a horizontal mesh region of 100 elements where the initial element width is 1/10000 the domain length.

(*Create a graded horizontal mesh segment*)
(*Initial element width is 1/10000 the domain length*)
seg = leftSegmentGrowth[xmax, 100, xmax/10000];
Print["Horizontal segment"]
rh = pointsToMesh@seg
(*Convert mesh region to element mesh*)
(*Extract Coords from horizontal region*)
crd = MeshCoordinates[rh];
(*Create element mesh*)
mesh = ToElementMesh[crd];
Print["ListPlot of exponential growth of element size"]
ListPlot[Transpose@mesh["Coordinates"]]

One can see the exponential growth of the element size as the element number increases.

Convert the mesh region into an element mesh and solve the PDE

I generally convert the MeshRegion to an 'ElementMesh'so that I can apply element and point markers if needed.

(*Solve PDE on graded mesh*)
{Hfun, Pfun} = 
  NDSolveValue[{eqsH, eqsP}, {H, P}, x ∈ mesh, {t, 0, tmax}, 
   Method -> {"MethodOfLines", 
     "SpatialDiscretization" -> {"FiniteElement"}}];
(*Animate Hfun solution*)
imgs = Plot[Hfun[x, #], x ∈ mesh, 
     PlotRange -> {0.9999999, 1.0018}] & /@ Subdivide[0, tmax, 120];
Print["Animation of Hfun solution"]
ListAnimate@imgs

Appendix: Anisotropic meshing examples

As I alluded to in the comment below, the bullet list below shows several examples where I have used anisotropic quad meshing to capture sharp interfaces that would otherwise be very computationally expensive. The code is functional, but not optimal and some of the functions have been modified over time. Use at your own risk

1. 2D-Stationary
   • Mathematica vs. MATLAB: why am I getting different results for PDE with non-constant boundary condition?
   • Improving mesh and NDSolve solution convergence
2. 2D-Transient
   • Controlling dynamic time step size in NDSolveValue
   • How to model diffusion through a membrane?
   • Mass Transport FEM Using Quad Mesh
   • NDSolve with equation system with unknown functions defined on different domains
3. 3D-Meshing
   • Create graded mesh
4. 3D-Stationary
   • How to Improve FEM Solution with NDSolve?
   • 3D FEM Vector Potential

If you have access to other tools, such as COMSOL, that have boundary layer functionality, you can import meshes via the FEMAddOns resource function. It will not work for 3D meshes that require additional element types like prisms and pyramids that are not currently supported in Mathematica's FEM.

Answer 2

This is the same approach that Tim describes in his solution, but implemented with the ToGradedMesh function added in version 13.0:

mesh = ToGradedMesh[
   Line[{{0}, {xmax}}], <|"Alignment" -> "Left", 
    "ElementCount" -> 100, "MinimalDistance" -> xmax/10000|>];
MeshRegion[mesh]

eqsH = {D[H[x, t], t] - dH D[H[x, t], x, x] == 
    NeumannValue[Ls Exp[a f Ea] P[x, t] - Ls Exp[-a f Ea] H[x, t], 
     x == 0], H[x, 0] == 1};
eqsP = {D[P[x, t], t] == 
    NeumannValue[-Ls Exp[a f Ea] P[x, t] + Ls Exp[-a f Ea] H[x, t], 
     x == 0], P[x, 0] == 1};

{Hfun, Pfun} = 
  NDSolveValue[{eqsH, eqsP}, {H, P}, x \[Element] mesh, {t, 0, tmax}, 
   Method -> {"MethodOfLines", "SpatialDiscretization" -> {"FiniteElement"}}];

imgs = Plot[Hfun[x, #], x \[Element] mesh, 
     PlotRange -> {{0, 0.015}, {0.99999, 1.002}}] & /@ Subdivide[0, 2, 120];
ListAnimate@imgs

Answer 3

What about this?

lst1 = Partition[
  Join[Table[0.01*i, {i, 0, 5}], Table[0.1*i, {i, 0, 15}]], 1];
lst2 = Table[{i, i + 1}, {i, 1, Length[lst1] - 1}];

<< NDSolve`FEM`

mesh2 = ToElementMesh["Coordinates" -> lst1, 
  "MeshElements" -> {LineElement[lst2]}]

(*  ElementMesh[{{0., 1.5}}, {LineElement["<" 21 ">"]}]  *)

Let us visualize it:

mesh2["Wireframe"["MeshElementIDStyle" -> Red]]

The red figures indicate the mesh elements. The place where they overlap is in fact the one where the mesh is 10 times denser (see the blown up image below):

Have fun!