NDSolve refuses to initialise when solving an integro-differential equation with custom pdetoode solver scheme

Question

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:

1. $\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;
2. 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:

1. Is my approach to this problem computationally sound?
2. Is there any severe inefficiency in my approach that I can simplify for the problem?
3. 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?

Answer 1

Aside from those simple mistakes already corrected in the question, there're 2 issues worth elaborating a bit, so let me write an answer.

1

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.

This is because pdetoode internally uses PossibleZeroQ to detect whether an equation lies on the boundary of the computational domain, but:

With[{lb = 0.01, rb = 50, points = 100}, 
 Array[# &, points, {lb, rb}][[-1]] - rb]
(* -7.10543*10^-15 *)
% // PossibleZeroQ
(* False *)

As we can see, the numeric error caused by MachinePrecision number is so large that PossibleZeroQ doesn't think the 2 numbers equal. Actually this is a common problem when using MachinePrecision number for indexing, you may check the following posts for more info:

Strange behavior when defining a symbol/function pointwise

Strange behavior of MemberQ, Position

Position function not always returning an answer even with no apparent problems

Can Someone Please Explain Internal`$SameQTolerance?

So, by changing definition of lb to lb = 1/100 or lb = 0.01`16, and adding Method -> {EquationSimplification -> Solve} or Method -> {EquationSimplification -> Residual} or SolveDelayed -> True or Method -> {EquationSimplification -> MassMatrix} to NDSolveValue, the equation system will be solved without difficulty:

showStatus[status_]:=LinkWrite[$ParentLink,
  SetNotebookStatusLine[
    FrontEnd`EvaluationNotebook[],ToString[status]]];
clearStatus[]:=showStatus[""];
clearStatus[]
jianshi[t_]:=EvaluationMonitor:>showStatus["t = "<>ToString[CForm[t]]]

tend = 200;

sollst = 
   NDSolveValue[{ode, odIC, odBC}, var, {t, 0, tend}, jianshi[t], 
    Method -> {EquationSimplification -> MassMatrix}]; // AbsoluteTiming
(* {4.63116, Null} *)

solfunc = rebuild[sollst, grid];

Manipulate[
 Plot[solfunc[t, r], {r, lb, rb}, PlotRange -> {0, 20}], {t, 0, tend}]

2

The other issue, though already circumvented by the code in Edits, deserves to be mentioned in an answer so it can be easily searched because it's caused by a surprsing and long-standing bug of NDSolve`FiniteDifferenceDerivative.

As mentioned above, the original Int sometimes give Indeterminate as the output, but:

grid = Range@6/10;
values = 0.1 ReplacePart[h /@ grid, {1} -> Indeterminate];
NDSolve`FiniteDifferenceDerivative[1, grid, values]

As we can see, when Indeterminate involves in, NDSolve`FiniteDifferenceDerivative arbitrarily changes the precision of values, including the those inside h. What's more surprising:

grid = Range@6/10;
values = 0.1 ReplacePart[h /@ grid, {1} -> Indeterminate];
code := NDSolve`FiniteDifferenceDerivative[0, grid, values]
code
(* Start a new cell here. *)
code

Cell division influences the output! These behaviors definitely look like bugs. Already reported to WRI.