Add a soil carbon decomposition module: disordered kinetics

experiments/standalone/Biogeochemistry/disordered_kinetics.jl

@@ -0,0 +1,90 @@
+import SciMLBase
+import ClimaComms
+ClimaComms.@import_required_backends
+import ClimaTimeSteppers as CTS
+using ClimaCore
+using ClimaLand
+using ClimaLand.Domains: Column
+using ClimaLand.Soil
+using ClimaLand.Soil.DisorderedKinetics
+using Dates
+
+import ClimaLand.Parameters as LP
+
+# Define simulation times
+t0 = Float64(0);
+tf = Float64(100);
+dt = Float64(0.01);
+
+FT = Float64;
+
+model_params = DisorderedKineticsModelParameters{FT}(
+    mu = FT(-1.0),
+    sigma = FT(2.0));
+zmax = FT(0);
+zmin = FT(-1);
+nelems = 1;
+# we use a vertical column with one layer to in the future interface with other variables like temperature and moisture that are depth dependent
+lsm_domain = Column(; zlim = (zmin, zmax), nelements = nelems);
+
+# Make biogeochemistry model args
+NPP = PrescribedSOCInputs{FT}(TimeVaryingInput((t) -> 500));
+model_sources = (LitterInput{FT}(),);
+
+model = DisorderedKineticsModel{FT}(;
+    parameters=model_params,
+    npools=100,
+    domain=lsm_domain,
+    sources=model_sources,
+    drivers=NPP
+);
+
+Y, p, coords = initialize(model);
+set_initial_cache! = make_set_initial_cache(model);
+
+
+
+
+function init_model!(Y, p, model)
+    N = NTuple{model.npools, FT}
+    function set_k(k::N) where {N}
+        k = collect(exp.(range(model.parameters.mu-8*model.parameters.sigma,model.parameters.mu+8*model.parameters.sigma,length=model.npools)));
+        return ntuple(x->k[x],model.npools)
+    end
+    p.soilC.ks .= set_k.(p.soilC.ks);
+end
+
+
+init_model!(Y, p, model);
+
+set_initial_cache!(p, Y, t0);
+disordered_kinetics_exp_tendency! = make_exp_tendency(model);
+
+timestepper = CTS.RK4();
+ode_algo = CTS.ExplicitAlgorithm(timestepper);
+
+saveat = collect(t0:FT(10 * dt):tf);
+sv = (;
+    t = Array{Float64}(undef, length(saveat)),
+    saveval = Array{NamedTuple}(undef, length(saveat)),
+);
+saving_cb = ClimaLand.NonInterpSavingCallback(sv, saveat);
+# updateat = deepcopy(saveat)
+# drivers = ClimaLand.get_drivers(model)
+# updatefunc = ClimaLand.make_update_drivers(drivers)
+# driver_cb = ClimaLand.DriverUpdateCallback(updateat, updatefunc)
+# cb = SciMLBase.CallbackSet(driver_cb, saving_cb)
+
+prob = SciMLBase.ODEProblem(
+    CTS.ClimaODEFunction(T_exp! = disordered_kinetics_exp_tendency!),
+    Y,
+    (t0, tf),
+    p,
+)
+# sol = SciMLBase.solve(prob, ode_algo; dt = dt, callback = cb)
+sol = SciMLBase.solve(prob, ode_algo; dt = dt, callback = saving_cb);
+
+# Check that simulation still has correct float type
+@assert eltype(sol.u[end].soilC) == FT;
+
+

experiments/yinon_script/code.jl

@@ -0,0 +1,67 @@
+# parameters: mu and sigma (mean and std distributions)
+# auxiliary: k is an array of decay rate (constant)
+# state variable: u (array of C with different decay rate)
+# input: input of carbon per area per unit time
+# we prescribe input for now, but it should be given by another
+# module prognostically, e.g., litter fall etc.
+
+using QuadGK
+using DifferentialEquations
+using StaticArrays
+using Distributions
+using DataFrames
+using CSV
+using BenchmarkTools
+# integral, error = quadgk(x -> cos(200x), 0, 1)
+
+# function dc_kt_dt(u, p, t)#, fixed=true)
+#     k,mu,sigma,input = p
+#     if fixed==true
+#         return SA[input*pdf(LogNormal(mu,sigma),k) - k*u]
+#     else
+#         ind = min(int(t+1),length(input))
+#         return SA[input[ind]*pdf(LogNormal(mu,sigma),k) - k*u]
+#     end
+
+
+# end
+
+function dc_kt_dt(u, p, t)#, fixed=true)
+    k,mu,sigma,input,fixed = p
+    if fixed== true
+        SA[input.*pdf(LogNormal(mu,sigma),k) .- k.*u[1]]
+    else
+        ind = min(floor(Int,t+1),length(input))
+        return SA[input[ind]*pdf(LogNormal(mu,sigma),k) - k*u[1]]
+    end
+
+    # SA[p[4].*pdf(LogNormal(p[2],p[3]),p[1]) .- p[1].*u[1]]
+end
+# prob = ODEProblem(dc_kt_dt, SA[0.], (0.,100.),[0.1,1.,2.,0.25,true])
+# solve(prob)
+
+
+function solve_ode(k,tau,age,input,fixed=true)
+    sigma = sqrt(log(age/tau))
+    mu = - log(sqrt(tau^3/age));
+    ps = (k,mu,sigma,input,fixed)
+    tspan = (0.,100.)
+    u0 = SA[0.]
+    prob = ODEProblem(dc_kt_dt, u0,tspan , ps)
+    solve(prob,saveat=1)
+end
+function run_diskin(tau,age,input,fixed=true,discrete=false)
+    if discrete==true
+        krange = exp.(range(-10,stop=10,length=100000));
+        dk = diff(vcat(0,krange));
+        res = sum(reduce(hcat,[solve_ode(k,tau,age,input,fixed).u.*d for (k,d) in zip(krange,dk)]),dims=2);
+    else
+        res, error = quadgk(x -> solve_ode(x,tau,age,input,fixed), 0, Inf)
+    end
+    res
+end
+NPP = Array(CSV.read("experiments/yinon_script/test_NPP.csv",DataFrame));
+
+cont = run_diskin(17,1000,0.25,true,false);
+@btime noncont = run_diskin(17,1000,0.25,true,true);
+cont'./noncont

src/shared_utilities/drivers.jl

@@ -18,6 +18,7 @@ export AbstractAtmosphericDrivers,
     PrescribedAtmosphere,
     PrescribedPrecipitation,
     PrescribedSoilOrganicCarbon,
+    PrescribedSOCInputs,
     CoupledAtmosphere,
     PrescribedRadiativeFluxes,
     CoupledRadiativeFluxes,
@@ -95,6 +96,23 @@ end
 PrescribedSoilOrganicCarbon{FT}(soc) where {FT} =
     PrescribedSoilOrganicCarbon{FT, typeof(soc)}(soc)
 
+
+"""
+     PrescribedSOCInputs{FT}
+
+A type for prescribing soil organic carbon inputs into our disordered kinetics model.
+$(DocStringExtensions.FIELDS)
+"""
+struct PrescribedSOCInputs{FT, SOC_input <: AbstractTimeVaryingInput} <:
+       AbstractClimaLandDrivers{FT}
+    "Soil organic carbon input, function of time and space: kg C/m^2/yr"
+    soc_input::SOC_input
+end
+
+PrescribedSOCInputs{FT}(soc_input) where {FT} =
+    PrescribedSoilOrganicCarbon{FT, typeof(soc_input)}(soc_input)
+
+
 """
     PrescribedAtmosphere{FT, CA, DT} <: AbstractAtmosphericDrivers{FT}