Is it possible to solve a differential equation with a user-defined variable mesh in NDSolve?

Question

For some differential equations, its solution may evolve to a cusp, for example, singular behavior. One may want to introduce additional mesh points near the cusp to accurately follow the solution throughout the formation of the local singularity.

General question: Can we define a variable mesh scheme in NDSolve to solve a nonlinear PDE about $h(x,t)$ with finite difference method (FDM), in which the number of mesh points near a cusp may be determined by the value of the local minimum or maximum, $h_{min}$ or $h_{max}$. Let us consider the case of a cusp having $h_{min}$. I noted that there are three additional processes should be added: (1) after each time step, a new set of mesh needs to be defined with additional points introduced; (2) the values of $h(x^{n+1}_j,t^{n+1})$ on the new mesh needs to be obtained from those $h(x^n_i,t^n)$ on the old mesh by interpolating; (3) the time step should also be adjusted according to increasing mesh resolution to ensure a consistent evolution in local mesh scales.

For example, consider the following PDE

$$h_t+(h^3h_{xxx}+h^{-1}h_x)_x=0,$$ which is defined on $x\in[0, L]$ with $L=8$, subject to the boundary conditions $h_x=h_{xxx}=0$ at the two end points of the interval. This equation will evolve a cusp at $x=0$ with a harmonic initial condition $h(x,0)=1-0.1\cos(2\pi x/L)$.

Pleas note that I can solve this equation using NDSolve with Method of Line and some common settings, which can be found easily on the site. However, to have more control on the numerical scheme and to obtain a more accurate solution especially near the singular point, I want to use the above mentioned variable mesh scheme, which begins with, say, $20$ points and increases to a much more dense mesh by $h_{min}=10^{-5}$. Specifically, is it possible to define a mesh in NDSolve with the following transformation $x(s)=h^2_{min}(t)\sinh(s)$, where $s(x,t)$ is the coordinate used in the numerical scheme. The mesh points are uniform in $s$ and thus the mesh spacing in $x$ is nearly uniform near $x=0$ but its spacing increases away from $x=0$. Please refer to the $\sinh$ curve.

After each time step, a new mesh is produced using the current value of $h_{min}(t^n)$ with additional mesh points introduced to keep the mesh spacing in $s$ roughly constant. The values of $h(x^{n+1}_j,t^{n+1})$ on the new mesh are obtained from those $h(x^n_i,t^n)$ on the old mesh by cubic interpolation, say. Also, the time step is chosen by

$$\Delta t=\Delta x_0^4/h^3_{min}\approx h_{min}^5 \Delta s^4,$$ where $\Delta x_0$ is the mesh spacing near $x=0$.

Answer 1

This is hard (if not impossible) to achieve inside NDSolve AFAIK. I sincerely hope I'm wrong. As mentioned in the tutorial The Numerical Method of Lines, particularly in Spatial Error Estimates section:

Ideally, the spatial error estimates would be monitored over time and the spatial mesh updated according to the evolution of the solution. The problem of grid adaptivity is difficult enough for a specific type of PDE and certainly has not been solved in any general way. Furthermore, techniques such as local refinement can be problematic with the method of lines since changing the number of mesh points requires a complete restart of the ODE methods. There are moving mesh techniques that appear promising for this approach, but at this point, NDSolve uses a static grid.

The tutorial above is about TensorProductGrid method. As to FiniteElement, as mentioned by user21 in this answer:

No, adaptive mesh refinement is not implemented in version 11.3 for PDE solving with the FEM. NIntegrate with the FEM uses an adaptive mesh refinement.

Additionally, your PDE is 4th order in $x$ direction, while currently FiniteElement method doesn't support PDE whose spatial derivative order exceeds 2.