6
$\begingroup$

I want to see the effect of mesh size on the solution curves while simulating convection-diffusion - point sink via NDSolve FEM.

In the following code (ref also related to this post):

I'd like to know how to visualize the mesh size distribution i.e mesh element size defined in mesh = ToElementMesh[region, "IncludePoints" -> includePoints, "MaxCellMeasure" -> 1];. Also, from what I read here, I understand MaxCellMeasure sets the upper limit of the mesh element size. I want to know if there is a similar option for defining the minimum size of mesh elements. e.g. MinCellMeasure argument ? I couldn't find this in the documentation.

Needs["NDSolve`FEM`"]
region = Line[{{0}, {100}}];
includePoints = {{10}};
mesh = ToElementMesh[region, "IncludePoints" -> includePoints, 
   "MaxCellMeasure" -> 1];
vars = {c[t, x], t, {x}};
RegularizedDeltaPoint[g_, X_List, Xs_List] := 
 Piecewise[{{Times @@ Thread[1/(4 g) (1 + Cos[\[Pi]/(2 g) (X - Xs)])],
     And @@ Thread[RealAbs[X - Xs] <= 2 g]}, {0, True}}]
Subscript[h, mesh] = Sqrt[Min[mesh["MeshElementMeasure"]]];
Subscript[gamma, reg] = Subscript[h, mesh]/2;
temp = RegularizedDeltaPoint[Subscript[gamma, reg], {x}, 
   includePoints[[1]]];
parameters = {kappa -> {{1000}}, v1 -> 100, 
   gamma -> Subscript[gamma, reg], Qp -> 3};
pde = {Derivative[1, 0][c][t, x] + 
      Inactive[
        Div][(-kappa).Inactive[Grad][
         c[t, x], {x}], {x}] + {v1}.Inactive[Grad][c[t, x], {x}] + 
      Qp*RegularizedDeltaPoint[gamma, {x}, {10}] == 0, 
    c[0, x] == 1} /. parameters;

tEnd = 3;
cfun = NDSolveValue[
   pde~Join~{DirichletCondition[c[t, x] == 4 (1 - Exp[-1000 t]) + 1, 
      x == 0]}, c, {t, 0, tEnd}, {x} \[Element] mesh];
With[{i = Flatten[{0}~Join~includePoints]}, 
 Plot[Evaluate[cfun[t, #] & /@ i], {t, 0, tEnd}, 
  PlotRange -> {{0, tEnd}, {0, 5.1}}, 
  PlotLegends -> (StringTemplate["c(t,``)"] /@ i)]]
ListPlot[Table[cfun[t, x], {x, 0, 100, 10}, {t, 0, tEnd, tEnd/100}], 
 Joined -> True, DataRange -> {0, tEnd}]

Suggestions will be really helpful.

$\endgroup$
8
$\begingroup$

There is no MinCellMeasure, but you can build your own mesh, thereby giving you complete control over the discretization. I will show an example of how to do this at the end of the post. With respect to mesh visualization, your current mesh is pretty finely discretized, making it difficult to view the distributions. Also, your parameters make it difficult to see your point source term. So, let us make some parameter changes (coarsen the mesh, lower the diffusivity and velocity, and crank up the point source term) to better visualize some of the features I think the OP wants to capture.

Needs["NDSolve`FEM`"]
region = Line[{{0}, {100}}];
includePoints = {{10}};
mesh = ToElementMesh[region, "IncludePoints" -> includePoints, 
   "MaxCellMeasure" -> 5];
vars = {c[t, x], t, {x}};
RegularizedDeltaPoint[g_, X_List, Xs_List] := 
 Piecewise[{{Times @@ Thread[1/(4 g) (1 + Cos[π/(2 g) (X - Xs)])],
     And @@ Thread[RealAbs[X - Xs] <= 2 g]}, {0, True}}]
hmesh = Sqrt[Min[mesh["MeshElementMeasure"]]];
gammareg = hmesh/2;
temp = RegularizedDeltaPoint[gammareg, {x}, includePoints[[1]]];
parameters = {kappa -> {{1}}, v1 -> 10, gamma -> gammareg, Qp -> 10};
pde = {Derivative[1, 0][c][t, x] + 
      Inactive[Div][(-kappa) . 
        Inactive[Grad][c[t, x], {x}], {x}] + {v1} . 
       Inactive[Grad][c[t, x], {x}] + 
      Qp*RegularizedDeltaPoint[gamma, {x}, {10}] == 0, c[0, x] == 1} /. 
   parameters;

tEnd = 3;
cfun = NDSolveValue[
   pde~Join~{DirichletCondition[c[t, x] == 4 (1 - Exp[-1000 t]) + 1, 
      x == 0]}, c, {t, 0, tEnd}, {x} ∈ mesh];
With[{i = Flatten[{0}~Join~includePoints]}, 
 Plot[Evaluate[cfun[t, #] & /@ i], {t, 0, tEnd}, 
  PlotRange -> {{0, tEnd}, {-0.5, 5.1}}, 
  PlotLegends -> (StringTemplate["c(t,``)"] /@ i), Frame -> True, 
  FrameLabel -> {{"Concentration", None}, {"Time", 
     "Concentrations @ 0 an 10 x units"}}]]
ListPlot[Table[cfun[t, x], {x, 0, 100, 10}, {t, 0, tEnd, tEnd/100}], 
 Joined -> True, DataRange -> {0, tEnd}, 
 PlotRange -> {{0, tEnd}, {-0.5, 5.1}}, 
 PlotLegends -> (StringTemplate["c(t,``)"] /@ Subdivide[100, 10]), 
 Frame -> True, 
 FrameLabel -> {{"Concentration", None}, {"Time", 
    "Concentrations every 10 x units"}}]
Plot[cfun[0.01, x], x ∈ mesh, 
 PlotRange -> {{9, 11}, {0.4, 1.05}}, Frame -> True, 
 FrameLabel -> {{"Concentration", None}, {"X", 
    "Concentration near x=10 for small time"}}]

Coarse solution

Because the mesh is so coarse, NDSolve complains about high Péclet numbers.

Mesh visualization

For a 1D mesh, you could use ListPlot on the mesh coordinates, or you can convert the mesh into a MeshRegion and use HighlightMesh to visualize the discretization.

ListPlot[Union@Flatten@mesh["Coordinates"], Frame -> True, 
 FrameLabel -> {{"Node Position", None}, {"Node ID", 
    "Node Position vs Node ID"}}]
HighlightMesh[MeshRegion@mesh, 0]

Coarse mesh discretization

Anisotropic meshing

In my answer here 240408, there are many links to examples where I have answered questions using an anisotropic meshing approach. This approach will allow you to capture very sharp features without much computational expense.

Helper functions

First, we will define some helper functions to create the anisotropic mesh.

(*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, 0.00000001, 
      100000000/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 an anisotropic mesh

The above helper functions allow build anisotropic mesh segments and glue them together. The following will create an anisotropic mesh with refinements at 0 and to either side of 10.

(*Define some parameters*)
x1 = 10; x2 = 90; nelm1 = 40; nelm2 = 50;
(*Create mesh segment with finer mesh at x=0 and 10*)
seg1 = reflectRight@leftSegmentGrowth[x1/2, nelm1, x1/200];
Print["Left of 10 domain"]
r1 = pointsToMesh@seg1
seg2 = leftSegmentGrowth[x2, nelm2, x1/200];
Print["Right of 10"]
r2 = pointsToMesh@seg2
Print["Full domain"]
regfull = pointsToMesh@extendMesh[seg1, seg2]
(*Extract Coords from region*)
crd = MeshCoordinates[regfull];
(*Create element mesh*)
mesh = ToElementMesh[crd];
Print["Coordinate positions on a log scale"]
ListLogPlot[Transpose@mesh["Coordinates"], PlotRange -> All, 
 Frame -> True, 
 FrameLabel -> {{"Node Position", None}, {"Node ID", 
    "Node Position vs Node ID"}}]

Anisotropic mesh

Solve on an anisotropic mesh

gammareg = hmesh/400;
temp = RegularizedDeltaPoint[gammareg, {x}, includePoints[[1]]];
parameters = {kappa -> {{1}}, v1 -> 10, gamma -> gammareg, Qp -> 10};
pde = {Derivative[1, 0][c][t, x] + 
      Inactive[Div][(-kappa) . 
        Inactive[Grad][c[t, x], {x}], {x}] + {v1} . 
       Inactive[Grad][c[t, x], {x}] + 
      Qp*RegularizedDeltaPoint[gamma, {x}, {10}] == 0, c[0, x] == 1} /. 
   parameters;

tEnd = 3;
cfun = NDSolveValue[
   pde~Join~{DirichletCondition[c[t, x] == 4 (1 - Exp[-1000 t]) + 1, 
      x == 0]}, c, {t, 0, tEnd}, {x} ∈ mesh];
With[{i = Flatten[{0}~Join~includePoints]}, 
 Plot[Evaluate[cfun[t, #] & /@ i], {t, 0, tEnd}, 
  PlotRange -> {{0, tEnd}, {-0.5, 5.1}}, 
  PlotLegends -> (StringTemplate["c(t,``)"] /@ i), Frame -> True, 
  FrameLabel -> {{"Concentration", None}, {"Time", 
     "Concentrations @ 0 an 10 x units"}}]]
ListPlot[Table[cfun[t, x], {x, 0, 100, 10}, {t, 0, tEnd, tEnd/100}], 
 Joined -> True, DataRange -> {0, tEnd}, 
 PlotRange -> {{0, tEnd}, {-0.5, 5.1}}, 
 PlotLegends -> (StringTemplate["c(t,``)"] /@ Subdivide[100, 10]), 
 Frame -> True, 
 FrameLabel -> {{"Concentration", None}, {"Time", 
    "Concentrations every 10 x units"}}]
Plot[cfun[0.01, x], x ∈ mesh, 
 PlotRange -> {{9, 11}, {0.4, 1.05}}, Frame -> True, 
 FrameLabel -> {{"Concentration", None}, {"X", 
    "Concentration near x=10 for small time"}}]

Solutions on an anisotropic mesh

These solutions behave much better than the initial coarse solution and it allows one to resolve a very fine point source without blowing up the model size.

$\endgroup$
2
  • $\begingroup$ Thanks a lot for the detailed answer. I'm sorry I am new to MMA and FEM. Please excuse me for a few more questions. Could you please explain the legend of the second plot in the post? I would like to know what the green and orange colored solution curves correspond to. $\endgroup$
    – Natasha
    Jun 28 at 6:19
  • 1
    $\begingroup$ @Natasha I added some FrameLabels and legends for clarity. $\endgroup$
    – Tim Laska
    Jun 28 at 12:44

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.