4

I have a 4th order nonlinear PDE with the following BCs and IC:

(* Boundary conditions *)   
(h^(1,0,0))[0,y,t]==0.0,
(h^(1,0,0))[L,y,t]==0.0,
(h^(0,1,0))[x,0,t]==0.0,
(h^(0,1,0))[x,L,t]==0.0,

(h^(3,0,0))[0,y,t]==0,
(h^(3,0,0))[L,y,t]==0,
(h^(0,3,0))[x,0,t]==0,
(h^(0,3,0))[x,L,t]==0,

(* Initial condition *)
h[x, y, 0] == 1 + (-0.05 Cos[2 π x/L] - 0.05 Sin[2 π x/L]) Cos[2 π y/L]

I realized that the first derivative and third derivative of my initial condition, h[x,y,t=0] wrt to x evaluated at x=0 and x=L are not 0 and hence pose an inconsistency (I get the ibcinc error).

I can't change my boundary conditions as the example problem here does because my boundary conditions represent static contact angles and they are either a zero or a number. And I don't want to change the ICs because this is the initial condition I have been using all throughout and I'd like consistency in comparison of my results.

What should I do about this? The entire code follows:

$HistoryLength = 0;
Needs["VectorAnalysis`"]
Needs["DifferentialEquations`InterpolatingFunctionAnatomy`"];
Clear[Eq0, EvapThickFilm, h, Bo, ε, K1, δ, Bi, m, r]
Eq0[h_, {Bo_, ε_, K1_, δ_, Bi_, m_, r_}] := D[h, t] + Div[-h^3 Bo Grad[h] + 
    h^3 Grad[Laplacian[h]] + (δ h^3)/(Bi h + K1)^3 Grad[h] + 
    m (h/(K1 + Bi h))^2 Grad[h]] + ε/(Bi h + K1) + 
    (r) D[D[(h^2/(K1 + Bi h)), x] h^3, x] == 0;
SetCoordinates[Cartesian[x, y, z]];
EvapThickFilm[Bo_, ε_, K1_, δ_, Bi_, m_, r_] := Eq0[h[x, y, t], {Bo, ε, K1, δ, Bi, m, r}];
TraditionalForm[EvapThickFilm[Bo, ε, K1, δ, Bi, m, r]];

L=2*92.389;TMax=2491.29*100;

Clear[Kvar];
Kvar[t_]:=Piecewise[{{1,t<=1},{2,t>1}}]
(*Ktemp=Array[0.001+0.001#^2&,13]*)
hSol=h/.NDSolve[{
(*S,G,E,K,D,VR,M*)
EvapThickFilm[0.003,0,1,0,1,0.025,0],
(*h[0,y,t]==h[L,y,t],
h[x,0,t]==h[x,L,t],*)

(* Boundary conditions *)   
Derivative[1,0,0][h][0,y,t]==0.0,
Derivative[1,0,0][h][L,y,t]==0.0,
Derivative[0,1,0][h][x,0,t]==0.0,
Derivative[0,1,0][h][x,L,t]==0.0,

Derivative[3,0,0][h][0,y,t]==0,
Derivative[3,0,0][h][L,y,t]==0,
Derivative[0,3,0][h][x,0,t]==0,
Derivative[0,3,0][h][x,L,t]==0,

(* Initial conditions *)
(*h[x,y,0] == 1.1+Cos[km x] Sin[2 km y] *)
h[x,y,0]==1+(-0.05 Cos[2 π x/L]-0.05 Sin[2 π x/L]) Cos[2 π y/L]
},
h,
{x, 0, L},
{y,0, L},
{t, 0, TMax},
Method->{"BDF","MaxDifferenceOrder"->1},MaxStepFraction->1/50
(*StepMonitor:> Print["time= ",t] *)
][[1]]

What I have tried so far:

  1. I tried using Piecewise but my Piecewise function needed more than 2 arguments so I failed.
  2. I tried changing the solution method to

    Method -> {"MethodOfLines", "SpatialDiscretization" -> {"TensorProductGrid", "MinPoints" -> 100}} This didn't yield anything either.

  3. I also tried setting "DifferentiateBoundaryConditions"->False, I still get the same error.

Is there a fix to this "inconsistent ic/bc" issue?

Update 5/31/2012

I changed the boundary conditions to reflect the initial condition (I changed the BCs to the first and third derivatives of the initial condition so as to make them consistent) but still, I get the same error.

dearN
  • 5,341
  • 3
  • 35
  • 74
  • 1
    your input is not valid M- code. Also, it were good if you could post the NDSolve command you are using. Concerning ICs and BCs, they need to be consistent otherwise the PDE can not be solved reliably. –  May 28 '12 at 15:23
  • 2
    Could you repost the input using plain text only? If h^(1,0,0)[...] is a derivative, use Derivative[1,0,0][h][...]. – Szabolcs May 28 '12 at 15:36
  • 2
    @ruebenko There seems to be a copying bug in the front end: evaluate Derivative[1, 0, 0][h][x, y, z], select the output, Copy As Input Text, paste back into the notebook, and you'll get something invalid. The culprits are those \* at the end of each Superscipt line. – Szabolcs May 28 '12 at 15:45
  • @Szabolcs Reposted the equation in plain text (you'll see it if scroll down to the end) – dearN May 28 '12 at 19:09
  • I have posted my entire code in plain text. – dearN May 28 '12 at 19:10
  • 1
    @DNA I've formatted your code, but I think you have a typo in your IC in the "plain text" code block. You have h[x,y,0]==h[x,y,0]==... (twice), whereas it is correct in the first paragraph. I suggest looking into it and correcting it – rm -rf May 28 '12 at 21:12
  • 4
    @DNA When posting large pieces of code it is a good practice to copy it from the pasted code here to a blank notebook and run it. You will find the errors much easier than we can – Dr. belisarius May 28 '12 at 21:39
  • Hi All. Will read your comments and post the outcomes soon. Thanks. @R.M, yes, I noticed the typo and corrected it but forgot to update it here! Just did! – dearN May 29 '12 at 02:20
  • @Szabolcs, filed as a bug. Thanks. –  May 29 '12 at 07:06
  • @DNA, the code still does not work. –  May 29 '12 at 07:09
  • @DNA I meant: please post code that can be pasted back into a notebook, and is syntactically correct. h^(1,0,0) is pretty (easy to read for a human), but the system won't be able to interpret it. I tried to fix up your code (see my edit). – Szabolcs May 29 '12 at 12:12
  • @Szabolcs I apologize. I could never understand how to paste code back and forth in mathematica and other applications... So did you copy the code as input text? I just tried that it seems the best bet... – dearN May 29 '12 at 17:55
  • @DNA I didn't copy as input text. Unfortunately it seems there's a bug in copying as input text (see my comment above about it). I had to enter the Derivative lines manually. – Szabolcs May 30 '12 at 10:19
  • @Szabolcs Check! – dearN May 30 '12 at 14:51
  • Updated question. – dearN May 30 '12 at 15:27
  • As to the ibcinc warning, you may want to read this post. In my opinion we can simply live with it, as long as we can make the inconsistent b.c. really play a role in the solution. – xzczd Sep 29 '16 at 08:31

1 Answers1

2

Your approach as of the 31st May update should get rid of the error. I did it by defining:

g[x_,y_]=1+(-0.05 Cos[2 π x/L]-0.05 Sin[2 π x/L]) Cos[2 π y/L]

and setting the boundary conditions like this:

Derivative[1,0,0][h][0,y,t] == Derivative[1,0][g][0,y]

etc...

Simon Woods
  • 84,945
  • 8
  • 175
  • 324
  • Well, that may get rid of the error. But however, as a side note, it deflates the physics of the problem but obviously, this isn't the forum for that. Thanks. – dearN Jun 01 '12 at 14:07