Add subtimestepping for the internal tracer calculation#3014
Conversation
The subtimestepping approach clearly works. An experiment based on backwater with different basins show there's no lag in terms of flow/level, but there is a significant one for tracer if you go to more basins. Subtimestepping nicely fixes this. |
| cumulative_in .= vertical_flux.drainage * dt | ||
| cumulative_in .+= vertical_flux.surface_runoff * dt | ||
| # Determine number of substeps needed based on a fraction of the residence time | ||
| safe = x -> isinf(x) ? 0 : x |
There was a problem hiding this comment.
Looks like this type unstable anonymous function should be made stable, and defined outside the callback.
| cumulative_in .+= vertical_flux.surface_runoff * dt | ||
| # Determine number of substeps needed based on a fraction of the residence time | ||
| safe = x -> isinf(x) ? 0 : x | ||
| @. nsubsteps = 2^clamp(floor(Int, safe(log2(dt / (concentration_state[:, Substance.ResidenceTime] * substep_ratio)))), 0, substep_depth) |
| # Debug: check continuity tracer after dilution step | ||
| c_cont = concentration_state[node_id.idx, Substance.Continuity] | ||
| if !iszero(storage_only_in) && abs(c_cont - 1.0) > 0.01 | ||
| @warn "Continuity drift after dilution" node_id substep c_cont storage_only_in cumulative_in = cumulative_in[node_id.idx] s_start mass_cont = mass[node_id.idx][Substance.Continuity] |
| # Take care of infinitely small masses, possibly becoming negative due to truncation. | ||
| for I in eachindex(mass_node) | ||
| if (-eps(Float64)) < mass_node[I] < (eps(Float64)) | ||
| mass_node[I] = 0.0 | ||
| end | ||
| end |
There was a problem hiding this comment.
This was already there, but I wonder if this will bite us when we start simulating substances that only appear in trace amounts. eps(Float64) == eps(1.0) == 2.220446049250313e-16 is not that small for a Float64, and closer to 0 the eps is much smaller, in the limit eps(0.0) == 5.0e-324.
There was a problem hiding this comment.
We're talking about mass here, and waterquality tends to be about microgram, or even milligram / L. So the limit for a basin of 1 m3 will already be in the -19, and for 1000 m3 it becomes -22, close to the order of a single molecule. Anyway, these quantities are dwarfed by the uncertainty of the calculation.
|
Have you measured performance impact? A possibility for 1024 sub time steps sounds impactfull. It may be worth while to try and use the analytical solution to: |
3 participants
Fixes #2438
concentrationto config (moving a solver option)max_subtimestepsto solver stats outputDraft for now, will remove the debug statements and do some validation, ideally using the HWS model to see a different response on the Rhine. I'm also considering to rename concentration to tracer in at least the config.
The implementation follows Delwaq scheme 24 with subtimesteps (also used for the new prototype). Based on the mean age of the water from the previous timestep in a basin (called ResidenceTime currently, half correct), the dt and a safety factor, we determine a number of required substeps. That is, the mean age should be less than
dt * substep_ratio(default of0.2) , otherwise we create subtimesteps. The number of subtimesteps can vary per basin, and come in powers of two, up to2^substep_depth(default of10). So the number of subtimesteps can vary between1and1024.We take the maximum number of timesteps, and use that as a number of iterations, only triggering the basins that require that iteration (basin with subtimestep of 1 will trigger once at the end, subts of 2 halfway and at the end, and a subts of 1024 will trigger each iteration). By using powers of two all the basins synchronize nicely without odd fractions, bookkeeping, while also enabling a nice depth (range from a day to a minute or two).