Problem Background
Recently I'm attempting to replicate the result of the following research paper on the Nonlocal description of Evaporating Drops. The equation of motion of a evaporating, spreading drop is described as
Where the kernel $K(r,r')$ is defined as such,
My current approach
To deal with the coupled system of integro-differential equations, I've decided to discretise the system with the useful pdetoode
into a massive system of coupled ODEs. By expressing the integral as a sum across spatial coordinates, I can simplify the problem to a coupled ODE problem.
As mentioned by this post, the high order of the differentiation leads to problem with conservation, necessitating numerical calculation of the respective intermediates. The equations are therefore defined as such:
(*Defining Variables*)
a = 0.001;(*pseudo hanmaker's constant*)
hf = 0.01;(*precursor film thickness*)
μ = 1;(*viscosity*)
β = 1;(*Evaporation Constant*)
(*Defining Equations*)
With[{h = h[r, t], P = P[r, t], Q = Q[r, t], J = J[r, t]},
EqnP = P == D[h, r, r] + D[h, r]/r + a^2/h^3;
EqnQ = Q == 1/3 h^3 r D[P, r];
EqnC = D[h, t] + (μ D[Q, r])/r == -J;
EqnJ = J == (β D[Int[r, t], r])/r;]
(*Defining Initial and Boundary Conditions*)
IC = {h[r, t] == 20 E^-(r/20)^2 + hf} /. t -> 0
BC = {{D[h[r, t], r] == 0, D[h[r, t], r, r, r] == 0, D[h[x, t], x] == 0,
D[h[x, t], x, x, x] == 0} /. {r -> lb, x -> rb}}
The integral, $\int_{0}^{\infty} K(r,r')\frac{\partial}{\partial r'} (\frac{h}{hf})^3 \,dr'$ is completed as such:
- $\frac{\partial}{\partial r'} (\frac{h}{hf})^3$ is discretised into a list of differences of $h(r,t)$, and the integrand is integrated via the midpoint rule;
- For $r'\approx r$, there is an integrable singularity for $K(r,r')$. We approximate $K(r,r')$ as $\frac{r}{\pi}\ln|r'-r|$ and integrate it symbolically to find the integral for this region.
(*Discretising the Integral Function*)
lb = 1/100; rb = 50;
difforder = 2;
points = 100;
unitStepExpand = Simplify`PWToUnitStep@PiecewiseExpand@# &;
step = (rb - lb)/points;
grid = Array[# &, points, {lb, rb}];
hgrid = h[#][t] & /@ grid;
hdgrid = fdd[Derivative[1], grid, hgrid];
Kernel[r_, i_] := (
2 r)/π unitStepExpand@
If[r > i, r (EllipticK[i/r] - EllipticE[i/r]),
i (EllipticK[r/i] - EllipticE[r/i])];
IntegratedSingularity[r_, i_] :=
r/π Sqrt[(r - i)^2] (Log[Sqrt[(r - i)^2]] - 1)
Integrand =
D[(hf/h[r, t])^3, r] /. Thread[{D[h[r, t], r], h[r, t]} -> {hdgrid, hgrid}];
Int[r_][t_] :=
Total@Table[
unitStepExpand@
If[Sqrt[(r - lb - (i + 0.5) step)^2] <=
3 step, (IntegratedSingularity[r, lb + (i + 1) step] -
IntegratedSingularity[r, lb + i step])*0.5 (Integrand[[i]] +
Integrand[[i + 1]]),
step*Kernel[r,
lb + (i + 0.5) step]*0.5 (Integrand[[i]] +
Integrand[[i + 1]])], {i, 1, points - 1}];
The rest of the system is then discretised as such:
(*Discretising the rest of the system*)
ptoofunc =
pdetoode[{h[r, t], P[r, t], Q[r, t], J[r, t], Int[r, t]}, t, grid,
difforder];
removeredundant = #[[3 ;; -3]] &;
odeadd = Map[ptoofunc, {EqnP, EqnQ, EqnJ}, {2}];
ode = Block[{P, Q, J}, Set @@@ odeadd; ptoofunc@EqnC] //
removeredundant;
odIC = ptoofunc@IC;
(*For some weird reason if you apply ptoofunc to a boundary condition at a finite value ie /
D[h[r,t],r]/.r->50==0,it would discretise all values of D[h[r,t],r] prior to 50 (ie /
49,48,47,46...) as 0 as well. I have to manually choose the correct boundary conditions from /
the given list*)
(*Edit: Resolved with fractional lb instead of MachinePrecision*)
With[{sf = 1}, odBC = diffbc[t, sf]@BC // ptoofunc];
(*Solving the System*)
var = h /@ grid;
Monitor[sollst =
NDSolveValue[{ode, odIC, odBC}, var, {t, 0, 1},
EvaluationMonitor :> (time = t)], time]
Current Situation
At the NDSolveValue
step, however, the code refuses to initialise and the "time" under the EvaluatioMonitor
does not change into any numbers. Leaving NDSolveValue
to run on its own eventually causes the computer to lag severely and eventually crash. In fact, NDSolveValue
does not stop evaluating even if it encountered an error, and refuses to abort after running for a certain amount of time.
Personally, I'm not sure whether this is an issue with computational power, or that there are certain parts of my approach that is erroneous but I am personally unable to see. I would therefore like to ask the following questions:
- Is my approach to this problem computationally sound?
- Is there any severe inefficiency in my approach that I can simplify for the problem?
- Are there any other potential issues with my code?
Thank everyone in advance for any help you could offer!
Edits
At the advice of @xzczd, Ive looked through the integration function again and realised that the integral would take on indeterminate values whenever $r=r'$. This is because
$$\int_{a}^{b}\frac{r}{\pi}\ln|r'-r|\,dr'=\frac{r}{\pi}(r-b)(\ln|r-b|-1)-\frac{r}{\pi}(r-a)(\ln|r-a|-1)$$
Notice that whenever $a=r$ or $b=r$, $\ln|a-r|=-\infty$ resulting in the indefinite value as given.
In fact in our case, we don't even have to consider the case of the integrable singularity at all. Since we are taking the trapezoidal rule, $r$ can only take on the values lb+step*points
while $r'$ can only take on the values lb+(step+0.5)*points
, so the singularity does not occur in our entire integration. Using purely the trapezoidal rule to integrate $K(r,r')$, and comparing the results to NIntegrate
:
The fit is honestly not bad.
The integral function discretisation is now redefined as such:
(*Discretising the Integral Function*)
lb = 0.01; rb = 50;
difforder = 2;
points = 100;
unitStepExpand = Simplify`PWToUnitStep@PiecewiseExpand@# &;
step = (rb - lb)/points;
grid = Array[# &, points, {lb, rb}];
hgrid = h[#][t] & /@ grid;
hdgrid = fdd[Derivative[1], grid, hgrid];
Kernel[r_, i_] := (
2 r)/π unitStepExpand@
If[r > i, r (EllipticK[i/r] - EllipticE[i/r]),
i (EllipticK[r/i] - EllipticE[r/i])];
Integrand =
D[(hf/h[r, t])^3, r] /. Thread[{D[h[r, t], r], h[r, t]} -> {hdgrid, hgrid}];
Int[r_][t_] :=
Total@Table[
step*Kernel[r,
lb + (i + 0.5) step]*0.5 (Integrand[[i]] +
Integrand[[i + 1]]), {i, 1, points - 1}];
However, the computational problem persists. When the full solver is ran, the computer begins lagging severely and eventually crashes completely by switching off.
A cursory check via Task Manager indicates that the memory usage of NDSolve
reaches nearly 98% to 99% of all available memory prior to the crash, potentially explaining the lagging behaviour of the computer.
In that case, is there any method for reducing the memory usage of NDSolve
when computing the expression? Perhaps is it possible to conduct numerical preprocessing of the equations before plugging them into NDSolve
, or are there any other possible methods?
MachinePrecision
number to definelb
, which causes error accumulation. Defininglb
as e.g.1/100
or0.01`16
will fix the problem. $\endgroup$Indeterminate
generated byInt
is causing problem inNDSolve`FiniteDifferenceDerivative
. (Looks like a bug to me. ) You need to properly modifyInt
to avoid generatingIndeterminate
as output ofInt
to circumvent the issue. $\endgroup$EqnC
. 5. Since you've useddiffbc
, you should not remove anything from the discretized i.c.. $\endgroup$