# Problem with NDSolveValue

I am trying to solve some 3D PDEs (for example, 3D diffusion equation) with step function or DiracDelta function as initial condition. I know those functions have discontinuity so I use some continuous function to approximate them (Tanh for step function, Gaussian for DiracDelta). The following code is an example of using a Gaussian function to approximate a DiracDelta function as initial condition:

\[CapitalOmega] = Cuboid[{0, 0, 0}, {1, 1, 1}];
sol = NDSolveValue[{D[u[t, x, y, z], t] -
Laplacian[u[t, x, y, z], {x, y, z}] == 0,
PeriodicBoundaryCondition[u[t, x, y, z], x == 0,
TranslationTransform[{1, 0, 0}]],
PeriodicBoundaryCondition[u[t, x, y, z], y == 0,
TranslationTransform[{0, 1, 0}]],
PeriodicBoundaryCondition[u[t, x, y, z], z == 0,
TranslationTransform[{0, 0, 1}]],
u[0, x, y, z] ==
1000000/\[Pi]^(3/2)
E^(-10000 ((x - 0.5)^2 + (y - 0.5)^2 + (z - 0.5)^2))},
u, {t, 0, 1}, {x, y, z} \[Element] \[CapitalOmega]]

Plot[ {1000000/\[Pi]^(3/2)
E^(-10000 ((x - 0.5)^2 + (x - 0.5)^2 + (x - 0.5)^2)),
sol[0, x, x, x]}, {x, 0, 1}, PlotRange -> All]


If we plot both of the initial condition and sol[0,x,y,z] in the direction of the line x=y=z, they are inconsistent. The sol[0,x,y,z] has some weird negative valleys around the peak and they are not exactly the same height either.

I wonder if this means that sol[0,x,y,z] doesn't represent the solution at the very first time point, or it does represent and the valleys are caused by numerical errors? If it is the latter, is it because of the Gaussian is so narrow (cause it is supposed to approximate a DiracDelta) and how to fix this?

Trying to solve this type of equations for a long time... Any help is appreciated!!

As @Tim Laska suggested, using the anisotropic mesh seems to be a promising approach. But a way to generate an anisotropic mesh centered at an arbitrary chosen point in the computational box is needed. I am a newbie in FEM so could anyone help to illustrate a way to do this?

Update: Another issue I discovered is that the conservation of the solution doesn't hold. As a test, I integrated the solution at t = 0, 0.5, 1, respectively, and I also integrated the RegularizedDeltaPoint used in the initial condition whose results are given below.

As we see, the error is about 25% which is huge. I wonder if anyone has an idea how to improve the conservation property? (Maybe increasing grid points helps, but it will require crazy amount of computation time... If the error is caused by the rough linear interpolation between grid points as seen in the plot, then I guess there might be a way to use higher order interpolation to interpolate the solution among those grid points to make it smoother? I tried to search for such way to do it but got none.)

• It seems to me unusual that DiracDelta appears inside initial condition. I would expect it at the right handside of the pde (Green's function) Jul 13 at 8:38
• It's based on some specific physical scenario where we want to find the probability distribution of a particle denoted by u at any time t with the particle confined at one spacial point r and its probability diffuse with time. Jul 13 at 11:32
• In this case I would expect a convolution with the diracdistribution... Jul 13 at 11:35
• Sorry, could you explain a bit more about how to do it by a convolution with diracdistribution? Jul 13 at 13:10
• The delta distribution can't be approximated by a continuous function, therefor I doubt your "initial conditions". One can indeed approximate the "dirac measure" ( Integral containing dirac) quite well. Jul 13 at 13:24

I will present an anisotropic meshing solution that will allow you to resolve a sharp normalized peak at the center of your cube. Many links applying this technique are in my answer to the question Future enhancements for the finite element method.

# Helper functions

I will use some of the following helper functions to construct an anisotropic mesh that highly resolves the cube's center.

(*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, 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]


# Mesh construction

The following code will create a 1D segment with a highly resolved mesh in the center that geometrically expands towards the endpoints. We will convert the 1D segment into a 3D hexahedral mesh using RegionProduct

(*Define some parameters*)
x1 = 1; nelm1 = 12;
(*Create mesh segment with finer mesh at x=0.5*)
seg1 = reflectRight@rightSegmentGrowth[x1/2, nelm1, x1/400];
Print["1D Discretization"]
r1 = pointsToMesh@seg1
regfull = RegionProduct[r1, r1, r1];
Print["Hexahedral mesh"]
HighlightMesh[regfull, 1]
(*Extract Coords from region*)
crd = MeshCoordinates[regfull];
(*Create element mesh*)
mesh = ToElementMesh[crd];
Print["Tetrahedral element mesh"]
mesh["Wireframe"]


# Solution

Before we attempt to solve we will introduce the concept of a regularized Delta function, RegularizedDeltaPoint, as described in the Heat transfer tutorial. This will provide a sharp peak without underflow errors in the asymptote.

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}}]
sol = NDSolveValue[{D[u[t, x, y, z], t] -
Laplacian[u[t, x, y, z], {x, y, z}] == 0,
PeriodicBoundaryCondition[u[t, x, y, z], x == 0,
TranslationTransform[{1, 0, 0}]],
PeriodicBoundaryCondition[u[t, x, y, z], y == 0,
TranslationTransform[{0, 1, 0}]],
PeriodicBoundaryCondition[u[t, x, y, z], z == 0,
TranslationTransform[{0, 0, 1}]],
u[0, x, y, z] ==
RegularizedDeltaPoint[0.01, {x, y, z}, {0.5, 0.5, 0.5}]},
u, {t, 0, 1}, {x, y, z} ∈ mesh];

Plot[{RegularizedDeltaPoint[0.01, {x, x, x}, {0.5, 0.5, 0.5}],
sol[0, x, x, x]}, {x, 0.45, 0.55}, PlotRange -> All]


As you can see, the solution matches the RegularizedDeltaPoint to the level of discretization quite well.

# Update in response to a comment about off-center refinement

You can use extendMesh to glue as many refined 1D segments together as you would like. Here is an example workflow to build a 1D segment with refinement at one-third of one of the dimensions of the unit cube.

nelm = 24;
xfrac = 1/3;
nelml = xfrac nelm // Round;
nelmr = nelm - nelml;
felm = 1/400;
segl = rightSegmentGrowth[xfrac, nelml, felm];
Print["Left-hand side of the segment"]
pointsToMesh@segl
segr = leftSegmentGrowth[1 - xfrac, nelmr, felm];
Print["Right-hand side of the segment"]
pointsToMesh@segr
Print["Full extended segment"]
pointsToMesh@extendMesh[segl, segr]


Now, we can wrap the workflow into a module like so:

Clear[fracUnitSegment]
fracUnitSegment[xfrac_, nelm_, fElm_] :=
Module[{nelml, nelmr, segl, segr},
nelml = xfrac nelm // Round;
nelmr = nelm - nelml;
segl = rightSegmentGrowth[xfrac, nelml, fElm];
segr = leftSegmentGrowth[1 - xfrac, nelmr, fElm];
pointsToMesh@extendMesh[segl, segr]]


Now we can use fracUnitSegment to create an off-center refined mesh:

rx = fracUnitSegment[1/4, 24, 1/400]
ry = fracUnitSegment[1/3, 24, 1/400]
rz = fracUnitSegment[3/4, 24, 1/400]
rp = RegionProduct[rx, ry, rz];
Print["Hexahedral mesh"]
HighlightMesh[rp, 1]
(*Extract Coords from region*)
crd = MeshCoordinates[rp];
(*Create element mesh*)
mesh = ToElementMesh[crd];
Print["Tetrahedral element mesh"]
mesh["Wireframe"]


# Update to comment about conservation

It will always be difficult to model a tiny feature. For example, if we use a regularization parameter of $$\gamma=0.01$$ and want to resolve the RegularizedDeltaPoint with 10 elements, then a uniform mesh on the unit cube would require 1 billion elements.

The anisotropic meshing approach will allow you to dramatically reduce the element count, but one still needs to conduct a mesh sensitivity analysis to determine the trade-offs between accuracy and simulation cost.

The following mesh may be a good compromise. It contains 27,000 elements, shows about a 6.5% error, and good conservation at different time steps.

(*Define some parameters*)
xf = 1/2; nelm = 32; gamma = 0.01;
(*Create mesh segment with finer mesh at x=0.5*)
rx = fracUnitSegment[xf, nelm, gamma/5]
rp = RegionProduct[rx, rx, rx];
Print["Hexahedral mesh"]
HighlightMesh[rp, 1]
(*Extract Coords from region*)
crd = MeshCoordinates[rp];
(*Create element mesh*)
(*grab hexa element incidents RegionProduct mesh*)
inc = Delete[0] /@ MeshCells[rp, 3];
mesh = ToElementMesh["Coordinates" -> crd,
"MeshElements" -> {HexahedronElement[inc]}];
Print["Hexahedral element mesh"]
mesh["Wireframe"]
sol = NDSolveValue[{D[u[t, x, y, z], t] -
Laplacian[u[t, x, y, z], {x, y, z}] == 0,
PeriodicBoundaryCondition[u[t, x, y, z], x == 0,
TranslationTransform[{1, 0, 0}]],
PeriodicBoundaryCondition[u[t, x, y, z], y == 0,
TranslationTransform[{0, 1, 0}]],
PeriodicBoundaryCondition[u[t, x, y, z], z == 0,
TranslationTransform[{0, 0, 1}]],
u[0, x, y, z] ==
RegularizedDeltaPoint[gamma, {x, y, z}, {0.5, 0.5, 0.5}]},
u, {t, 0, 1}, {x, y, z} ∈ mesh];

Plot[{RegularizedDeltaPoint[gamma, {x, x, x}, {0.5, 0.5, 0.5}],
sol[0, x, x, x]}, {x, 0.5 - 3 gamma, 0.5 + 3 gamma},
PlotRange -> All]


Note that there is some sensitivity to the integration method used. In particular, using the mesh specification, the initial time point appears to be off. So you may be better off split simply defining the domain in NIntegrate.

NIntegrate[sol[0, x, y, z], {x, 0, 1}, {y, 0, 1}, {z, 0, 1}]
NIntegrate[sol[1/2, x, y, z], {x, 0, 1}, {y, 0, 1}, {z, 0, 1}]
NIntegrate[sol[1, x, y, z], {x, 0, 1}, {y, 0, 1}, {z, 0, 1}]
NIntegrate[sol[0, x, y, z], {x, y, z} ∈ mesh]
NIntegrate[sol[1/2, x, y, z], {x, y, z} ∈ mesh]
NIntegrate[sol[1, x, y, z], {x, y, z} ∈ mesh]
(* 1.0644 *)
(* 1.06439 *)
(* 1.06439 *)
(* 1.42413 *)
(* 1.06439 *)
(* 1.06439 *)

• Thanks! This seems to be a promising way to tackle the problem, which I will definitely try. However, could you show how to extend it to the case where the peak of initial Delta Function is at arbitrary position in the box instead of just at the center (so that the mesh is the most condensed at an arbitrarily chosen position)? Jul 14 at 13:31
• @Dennis, you are welcome! See updated answer. Jul 14 at 18:08
• Thanks for the feedback! It's very helpful! I will try to solve the equation and see if I can get expected behavior out of it. Jul 14 at 18:55
• Hey Tim Laska: I tried out the method you offered and it turned out that there is an issue with the conservation of the solution. I described it in details as an update in the question. Do you have an idea for it? Jul 15 at 19:49