Add option to load all data at once.

Project.toml

@@ -1,7 +1,7 @@
 name = "EarthSciData"
 uuid = "a293c155-435f-439d-9c11-a083b6b47337"
 authors = ["EarthSciML Authors"]
-version = "0.9.4"
+version = "0.10.0"
 
 [deps]
 DataInterpolations = "82cc6244-b520-54b8-b5a6-8a565e85f1d0"
@@ -26,15 +26,15 @@ Scratch = "6c6a2e73-6563-6170-7368-637461726353"
 Symbolics = "0c5d862f-8b57-4792-8d23-62f2024744c7"
 
 [compat]
-AllocCheck = "0.1"
+AllocCheck = "0.2"
 DataInterpolations = "3, 4, 5, 6"
-DiffEqCallbacks = "2, 3"
+DiffEqCallbacks = "2, 3, 4"
 DifferentialEquations = "7"
 DocStringExtensions = "0.9"
 DomainSets = "0.7"
 Downloads = "1"
 DynamicQuantities = "0.13, 1"
-EarthSciMLBase = "0.16"
+EarthSciMLBase = "0.18"
 Interpolations = "0.15"
 LaTeXStrings = "1"
 Latexify = "0.11, 0.12, 0.13, 0.14, 0.15, 0.16"

bench/runbenchmarks.jl

@@ -44,7 +44,7 @@ function EarthSciMLBase.couple2(sys::SysCoupler, emis::EarthSciData.NEI2016Month
     sys, emis = sys.sys, emis.sys
     operator_compose(sys, emis)
 end
-function nei_simulator()
+function nei_simulator(st)
     @parameters lat = deg2rad(40.0) lon = deg2rad(-97.0) lev = 0.0
     @variables ACET(t) = 0.0 [unit = u"kg*m^-3"]
     @constants c = 1000 [unit = u"s"]
@@ -64,22 +64,26 @@ function nei_simulator()
             lon ∈ Interval(deg2rad(-115), deg2rad(-68.75)),
             lat ∈ Interval(deg2rad(25), deg2rad(53.7)),
             lev ∈ Interval(1, 2)
-        ))
+        );
+        grid_spacing=[deg2rad(1.25), deg2rad(1.2), 1])
 
     csys = couple(sys, emis, domain)
 
-    Simulator(csys, [deg2rad(1.25), deg2rad(1.2), 1])
+    ODEProblem(csys, st)
 end
 
 suite["NEI Simulator"] = BenchmarkGroup()
-sim = nei_simulator()
 st = SimulatorStrangSerial(Tsit5(), Euler(), 100.0)
-suite["NEI Simulator"]["Serial"] = @benchmarkable run!($sim, $st)
+sim = nei_simulator(st)
+suite["NEI Simulator"]["Serial"] = @benchmarkable solve($sim, dt=100.0,
+    save_on=false, save_start=false, save_end=false, initialize_save=false)
 
 st = SimulatorStrangThreads(Tsit5(), Euler(), 100.0)
-suite["NEI Simulator"]["Threads"] = @benchmarkable run!($sim, $st)
+sim = nei_simulator(st)
+suite["NEI Simulator"]["Threads"] = @benchmarkable solve($sim, dt=100.0,
+    save_on=false, save_start=false, save_end=false, initialize_save=false)
 
-tune!(suite, verbose = true)
-results = run(suite, verbose = true)
+tune!(suite, verbose=true)
+results = run(suite, verbose=true)
 
 BenchmarkTools.save("output.json", median(results))

docs/Project.toml

@@ -16,7 +16,7 @@ DifferentialEquations = "7"
 Documenter = "1"
 DomainSets = "0.7"
 DynamicQuantities = "0.13, 1"
-EarthSciMLBase = "0.16"
+EarthSciMLBase = "0.18"
 MethodOfLines = "0.11"
 ModelingToolkit = "9"
 Plots = "1"

docs/src/geosfp.md

@@ -6,21 +6,19 @@ First, let's initialize some packages and set up the [GEOS-FP](@ref GEOSFP) equa
 
 ```@example geosfp
 using EarthSciData, EarthSciMLBase
-using DomainSets, ModelingToolkit, MethodOfLines, DifferentialEquations
+using ModelingToolkit, MethodOfLines, DifferentialEquations
 using ModelingToolkit: t, D
 using Dates, Plots, DataFrames
 using DynamicQuantities
 using DynamicQuantities: dimension
 
 # Set up system
-@parameters(
-    lev, 
-    lon, [unit=u"rad"], 
-    lat, [unit=u"rad"],
-)
-geosfp, geosfp_updater = GEOSFP("4x5")
+domain = DomainInfo(DateTime(2022, 1, 1), DateTime(2022, 1, 3);
+    latrange=deg2rad(-85.0f0):deg2rad(2):deg2rad(85.0f0),
+    lonrange=deg2rad(-180.0f0):deg2rad(2.5):deg2rad(175.0f0),
+    levelrange=1:11, dtype=Float32)
 
-geosfp
+geosfp = GEOSFP("4x5", domain)
 ```
 
 Note that the [`GEOSFP`](@ref) function returns to things, an equation system and an object that can used to update the time in the underlying data loaders. 
@@ -42,6 +40,8 @@ To fix this, we create another equation system that is an ODE.
 
 ```@example geosfp
 function Example()
+    @parameters lat=0.0 [unit=u"rad"]
+    @parameters lon=0.0 [unit=u"rad"]
     @variables c(t) = 5.0 [unit=u"s"]
     ODESystem([D(c) ~ sin(lat * 6) + sin(lon * 6)], t, name=:Docs₊Example)
 end
@@ -51,21 +51,10 @@ examplesys = Example()
 Now, let's couple these two systems together, and also add in advection and some information about the domain:
 
 ```@example geosfp
-domain = DomainInfo(
-    partialderivatives_δxyδlonlat,
-    constIC(0.0, t ∈ Interval(Dates.datetime2unix(DateTime(2022, 1, 1)), Dates.datetime2unix(DateTime(2022, 1, 3)))),
-    zerogradBC(lat ∈ Interval(deg2rad(-80.0f0), deg2rad(80.0f0))),
-    periodicBC(lon ∈ Interval(deg2rad(-180.0f0), deg2rad(180.0f0))),
-    zerogradBC(lev ∈ Interval(1.0f0, 11.0f0)),
-)
-
-composed_sys = couple(examplesys, domain, geosfp, geosfp_updater)
+composed_sys = couple(examplesys, domain, geosfp)
 pde_sys = convert(PDESystem, composed_sys)
 ```
 
-You can see above that we add both `geosfp` and `geosfp_updater` to our coupled system. 
-If we didn't include the updater, the resulting model would not give the correct results.
-
 Now, finally, we can run the simulation and plot the GEOS-FP wind fields in the result:
 
 (The code below is commented out because it is very slow right now. A faster solution is coming soon!)

docs/src/nei2016.md

@@ -20,11 +20,21 @@ export JULIA_NO_VERIFY_HOSTS=gaftp.epa.gov
 This is what its equation system looks like:
 
 ```@example nei2016
-using EarthSciData, ModelingToolkit, DynamicQuantities, DataFrames
+using EarthSciData, EarthSciMLBase
+using ModelingToolkit, DynamicQuantities, DataFrames
 using ModelingToolkit: t
 using DynamicQuantities: dimension
-@parameters lat, [unit=u"rad"], lon, [unit=u"rad"], lev [unit=u"rad"]
-emis, emis_updater = NEI2016MonthlyEmis("mrggrid_withbeis_withrwc", lon, lat, lev)
+using Dates
+
+domain = DomainInfo(
+        DateTime(2016, 5, 1), DateTime(2016, 5, 2);
+        lonrange=deg2rad(-115):deg2rad(2.5):deg2rad(-68.75),
+        latrange=deg2rad(25):deg2rad(2):deg2rad(53.7),
+        levelrange=1:10, 
+        dtype=Float32,
+)
+
+emis = NEI2016MonthlyEmis("mrggrid_withbeis_withrwc", domain)
 ```
 
 ## Variables

src/EarthSciData.jl

@@ -12,7 +12,6 @@ using Scratch
 # General utilities
 include("load.jl")
 include("utils.jl")
-include("update_callback.jl")
 
 # Specific data sets
 include("netcdf.jl")

src/geosfp.jl

@@ -219,16 +219,13 @@ Domain options (as of 2022-01-30):
 - NATIVE
 - c720
 
-`coord_defaults` can be used to provide default values for the coordinates of the
-domain. For example if we want to perform a 2D simulation with a vertical dimension,
-we can set `coord_defaults = Dict(:lev => 1)`.
+The native data type for this dataset is Float32.
 
-`dtype` represents the desired data type of the interpolated values. The native data type
-for this dataset is Float32.
+`stream_data` specifies whether the data should be streamed in as needed or loaded all at once.
 
-See http://geoschemdata.wustl.edu/ExtData/ for current options.
+See http://geoschemdata.wustl.edu/ExtData/ for current data domain options.
 """
-function GEOSFP(domain::AbstractString; coord_defaults=Dict{Symbol,Number}(), spatial_ref="+proj=longlat +datum=WGS84 +no_defs", dtype=Float32, name=:GEOSFP, kwargs...)
+function GEOSFP(domain::AbstractString, domaininfo::DomainInfo; name=:GEOSFP, stream_data=true)
     filesets = Dict{String,GEOSFPFileSet}(
         "A1" => GEOSFPFileSet(domain, "A1"),
         "A3cld" => GEOSFPFileSet(domain, "A3cld"),
@@ -237,48 +234,40 @@ function GEOSFP(domain::AbstractString; coord_defaults=Dict{Symbol,Number}(), sp
         "A3mstE" => GEOSFPFileSet(domain, "A3mstE"),
         "I3" => GEOSFPFileSet(domain, "I3"))
 
-    sample_time = DateTime(2022, 5, 1) # Dummy time to get variable names and dimensions from data.
-    eqs = []
-    vars = []
-    itps = []
+    starttime, endtime = EarthSciMLBase.tspan_datetime(domaininfo)
+    pvdict = Dict([Symbol(v) => v for v in EarthSciMLBase.pvars(domaininfo)]...)
+    eqs = Equation[]
+    vars = Num[]
     for (filename, fs) in filesets
-        for varname ∈ varnames(fs, sample_time)
-            itp = DataSetInterpolator{dtype}(fs, varname, sample_time; spatial_ref, kwargs...)
-            dims = dimnames(itp, sample_time)
+        for varname ∈ varnames(fs, starttime)
+            dt = EarthSciMLBase.dtype(domaininfo)
+            itp = DataSetInterpolator{dt}(fs, varname, starttime, endtime,
+                domaininfo.spatial_ref; stream_data=stream_data)
+            dims = dimnames(itp)
             coords = Num[]
-            for dim ∈ dims
+            for dim in dims
                 d = Symbol(dim)
-                if d ∈ keys(coord_defaults) # Set a default value for this coordinate.
-                    v = (@parameters $d = coord_defaults[d])[1]
-                else # No default value.
-                    v = (@parameters $d)[1]
-                end
-                push!(coords, v)
+                @assert d ∈ keys(pvdict) "GEOSFP coordinate $d not found in domaininfo coordinates ($(pvs))."
+                push!(coords, pvdict[d])
             end
-            eq = create_interp_equation(itp, filename, t, sample_time, coords)
+            eq = create_interp_equation(itp, filename, t, coords)
             push!(eqs, eq)
             push!(vars, eq.lhs)
-            push!(itps, itp)
         end
     end
 
     # Implement hybrid grid pressure: https://wiki.seas.harvard.edu/geos-chem/index.php/GEOS-Chem_vertical_grids
     @constants P_unit = 1.0 [unit = u"Pa", description = "Unit pressure"]
     @variables P(t) [unit = u"Pa", description = "Pressure"]
     @variables I3₊PS(t) [unit = u"Pa", description = "Pressure at the surface"]
-    d = :lev
-    if d ∈ keys(coord_defaults) # Set a default value for this coordinate.
-        lev = (@parameters $d = coord_defaults[d])[1]
-    else # No default value.
-        lev = (@parameters $d)[1]
-    end
+    @assert :lev in keys(pvdict) "GEOSFP coordinate :lev not found in domaininfo coordinates ($(pvs))."
+    lev = pvdict[:lev]
     pressure_eq = P ~ P_unit * Ap(lev) + Bp(lev) * I3₊PS
     push!(eqs, pressure_eq)
 
-    sys = ODESystem(eqs, t, name=name,
+    sys = ODESystem(eqs, t, vars, [pvdict[:lon], pvdict[:lat], lev], name=name,
         metadata=Dict(:coupletype => GEOSFPCoupler))
-    cb = UpdateCallbackCreator(sys, vars, itps)
-    return sys, cb
+    return sys
 end
 
 function EarthSciMLBase.couple2(mw::EarthSciMLBase.MeanWindCoupler, g::GEOSFPCoupler)