EarthSciML Example

Example

  1. Figure 1
  2. Figure 2
  3. Figure 3
  4. Supplemental code to skip unit enforcement
  5. Figure 4
  6. Supplemental plotting code

Figure 1

using EarthSciData, EarthSciMLBase, GasChem,
    DomainSets, ModelingToolkit, MethodOfLines, 
    DifferentialEquations, Dates, Distributions,
    Latexify, Plots, SciMLBase,
    CairoMakie, GraphMakie, MetaGraphsNext

@parameters t lev lon lat
model = SuperFast(t)
ls = latexify(model.rxn_sys) # TODO: Change to model.rxn_sys

This is what ls looks like:

O3+OH→eT1T2.46e10k1HO2+O2HO2+O3→eT2T2.46e10k22O2+OHHO2+OH→eT3T2.46e10k3HO2+O2NO+O3→eT4T2.46e10k4NO2+O2HO2+NO→eT5T2.46e10k5NO2+OHCH4+OH→eT6T2.46e10k6CH3O2+H2OCH2O+OH→eT7T2.46e10k7CO+H2O+HO2CH3O2+HO2→eT8T2.46e10k8CH3OOH+O2CH3OOH+OH→eT9T2.46e10k9CH3O2+H2OCH3OOH+OH→eT10T2.46e10k10CH3O+H2O+OHCH3O2+NO→eT11T2.46e10k11CH2O+HO2+NO220CH3O2→eT12T2.46e10k1220CH2O+8H2ODMS+OH→eT13T2.46e10k13SO2ISOP+OH→eT14T2.46e10k142CH3O22ISOP+2OH→eT15T2.46e10k152ISOP+OHISOP+O3→eT16T2.46e10k160.87CH2O+1.86CH3O2+0.06HO2+0.05COO3→1.0000000000000001e−21jO31DO1d+O2O1d+H2O→100j2OH2OHH2O2→jH2O22OHNO2→jNO2NO+O3CH2O→jCH2OaCO+2HO2CH3OOH→jCH3OOHCH2O+HO2+OH2HO2→eT17T2.46e10k17H2O2+O2OH+H2O2→2.46e10k18H2O+HO2OH+CO→2.46e10k19HO2\begin{aligned} \mathrm{O3} + \mathrm{OH} &\xrightarrow{e^{\frac{T1}{T}} 2.46e10 k1} \mathrm{HO2} + \mathrm{O2} \\ \mathrm{HO2} + \mathrm{O3} &\xrightarrow{e^{\frac{T2}{T}} 2.46e10 k2} 2 \mathrm{O2} + \mathrm{OH} \\ \mathrm{HO2} + \mathrm{OH} &\xrightarrow{e^{\frac{T3}{T}} 2.46e10 k3} \mathrm{HO2} + \mathrm{O2} \\ \mathrm{NO} + \mathrm{O3} &\xrightarrow{e^{\frac{T4}{T}} 2.46e10 k4} \mathrm{NO2} + \mathrm{O2} \\ \mathrm{HO2} + \mathrm{NO} &\xrightarrow{e^{\frac{T5}{T}} 2.46e10 k5} \mathrm{NO2} + \mathrm{OH} \\ \mathrm{CH4} + \mathrm{OH} &\xrightarrow{e^{\frac{T6}{T}} 2.46e10 k6} \mathrm{CH3O2} + \mathrm{H2O} \\ \mathrm{CH2O} + \mathrm{OH} &\xrightarrow{e^{\frac{T7}{T}} 2.46e10 k7} \mathrm{CO} + \mathrm{H2O} + \mathrm{HO2} \\ \mathrm{CH3O2} + \mathrm{HO2} &\xrightarrow{e^{\frac{T8}{T}} 2.46e10 k8} \mathrm{CH3OOH} + \mathrm{O2} \\ \mathrm{CH3OOH} + \mathrm{OH} &\xrightarrow{e^{\frac{T9}{T}} 2.46e10 k9} \mathrm{CH3O2} + \mathrm{H2O} \\ \mathrm{CH3OOH} + \mathrm{OH} &\xrightarrow{e^{\frac{T10}{T}} 2.46e10 k10} \mathrm{CH3O} + \mathrm{H2O} + \mathrm{OH} \\ \mathrm{CH3O2} + \mathrm{NO} &\xrightarrow{e^{\frac{T11}{T}} 2.46e10 k11} \mathrm{CH2O} + \mathrm{HO2} + \mathrm{NO2} \\ 20 \mathrm{CH3O2} &\xrightarrow{e^{\frac{T12}{T}} 2.46e10 k12} 20 \mathrm{CH2O} + 8 \mathrm{H2O} \\ \mathrm{DMS} + \mathrm{OH} &\xrightarrow{e^{\frac{T13}{T}} 2.46e10 k13} \mathrm{SO2} \\ \mathrm{ISOP} + \mathrm{OH} &\xrightarrow{e^{\frac{T14}{T}} 2.46e10 k14} 2 \mathrm{CH3O2} \\ 2 \mathrm{ISOP} + 2 \mathrm{OH} &\xrightarrow{e^{\frac{T15}{T}} 2.46e10 k15} 2 \mathrm{ISOP} + \mathrm{OH} \\ \mathrm{ISOP} + \mathrm{O3} &\xrightarrow{e^{\frac{T16}{T}} 2.46e10 k16} 0.87 \mathrm{CH2O} + 1.86 \mathrm{CH3O2} + 0.06 \mathrm{HO2} + 0.05 \mathrm{CO} \\ \mathrm{O3} &\xrightarrow{1.0000000000000001e-21 jO31D} \mathrm{O1d} + \mathrm{O2} \\ \mathrm{O1d} + \mathrm{H2O} &\xrightarrow{100 j2OH} 2 \mathrm{OH} \\ \mathrm{H2O2} &\xrightarrow{jH2O2} 2 \mathrm{OH} \\ \mathrm{NO2} &\xrightarrow{jNO2} \mathrm{NO} + \mathrm{O3} \\ \mathrm{CH2O} &\xrightarrow{jCH2Oa} \mathrm{CO} + 2 \mathrm{HO2} \\ \mathrm{CH3OOH} &\xrightarrow{jCH3OOH} \mathrm{CH2O} + \mathrm{HO2} + \mathrm{OH} \\ 2 \mathrm{HO2} &\xrightarrow{e^{\frac{T17}{T}} 2.46e10 k17} \mathrm{H2O2} + \mathrm{O2} \\ \mathrm{OH} + \mathrm{H2O2} &\xrightarrow{2.46e10 k18} \mathrm{H2O} + \mathrm{HO2} \\ \mathrm{OH} + \mathrm{CO} &\xrightarrow{2.46e10 k19} \mathrm{HO2} \end{aligned}

Figure 2

# Set start and end times.
start, finish = Dates.datetime2unix.([
    Dates.DateTime(2022, 5, 1),
    Dates.DateTime(2022, 5, 3),
])

# Create a model by combining SuperFast chemistry and FastJX photolysis.
model_ode = SuperFast(t) + FastJX(t)

# Run the simulation.
prob = ODEProblem(structural_simplify(get_mtk(model_ode)), [], (start, finish), [])
sol = solve(prob, TRBDF2(), saveat=1800.0)

# Plot result.
ticks = Dates.datetime2unix.([Dates.DateTime(2022, 5, 1), Dates.DateTime(2022, 5, 2), 
    Dates.DateTime(2022, 5, 3)])
tickstrings = Dates.format.(Dates.unix2datetime.(ticks), "m-d-yy")
Plots.plot(sol,ylims=(0,30),xlabel="Date", ylabel="Concentration (ppb)", 
    ylabelfontsize=9, xlabelfontsize=9, 
    xticks=(ticks, tickstrings), legend=:outertopright, size=(500, 310))

Figure 3

struct Emissions <: EarthSciMLODESystem
    sys
    function Emissions(t, μ_lon, σ)
        @variables NO(t) ISOP(t)
        dist = MvNormal([start,μ_lon],[3600.0,σ])
        @parameters lon
        D = Differential(t)
        new(ODESystem([
            D(NO) ~ pdf(dist, [t, lon]) * 50, 
            D(ISOP) ~ pdf(dist, [t, lon]) * 50,
        ], t, name=:emissions))
    end
end
Base.:(+)(e::Emissions, b::SuperFast) = operator_compose(b, e)
Base.:(+)(b::SuperFast, e::Emissions) = e + b

Supplemental code to skip unit enforcement

function ModelingToolkit.check_units(eqs...) # Skip unit enforcement for now
    nothing
    #ModelingToolkit.validate(eqs...) || @info "Some equations had invalid units. See warnings for details."
end

# Skip returning the observed variables (i.e. variables that are not state variables)
# because it is currently extremely slow to do so for this demo. 
SciMLBase.observed(sol::SciMLBase.AbstractTimeseriesSolution, sym, i::Colon) = zeros(Float64, length(sol.t))

Figure 4

domain = DomainInfo(
    partialderivatives_lonlat2xymeters,
    constIC(1.0f0, t ∈ Interval(start, finish)),
    zerogradBC(lon ∈ Interval(-130f0, -100f0)))

geos = GEOSFP("0.25x0.3125_NA", t; coord_defaults = Dict(:lat => 34.0, :lev => 1))

model = SuperFast(t) + FastJX(t) + domain +
    Emissions(t, -118.2, 0.2)+ Advection() + geos

g = graph(model)

f, ax, p = graphplot(g; ilabels=labels(g))
hidedecorations!(ax); hidespines!(ax); ax.aspect = DataAspect()

MethodError: no method matching EarthSciData.GEOSFP(::String, ::Symbolics.Num; coord_defaults::Dict{Symbol, Real})

// Image matching '/assets/tutorials/example/code/graph' not found. //

discretization = MOLFiniteDifference([lon => 50], t, approx_order=2)
prob = discretize(get_mtk(model), discretization)
sol = solve(prob, TRBDF2(), saveat=3600.0)
MethodError: no method matching discretize(::ModelingToolkit.ODESystem, ::MethodOfLines.MOLFiniteDifference{MethodOfLines.CenterAlignedGrid, MethodOfLines.ScalarizedDiscretization})

Closest candidates are:
  discretize(!Matched::ModelingToolkit.PDESystem, ::PDEBase.AbstractEquationSystemDiscretization; analytic, kwargs...)
   @ PDEBase ~/.julia/packages/PDEBase/Aqj4G/src/discretization_state.jl:55

Supplemental plotting code

discrete_lon = sol[lon]
discrete_t = sol[t]

@variables superfast₊O3(..) superfast₊NO(..) superfast₊ISOP(..) superfast₊NO2(..)
sol_isop = sol[superfast₊ISOP(t, lon)]
sol_o3 = sol[superfast₊O3(t, lon)]
sol_no = sol[superfast₊NO(t, lon)]
sol_no2 = sol[superfast₊NO2(t, lon)]

using LaTeXStrings
Plots.plot(
    Plots.heatmap(discrete_t, discrete_lon[1:end], sol_isop[:, 1:end]', 
        xticks=(ticks, tickstrings), 
        ylabel="Longitude (°)", title="Isoprene", titlefontsize=12,
        margin=0Plots.mm,
        size=(600, 400), dpi=300),
    Plots.heatmap(discrete_t, discrete_lon[1:end], sol_no[:, 1:end]', 
        xticks=(ticks, tickstrings), title="NO", titlefontsize=12),
    Plots.heatmap(discrete_t, discrete_lon[1:end], sol_no2[:, 1:end]', 
        xticks=(ticks, tickstrings), xlabel="Date", title=L"\textrm{NO_2}"),
    Plots.heatmap(discrete_t, discrete_lon[1:end], sol_o3[:, 1:end]',colorbar_title="Concentration (ppb)",
        xticks=(ticks, tickstrings), title=L"\textrm{O_3}", titlefontsize=12),
)
Plots.savefig(joinpath(@OUTPUT, "pde.svg"))
ArgumentError: lon is neither an observed nor a state variable.

// Image matching '/assets/tutorials/example/code/pde' not found. //

CC BY-SA 4.0 EarthSciML Authors and Contributors. Last modified: February 18, 2024. Website built with Franklin.jl and the Julia programming language.