2D nonlinear PDE describing leveling process of fluid

Question

I already asked a question regarding how to solve a nonlinear pde in mathematica which was answered nicely.

https://mathematica.stackexchange.com/questions/134145/solving-a-non-linear-pde

Actually this was a 1-dimensional form of a general 2D problem that I was trying to solve with matlab, but first I wanted to have a clue about the solution using mathematica. Now I think I have to restrict myself into solving it in mathematica. Here is the 2D form of the problem with its initial and boundary conditions:

$$\frac{\partial h}{\partial t} = -0.01 \bigg(\frac{\partial}{\partial x}\Big(h^3 \big(\frac{\partial^3h}{\partial x^3}+\frac{\partial ^3h}{\partial x\partial y^2}\big)\Big)+ \frac{\partial }{\partial y}\Big(h^3 \big(\frac{\partial ^3h}{\partial y^3}+\frac{\partial ^3h}{\partial y\partial x^2}\big)\Big)\bigg)$$ $$\frac{\partial h}{\partial x}=0,\ \ \frac{\partial^3h}{\partial x^3}+\frac{\partial ^3h}{\partial x\partial y^2}=0\ \ \text{when}\ \ x=\pm1$$ $$\frac{\partial h}{\partial y}=0,\ \ \frac{\partial ^3h}{\partial y^3}+\frac{\partial ^3h}{\partial y\partial x^2}=0\ \ \text{when}\ \ y=\pm1$$ $$h(0,x,y)=1+\cos(\pi x) \cos(\pi y)$$ I tried to do the same as the 1D case, however it gives me an error with the boundary conditions. I tried this:

t0 = 6;
BCLx1 = Derivative[0, 1, 0][u][t, -1, y];

BCRx1 = Derivative[0, 1, 0][u][t, 1, y];

BCLx3 = Derivative[0, 3, 0][u][t, x, -1]+Derivative[0, 1, 2][u][t, -1, y];

BCRx3 = Derivative[0, 3, 0][u][t, x, 1]+Derivative[0, 1, 2][u][t, 1, y];

BCLy1 = Derivative[0, 0, 1][u][t, x, -1];

BCRy1 = Derivative[0, 0, 1][u][t, x, 1];

BCLy3 = Derivative[0, 0, 3][u][t, x, -1]+Derivative[0, 2, 1][u][t, x, -1];

BCRy3 = Derivative[0, 0, 3][u][t, x, 1]+Derivative[0, 2, 1][u][t, x, 1];

sol = NDSolve[{D[u[t, x, y], 
 t] == -0.0192*
               (D[u[t, x, y]^3*(D[u[t, x, y], {x, 3}] + 
                D[D[u[t, x, y], {y, 2}], x]), x] + 
                D[u[t, x, y]^3*(D[u[t, x, y], {y, 3}] + 
                D[D[u[t, x, y], {x, 2}], y]), y]), 
          u[0, x, y] == Cos[π y] Cos[π x] + 1,
          BCLx1 == 0, 
          BCRx1 == 0, 
          BCLx3 == 0, 
          BCRx3 == 0, 
          BCLy1 == 0, 
          BCRy1 == 0, 
          BCLy3 == 0, 
          BCRy3 == 0}, 
      u, {t, 0, t0}, {x, -1, 1}, {y, -1, 1}, Method -> {"MethodOfLines", 
"SpatialDiscretization" -> {"TensorProductGrid", "MinPoints" -> 100}}, 
PrecisionGoal -> 2]; Plot3D[{u[t0, x, y] /. sol}, {y, -1, 1}, {x, -1, 1}, PlotRange -> All]

It states that the boundary conditions include non-normal derivatives! but they are all meaningful BCs based on the physics of the problem.

Any help in this regard is highly appreciated.

Answer 1

OK, even after the incorrect translation for BCRx3 and BCLx3 is amended, NDSolve can't handle this kind of boundary conditions, so let's discretize spatial dimensions outside of NDSolve, with the help of pdetoode:

t0 = 6;
BCLx1 = Derivative[0, 1, 0][u][t, -1, y];
BCRx1 = Derivative[0, 1, 0][u][t, 1, y];
BCLx3 = Derivative[0, 3, 0][u][t, -1, y] + Derivative[0, 1, 2][u][t, -1, y];
BCRx3 = Derivative[0, 3, 0][u][t, 1, y] + Derivative[0, 1, 2][u][t, 1, y];
BCLy1 = Derivative[0, 0, 1][u][t, x, -1];
BCRy1 = Derivative[0, 0, 1][u][t, x, 1];
BCLy3 = Derivative[0, 0, 3][u][t, x, -1] + Derivative[0, 2, 1][u][t, x, -1];
BCRy3 = Derivative[0, 0, 3][u][t, x, 1] + Derivative[0, 2, 1][u][t, x, 1];
eqn = D[u[t, x, y], 
    t] == -0.0192*(D[u[t, x, y]^3*(D[u[t, x, y], {x, 3}] + D[D[u[t, x, y], {y, 2}], x]), 
       x] + D[u[t, x, y]^3*(D[u[t, x, y], {y, 3}] + D[D[u[t, x, y], {x, 2}], y]), y]);
ic = u[0, x, y] == Cos[π y] Cos[π x] + 1;
bc = {BCLx1 == 0, BCRx1 == 0, BCLx3 == 0, BCRx3 == 0, BCLy1 == 0, BCRy1 == 0, BCLy3 == 0,
    BCRy3 == 0};

xdomain = ydomain = {-1, 1};
xpoints = ypoints = 30;
xgrid = Array[# &, xpoints, xdomain];
ygrid = Array[# &, ypoints, ydomain];
difforder = 4;
var = Outer[u, xgrid, ygrid];
(* Definition of pdetoode isn't included in this code piece, 
   please find it in the link above. *)
ptoo = pdetoode[u[t, x, y], t, {xgrid, ygrid}, difforder];

del = Delete[#, {{1}, {2}, {-2}, {-1}}] &;
ode = del /@ del@ptoo@eqn;
odeic = del /@ del@ptoo@ic;
odebc = MapAt[del, ptoo@bc, List /@ Range@4];

sollst = NDSolveValue[{ode, odeic, odebc}, var, {t, 0, t0}, 
                      SolveDelayed -> True]; // AbsoluteTiming
(* {8.970609, Null} *)

sol = rebuild[sollst, {xgrid, ygrid}, -1];

Animate[
 Plot3D[sol[x, y, t], {x, xdomain} // Flatten // Evaluate, {y, ydomain} // Flatten // 
   Evaluate, PlotRange -> {-1, 2}, PerformanceGoal -> "Quality"], {t, 0, 1/2}]

SolveDelayed -> True option is necessary at least for v10.2, its color is red, but don't worry.