Adaptive time-stepping

[aedolfi]

Due to the non-disipative nature of the shallow water equation, the time scale soley depends on input parameters containing a dimension of time ($\gamma$ and $\mu$), but not on the relaxation time $\tau$. That allows for cheap addaptive time-stepping like so:

Write a new presssure and friction force function

function filmpressure_dt!(state::LBM_state_1D, sys::SysConst_1D, dt::Float64)
    # println("This is the pressure I am using!")
    hip, him = viewneighbors_1D(state.dgrad)
    # Straight elements j+1, i+1, i-1, j-1
    circshift!(hip, state.height, 1)
    circshift!(him, state.height, -1)

        if sys.n==9 && sys.m==3
                state.pressure .= -sys.γ*dt*dt .* (1 .- cospi.(sys.θ)) .* (sys.n - 1) .* (sys.m - 1) ./ ((sys.n - sys.m) * sys.hmin) .* Swalbe.fast_disj_93.(sys.hmin./(state.height .+ sys.hcrit))
        elseif sys.n==3 && sys.m==2
                state.pressure .= -sys.γ*dt*dt .* (1 .- cospi.(sys.θ)) .* (sys.n - 1) .* (sys.m - 1) ./ ((sys.n - sys.m) * sys.hmin) .* Swalbe.fast_disj_32.(sys.hmin./(state.height .+ sys.hcrit))
        else
                throw(DomainError((sys.n,sys.m), "This disjoining pressure is not implemented, Options currently are (n,m)=(9,3) or (n,m)=(3,2). Use those or implement a new option."))
        end
        state.pressure .-= sys.γ*dt*dt .* (hip .- 2 .* state.height .+ him)
    return nothing
end

and

function slippage_dt!(state::LBM_state_1D, sys::SysConst_1D, dt::Float64)
    state.slip .= (6*sys.μ*dt .* state.height .* state.vel) ./ (2 .* state.height.^2 .+ 6*sys.δ .* state.height .+ 3*sys.δ^2 )
    return nothing
end

With those function a minimalistic siumlation run has to look like this:

using Pkg

Pkg.activate(".")

using Swalbe


function run(;
    h=2,
    eps=0.001,
    Tmax=Int(1e8),
    dumps=100,
    tdump=Int(floor(Tmax/dumps)),
    L=2^8,
    delta=0.0,
    n=9, m=3,
    gamma=0.01,
    theta=1/9,
    data=dataall,
    plot=true,
    run=true,
    hmin=0.2,
    mu=1/6,
    dt_init=1.0,
    F_max=1e-4,
    relax_F=1e-3,
    dt_max=100
    )

    #Setup system
    sys = Swalbe.SysConst_1D(L=L, Tmax=Tmax, tdump=tdump, γ=gamma, δ=delta, n=n, m=m, hmin=hmin, θ=theta,  μ=mu)
    state = Swalbe.Sys(sys)
    state.height .= Swalbe.randominterface!(state.height, h, eps, L/2-L/10)
    println("Starting LBM time loop")
    t=0.0
    dt=dt_init
    i=0
    while t<=sys.Tmax
        if (t % sys.tdump + dt) >= sys.tdump
            #Store data
           
            #Plot data
            if plot
             
            end
            #status report
            println("Time step t=$(round(t,sigdigits=3)), mass m=$(round(sum(state.height),sigdigits=5)), diff d=$(round(maximum(state.height)-minimum(state.height), sigdigits=7)),dt=$(round(dt, sigdigits=3))")
            i+=1
        end
        if run
            Swalbe.filmpressure_dt!(state, sys, dt)
            Swalbe.h∇p!(state)
            Swalbe.slippage_dt!(state,sys, dt)
            state.F .= -state.h∇p  .-  state.slip
            Swalbe.equilibrium!(state)
            Swalbe.BGKandStream!(state)
            Swalbe.moments!(state)
            dt = min(dt*(1+relax_F*(F_max - abs(maximum(state.F)))/F_max),dt_max)
        end
        t+=dt
    end
end


run()

The idea is that numerical instability always comes from to high forcing, while to small forcing means that the code could be run faster. Thus setting a optimal force, that is numerically stable but not to small, and going up or down in the time step when the force becomes to small or to high will let the code run with the biggest possible time step without leading to too high forces. Then we have to multiply $\gamma$ by $dt^2$ and $\mu$ by $dt$ for dimensional reasons.

Tests

The insight is pretty new so I havn't tested much, but equilibrium properties of films should not be effected by $\gamma$ and $\mu$. Concerning time evolution the LSA of spinoidal rupture is fulfilled on small wavenumbers, there is just a lot more noise on the high wavenumbers. That was to be expected as we are permanently running at boaderline high forcing:

top left: film height, top right: FFT of filmheight vs LSA, bottom left: friction force, bottom right: pressure

After Film rupture of course the LSA is not expected to fit the data anymore, but we can see how the firction force is kept at a level that is close to creating numerical problems but is still tolerable

Time steps

Here is a plot of $dt$ indicating a speed up of factor almost 15 before film-rupture and 1.5 after film rupture for the parameters I have chosen

F_max

It seems like $F_{\max}=10^{-4}$ seems to be the most that is feasible. Above the code crashes or produces unphysical output. But playing arround with paramaters one can try to figure out what works and what doesn't.

[aedolfi]

@Zitzeronion

[Zitzeronion]

That looks like a very powerful feature for future versions.

Of course, we have to make sure that it works as intended. Would you mind (@aedolfi) compiling a PR with a short section of this in the docs?

Thanks!

[aedolfi]

I put it on the agenda. But I will take a while to find time.