Use DSS for land aux state

[juliasloan25]

We're seeing problems with water conservation in the coupler when snow is included in a simulation. We see that water is conserved when input evaporation is a uniform field, but large errors arise when evaporation varies randomly across space. In the land model, DSS is currently applied to the prognostic variables, but not to the auxiliary state p. We think that the incorrect water content calculations may be due to this lack of DSS.

Part of #404

Update: As of 9/19/23, we think the water conservation errors arose because dss is being applied to the state vector via the step! call, but it is never applied to the change in water content on the RHS (diff_PE below). If we modify the MWE to call dss!(diff_PE) just before the test comparison, we see similarly small errors to those in the uniform case: (-2.842170943040401e-14, 5.684341886080802e-14). Note that these extrema values are identical whether we call dss!(diff_PE) before or after step!.

MWE

• The MWE below can be run from the coupler to recreate this issue: Results:

  # uniform:
  # julia> extrema( @. diff_PE + diff_WL) = (-5.684341886080802e-14, 2.842170943040401e-14)
  # random:
  # julia> extrema( @. diff_PE + diff_WL) = (-127.08222748551799, 122.94495463691425)

Script:

import SciMLBase: step!, ODEProblem, solve, SSPRK33, savevalues!, init
using ClimaCore: Fields, Spaces
using ClimaCoupler: TestHelper
import Thermodynamics as TD
using ClimaComms
using Dates

FT = Float64
Δt_cpl = FT(200)

space = TestHelper.create_space(FT)
include("components/land/bucket_init.jl")
include("components/land/bucket_utils.jl")

test_sim = bucket_init(
    FT,
    FT.((0, 432000)),
    "sphere",
    "map_static",
    ClimaComms.SingletonCommsContext(),
    ".";
    dt = Δt_cpl,
    space = space,
    saveat = FT(3600),
    area_fraction = ones(space),
    date_ref = DateTime("19790101", dateformat"yyyymmdd"),
    t_start = FT(0),
);
test_sim.integrator.u.bucket.W .= FT(0.482052) 
test_sim.integrator.u.bucket.Ws .= FT(0)
test_sim.integrator.u.bucket.σS .= FT(-4.45e-322) 
test_sim.integrator.u.bucket.T .= FT(275.0) 

test_sim.integrator.p.bucket.P_liq .= FT(0.0) 
test_sim.integrator.p.bucket.P_snow .= FT(0.0) 
uniform = ones(space)
random = zeros(space)
parent(random) .= rand(FT, (4,4,1,96))
#evap = test_sim.integrator.p.bucket.evaporation .= random
evap = test_sim.integrator.p.bucket.evaporation .= uniform
precip = test_sim.integrator.p.bucket.P_liq .+ test_sim.integrator.p.bucket.P_snow

test_sim.integrator.p.bucket.R_n .= FT(-186.145)
test_sim.integrator.p.bucket.T_sfc .= FT(305.55)
test_sim.integrator.p.bucket.q_sfc .= FT(0.0286255)
test_sim.integrator.p.bucket.ρ_sfc .= FT(1.20363)
test_sim.integrator.p.bucket.α_sfc .= FT(0.0749572)
test_sim.integrator.p.bucket.turbulent_energy_flux .= FT(161.888)

WL_0 = deepcopy(@. test_sim.integrator.u.bucket.W + test_sim.integrator.u.bucket.Ws + test_sim.integrator.u.bucket.σS)

step!(test_sim.integrator, Δt_cpl, true)

WL = @. test_sim.integrator.u.bucket.W + test_sim.integrator.u.bucket.Ws + test_sim.integrator.u.bucket.σS
diff_WL = @. (WL - WL_0 )
diff_PE = @. (precip + evap) * Δt_cpl # bucket evap needs to be m3/m2 /s

### this dss reduces the error in the random evaporation case
dss_buffer = Spaces.create_dss_buffer(ClimaCore.Fields.zeros(space))
Spaces.weighted_dss!(diff_PE, dss_buffer)

extrema( @. diff_PE + diff_WL)

Alternative solutions

• The Bucket model is column-wise (i.e. no horizontal exchanges), which also may contribute to incorrect results when running a simulation with non-uniform evaporation

[juliasloan25]

@LenkaNovak @kmdeck tagging for visibility

[juliasloan25]

@simonbyrne Does this potential explanation of the bug make sense to you? Are we meant to apply DSS to both the prognostic variables and aux/cache state?

[simonbyrne]

The requirement is that any tendencies are continuous: that is, at any two colocated points, you get the same result.

For horizontal tendencies, we do this by calling DSS, for vertical tendencies and sources we do this by ensuring all inputs are continuous (so the result will also be).

If you're filling a field with random values, there is no guarantee of continuity, so any tendency computed from this may not be. So I would suggest DSS-ing the field after calling rand to ensure it is continuous.

[juliasloan25]

Closing because we don't actually need to implement this