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"}}]
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]
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"}}]
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"}}]
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.