From 8adc2303e359c9a9aa4024d86efd631348c39f9e Mon Sep 17 00:00:00 2001 From: Daines Date: Tue, 21 Jun 2022 18:44:07 +0100 Subject: [PATCH] updated documentation --- Project.toml | 2 +- docs/make.jl | 3 +- docs/src/HOWTOshowmodelandoutput.md | 49 ++++++++-- docs/src/PALEOmodel.md | 138 ---------------------------- docs/src/PALEOmodelOutput.md | 112 ++++++++++++++++++++++ docs/src/PALEOmodelSolvers.md | 95 +++++++++++++++++++ src/FieldRecord.jl | 4 +- src/ODE.jl | 38 ++++++-- src/ODEfixed.jl | 28 ++++-- src/OutputWriters.jl | 128 +++++++++++++++++--------- src/PALEOmodel.jl | 1 - src/PlotRecipes.jl | 40 +++++--- src/Run.jl | 37 ++++++-- src/SteadyState.jl | 48 +++++++--- src/SteadyStateKinsol.jl | 10 +- src/ThreadBarriers.jl | 2 +- 16 files changed, 493 insertions(+), 242 deletions(-) delete mode 100644 docs/src/PALEOmodel.md create mode 100644 docs/src/PALEOmodelOutput.md create mode 100644 docs/src/PALEOmodelSolvers.md diff --git a/Project.toml b/Project.toml index a0be3c2..c047e5d 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "PALEOmodel" uuid = "bf7b4fbe-ccb1-42c5-83c2-e6e9378b660c" authors = ["Stuart Daines "] -version = "0.14.6" +version = "0.14.7" [deps] DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0" diff --git a/docs/make.jl b/docs/make.jl index f19639f..067311a 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -17,7 +17,8 @@ makedocs(bib, sitename="PALEOmodel Documentation", "HOWTOsmallnegativevalues.md", ], "Reference" => [ - "PALEOmodel.md", + "PALEOmodelSolvers.md", + "PALEOmodelOutput.md", ], "References.md", "indexpage.md", diff --git a/docs/src/HOWTOshowmodelandoutput.md b/docs/src/HOWTOshowmodelandoutput.md index d700e0f..97c2757 100644 --- a/docs/src/HOWTOshowmodelandoutput.md +++ b/docs/src/HOWTOshowmodelandoutput.md @@ -70,11 +70,11 @@ To show linkage of a ReactionMethod Variable with localname "pO2PAL", Reaction " ## Display model output -Model output is stored in a [`PALEOmodel.OutputWriters.OutputMemory`](@ref) object, which is +Model output is stored in a [`PALEOmodel.AbstractOutputWriter`](@ref) object, which is available as `run.output`, ie the `output` field of the default [`PALEOmodel.Run`](@ref) instance created by the `COPSE_reloaded_reloaded.jl` script. -[`PALEOmodel.OutputWriters.OutputMemory`](@ref) stores model output by Domain: +The default [`PALEOmodel.OutputWriters.OutputMemory`](@ref) stores model output in memory, by Domain: julia> run.output # shows Domains @@ -111,17 +111,54 @@ Raw data arrays can also be accessed as Julia Vectors using `get_data`: ## Plot model output -The output can be plotted using the Julia Plots.jl package. Plot recipes are defined for [`PALEOmodel.FieldArray`](@ref), -so output data can be plotted directly: +The output can be plotted using the Julia Plots.jl package, see [Plotting output](@ref). Plot recipes are defined for [`PALEOmodel.FieldArray`](@ref), so output data can be plotted directly using the `plot` command: julia> using Plots julia> plot(run.output, "atm.pCO2atm") # plot output variable as a single command julia> plot(pCO2atm) # a PALEOmodel.FieldArray can be plotted julia> plot!(tmodel_raw, pCO2atm_raw, label="some raw data") # overlay data from standard Julia Vectors -## Spatial output +## Spatial or wavelength-dependent output -TODO - key point is that [`PALEOmodel.FieldArray`](@ref) includes coordinates to plot column and image data. +To analyze spatial or eg wavelength-dependent output (eg time series from a 1D column or 3D general circulation model, or quantities that are a function of wavelength or frequency), [`PALEOmodel.get_array`](@ref) takes additional arguments to take 1D or 2D slices from the spatial, spectral and timeseries data. The [`PALEOmodel.FieldArray`](@ref) returned includes coordinates to plot column (1D) and heatmap (2D) data. + +### Examples for a column-based model + +Visualisation of spatial and wavelength-dependent output from the PALEOdev.jl ozone photochemistry example (a single 1D atmospheric column): + +#### 1D column data + julia> plot(title="O3 mixing ratio", output, "atm.O3_mr", (tmodel=[0.0, 0.1, 1.0, 10.0, 100.0, 1000.0], column=1), + swap_xy=true, xaxis=:log, labelattribute=:filter_records) # plots O3 vs height + +Here the `labelattribute=:filter_records` keyword argument is used to generate plot labels from the `:filter_records` FieldArray attribute, which contains the `tmodel` values used to select the timeseries records. The plot recipe expands +the Vector-valued `tmodel` argument to overlay a sequence of plots. + +This is equivalent to first creating and then plotting a sequence of `FieldArray` objects: + + julia> O3_mr = PALEOmodel.get_array(run.output, "atm.O3_mr", tmodel=0.0, column=1) + julia> plot(title="O3 mixing ratio", O3_mr, swap_xy=true, xaxis=:log, labelattribute=:filter_records) + julia> O3_mr = PALEOmodel.get_array(run.output, "atm.O3_mr", tmodel=0.1, column=1) + julia> plot!(O3_mr, swap_xy=true, labelattribute=:filter_records) + +#### Wavelength-dependent data + julia> plot(title="direct transmittance", output, ["atm.direct_trans"], (tmodel=1e12, column=1, cell=[1, 80]), + ylabel="fraction", labelattribute=:filter_region) # plots vs wavelength + +Here `tmodel=1e12` selects the last model time output, and `column=1, cell=[1, 80]` selects the top and bottom cells within the first (only) 1D column. The `labelattribute=:filter_region` keyword argument is used to generate plot labels from the `:filter_region` FieldArray attribute, which contains the `column` and `cell` values used to select the spatial region. + +### Examples for a 3D GCM-based model + +Visualisation of spatial output from the 3D GENIE transport-matrix example (PALEOdev.jl repository) + +### Horizontal slices across levels + julia> heatmap(run.output, "ocean.O2_conc", (tmodel=1e12, k=1), swap_xy=true) + +Here `k=1` selects a horizontal level corresponding to model grid cells with index k=1, which is the ocean surface in the GENIE grid. + +### Vertical section at constant longitude + julia> heatmap(run.output, "ocean.O2_conc", (tmodel=1e12, i=10), swap_xy=true, mult_y_coord=-1.0) + +Here `i=10` selects a section at longitude corresponding to model grid cells with index i=10. ## Save and load output diff --git a/docs/src/PALEOmodel.md b/docs/src/PALEOmodel.md deleted file mode 100644 index 5a39825..0000000 --- a/docs/src/PALEOmodel.md +++ /dev/null @@ -1,138 +0,0 @@ - - -# PALEOmodel - -```@meta -CurrentModule = PALEOmodel -``` -```@docs -Run -``` -## Create and initialize -```@docs -initialize! -``` - - -## Integrate -```@meta -CurrentModule = PALEOmodel.ODE -``` -### DifferentialEquations solvers -```@docs -integrate -integrateDAE -integrateForwardDiff -integrateDAEForwardDiff -``` -### Fixed timestep solvers -```@meta -CurrentModule = PALEOmodel.ODEfixed -``` -```@docs -integrateEuler -integrateSplitEuler -integrateEulerthreads -integrateSplitEulerthreads -``` -```@meta -CurrentModule = PALEOmodel.ODELocalIMEX -``` -```@docs -integrateLocalIMEXEuler -``` -#### Fixed timestep wrappers - -```@meta -CurrentModule = PALEOmodel.ODEfixed -``` -```@docs -integrateFixed -integrateFixedthreads -``` - -### Steady-state solvers (Julia NLsolve based) -```@meta -CurrentModule = PALEOmodel.SteadyState -``` -```@docs -steadystate -steadystateForwardDiff -steadystate_ptc -steadystate_ptcForwardDiff -``` - -### Steady-state solvers (Sundials Kinsol based): -```@meta -CurrentModule = PALEOmodel.SteadyStateKinsol -``` -```@docs -steadystate_ptc -``` -```@meta -CurrentModule = PALEOmodel -``` -```@docs -Kinsol -Kinsol.kin_create -Kinsol.kin_solve -``` - -## Field Array - -```@meta -CurrentModule = PALEOmodel -``` -[`FieldArray`](@ref) provides a generic array type with named dimensions `PALEOboxes.NamedDimension` and optional coordinates `PALEOboxes.FixedCoord` for processing of model output. - -```@docs -FieldArray -get_array -``` - -```@docs -FieldRecord -``` - -## Output -```@meta -CurrentModule = PALEOmodel -``` -```@docs -OutputWriters -``` -```@meta -CurrentModule = PALEOmodel.OutputWriters -``` -```@docs -OutputMemory -OutputMemoryDomain -save_jld2 -load_jld2! -initialize! -add_record! -``` - -## Plot output - -```@meta -CurrentModule = PALEOmodel -``` -```@docs -RecipesBase.apply_recipe(::Dict{Symbol, Any}, fa::FieldArray) -RecipesBase.apply_recipe(::Dict{Symbol, Any}, output::AbstractOutputWriter, vars::Union{AbstractString, Vector{<:AbstractString}}, selectargs::NamedTuple) -RecipesBase.apply_recipe(::Dict{Symbol, Any}, outputs::Vector{<:AbstractOutputWriter}, vars::Union{AbstractString, Vector{<:AbstractString}}, selectargs::NamedTuple) -RecipesBase.apply_recipe(::Dict{Symbol, Any}, fr::FieldRecord, selectargs::NamedTuple) -RecipesBase.apply_recipe(::Dict{Symbol, Any}, fas::Vector{<:FieldArray}) - -PlotPager -Plot.test_heatmap_edges - -``` -## Analyze reaction network -```@meta -CurrentModule = PALEOmodel -``` -```@docs -ReactionNetwork -``` diff --git a/docs/src/PALEOmodelOutput.md b/docs/src/PALEOmodelOutput.md new file mode 100644 index 0000000..373338f --- /dev/null +++ b/docs/src/PALEOmodelOutput.md @@ -0,0 +1,112 @@ + + +# PALEOmodel output + +## Run + +```@meta +CurrentModule = PALEOmodel +``` +```@docs +Run +``` + +## Output +PALEO model output is accumulated into a container such as an [OutputMemory](@ref) instance that implements the [AbstractOutputWriter interface](@ref). + +### AbstractOutputWriter interface + +```@meta +CurrentModule = PALEOmodel +``` +```@docs +AbstractOutputWriter +``` + +```@meta +CurrentModule = PALEOmodel.OutputWriters +``` +#### Writing output +```@docs +initialize! +add_record! +``` +#### Modifying output +```@docs +PB.add_field!(output::PALEOmodel.AbstractOutputWriter, fr::PALEOmodel.FieldRecord) +``` +#### Querying output +```@docs +PB.get_table(output::PALEOmodel.AbstractOutputWriter, domainname::AbstractString) +PB.show_variables(output::PALEOmodel.AbstractOutputWriter) +PB.has_variable(output::PALEOmodel.AbstractOutputWriter, varname::AbstractString) +``` + +#### Accessing output data +```@docs +PALEOmodel.get_array(output::PALEOmodel.AbstractOutputWriter, varname::AbstractString; kwargs...) +PB.get_field(output::PALEOmodel.AbstractOutputWriter, varname::AbstractString) +PB.get_data(output::PALEOmodel.AbstractOutputWriter, varname::AbstractString; records=nothing) +PB.get_mesh(output::PALEOmodel.AbstractOutputWriter, domainname::AbstractString) +``` +```@meta +CurrentModule = PALEOmodel +``` +```@docs +FieldRecord +``` + +### OutputMemory + +```@meta +CurrentModule = PALEOmodel.OutputWriters +``` +```@docs +OutputMemory +OutputMemoryDomain +save_jld2 +load_jld2! +``` + +## Field Array + +```@meta +CurrentModule = PALEOmodel +``` +[`FieldArray`](@ref) provides a generic array type with named dimensions `PALEOboxes.NamedDimension` and optional coordinates `PALEOboxes.FixedCoord` for processing of model output. + +```@docs +FieldArray +get_array(f::PB.Field{D, PB.ScalarSpace, V, N, M}; attributes=nothing) where {D, V, N, M} +``` + +## Plotting output + +### Plot recipes +Plotting using the Julia [Plots.jl](https://github.com/JuliaPlots/Plots.jl) package is implemented by [plot recipes](https://docs.juliaplots.org/latest/recipes/) that enable plotting of PALEO data types using the `plot` command. + +The general principles are that: +- Data is extracted from model output into [`FieldArray`](@ref)s with attached coordinates +- Vector-valued arguments are "broadcast" to allow multiple line plot series to be overlaid in a single plot panel + +```@meta +CurrentModule = PALEOmodel +``` +```@docs +RecipesBase.apply_recipe(::Dict{Symbol, Any}, output::AbstractOutputWriter, vars::Union{AbstractString, Vector{<:AbstractString}}, selectargs::NamedTuple) +RecipesBase.apply_recipe(::Dict{Symbol, Any}, fr::FieldRecord, selectargs::NamedTuple) +RecipesBase.apply_recipe(::Dict{Symbol, Any}, fa::FieldArray) +``` +### Assembling multi-plot panels +```@docs +PlotPager +DefaultPlotPager +``` + +## Analyze reaction network +```@meta +CurrentModule = PALEOmodel +``` +```@docs +ReactionNetwork +``` diff --git a/docs/src/PALEOmodelSolvers.md b/docs/src/PALEOmodelSolvers.md new file mode 100644 index 0000000..ca20e51 --- /dev/null +++ b/docs/src/PALEOmodelSolvers.md @@ -0,0 +1,95 @@ + + +# PALEOmodel solvers + + +## Initialization +```@meta +CurrentModule = PALEOmodel +``` +```@docs +initialize! +``` + +## [DifferentialEquations](https://github.com/SciML/DifferentialEquations.jl) solvers + +Wrappers for the Julia [DifferentialEquations](https://github.com/SciML/DifferentialEquations.jl) package ODE and DAE solvers. These are usually appropriate for smaller biogeochemical models. + +NB: see [Managing small and negative values](@ref) for best practices and common issues when using ODE or DAE solvers. + +```@meta +CurrentModule = PALEOmodel.ODE +``` +```@docs +integrate +integrateDAE +``` +## Fixed timestep solvers +```@meta +CurrentModule = PALEOmodel.ODEfixed +``` +PALEO native fixed-timestep, first-order Euler integrators, with split-explicit and multi-threaded options. +These are usually appropriate for larger biogeochemical models (eg ocean models using GCM transport matrices). + +The low-level timestepping is provided by [`integrateFixed`](@ref) and [`integrateFixedthreads`](@ref), +with higher-level wrappers for common options provided by [`integrateEuler`](@ref) etc. + +### High-level wrappers +```@meta +CurrentModule = PALEOmodel.ODEfixed +``` +```@docs +integrateEuler +integrateSplitEuler +integrateEulerthreads +integrateSplitEulerthreads +``` +```@meta +CurrentModule = PALEOmodel.ODELocalIMEX +``` +```@docs +integrateLocalIMEXEuler +``` +### Low-level timesteppers +```@meta +CurrentModule = PALEOmodel.ODEfixed +``` +```@docs +integrateFixed +integrateFixedthreads +``` + +### Thread barriers +```@meta +CurrentModule = PALEOmodel.ThreadBarriers +``` +```@docs +ThreadBarrierAtomic +ThreadBarrierCond +``` + +## Steady-state solvers (Julia [`NLsolve`](https://github.com/JuliaNLSolvers/NLsolve.jl) based) +```@meta +CurrentModule = PALEOmodel.SteadyState +``` +```@docs +steadystate +steadystate_ptc +``` + +## Steady-state solvers (Sundials Kinsol based): +```@meta +CurrentModule = PALEOmodel.SteadyStateKinsol +``` +```@docs +steadystate_ptc +``` +```@meta +CurrentModule = PALEOmodel +``` +```@docs +Kinsol +Kinsol.kin_create +Kinsol.kin_solve +``` + diff --git a/src/FieldRecord.jl b/src/FieldRecord.jl index e45ecf2..f6e9d25 100644 --- a/src/FieldRecord.jl +++ b/src/FieldRecord.jl @@ -5,7 +5,7 @@ A series of records each containing a `PALEOboxes.Field`. # Implementation Fields with array values are stored in `records` as a Vector of arrays. -Fields with single values (field_single_element true) are stored as a Vector of eltype(Field.values). +Fields with single values (`field_single_element` true) are stored as a Vector of `eltype(Field.values)`. """ struct FieldRecord{D <: PB.AbstractData, S <: PB.AbstractSpace, V, N, M, R} records::Vector{R} @@ -137,7 +137,7 @@ spatial region defined by `selectargs`. `recordarg` can be one of: - `records=r::Int` to select a single record, or `records = first:last` to select a range. -- ``: (where eg =tmodel), `=t::Float64 to select a single record +- ``: (where eg `=tmodel`), `=t::Float64` to select a single record with the nearest value of `fr.coords_record`, or `=(first::Float64, last::Float64)` (a Tuple) to select a range starting at the record with the nearest value of `fr.coords_record` before `first` and ending at the nearest record after `last`. diff --git a/src/ODE.jl b/src/ODE.jl index 1b88598..450e57d 100644 --- a/src/ODE.jl +++ b/src/ODE.jl @@ -146,22 +146,31 @@ end """ integrate(run, initial_state, modeldata, tspan [; kwargs...] ) + integrateForwardDiff(run, initial_state, modeldata, tspan [;kwargs...]) Integrate run.model as an ODE or as a DAE with constant mass matrix, and write to `outputwriter` -# Arguments +Provides a wrapper around the Julia SciML [DifferentialEquations](https://github.com/SciML/DifferentialEquations.jl) +package ODE solvers, with PALEO-specific additional setup. Keyword arguments `alg` and `solvekwargs` are passed through to the +`DifferentialEquations` `solve` method. + +`integrateForwardDiff` sets keyword arguments `jac_ad=:ForwardDiffSparse`, `alg=Sundials.CVODE_BDF(linear_solver=:KLU)` +to use the Julia [ForwardDiff](https://github.com/JuliaDiff/ForwardDiff.jl) package to provide the Jacobian with +forward-mode automatic differentiation and automatic sparsity detection. +# Arguments - `run::Run`: struct with `model::PB.Model` to integrate and `output` field - `initial_state::AbstractVector`: initial state vector -- `modeldata::Modeldata`: ModelData struct with appropriate element type for forward model +- `modeldata::Modeldata`: ModelData struct - `tspan`: (tstart, tstop) integration start and stop times # Keywords - -- `alg=Sundials.CVODE_BDF()`: ODE algorithm to use -- `outputwriter=run.output`: PALEOmodel.AbstractOutputWriter instance to hold output +- `alg=Sundials.CVODE_BDF()`: ODE algorithm to use, passed through to DifferentialEquations.jl `solve` method. + The default is appropriate for a stiff system of equations (common in biogeochemical models), + see for other options. - `solvekwargs=NamedTuple()`: NamedTuple of keyword arguments passed through to DifferentialEquations.jl `solve` (eg to set `abstol`, `reltol`, `saveat`, see ) +- `outputwriter=run.output`: PALEOmodel.AbstractOutputWriter instance to hold output - `jac_ad=:NoJacobian`: Jacobian to generate and use (:NoJacobian, :ForwardDiffSparse, :ForwardDiff) - `jac_ad_t_sparsity=tspan[1]`: model time at which to calculate Jacobian sparsity pattern - `steadystate=false`: true to use `DifferentialEquations.jl` `SteadyStateProblem` @@ -231,10 +240,23 @@ end """ - integrateDAE(run, initial_state, modeldata, tspan;alg=IDA()) + integrateDAE(run, initial_state, modeldata, tspan [; kwargs...]) + integrateDAEForwardDiff(run, initial_state, modeldata, tspan [; kwargs...]) + +Integrate `run.model` as a DAE and copy output to `outputwriter`. + +Provides a wrapper around the Julia SciML [DifferentialEquations](https://github.com/SciML/DifferentialEquations.jl) +package DAE solvers, with PALEO-specific additional setup. Keyword arguments `alg` and `solvekwargs` are passed through to the +`DifferentialEquations` `solve` method. + +`integrateDAEForwardDiff` sets keyword arguments `jac_ad=:ForwardDiffSparse`, `alg=Sundials.CVODE_BDF(linear_solver=:KLU)` +to use the Julia [ForwardDiff](https://github.com/JuliaDiff/ForwardDiff.jl) package to provide the Jacobian with +forward-mode automatic differentiation and automatic sparsity detection. -Integrate run.model as a DAE and copy output to `outputwriter`. -Arguments as [`integrate`](@ref). +# Keywords +As [`integrate`](@ref), with defaults: +- `alg=Sundials.IDA()` (`integrateDAE`) +- `alg=Sundials.IDA(linear_solver=:KLU)` (`integrateDAEForwardDiff`) """ function integrateDAE( run, initial_state, modeldata, tspan; diff --git a/src/ODEfixed.jl b/src/ODEfixed.jl index 7ade509..fc15e72 100644 --- a/src/ODEfixed.jl +++ b/src/ODEfixed.jl @@ -12,18 +12,22 @@ import PALEOmodel ########################################################################### """ - integrateEuler(run, initial_state, modeldata, tspan, Δt; [,outputwriter]) + integrateEuler(run, initial_state, modeldata, tspan, Δt [; kwargs]) -Integrate run.model using first-order Euler with fixed timestep. +Integrate `run.model` from `initial_state` using first-order Euler with fixed timestep. + +Calls [`integrateFixed`](@ref) # Arguments - `run::Run`: struct with `model::PB.Model` to integrate and `output` field - `initial_state::AbstractVector`: initial state vector - `modeldata::Modeldata`: ModelData struct with appropriate element type for forward model - `tspan`: (tstart, toutput1, toutput2, ..., tstop) integration start, output, stop times -- `Δt`: (yr) fixed timestep -- [`outputwriter`: `PALEOmodel.AbstractOutputWriter` instance to write model output to] -- [`report_interval`: number of timesteps between progress update to console] +- `Δt`: (yr) fixed timestep + +# Keywords +- `outputwriter::PALEOmodel.AbstractOutputWriter=run.output`: container to write model output to +- `report_interval=1000`: number of timesteps between progress update to console """ function integrateEuler( run, initial_state, modeldata, tspan, Δt; @@ -79,8 +83,8 @@ NB: the combined time derivative is written to `outputwriter`. # Keywords - `cellranges_outer`: Vector of `CellRange` with `operatorID` defining `f_outer`. - `cellranges_inner`: Vector of `CellRange` with `operatorID` defining `f_inner`. -- [`outputwriter`: `PALEOmodel.AbstractOutputWriter` instance to write model output to] -- [`report_interval`: number of outer timesteps between progress update to console] +- `outputwriter::PALEOmodel.AbstractOutputWriter=run.output`: container to write model output to +- `report_interval=1000`: number of outer timesteps between progress update to console """ function integrateSplitEuler( run, initial_state, modeldata, tspan, Δt_outer, n_inner; @@ -144,7 +148,11 @@ Integrate run.model using first-order Euler with fixed timestep `Δt`, with tili - `modeldata::Modeldata`: ModelData struct with appropriate element type for forward model - `cellranges::Vector{Vector{AbstractCellRange}}`: Vector of Vector-of-cellranges, one per thread (so length(cellranges) == Threads.nthreads). - `tspan`: (tstart, toutput1, toutput2, ..., tstop) integration start, output, stop times -- `Δt`: (yr) fixed outer timestep +- `Δt`: (yr) fixed outer timestep + +# Keywords +- `outputwriter::PALEOmodel.AbstractOutputWriter=run.output`: container to write model output to +- `report_interval=1000`: number of outer timesteps between progress update to console """ function integrateEulerthreads( run, initial_state, modeldata, cellranges, tspan, Δt; @@ -213,8 +221,12 @@ Uses `Threads.nthreads` threads and tiling described by `cellranges_inner` and ` - `tspan`: (tstart, toutput1, toutput2, ..., tstop) integration start, output, stop times - `Δt_outer`: (yr) fixed outer timestep - `n_inner`: number of inner timesteps per outer timestep (0 for non-split solver) + +# Keywords - `cellranges_outer::Vector{Vector{AbstractCellRange}}`: Vector of list-of-cellranges, one list per thread (so length(cellranges) == Threads.nthreads), with `operatorID` defining `f_outer`. - `cellranges_inner::Vector{Vector{AbstractCellRange}}`: As `cellranges_outer`, with `operatorID` defining `f_inner`. +- `outputwriter::PALEOmodel.AbstractOutputWriter=run.output`: container to write model output to +- `report_interval=1000`: number of outer timesteps between progress update to console """ function integrateSplitEulerthreads( run, initial_state, modeldata, tspan, Δt_outer, n_inner::Int; diff --git a/src/OutputWriters.jl b/src/OutputWriters.jl index 378a2f8..ee19ad5 100644 --- a/src/OutputWriters.jl +++ b/src/OutputWriters.jl @@ -20,6 +20,64 @@ using Infiltrator # Julia debugger # AbstractOutputWriter interface ################################### +""" + AbstractOutputWriter + +Interface implemented by containers for PALEO model output. + +Implementations should define methods for: + +# Writing output +- [`initialize!`](@ref) +- [`add_record!`](@ref) + +# Modifying output +- [`PB.add_field!`](@ref) + +# Querying output +- [`PB.get_table`](@ref) +- [`PB.show_variables`](@ref) +- [`PB.has_variable`](@ref) + +# Accessing output data +- [`PALEOmodel.get_array`](@ref) +- [`PB.get_field`](@ref) +- [`PB.get_mesh`](@ref) +- [`PB.get_data`](@ref) + +""" +PALEOmodel.AbstractOutputWriter + +""" + initialize!(output::PALEOmodel.AbstractOutputWriter, model, modeldata, nrecords [;rec_coord=:tmodel]) + +Initialize from a PALEOboxes::Model, reserving memory for an assumed output dataset of `nrecords`. + +The default for `rec_coord` is `:tmodel`, for a sequence of records following the time evolution +of the model. +""" +function initialize!( + output::PALEOmodel.AbstractOutputWriter, model::PB.Model, modeldata::PB.ModelData, nrecords; + rec_coord::Symbol=:tmodel +) +end + +""" + add_record!(output::PALEOmodel.AbstractOutputWriter, model, modeldata, rec_coord) + +Add an output record for current state of `model` at record coordinate `rec_coord`. +The usual case (set by [`initialize!`](@ref)) is that the record coordinate is model time `tmodel`. +""" +function add_record!(output::PALEOmodel.AbstractOutputWriter, model, modeldata, rec_coord) end + +""" + add_field!(output::PALEOmodel.AbstractOutputWriter, fr::PALEOmodel.FieldRecord) + +Add [`PALEOmodel.FieldRecord`](@ref) `fr` to `output`, with Domain and name defined by `fr.attributes[:var_domain]` and +`fr.attributes[:var_name]`. +""" +function PB.add_field!(output::PALEOmodel.AbstractOutputWriter, fr::PALEOmodel.FieldRecord) end + """ get_table(output::PALEOmodel.AbstractOutputWriter, domainname::AbstractString) -> Table get_table(output::PALEOmodel.AbstractOutputWriter, varnames::Vector{<:AbstractString}) -> Table @@ -35,13 +93,27 @@ True if model `output` contains Variable `varname`. """ function PB.has_variable(output::PALEOmodel.AbstractOutputWriter, varname::AbstractString) end + +""" + show_variables(output::PALEOmodel.AbstractOutputWriter; kwargs...) -> Table + show_variables(output::PALEOmodel.AbstractOutputWriter, domainname::AbstractString; kwargs...) -> Table + +# Keywords: +- `attributes=[:units, :vfunction, :space, :field_data, :description]`: Variable attributes to include +- `filter = attrb->true`: function to filter by Variable attributes. + Example: `filter=attrb->attrb[:vfunction]!=PB.VF_Undefined` to show state Variables and derivatives. + +# Examples: + julia> vscodedisplay(PB.show_variables(run.output)) # show all output Variables in VS code table viewer + julia> vscodedisplay(PB.show_variables(run.output, ["atm.pCO2PAL", "fluxOceanBurial.flux_P_total"])) # show subset of output Variables in VS code table viewer +""" function PB.show_variables(output::PALEOmodel.AbstractOutputWriter) end """ get_array(output::PALEOmodel.AbstractOutputWriter, varname::AbstractString; kwargs...) -> FieldArray -Return a FieldArray containing data values and any attached coordinates, for the -FieldRecord for `varname`, and records and spatial region defined by `kwargs` +Return a [`PALEOmodel.FieldArray`](@ref) containing data values and any attached coordinates, for the +[`PALEOmodel.FieldRecord`](@ref) for `varname`, and records and spatial region defined by `kwargs` Equivalent to `PALEOmodel.get_array(PB.get_field(output, varname), kwargs...)` """ @@ -57,26 +129,24 @@ end """ get_field(output::PALEOmodel.AbstractOutputWriter, varname::AbstractString) -> FieldRecord -Return the `FieldRecord` for `varname`. +Return the [`PALEOmodel.FieldRecord`](@ref) for `varname`. """ function PB.get_field(output::PALEOmodel.AbstractOutputWriter, varname::AbstractString) end +""" + get_data(output::PALEOmodel.AbstractOutputWriter, varname; records=nothing) -> values + +Get Variable `varname` raw data array, optionally restricting to `records` +""" function PB.get_data(output::PALEOmodel.AbstractOutputWriter, varname::AbstractString; records=nothing) end """ get_mesh(output::PALEOmodel.AbstractOutputWriter, domainname::AbstractString) -> grid::Union{AbstractMesh, Nothing} -Return `grid`` for `output` Domain `domainname`. +Return `grid` for `output` Domain `domainname`. """ function PB.get_mesh(output::PALEOmodel.AbstractOutputWriter, domainname::AbstractString) end -""" - add_field!(output::PALEOmodel.AbstractOutputWriter, fr::PALEOmodel.FieldRecord) - -Add `FieldRecord` `fr` to `output`, with Domain and name defined by `fr.attributes[:var_domain]` and -`fr.attributes[:var_name]`. -""" -function PB.add_field!(output::PALEOmodel.AbstractOutputWriter, fr::PALEOmodel.FieldRecord) end ########################## # OutputMemoryDomain @@ -371,8 +441,13 @@ end """ OutputMemory -In-memory model output, organized by model Domains. +In-memory container for model output, organized by model Domains. + +Implements the [`PALEOmodel.AbstractOutputWriter`](@ref) interface, with additional methods +[`save_jld2`](@ref) and [`load_jld2!`](@ref) to save and load data. + +# Implementation Field `domains::Dict{String, OutputMemoryDomain}` contains per-Domain model output. """ struct OutputMemory <: PALEOmodel.AbstractOutputWriter @@ -437,19 +512,6 @@ function PB.get_mesh(output::OutputMemory, domainname::AbstractString) return output.domains[domainname].grid end -""" - show_variables(output::OutputMemory; kwargs...) -> Table - show_variables(output::OutputMemory, domainname::AbstractString; kwargs...) -> Table - -# Keywords: -- `attributes=[:units, :vfunction, :space, :field_data, :description]`: Variable attributes to include -- `filter = attrb->true`: function to filter by Variable attributes. - Example: `filter=attrb->attrb[:vfunction]!=PB.VF_Undefined` to show state Variables and derivatives. - -# Examples: - julia> vscodedisplay(PB.show_variables(run.output)) # show all output Variables in VS code table viewer - julia> vscodedisplay(PB.show_variables(run.output, ["atm.pCO2PAL", "fluxOceanBurial.flux_P_total"])) # show subset of output Variables in VS code table viewer -""" function PB.show_variables( output::OutputMemory, domainname::AbstractString; kwargs... ) @@ -555,14 +617,7 @@ function Base.append!(output1::OutputMemory, output2::OutputMemory) return output1 end -""" - initialize!(output::OutputMemory, model, modeldata, nrecords [;rec_coord=:tmodel]) - -Initialize from a PALEOboxes::Model, reserving memory for an assumed output dataset of `nrecords`. -The default for `rec_coord` is `:tmodel`, for a sequence of records following the time evolution -of the model. -""" function initialize!( output::OutputMemory, model::PB.Model, modeldata::PB.ModelData, nrecords; rec_coord::Symbol=:tmodel @@ -580,12 +635,7 @@ function initialize!( return nothing end -""" - add_record!(output::OutputMemory, model, modeldata, rec_coord) -Add an output record for current state of `model` at record coordinate `rec_coord`. -The usual case (set by [`initialize!`](@ref)) is that the record coordinate is model time `tmodel`. -""" function add_record!(output::OutputMemory, model, modeldata, rec_coord) for dom in model.domains @@ -625,11 +675,7 @@ function PB.add_field!(output::OutputMemory, fr::PALEOmodel.FieldRecord) return PB.add_field!(output.domains[domainname], fr) end -""" - get_data(output::OutputMemory, varname; records=nothing) -> values -Get Variable `varname` raw data array, optionally restricting to `records` -""" PB.get_data(output::OutputMemory, varname::AbstractString; records=nothing) = _get_outputvar(output, nothing, varname, records) diff --git a/src/PALEOmodel.jl b/src/PALEOmodel.jl index 290925a..3a26d28 100644 --- a/src/PALEOmodel.jl +++ b/src/PALEOmodel.jl @@ -10,7 +10,6 @@ import SparsityTracing PB.value_ad(x::SparsityTracing.ADval) = SparsityTracing.value(x) PB.value_ad(x::ForwardDiff.Dual) = ForwardDiff.value(x) -"base type for OutputWriters" abstract type AbstractOutputWriter end diff --git a/src/PlotRecipes.jl b/src/PlotRecipes.jl index e253593..a2ee3d6 100644 --- a/src/PlotRecipes.jl +++ b/src/PlotRecipes.jl @@ -12,10 +12,10 @@ import Infiltrator Accumulate plots into subplots. `layout` is supplied to `Plots.jl` `layout` keyword, may be an Int or a Tuple (ny, nx), -see https://docs.juliaplots.org/latest/ +see . -Optional `kwargs` provides keyword arguments supplied to `plot` -(eg (legend_background_color=nothing, ) to set all subplot legends to transparent backgrounds) +Optional `kwargs::NamedTuple` provides additional keyword arguments passed through to `plot` +(eg `(legend_background_color=nothing, )` to set all subplot legends to transparent backgrounds) # Usage @@ -78,7 +78,7 @@ function (pp::PlotPager)(p, ps... ) end """ - DefaultPlotPager + DefaultPlotPager() Dummy version of [`PlotPager`](@ref) that just calls `display` for each plot. """ @@ -101,7 +101,7 @@ function (pp::DefaultPlotPager)(p, ps... ) end """ - NoPlotPager + NoPlotPager() Dummy version of [`PlotPager`](@ref) that doesn't display the plot """ @@ -116,12 +116,18 @@ end """ plot(output::AbstractOutputWriter, vars::Union{AbstractString, Vector{<:AbstractString}}, selectargs::NamedTuple=NamedTuple()) + heatmap(output::AbstractOutputWriter, vars::Union{AbstractString, Vector{<:AbstractString}}, selectargs::NamedTuple=NamedTuple()) + plot(outputs::Vector{<:AbstractOutputWriter}, vars::Union{AbstractString, Vector{<:AbstractString}}, selectargs::NamedTuple=NamedTuple()) + Plot recipe that calls `PB.get_field(output, var)`, and passes on to `plot(fr::FieldRecord, selectargs)` +(see [`RecipesBase.apply_recipe(::Dict{Symbol, Any}, fr::FieldRecord, selectargs::NamedTuple)`](@ref)) + +Vector-valued `outputs` or `vars` are "broadcast" to create a plot series for each element. +A `labelprefix` (index in `outputs` Vector) is added to identify each plot series produced. -Additional features: -- if `vars` is a Vector, create a plot series for each element. -- if `var` is of form `..`, then set `structfield` to take a single field from a `struct` Variable +If `var` is of form `..`, then sets the `structfield` keyword argument +to take a single field from a `struct` Variable. """ RecipesBase.@recipe function f( output::AbstractOutputWriter, @@ -174,10 +180,14 @@ end """ plot(fr::FieldRecord, selectargs::NamedTuple) + heatmap(fr::FieldRecord, selectargs::NamedTuple) -Plot recipe that calls `get_array(fr; selectargs...)` and passes on to `plot(fa::FieldArray)`. +Plot recipe to plot a single [`FieldRecord`](@ref) -Vector-valued fields in `selectargs` are broadcasted (generating a separate plot series for each combination) +Calls `get_array(fr; selectargs...)` and passes on to `plot(fa::FieldArray)` +(see [`RecipesBase.apply_recipe(::Dict{Symbol, Any}, fa::FieldArray)`](@ref)). + +Vector-valued fields in `selectargs` are "broadcasted" (generating a separate plot series for each combination) """ RecipesBase.@recipe function f(fr::FieldRecord, selectargs::NamedTuple) @@ -194,9 +204,15 @@ end """ plot(fa::FieldArray; kwargs...) + heatmap(fa::FieldArray; kwargs...) + plot(fas::Vector{<:FieldArray}; kwargs...) + +Plot recipe that plots a single [`FieldArray`] or Vector of [`FieldArray`]. + +If `fa` has a single dimension, this is suitable for a line-like plot, if two dimensions, a heatmap. -Plot recipe that plots `fa`. If `fa` has a single dimension, this is suitable for a line-like plot, -if two dimensions, a heatmap. +If `fas::Vector` is supplied, this is "broadcast" generating one plot series for each element, +with the Vector index prepended to `labelprefix` to identify the plot series (unless overridden by `labellist` or `labelattribute`) # Keywords - `swap_xy::Bool=false`: true to swap x and y axes diff --git a/src/Run.jl b/src/Run.jl index c8a4eb7..9533e1b 100644 --- a/src/Run.jl +++ b/src/Run.jl @@ -15,30 +15,49 @@ Base.@kwdef mutable struct Run end -"compact form" function Base.show(io::IO, run::Run) - print(io, "Run(model=", run.model,")") + print(io, "Run(model='", run.model, "', output='", run.output, "')") end -"multiline form" function Base.show(io::IO, ::MIME"text/plain", run::Run) println(io, "PALEOmodel.Run") println(io, " model='", run.model,"'") + println(io, " output='", run.output,"'") end initialize!(run::Run; kwargs...) = initialize!(run.model; kwargs...) """ - initialize!(model::PB.Model; kwargs...) -> (initial_state, modeldata) + initialize!(model::PB.Model; kwargs...) -> (initial_state::Vector, modeldata::PB.ModelData) + initialize!(run::Run; kwargs...) -> (initial_state::Vector, modeldata::PB.ModelData) -Initialize `model` and return `initial_state` Vector and `modeldata` struct +Initialize `model` or `run.model` and return: +- an `initial_state` Vector +- a `modeldata` struct containing the model arrays. -# Keywords: -- `eltype::Type=Float64`: default data type to use for model arrays -- `eltypemap=Dict{String, DataType}`: Dict of data types to look up Variable :datatype attribute +# Initialising the state vector +With default arguments, the `model` state variables are initialised to the values defined in the `.yaml` +configuration file used to create the `model`. + +The optional `pickup_output` argument can be used to provide an `OutputWriter` instance with pickup data to initialise from. +This is applied afer the default initialisation, hence can be used to (re)initialise a subset of model state variables. + +# `DataType`s for model arrays +With default arguments, the `model` arrays use `Float64` as the element type. The `eltype` keyword argument can be used to +specify a different Julia `DataType`, eg for use with automatic differentiation. Per-Variable `DataType` can be specified by using +the `:datatype` Variable attribute to specify `String`-valued tags, in combination with the `eltypemap` keyword argument to +provide a `Dict` of tag names => `DataType`s. + +# Thread safety +A thread-safe model can be created with `threadsafe=true` (to create Atomic Variables for those Variables with attribute `:atomic==true`), +and supplying `method_barrier` (a thread barrier to add to `ReactionMethod` dispatch lists between dependency-free groups) + +# Keyword summary - `pickup_output=nothing`: OutputWriter with pickup data to initialise from -- `threadsafe=false`: true to create thread safe Atomic Variables where :atomic attribute = true +- `eltype::Type=Float64`: default data type to use for model arrays +- `eltypemap=Dict{String, DataType}()`: Dict of data types to look up Variable :datatype attribute +- `threadsafe=false`: true to create thread safe Atomic Variables where Variable attribute `:atomic==true` - `method_barrier=nothing`: thread barrier to add to dispatch lists if `threadsafe==true` - `expect_hostdep_varnames=["global.tforce"]`: non-state-Variable host-dependent Variable names expected """ diff --git a/src/SteadyState.jl b/src/SteadyState.jl index f2c65c0..1c9490b 100644 --- a/src/SteadyState.jl +++ b/src/SteadyState.jl @@ -20,10 +20,14 @@ using SparseDiffTools """ steadystate(run, initial_state, modeldata, tss [; kwargs...] ) + steadystateForwardDiff(run, initial_state, modeldata, tss [; kwargs...] ) Find steady-state solution (using [NLsolve.jl](https://github.com/JuliaNLSolvers/NLsolve.jl) package) and write to `outputwriter` (two records are written, for `initial_state` and the steady-state solution). +`steadystateForwardDiff` has default keyword argument `jac_ad=:ForwardDiffSparse` to use automatic differentiation +for sparse Jacobian. + # Arguments - `run::Run`: struct with `model::PB.Model` to integrate and `output` field - `initial_state::AbstractVector`: initial state vector @@ -31,7 +35,7 @@ and write to `outputwriter` (two records are written, for `initial_state` and th - `tss`: (yr) model tforce time for steady state solution # Optional Keywords -- `outputwriter=run.output`: PALEOmodel.AbstractOutputWriter instance to hold output +- `outputwriter::PALEOmodel.AbstractOutputWriter=run.output`: container to hold output - `initial_time=-1.0`: tmodel to write for first output record - `solvekwargs=NamedTuple()`: NamedTuple of keyword arguments passed through to [NLsolve.jl](https://github.com/JuliaNLSolvers/NLsolve.jl) (eg to set `method`, `ftol`, `iteration`, `show_trace`, `store_trace`). @@ -182,12 +186,16 @@ end ############################################################## """ - steadystate_ptc(run, initial_state, modeldata, tss, deltat_initial, tss_max [; kwargs...]) - + steadystate_ptc(run, initial_state, modeldata, tspan, deltat_initial; kwargs...) + steadystate_ptcForwardDiff(run, initial_state, modeldata, tspan, deltat_initial; kwargs...) + Find steady-state solution and write to `outputwriter`, using naive pseudo-transient-continuation -with first order implicit Euler pseudo-timesteps from `tss` to `tss_max` +with first order implicit Euler pseudo-timesteps from `tspan[1]` to `tspan[2]` and [NLsolve.jl](https://github.com/JuliaNLSolvers/NLsolve.jl) as the non-linear solver. +`steadystate_ptcForwardDiff` has keyword argument default `jac_ad=:ForwardDiffSparse` to use automatic differentiation for +sparse Jacobian. + Each pseudo-timestep solves the nonlinear system S(t+Δt) = S(t) + Δt dS/dt(t+Δt) for S(t+Δt), using a variant of Newton's method. @@ -198,8 +206,16 @@ until pseudo-time `tss_max` is reached. If an iteration fails, Δt is divided by NB: this is a _very_ naive initial implementation, there is currently no reliable error control to adapt pseudo-timesteps to the rate of convergence, so requires some trial-and-error to set an appropiate `deltat_fac` for each problem. +# Arguments +- `run::Run`: struct with `model::PB.Model` to integrate and `output` field +- `initial_state::AbstractVector`: initial state vector +- `modeldata::Modeldata`: ModelData struct with appropriate element type for forward model +- `tspan`: Vector or Tuple with `(initial_time, final_time)` +- `deltat_initial`: initial pseudo-timestep + # Keywords - `deltat_fac=2.0`: factor to increase pseudo-timestep on success +- `tss_output=[]`: Vector of model times at which to save output (empty Vector to save all output timesteps) - `outputwriter=run.output`: output destination - `solvekwargs=NamedTuple()`: arguments to pass through to NLsolve - `jac_ad=:NoJacobian`: AD Jacobian to use @@ -208,14 +224,11 @@ to the rate of convergence, so requires some trial-and-error to set an appropiat (eg to restrict to an approximate Jacobian) - `enforce_noneg=false`: fail pseudo-timesteps that generate negative values for state variables. - `use_norm=false`: true to apply PALEO norm_value to state variables -- `verbose=false`: true to detailed output +- `verbose=false`: true for detailed output - `BLAS_num_threads=1`: restrict threads used by Julia BLAS (likely irrelevant if using sparse Jacobian?) - -See [`steadystate`](@ref) for more details of arguments. - """ function steadystate_ptc( - run, initial_state, modeldata, tss::Float64, deltat_initial::Float64, tss_max::Float64; + run, initial_state, modeldata, tspan, deltat_initial::Float64; deltat_fac=2.0, tss_output=[], outputwriter=run.output, @@ -231,6 +244,9 @@ function steadystate_ptc( nevals = 0 + # start, end times + tss, tss_max = tspan + # workaround Julia BLAS default (mis)configuration that defaults to multi-threaded LinearAlgebra.BLAS.set_num_threads(BLAS_num_threads) @info "steadystate_ptc: using BLAS with $(LinearAlgebra.BLAS.get_num_threads()) threads" @@ -525,15 +541,14 @@ function steadystate_ptc( return nothing end -"[`steadystate_ptc`](@ref) with argument defaults to use ForwardDiff AD Jacobian" function steadystate_ptcForwardDiff( - run, initial_state, modeldata, tss, deltat_initial, tss_max; + run, initial_state, modeldata, tspan, deltat_initial::Float64; jac_ad=:ForwardDiffSparse, kwargs... ) return steadystate_ptc( - run, initial_state, modeldata, tss, deltat_initial, tss_max; + run, initial_state, modeldata, tspan, deltat_initial; # (:jac_ad=>:ForwardDiffSparse, kwargs...)... jac_ad=jac_ad, kwargs... @@ -541,4 +556,13 @@ function steadystate_ptcForwardDiff( end +steadystate_ptc( + run, initial_state, modeldata, tss::Float64, deltat_initial::Float64, tss_max::Float64; kwargs... +) = steadystate_ptc(run, initial_state, modeldata, (tss, tss_max), deltat_initial; kwargs...) + +steadystate_ptcForwardDiff( + run, initial_state, modeldata, tss::Float64, deltat_initial::Float64, tss_max::Float64; kwargs... +) = steadystate_ptcForwardDiff(run, initial_state, modeldata, (tss, tss_max), deltat_initial; kwargs...) + + end # module diff --git a/src/SteadyStateKinsol.jl b/src/SteadyStateKinsol.jl index aed25b1..795dfec 100644 --- a/src/SteadyStateKinsol.jl +++ b/src/SteadyStateKinsol.jl @@ -14,7 +14,7 @@ using SparseDiffTools """ - steadystate_ptc(run, initial_state, modeldata, tss, deltat_initial, tss_max; + steadystate_ptc(run, initial_state, modeldata, tspan, deltat_initial; [,deltat_fac=2.0] [,tss_output] [,outputwriter] [,createkwargs] [,solvekwargs] [,jac_cellranges] [, use_directional_ad] [, directional_ad_eltypestomap] [,verbose] [, BLAS_num_threads] ) @@ -44,7 +44,7 @@ the default finite difference approximation). that should be mapped to the AD directional derivative datatype hence included in the AD directional derivative. """ function steadystate_ptc( - run, initial_state, modeldata, tss, deltat_initial, tss_max; + run, initial_state, modeldata, tspan, deltat_initial::Float64; deltat_fac=2.0, tss_output=[], outputwriter=run.output, @@ -57,6 +57,8 @@ function steadystate_ptc( verbose=false, BLAS_num_threads=1 ) + # start, end times + tss, tss_max = tspan # workaround Julia BLAS default (mis)configuration that defaults to multi-threaded LinearAlgebra.BLAS.set_num_threads(BLAS_num_threads) @@ -312,4 +314,8 @@ function steadystate_ptc( return nothing end +steadystate_ptc( + run, initial_state, modeldata, tss::Float64, deltat_initial::Float64, tss_max::Float64; kwargs... +) = steadystate_ptc(run, initial_state, modeldata, (tss, tss_max), deltat_initial; kwargs...) + end # module diff --git a/src/ThreadBarriers.jl b/src/ThreadBarriers.jl index 386db26..5c3dc4e 100644 --- a/src/ThreadBarriers.jl +++ b/src/ThreadBarriers.jl @@ -142,4 +142,4 @@ function wait_barrier(barrier::ThreadBarrierCond) #, cond) return nothing end -end +end # module