5
$\begingroup$

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

enter image description here

Where the kernel $K(r,r')$ is defined as such,

enter image description here

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:

enter image description here

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?

$\endgroup$
12
  • $\begingroup$ 1. It's not a good idea to create sock puppet in this site. 2. "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." It's because you've used MachinePrecision number to define lb, which causes error accumulation. Defining lb as e.g. 1/100 or 0.01`16 will fix the problem. $\endgroup$
    – xzczd
    Dec 17, 2022 at 3:31
  • $\begingroup$ 1. Apologies! Honestly didn't know that it is bad to post under multiple accounts, will note that in the future! (However, since this question is already posted I'll just leave this here to avoid confusion) 2. Thanks for the advice - I'll try that out! $\endgroup$
    – FLP
    Dec 17, 2022 at 4:02
  • $\begingroup$ You can merge the accounts: mathematica.stackexchange.com/help/merging-accounts 3. The Indeterminate generated by Int is causing problem in NDSolve`FiniteDifferenceDerivative. (Looks like a bug to me. ) You need to properly modify Int to avoid generating Indeterminate as output of Int to circumvent the issue. $\endgroup$
    – xzczd
    Dec 17, 2022 at 4:10
  • $\begingroup$ 4. There's a simple mistake in definition of EqnC. 5. Since you've used diffbc, you should not remove anything from the discretized i.c.. $\endgroup$
    – xzczd
    Dec 17, 2022 at 4:32
  • $\begingroup$ 4.,5. Thanks and edited! $\endgroup$
    – FLP
    Dec 17, 2022 at 4:36

1 Answer 1

3
$\begingroup$

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

enter image description here

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]

enter image description here

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

enter image description here

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

$\endgroup$
8
  • $\begingroup$ already read this link you suggested, I tried a piecewice IC, IC= {h[r, t] == If[r >= lb && r <= 1, 3/5 (1 - r^2) + hf, hf]} /. t -> 0, with an additional modification odIC = ptoofunc@unitStepExpand@IC; NDSolveValue gives an error <The number of constraints (198) (ics) is not equal to the total differential order of the system plus the number of discrete variables (100)>. (to continue) $\endgroup$
    – lxy
    Mar 20, 2023 at 14:42
  • $\begingroup$ @jsxs Check carefully about how I use unitStepExpand. If you still feel confused, observe the output of PiecewiseExpand@IC and think about what trouble will arise then. $\endgroup$
    – xzczd
    Mar 20, 2023 at 15:11
  • $\begingroup$ Oops, I should have extracted the i.c. in the form of an If function with IC[[1, 2]]. I have modified the error, but when executing odIC = ptoofunc@unitStepExpand@IC[[1, 2]], i got another error >"Part 1 of {} does not exist". I have also checked unitStepExpand@IC[[1, 2]], which looks correct. Thank you for the help! $\endgroup$
    – lxy
    Mar 21, 2023 at 1:26
  • $\begingroup$ @jsxs …It should be IC= {h[r, t] == unitStepExpand@If[r >= lb && r <= 1, 3/5 (1 - r^2) + hf, hf]} /. t -> 0. $\endgroup$
    – xzczd
    Mar 21, 2023 at 2:56
  • $\begingroup$ @jsxs by the way if you are still interested in this simulation I think the original code has some additional problems as well that was detected later on... Didn't really include it into the edits (the post is already long enough) but may post another answer in the near future to include all necessary edits. $\endgroup$
    – FLP
    Apr 6, 2023 at 7:13

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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