Divergence of free boundary in a liquid film problem

These days I’ve done a code to simulate a liquid film with a periodic free surface.
After some steps, instability is observed on the upstream of the peak on the free surface. And it will evolution until the blowing up of this simulation in the end.




I expect to observe a traveling wave solution.

Here is my script:
liquidfilm.edp (5.7 KB)

In the present script, curvature is computed by:

    for[i,v:dxh[]]{
        v=(hy[][(i+1)%sizeL]-hy[][(i-1+sizeL)%sizeL])/(2*hc);//centre scheme
    }
    //compute dh/dx on sides
    dxh[][0]=(hy[][1]-hy[][sizeL-2])/(2*hc);
    dxh[][sizeL-1]=dxh[][0];

    for[i,v:c[]]{
        v=(hy[][(i-1+sizeL)%sizeL]-2*hy[][i]+hy[][(i+1)%sizeL])/(hc^2);
    }
    c[][0]=(hy[][(sizeL-2)%sizeL]-2*hy[][0]+hy[][1])/(hc^2);
    c[][sizeL-1]=c[][0];
    for[i,v:c[]]{
        v*=1.0/((1+dxh[][i]^2)^(1.5));
    }

The mesh size on the border seems to be too large to resolve the boundary, when the flow is not smooth.
So I think mesh size Lx/nx=0.1 maybe too large, but it is expensive to do simulation on a fine mesh. Now I am trying to get a more smooth curvature, from which I do not need a fine mesh on the whole domain but only refine the free boundary.

Any comments are welcome.

Best regards,
Haocheng Weng