8
$\begingroup$

PROBLEM STATEMENT

Recently, I was trying to verify the solution of a set of thin film spreading equations given by this paper.

enter image description here

Where $$p1=-(h_{1xx}+\sigma h_{2xx})$$

$$p2=-\sigma h_{2xx}-\Pi[h2-h1]$$

$$\Pi[h]=\frac{k}{h_f}[(\frac{h_f}{h})^4-(\frac{h_f}{h})^3]$$

However, given the high order and non-linear nature of the equations, Mathematica keeps returning me enormous spatial errors even under extremely fine meshes of 1000 points and a difforder of 5 in the finite difference scheme.

NDSolveValue::eerr: Warning: scaled local spatial error estimate of 68548.44786079538` at t = 0.03579162478117444` in the direction of independent variable x is much greater than the prescribed error tolerance. Grid spacing with 1000 points may be too large to achieve the desired accuracy or precision. A singularity may have formed or a smaller grid spacing can be specified using the MaxStepSize or MinPoints method options.

The enormous error, in turn, leads to non-physical behaviours in the solution field, which hinders me from solving it under reasonable time-scales.

enter image description here

While increasing both the mesh size and the difference order modestly reduce the error, both step increase the solve time significantly. (it now takes more than 2 hours to solve for 0.1s of spreading time for 5000 points and a difference order of 11) Since I am planning to solve an even more complex version of this equation, is there any way to reduce the error of the system through alternative ways?

CODE

My current attempt is given in the code below.

(*Setting Known Functions*)

k = 1; hf = 0.01;
Π[h_] := k/hf ((hf/h)^4 - (hf/h)^3)

G = 1; μ = 1; λ = 5; hini = 5;

With[{h1 = h1[x, t], h2 = h2[x, t]},
  P1[x_, t_] := -D[h1, x, x] - G D[h2, x, x];
  P2[x_, t_] := -G D[h2, x, x] - Π[h2 - h1]];

(*Equations*)

With[{h1 = h1[x, t], h2 = h2[x, t], P1 = P1[x, t], P2 = P2[x, t]},
  Q = h1 (h1 - h2) D[P2, x] - D[P1, x]/2 h1^2 - 
    D[P2, x]/μ h1 (h1/2 - h2);
  eqn1 = {D[h1, t] == 
     D[(h1^3/3 D[P1, x] - h1^2/2 D[P2, x] (h1 - h2)), x]};
  eqn2 = {D[(h2 - h1), t] == 
     D[(D[P2, x]/μ (h2^3/3 + h1^3/6 - (h1^2 h2)/2) - Q (h2 - h1)),
       x]};
  
  bc = {{D[h1, x] == 0, D[h1, x, x, x] == 0, D[h2, x] == 0, 
      D[h2, x, x, x] == 0} /. 
     x -> -10, {D[h1, x] == 0, D[h1, x, x, x] == 0, D[h2, x] == 0, 
      D[h2, x, x, x] == 0} /. x -> 10};
  ic = {h1 == hini, 
     h2 == hf + hini + 
       1/3 (λ/π)^(1/2) E^(-λ x^2)} /. t -> 0];

(*Solving the equations*)
timer = 0.100;
Scheme[n : _Integer | {_Integer ..}, o_] := {"MethodOfLines", 
   "SpatialDiscretization" -> {"TensorProductGrid", "MaxPoints" -> n, 
     "MinPoints" -> n, "DifferenceOrder" -> o}};

time = 0; points = 5000; difforder = 11;

Monitor[{solh1, solh2} = 
  NDSolveValue[{eqn1, eqn2, ic, bc}, {h1, h2}, {t, 0, timer}, {x, -10,
     10}, Method -> Scheme[points, difforder], 
   StepMonitor :> (time = t)], time]
$\endgroup$
2
  • $\begingroup$ Your problem is x. $\endgroup$
    – John Manko
    Jul 7 at 15:41
  • $\begingroup$ @JohnManko Er… what do you mean by x? $\endgroup$
    – xzczd
    Jul 8 at 3:43

1 Answer 1

11
$\begingroup$

Yet another problem related to the hyperbolic conservation law. Based on the experience obtained in e.g. here, let's avoid symbolic expansion of D. I'll use pdetoode for the task.

k = 1; hf = 0.01;
Π[h_] := k/hf ((hf/h)^4 - (hf/h)^3)
G = 1; μ = 1; λ = 5; hini = 5;

With[{h1 = h1[x, t], h2 = h2[x, t]}, P1[x_, t_] := -D[h1, x, x] - G D[h2, x, x];
                                     P2[x_, t_] := -G D[h2, x, x] - Π[h2 - h1]];

With[{h1 = h1[x, t], h2 = h2[x, t], P1 = P1[x, t], P2 = P2[x, t], mid1 = mid1[x, t], 
   mid2 = mid2[x, t], midP2 = midP2[x, t], midP1 = midP1[x, t]}, 
  Q = h1 (h1 - h2) D[midP2, x] - D[midP1, x]/2 h1^2 - D[midP2, x]/μ h1 (h1/2 - h2);
  eqnaddP1 = midP1 == P1;
  eqnaddP2 = midP2 == P2;
  eqnadd1 = mid1 == (h1^3/3 D[midP1, x] - h1^2/2 D[midP2, x] (h1 - h2));
  eqnadd2 = mid2 == (D[midP2, x]/μ (h2^3/3 + h1^3/6 - (h1^2 h2)/2) - Q (h2 - h1));
  eqn1 = {D[h1, t] == D[mid1, x]};
  eqn2 = {D[(h2), t] == D[mid2 + mid1, x]};
  ic = {h1 == hini, h2 == hf + hini + 1/3 (λ/π)^(1/2) E^(-λ x^2)} /. t -> 0];

I've removed the original bc. I'll use periodic b.c.s instead in subsequent coding, because

  1. Considering the physics background, periodic b.c. should be a good enough approximation for boundary at infinity in this case.

  2. Periodic b.c. is convenient to set in pdetoode via last argument.

points = 100; difforder = 2;
domain = {-10, 10}; grid = Array[# &, points, domain];
(* Definition of pdetoode isn't included in this post,
   please find it in the link above. *)
ptoofunc = 
  pdetoode[{h1, h2, mid1, mid2, midP1, midP2}[x, t], t, grid, difforder, True];
odeadd = Map[ptoofunc, {eqnadd1, eqnadd2, eqnaddP1, eqnaddP2}, {2}];
ode = Block[{mid1, mid2, midP1, midP2}, Set @@@ odeadd; ptoofunc@{eqn1, eqn2}];
odeic = ptoofunc@ic;
var = Outer[#[#2] &, {h1, h2}, grid];
timer = 0.1;
sollst = NDSolveValue[{ode, odeic}, var, {t, 0, timer}]; // AbsoluteTiming
(* {1.04179, Null} *)
{h1func, h2func} = rebuild[#, grid, 2] & /@ sollst;

Manipulate[
 Plot[{h1func[x, t], h2func[x, t]}, {x, ##}, PlotRange -> {4.5, 5.5}] & @@ 
  domain, {t, 0, timer}]

enter image description here

$\endgroup$
2
  • $\begingroup$ Hello! I have recently been playing around with the code and is wondering if it is possible to allow some of the boundary conditions to be defined while letting others remain periodic? Thanks! $\endgroup$
    – FLP
    Jul 10 at 19:33
  • $\begingroup$ @flp Are you sure this is physically correct? Anyway, it's possible, but you cannot make use of the last argument of pdetoode in this case. You need to write down the periodic b.c. in traditional form e.g. u[a, t]==u[b, t], Derivative[1,0][u][a,t]==Derivative[1,0][u][b,t],… then use pdetoode to discretize it. Notice pdetoode cannot directly handle traditional periodic b.c. so you need to discretize it in a manner similar to odeadd above. $\endgroup$
    – xzczd
    Jul 11 at 2:16

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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