Update coupler stepping logic so that fluxes are exchanged as often as possible

[daverumph]

Purpose

Update coupler stepping logic so that fluxes are exchanged as often as possible
Closes #1756

Content

Before this PR, the coupler time step (dt) was constrained to be the same as the dt of
the largest component model. Consequently, fluxes were exchanged less often than
optimal when one model (e.g. ocean) stepped much more slowly than others. With
this change, the coupler can, and should, step at the rate of the fastest surface model.

• I have read and checked the items on the review checklist.

[juliasloan25]

cc @kmdeck @szy

With this change, when we run with dt_land < dt_ocean (for example), the coupler will re-compute fluxes and send them to the atmosphere at dt_land frequency. This should be more correct since the land model computes implicit fluxes within its step - now the atmosphere will see updated fluxes at the same time the land does.

Note that we'll still recompute fluxes over ocean and sea ice at every coupling timestep (dt_cpl). Since the ocean/sea ice states haven't changed, this is extra computation (and maybe introduces the same issue we previously had for the land surface, since the atmosphere will get updated fluxes when the ocean/sea ice don't). Maybe we want to do something like: when an ImplicitFluxSimulation steps, update the coupler and the atmosphere with the fluxes over that component, but don't recompute fluxes over other surfaces until we hit the end of an ExplicitFluxSimulation timestep

[juliasloan25]

We should add an assertion that: if an ImplicitFluxSimulation is present in cs.model_sims, dt_cpl == the dt of this simulation. If none is present, we want to set dt_cpl = max(smallest surface dt, atmos dt) - I think this is the generalization of how we've been trying to set dt_cpl.

[juliasloan25]

I'm a little concerned about the Float subtraction here. Instead of computing Δt as a float, could we instead compute the number of steps as an integer and use that to take the correct number of steps? And same for the other methods too

function Interfacer.step!(sim::ClimaLandSimulation, t::Float64)
    model_t = Float64(sim.integrator.t)
    model_dt = Float64(sim.integrator.dt)
    n_steps = round(Int, (t - model_t) / model_dt)
    for _ in 1:n_steps
        Interfacer.step!(sim.integrator)
    end
end

@ph-kev do you agree this should avoid the catastrophic cancellation?

[ph-kev]

I could be missing something, but aren't t and model_t relatively the same scale? I still think you would get catastrophic cancellation in this case. I haven't looked too deeply at what this function is doing, but Δt is only used for comparing against other values. For example, isapprox(Δt, model_dt) is the same as isapprox(t - model_t, model_dt) which is the same as isapprox(t, 2 * model_dt) which avoids the subtraction.

[juliasloan25]

The rounding should counteract the catastrophic cancellation - e.g. if (t - model_d) / model_dt gives 0.999997 instead of 1.0 (or even 0.9), it will get rounded to 1.

I think it's better to use concrete number of steps here rather than any isapprox (though I agree what you're suggesting is better), so that we reduce the uncertainty in how many steps will be taken.

[daverumph]

To avoid the catastrophic situation, use ITime/DateTime. That's why we're putting in support for it.
Is there a consensus on what we should do here?

[daverumph]

I don't see a viable way around this without changing the call signature. The problem is we're handed t, which is exactly what is getting larger and larger, and needs to be compared to the model_t to determine if it's time to step the model. Alternatively, we could have the coupler figure out when to step each model strictly on the basis of the ratio of dt's, but that would require an interface change I think. I suppose step! could ignore the t argument and just keep internal to the model and use something like model_dt_since_last_step, which would keep resetting when the model actually stepped. Each model would need to know dt_cpl though to know how much to advance model_dt_since_last_step.

I continue to think that using Float64 clocks is inherently flawed for long runs and we already have the solution of using ITime/DateTime for that case, so we shouldn't spend more time trying to make the Float64 case better. We should just document that it should not be used for long simulations.

[ph-kev]

I thought Julia's approach was reasonable by computing the number of steps to take to help with catastrophic cancellation. Is there any blockers with that approach? I agreed with for long runs, you would want to use ITime, but I don't think the Float64 case be neglected.

[daverumph]

I don't think Julia's method works correctly. The land and coupler dts are expected to be the same, so the problem doesn't surface, but I'm assuming we want to make the same change for other component models, where t increments are smaller than the model_dt. Take the ocean model as an example, with model_dt == 1800.0, and dt_cpl == 360.0.
With model_t == 0.0, i.e. just getting started, step! will be called sequentially with t == 360.0, 720.0, 1080.0, 1440.0, and finally 1800.0 as the coupler runs the faster-stepping land and atmos models. But the round() will cause the ocean model to step at t == 1080.0, because it starts rounding up at the halfway point, 900.0, not at t == 1800.0 as we want.