diff --git a/dev/.documenter-siteinfo.json b/dev/.documenter-siteinfo.json index 34799649..f7dee889 100644 --- a/dev/.documenter-siteinfo.json +++ b/dev/.documenter-siteinfo.json @@ -1 +1 @@ -{"documenter":{"julia_version":"1.11.1","generation_timestamp":"2024-10-31T21:37:22","documenter_version":"1.7.0"}} \ No newline at end of file +{"documenter":{"julia_version":"1.11.1","generation_timestamp":"2024-11-03T23:26:16","documenter_version":"1.7.0"}} \ No newline at end of file diff --git a/dev/advection/09a42f74.gif b/dev/advection/09a42f74.gif deleted file mode 100644 index 2057b5e5..00000000 Binary files a/dev/advection/09a42f74.gif and /dev/null differ diff --git a/dev/advection/e2bfbca2.gif b/dev/advection/e2bfbca2.gif new file mode 100644 index 00000000..c10e3e03 Binary files /dev/null and b/dev/advection/e2bfbca2.gif differ diff --git a/dev/advection/index.html b/dev/advection/index.html index be9d3aa9..a3e46205 100644 --- a/dev/advection/index.html +++ b/dev/advection/index.html @@ -17,8 +17,8 @@ @parameters x domain = DomainInfo(constIC(0.0, t ∈ Interval(0, 1.0)), constBC(1.0, x ∈ Interval(0, 1.0)))

Now we convert add advection to each of the state variables. We're also adding a constant wind (ConstantWind) in the x-direction, with a speed of 1.0.

sys_advection = couple(ExampleSys(), domain, ConstantWind(t, 1.0), Advection())
 sys_mtk = convert(PDESystem, sys_advection)

\[ \begin{align} -\frac{\mathrm{d}}{\mathrm{d}t} \mathtt{ExampleSys.y}\left( t, x \right) &= \mathtt{ExampleSys.p} - \mathtt{MeanWind.v\_x}\left( t, x \right) \frac{\mathrm{d}}{\mathrm{d}x} \mathtt{ExampleSys.y}\left( t, x \right) \\ \mathtt{MeanWind.v\_x}\left( t, x \right) &= \mathtt{ConstantWind.v\_1}\left( t, x \right) \\ +\frac{\mathrm{d}}{\mathrm{d}t} \mathtt{ExampleSys.y}\left( t, x \right) &= \mathtt{ExampleSys.p} - \mathtt{MeanWind.v\_x}\left( t, x \right) \frac{\mathrm{d}}{\mathrm{d}x} \mathtt{ExampleSys.y}\left( t, x \right) \\ \mathtt{ConstantWind.v\_1}\left( t, x \right) &= \mathtt{ConstantWind.c\_v1} \end{align} \]

Finally, we can discretize the system and solve it:

using MethodOfLines, DifferentialEquations, Plots
@@ -34,4 +34,4 @@
 anim = @animate for k in 1:length(discrete_t)
     plot(discrete_x, soly[k, 1:end], title="t=\$(discrete_t[k])", ylim=(0,2.5), lab=:none)
 end
-gif(anim, fps = 8)
Example block output +gif(anim, fps = 8)Example block output diff --git a/dev/api/index.html b/dev/api/index.html index dfd873d5..8120871d 100644 --- a/dev/api/index.html +++ b/dev/api/index.html @@ -1,15 +1,15 @@ -API Reference · EarthSciMLBase.jl
GitHub

API Index

API Documentation

EarthSciMLBase.CoupledSystem — Type

A system for composing together other systems using the couple function.

  • systems: Model components to be composed together

  • domaininfo: Initial and boundary conditions and other domain information

  • pdefunctions: A vector of functions where each function takes as an argument the resulting PDESystem after DomainInfo is added to this system, and returns a transformed PDESystem.

  • ops: A vector of operators to run during simulations.

  • callbacks: A vector of callbacks to run during simulations.

  • init_callbacks: Objects x with an init_callback(x, Simulator)::DECallback method.

Things that can be added to a CoupledSystem: * ModelingToolkit.ODESystems. If the ODESystem has a field in the metadata called :coupletype (e.g. ModelingToolkit.get_metadata(sys)[:coupletype] returns a struct type with a single field called sys) then that type will be used to check for methods of EarthSciMLBase.couple that use that type. * Operators * DomainInfos * Callbacks * Types X that implement a EarthSciMLBase.init_callback(::X, sys::CoupledSystem, sys_mtk, obs_eqs, domain::DomainInfo)::DECallback method * Other CoupledSystems * Types X that implement a EarthSciMLBase.couple2(::X, ::CoupledSystem) or EarthSciMLBase.couple2(::CoupledSystem, ::X) method. * Tuples or AbstractVectors of any of the things above.

source
EarthSciMLBase.DomainInfo — Type

Domain information for a ModelingToolkit.jl PDESystem. It can be used with the + operator to add initial and boundary conditions and coordinate transforms to a ModelingToolkit.jl ODESystem or Catalyst.jl ReactionSystem.

NOTE: The independent variable (usually time) must be first in the list of initial and boundary conditions.

  • partial_derivative_funcs: Function that returns spatial derivatives of the partially-independent variables, optionally performing a coordinate transformation first.

    Current function options in this package are:

    • partialderivatives_δxyδlonlat: Returns partial derivatives after transforming any variables named lat and lon

    from degrees to cartesian meters, assuming a spherical Earth.

    Other packages may implement additional functions. They are encouraged to use function names starting with partialderivatives_.

  • grid_spacing: The discretization intervals for the partial independent variables.

  • icbc: The sets of initial and/or boundary conditions.

  • spatial_ref: The spatial reference system for the domain.

  • time_offset: The time offset for the domain.

source
EarthSciMLBase.Operator — Type

Operators are objects that modify the current state of a Simulator system. Each operator should be define a function with the signature:

`EarthSciMLBase.get_scimlop(::Operator, csys::CoupledSystem, mtk_sys, domain::DomainInfo, obs_functions, coordinate_transform_functions, u0, p)::AbstractSciMLOperator`

which should return a SciMLOperators.AbstractSciMLOperator. Refer to the SciMLOperators.jl documentation for more information on how to define operators.

source
EarthSciMLBase.SolverIMEX — Type

A solver strategy based on implicit-explicit (IMEX) time integration. See here for additional information.

kwargs for ODEProblem constructor:

  • stiff_scimlop: Whether the stiff ODE function should be implemented as a SciMLOperator.
  • stiff_sparse: Whether the stiff ODE function should use a sparse Jacobian.
  • stiff_jac: Whether the stiff ODE function should use an analytical Jacobian.
  • stiffjacscimlop: Whether the stiff ODE function Jacobian should be implemented as a SciMLOperator. (Ignored if stiff_jac==false.)
  • stiff_tgrad: Whether the stiff ODE function should use an analytical time gradient.
  • u0: initial condtions; if "nothing", default values will be used.
  • p: parameters; if "nothing", default values will be used.
source
EarthSciMLBase.SolverStrang — Type

A simulator strategy based on Strang splitting. Choose either SimulatorStrangThreads or SimulatorStrangSerial to run the simulation.

kwargs for ODEProblem constructor:

  • u0: initial condtions; if "nothing", default values will be used.
  • p: parameters; if "nothing", default values will be used.
  • nonstiff_params: parameters for the non-stiff ODE system.
  • name: name of the system.
source
EarthSciMLBase.SolverStrangSerial — Type
# Specify the stiff ODE solver algorithm.
+API Reference · EarthSciMLBase.jl
GitHub
API Index
API Documentation
EarthSciMLBase.CoupledSystem — Type
A system for composing together other systems using the couple function.
  • systems: Model components to be composed together
  • domaininfo: Initial and boundary conditions and other domain information
  • pdefunctions: A vector of functions where each function takes as an argument the resulting PDESystem after DomainInfo is added to this system, and returns a transformed PDESystem.
  • ops: A vector of operators to run during simulations.
  • callbacks: A vector of callbacks to run during simulations.
  • init_callbacks: Objects x with an init_callback(x, Simulator)::DECallback method.
Things that can be added to a CoupledSystem:     * ModelingToolkit.ODESystems. If the ODESystem has a field in the metadata called         :coupletype (e.g. ModelingToolkit.get_metadata(sys)[:coupletype] returns a struct type         with a single field called sys)         then that type will be used to check for methods of EarthSciMLBase.couple that use that type.     * Operators     * DomainInfos     * Callbacks     * Types X that implement a EarthSciMLBase.init_callback(::X, sys::CoupledSystem, sys_mtk, obs_eqs, domain::DomainInfo)::DECallback method     * Other CoupledSystems     * Types X that implement a EarthSciMLBase.couple2(::X, ::CoupledSystem) or EarthSciMLBase.couple2(::CoupledSystem, ::X) method.     * Tuples or AbstractVectors of any of the things above.
source
EarthSciMLBase.DomainInfo — Type
Domain information for a ModelingToolkit.jl PDESystem. It can be used with the + operator to add initial and boundary conditions and coordinate transforms to a ModelingToolkit.jl ODESystem or Catalyst.jl ReactionSystem.
NOTE: The independent variable (usually time) must be first in the list of initial and boundary conditions.
  • partial_derivative_funcs: Function that returns spatial derivatives of the partially-independent variables, optionally performing a coordinate transformation first.
    Current function options in this package are:
    • partialderivatives_δxyδlonlat: Returns partial derivatives after transforming any variables named lat and lon
    from degrees to cartesian meters, assuming a spherical Earth.
    Other packages may implement additional functions. They are encouraged to use function names starting with partialderivatives_.
  • grid_spacing: The discretization intervals for the partial independent variables.
  • icbc: The sets of initial and/or boundary conditions.
  • spatial_ref: The spatial reference system for the domain.
  • time_offset: The time offset for the domain.
source
EarthSciMLBase.Operator — Type
Operators are objects that modify the current state of a Simulator system. Each operator should be define a function with the signature:
`EarthSciMLBase.get_scimlop(::Operator, csys::CoupledSystem, mtk_sys, domain::DomainInfo, obs_functions, coordinate_transform_functions, u0, p)::AbstractSciMLOperator`
which should return a SciMLOperators.AbstractSciMLOperator. Refer to the SciMLOperators.jl documentation for more information on how to define operators.
source
EarthSciMLBase.SolverIMEX — Type
A solver strategy based on implicit-explicit (IMEX) time integration. See here for additional information.
kwargs:
  • stiff_scimlop: Whether the stiff ODE function should be implemented as a SciMLOperator.
  • stiff_sparse: Whether the stiff ODE function should use a sparse Jacobian.
  • stiff_jac: Whether the stiff ODE function should use an analytical Jacobian.
  • stiffjacscimlop: Whether the stiff ODE function Jacobian should be implemented as a SciMLOperator. (Ignored if stiff_jac==false.)
  • stiff_tgrad: Whether the stiff ODE function should use an analytical time gradient.
Additional kwargs for ODEProblem constructor:
  • u0: initial condtions; if "nothing", default values will be used.
  • p: parameters; if "nothing", default values will be used.
  • name: name of the model.
source
EarthSciMLBase.SolverStrang — Type
A simulator strategy based on Strang splitting. Choose either SimulatorStrangThreads or SimulatorStrangSerial to run the simulation.
kwargs for ODEProblem constructor:
  • u0: initial condtions; if "nothing", default values will be used.
  • p: parameters; if "nothing", default values will be used.
  • nonstiff_params: parameters for the non-stiff ODE system.
  • name: name of the system.
source
EarthSciMLBase.SolverStrangSerial — Type
# Specify the stiff ODE solver algorithm.
 # `timestep` is the length of time for each splitting step.
-SimulatorStrangSerial(stiffalg, timestep; kwargs...)
Perform a simulation using Strang splitting, where the MTK system is assumed to be stiff and the operators are assumed to be non-stiff. The solution will be calculated in serial.
  • stiffalg: Stiff solver algorithm to use (see https://docs.sciml.ai/DiffEqDocs/stable/solvers/ode_solve/)
  • timestep: Length of each splitting time step
  • stiff_kwargs: Keyword arguments for the stiff ODEProblem constructor and solver.
source
EarthSciMLBase.SolverStrangThreads — Type
# Specify the number of threads and the stiff ODE solver algorithm.
+SimulatorStrangSerial(stiffalg, timestep; kwargs...)
Perform a simulation using Strang splitting, where the MTK system is assumed to be stiff and the operators are assumed to be non-stiff. The solution will be calculated in serial.
Additional kwargs for ODEProblem constructor:
  • u0: initial condtions; if "nothing", default values will be used.
  • p: parameters; if "nothing", default values will be used.
  • nonstiff_params: parameters for the Operators.
  • stiffalg: Stiff solver algorithm to use (see https://docs.sciml.ai/DiffEqDocs/stable/solvers/ode_solve/)
  • timestep: Length of each splitting time step
  • stiff_kwargs: Keyword arguments for the stiff ODEProblem constructor and solver.
source
EarthSciMLBase.SolverStrangThreads — Type
# Specify the number of threads and the stiff ODE solver algorithm.
 # `timestep` is the length of time for each splitting step.
 SimulatorStrangThreads(threads, stiffalg, timestep; kwargs...)
 # Use the default number of threads.
-SimulatorStrangThreads(stiffalg, timestep; kwargs...)
Perform a simulation using Strang splitting, where the MTK system is assumed to be stiff and the operators are assumed to be non-stiff. The solution of the stiff ODE system is parallelized across grid cells using the specified number of threads.
  • threads: Number of threads to use
  • stiffalg: Stiff solver algorithm to use (see https://docs.sciml.ai/DiffEqDocs/stable/solvers/ode_solve/)
  • timestep: Length of each splitting time step
  • stiff_kwargs: Keyword arguments for the stiff ODEProblem constructor and solver.
source
EarthSciMLBase.constBC — Type
Construct constant boundary conditions equal to the value specified by val.
  • val: The value of the constant boundary conditions.
  • partialdomains: The partial domains, e.g. [x ∈ Interval(x_min, x_max), y ∈ Interval(y_min, y_max)].
source
EarthSciMLBase.constIC — Type
Construct constant initial conditions equal to the value specified by val.
  • val: The value of the constant initial conditions.
  • indepdomain: The independent domain, e.g. t ∈ Interval(t_min, t_max).
source
EarthSciMLBase.periodicBC — Type
Construct periodic boundary conditions for the given partialdomains. Periodic boundary conditions are defined as when the value at one side of the domain is set equal to the value at the other side, so that the domain "wraps around" from one side to the other.
  • partialdomains: The partial domains, e.g. [x ∈ Interval(x_min, x_max), y ∈ Interval(y_min, y_max)].
source
EarthSciMLBase.zerogradBC — Type
Construct zero-gradient boundary conditions for the given partialdomains.
  • partialdomains: The partial domains, e.g. [x ∈ Interval(x_min, x_max), y ∈ Interval(y_min, y_max)].
source
Base.convert — Method
convert(, sys; name, simplify, kwargs...)
-
Get the ODE ModelingToolkit ODESystem representation of a CoupledSystem.
kwargs:
  • name: The desired name for the resulting ODESystem
  • simplify: if true, the observed variables that are not needed to specify the state variables   will be pruned and returned as a second return value after the ODESystem, which will be   structurally simplified.
source
EarthSciMLBase.ConstantWind — Method
ConstantWind(t, vals; name)
-
Construct a constant wind velocity model component with the given wind speed(s), which should include units. For example, ConstantWind(t, 1u"m/s", 2u"m/s").
source
EarthSciMLBase.MeanWind — Method
MeanWind(t, domain)
-
A model component that represents the mean wind velocity, where pvars is the partial dependent variables for the domain.
source
EarthSciMLBase.add_dims — Method
add_dims(expression, vars, dims)
+SimulatorStrangThreads(stiffalg, timestep; kwargs...)
Perform a simulation using Strang splitting, where the MTK system is assumed to be stiff and the operators are assumed to be non-stiff. The solution of the stiff ODE system is parallelized across grid cells using the specified number of threads.
  • threads: Number of threads to use
  • stiffalg: Stiff solver algorithm to use (see https://docs.sciml.ai/DiffEqDocs/stable/solvers/ode_solve/)
  • timestep: Length of each splitting time step
  • stiff_kwargs: Keyword arguments for the stiff ODEProblem constructor and solver.
source
EarthSciMLBase.constBC — Type
Construct constant boundary conditions equal to the value specified by val.
  • val: The value of the constant boundary conditions.
  • partialdomains: The partial domains, e.g. [x ∈ Interval(x_min, x_max), y ∈ Interval(y_min, y_max)].
source
EarthSciMLBase.constIC — Type
Construct constant initial conditions equal to the value specified by val.
  • val: The value of the constant initial conditions.
  • indepdomain: The independent domain, e.g. t ∈ Interval(t_min, t_max).
source
EarthSciMLBase.periodicBC — Type
Construct periodic boundary conditions for the given partialdomains. Periodic boundary conditions are defined as when the value at one side of the domain is set equal to the value at the other side, so that the domain "wraps around" from one side to the other.
  • partialdomains: The partial domains, e.g. [x ∈ Interval(x_min, x_max), y ∈ Interval(y_min, y_max)].
source
EarthSciMLBase.zerogradBC — Type
Construct zero-gradient boundary conditions for the given partialdomains.
  • partialdomains: The partial domains, e.g. [x ∈ Interval(x_min, x_max), y ∈ Interval(y_min, y_max)].
source
Base.convert — Method
convert(, sys; name, simplify, kwargs...)
+
Get the ODE ModelingToolkit ODESystem representation of a CoupledSystem.
kwargs:
  • name: The desired name for the resulting ODESystem
  • simplify: if true, the observed variables that are not needed to specify the state variables   will be pruned and returned as a second return value after the ODESystem, which will be   structurally simplified.
source
EarthSciMLBase.ConstantWind — Method
ConstantWind(t, vals; name)
+
Construct a constant wind velocity model component with the given wind speed(s), which should include units. For example, ConstantWind(t, 1u"m/s", 2u"m/s").
source
EarthSciMLBase.MeanWind — Method
MeanWind(t, domain)
+
A model component that represents the mean wind velocity, where pvars is the partial dependent variables for the domain.
source
EarthSciMLBase.add_dims — Method
add_dims(expression, vars, dims)
 add_dims(equation, vars, dims)
Add the given dimensions to each variable in vars in the given expression or equation.  Each variable in vars must be unidimensional, i.e. defined like @variables u(t) rather than @variables u(..).
Example:
using EarthSciMLBase, ModelingToolkit
 
 @parameters x y k t
@@ -18,39 +18,41 @@
 EarthSciMLBase.add_dims(exp, [u, q], [x, y, t])
 
 # output
-1 + 2u(x, y, t) + 3k*q(x, y, t)
source
EarthSciMLBase.add_scope — Method
add_scope(sys, v, iv)
-
Add a system scope to a variable name, for example so that x in system sys1 becomes sys1₊x. iv is the independent variable.
source
EarthSciMLBase.add_scope — Method
add_scope(sys, v, iv)
+
Add a system scope to a variable name, for example so that x in system sys1 becomes sys1₊x. iv is the independent variable.
source
EarthSciMLBase.copy_with_change — Method
copy_with_change(
     sys;
     eqs,
     name,
+    unknowns,
+    parameters,
     metadata,
     continuous_events,
     discrete_events
 )
-
Create a copy of an ODESystem with the given changes.
source
EarthSciMLBase.couple — Method
couple(systems...) -> CoupledSystem
-
Couple multiple ModelingToolkit systems together.
The systems that are arguments to this system can be of type ModelingToolkit.AbstractSystem, CoupledSystem, DomainInfo, or any type T that has a method couple(::CoupledSystem, ::T)::CoupledSystem or a method couple(::T, ::CoupledSystem)::CoupledSystem defined for it.
source
EarthSciMLBase.couple — Method
couple(systems...) -> CoupledSystem
+
Couple multiple ModelingToolkit systems together.
The systems that are arguments to this system can be of type ModelingToolkit.AbstractSystem, CoupledSystem, DomainInfo, or any type T that has a method couple(::CoupledSystem, ::T)::CoupledSystem or a method couple(::T, ::CoupledSystem)::CoupledSystem defined for it.
source
EarthSciMLBase.couple2 — Method
couple2()
 
Perform bi-directional coupling for two equation systems.
To specify couplings for system pairs, create methods for this function with the signature:
EarthSciMLBase.couple2(a::ACoupler, b::BCoupler)::ConnectorSystem
where ACoupler and BCoupler are :coupletypes defined like this:
struct ACoupler sys end
-@named asys = ODESystem([], t, metadata=Dict(:coupletype=>ACoupler))
source
EarthSciMLBase.dims — Method
dims(
     icbc::EarthSciMLBase.ICcomponent
 ) -> Vector{Symbolics.Num}
-
Returns the dimensions of the independent and partial domains associated with these initial or boundary conditions.
source
EarthSciMLBase.domains — Method
domains(icbc::EarthSciMLBase.ICcomponent) -> Vector
-
Returns the domains associated with these initial or boundary conditions.
source
EarthSciMLBase.dtype — Method
dtype(_)
-
Return the data type of the state variables for this domain, based on the data types of the boundary conditions domain intervals.
source
EarthSciMLBase.get_dv — Method
Return the dependent variable, which is the first argument of the term, unless the term is a time derivative, in which case the dependent variable is the argument of the time derivative.
source
EarthSciMLBase.get_needed_vars — Method
get_needed_vars(sys)
-
Return the indexes of the system variables that the state variables of the final simplified system depend on. This should be done before running structural_simplify on the system.
source
EarthSciMLBase.grid — Method
grid(d)
-
Return the ranges representing the discretization of the partial independent variables for this domain, based on the discretization intervals given in Δs.
source
EarthSciMLBase.icbc — Method
icbc(di, states)
-
Return a vector of equations that define the initial and boundary conditions for the given state variables.
source
EarthSciMLBase.init_callback — Method
Types that implement an:
init_callback(x, sys::CoupledSystem, obs_eqs, domain::DomainInfo)::DECallback
method can also be coupled into a CoupledSystem. The init_callback function will be run before the simulator is run to get the callback.
source
EarthSciMLBase.ivar — Method
ivar(di::DomainInfo) -> Any
-
Return the independent variable associated with these initial and boundary conditions.
source
EarthSciMLBase.observed_expression — Method
observed_expression(eqs, x)
-
Return an expression for the observed value of a variable x after substituting in the constants observed values of other variables. extra_eqs is a list of additional equations to use in the substitution.
source
EarthSciMLBase.observed_function — Method
observed_function(eqs, x, coords)
-
Return a function to  for the observed value of a variable x based on the input arguments in coords. extra_eqs is a list of additional equations to use to determine the value of x.
source
EarthSciMLBase.operator_compose — Function
operator_compose(a, b)
+
Returns the dimensions of the independent and partial domains associated with these initial or boundary conditions.
source
EarthSciMLBase.domains — Method
domains(icbc::EarthSciMLBase.ICcomponent) -> Vector
+
Returns the domains associated with these initial or boundary conditions.
source
EarthSciMLBase.dtype — Method
dtype(_)
+
Return the data type of the state variables for this domain, based on the data types of the boundary conditions domain intervals.
source
EarthSciMLBase.get_dv — Method
Return the dependent variable, which is the first argument of the term, unless the term is a time derivative, in which case the dependent variable is the argument of the time derivative.
source
EarthSciMLBase.get_needed_vars — Method
get_needed_vars(sys)
+
Return the indexes of the system variables that the state variables of the final simplified system depend on. This should be done before running structural_simplify on the system.
source
EarthSciMLBase.grid — Method
grid(d)
+
Return the ranges representing the discretization of the partial independent variables for this domain, based on the discretization intervals given in Δs.
source
EarthSciMLBase.icbc — Method
icbc(di, states)
+
Return a vector of equations that define the initial and boundary conditions for the given state variables.
source
EarthSciMLBase.init_callback — Method
Types that implement an:
init_callback(x, sys::CoupledSystem, obs_eqs, domain::DomainInfo)::DECallback
method can also be coupled into a CoupledSystem. The init_callback function will be run before the simulator is run to get the callback.
source
EarthSciMLBase.ivar — Method
ivar(di::DomainInfo) -> Any
+
Return the independent variable associated with these initial and boundary conditions.
source
EarthSciMLBase.observed_expression — Method
observed_expression(eqs, x)
+
Return an expression for the observed value of a variable x after substituting in the constants observed values of other variables. extra_eqs is a list of additional equations to use in the substitution.
source
EarthSciMLBase.observed_function — Method
observed_function(eqs, x, coords)
+
Return a function to  for the observed value of a variable x based on the input arguments in coords. extra_eqs is a list of additional equations to use to determine the value of x.
source
EarthSciMLBase.operator_compose — Function
operator_compose(a, b)
 operator_compose(a, b, translate)
-
Compose to systems of equations together by adding the right-hand side terms together of equations that have matching left-hand sides. The left hand sides of two equations will be considered matching if:
  1. They are both time derivatives of the same variable.
  2. The first one is a time derivative of a variable and the second one is the variable itself.
  3. There is an entry in the optional translate dictionary that maps the dependent variable in the first system to the dependent variable in the second system, e.g. Dict(sys1.sys.x => sys2.sys.y).
  4. There is an entry in the optional translate dictionary that maps the dependent variable in the first system to the dependent variable in the second system, with a conversion factor, e.g. Dict(sys1.sys.x => sys2.sys.y => 6).
source
EarthSciMLBase.param_to_var — Method
Replace the parameter p in the system sys with a new variable that has the same name, units, and description as p.
param_to_var(sys, ps)
-
This can be useful to replace a parameter that does not change in time in a model component with one specified by another system that does change in time (or space). For example, the code below specifies a first-order loss equation, and then changes the temperature (which determines the loss rate) with a temperature value that varies in time. ```
source
EarthSciMLBase.partialderivatives_δxyδlonlat — Method
partialderivatives_δxyδlonlat(pvars; default_lat)
-
Return partial derivative operator transform factors corresponding  for the given partial-independent variables after converting variables named lon and lat from degrees to x and y meters,  assuming they represent longitude and latitude on a spherical Earth.
source
EarthSciMLBase.prune_observed — Method
prune_observed(sys)
-
Remove equations from an ODESystem where the variable in the LHS is not present in any of the equations for the state variables. This can be used to remove computationally intensive equations that are not used in the final model.
source
EarthSciMLBase.pvars — Method
pvars(di::DomainInfo) -> Any
-
Return the partial independent variables associated with these initial and boundary conditions.
source
EarthSciMLBase.steplength — Method
steplength(timesteps)
-
Return the time step length common to all of the given timesteps. Throw an error if not all timesteps are the same length.
source

+

Compose to systems of equations together by adding the right-hand side terms together of equations that have matching left-hand sides. The left hand sides of two equations will be considered matching if:

  1. They are both time derivatives of the same variable.
  2. The first one is a time derivative of a variable and the second one is the variable itself.
  3. There is an entry in the optional translate dictionary that maps the dependent variable in the first system to the dependent variable in the second system, e.g. Dict(sys1.sys.x => sys2.sys.y).
  4. There is an entry in the optional translate dictionary that maps the dependent variable in the first system to the dependent variable in the second system, with a conversion factor, e.g. Dict(sys1.sys.x => sys2.sys.y => 6).
source
EarthSciMLBase.param_to_var — Method

Replace the parameter p in the system sys with a new variable that has the same name, units, and description as p.

param_to_var(sys, ps)
+

This can be useful to replace a parameter that does not change in time in a model component with one specified by another system that does change in time (or space). For example, the code below specifies a first-order loss equation, and then changes the temperature (which determines the loss rate) with a temperature value that varies in time. ```

source
EarthSciMLBase.partialderivatives_δxyδlonlat — Method
partialderivatives_δxyδlonlat(pvars; default_lat)
+

Return partial derivative operator transform factors corresponding for the given partial-independent variables after converting variables named lon and lat from degrees to x and y meters, assuming they represent longitude and latitude on a spherical Earth.

source
EarthSciMLBase.prune_observed — Method
prune_observed(sys)
+

Remove equations from an ODESystem where the variable in the LHS is not present in any of the equations for the state variables. This can be used to remove computationally intensive equations that are not used in the final model.

source
EarthSciMLBase.pvars — Method
pvars(di::DomainInfo) -> Any
+

Return the partial independent variables associated with these initial and boundary conditions.

source
EarthSciMLBase.steplength — Method
steplength(timesteps)
+

Return the time step length common to all of the given timesteps. Throw an error if not all timesteps are the same length.

source
diff --git a/dev/composition/index.html b/dev/composition/index.html index aabf2581..874c4426 100644 --- a/dev/composition/index.html +++ b/dev/composition/index.html @@ -54,10 +54,10 @@ )) end

Finally, we can compose the model components together to create our complete model. To do so, we just initialize each model component and add the components together using the couple function. We can use the convert function to convert the composed model to a ModelingToolkit ODESystem so we can see what the combined equations look like.

model = couple(Photolysis(), Chemistry(), Emissions())
 convert(ODESystem, model)

\[ \begin{align} -\frac{\mathrm{d} \mathtt{Chemistry.NO2}\left( t \right)}{\mathrm{d}t} &= \mathtt{Chemistry.Emissions\_NO2}\left( t \right) - \mathtt{Chemistry.jNO2}\left( t \right) \mathtt{Chemistry.NO2}\left( t \right) \\ \mathtt{Chemistry.jNO2}\left( t \right) &= \mathtt{Photolysis.j\_NO2}\left( t \right) \\ \mathtt{Chemistry.Emissions\_NO2}\left( t \right) &= \mathtt{Emissions.NO2}\left( t \right) \\ \mathtt{Photolysis.j\_NO2}\left( t \right) &= max\left( \sin\left( \frac{1}{86400} t \right), 0 \right) \\ +\frac{\mathrm{d} \mathtt{Chemistry.NO2}\left( t \right)}{\mathrm{d}t} &= \mathtt{Chemistry.Emissions\_NO2}\left( t \right) - \mathtt{Chemistry.jNO2}\left( t \right) \mathtt{Chemistry.NO2}\left( t \right) \\ \mathtt{Emissions.NO2}\left( t \right) &= \mathtt{Emissions.emis} \end{align} \]

Finally, we can use the graph function to visualize the model components and their connections.

using MetaGraphsNext
@@ -68,4 +68,4 @@
 f, ax, p = graphplot(g; ilabels=labels(g))
 hidedecorations!(ax); hidespines!(ax); ax.aspect = DataAspect()
 
-f
Example block output +fExample block output diff --git a/dev/coord_transforms/index.html b/dev/coord_transforms/index.html index d672aae6..c8d9c1fb 100644 --- a/dev/coord_transforms/index.html +++ b/dev/coord_transforms/index.html @@ -36,4 +36,4 @@ \end{array} \right] \end{equation} - \]

You can see an example of how this is implemented in the source code for the Advection model component.

Additional transformation functions may be defined in other packages. We recommend that the names of these functions start with partialderivatives_ to make it clear that they are intended to be used in this context.

+ \]

You can see an example of how this is implemented in the source code for the Advection model component.

Additional transformation functions may be defined in other packages. We recommend that the names of these functions start with partialderivatives_ to make it clear that they are intended to be used in this context.

diff --git a/dev/example_all_together/76a8d88f.svg b/dev/example_all_together/76a8d88f.svg new file mode 100644 index 00000000..a874ff4a --- /dev/null +++ b/dev/example_all_together/76a8d88f.svg @@ -0,0 +1,50 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/dev/example_all_together/efbbc156.svg b/dev/example_all_together/efbbc156.svg deleted file mode 100644 index d53b199e..00000000 --- a/dev/example_all_together/efbbc156.svg +++ /dev/null @@ -1,50 +0,0 @@ - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - diff --git a/dev/example_all_together/index.html b/dev/example_all_together/index.html index 7ba19bef..0c16d048 100644 --- a/dev/example_all_together/index.html +++ b/dev/example_all_together/index.html @@ -44,12 +44,12 @@ sys = couple(sys1, sys2) convert(ODESystem, sys)

\[ \begin{align} -\frac{\mathrm{d} \mathtt{Sys1.c_1}\left( t \right)}{\mathrm{d}t} &= \mathtt{Sys1.Sys2\_ddt\_c_{1}ˍt}\left( t \right) - 2 \mathtt{Sys1.c_1}\left( t \right) \\ -\frac{\mathrm{d} \mathtt{Sys1.c_2}\left( t \right)}{\mathrm{d}t} &= \mathtt{Sys1.Sys2\_ddt\_c_{2}ˍt}\left( t \right) + 2 \mathtt{Sys1.c_1}\left( t \right) \\ \mathtt{Sys1.Sys2\_ddt\_c_{1}ˍt}\left( t \right) &= \mathtt{Sys2.Sys2\_ddt\_c_{1}ˍt}\left( t \right) \\ \mathtt{Sys1.c_1}\left( t \right) &= \mathtt{Sys2.c_1}\left( t \right) \\ \mathtt{Sys1.Sys2\_ddt\_c_{2}ˍt}\left( t \right) &= \mathtt{Sys2.Sys2\_ddt\_c_{2}ˍt}\left( t \right) \\ \mathtt{Sys1.c_2}\left( t \right) &= \mathtt{Sys2.c_2}\left( t \right) \\ +\frac{\mathrm{d} \mathtt{Sys1.c_1}\left( t \right)}{\mathrm{d}t} &= \mathtt{Sys1.Sys2\_ddt\_c_{1}ˍt}\left( t \right) - 2 \mathtt{Sys1.c_1}\left( t \right) \\ +\frac{\mathrm{d} \mathtt{Sys1.c_2}\left( t \right)}{\mathrm{d}t} &= \mathtt{Sys1.Sys2\_ddt\_c_{2}ˍt}\left( t \right) + 2 \mathtt{Sys1.c_1}\left( t \right) \\ \mathtt{Sys2.Sys2\_ddt\_c_{1}ˍt}\left( t \right) &= - \mathtt{Sys2.p{_1}} \\ \mathtt{Sys2.Sys2\_ddt\_c_{2}ˍt}\left( t \right) &= \mathtt{Sys2.p{_2}} \end{align} @@ -59,15 +59,15 @@ \end{align} \]

observed(simplified_sys)

\[ \begin{align} \mathtt{Sys2.Sys2\_ddt\_c_{1}ˍt}\left( t \right) &= - \mathtt{Sys2.p{_1}} \\ -\mathtt{Sys2.Sys2\_ddt\_c_{2}ˍt}\left( t \right) &= \mathtt{Sys2.p{_2}} \\ \mathtt{Sys2.c_1}\left( t \right) &= \mathtt{Sys1.c_1}\left( t \right) \\ +\mathtt{Sys2.Sys2\_ddt\_c_{2}ˍt}\left( t \right) &= \mathtt{Sys2.p{_2}} \\ \mathtt{Sys2.c_2}\left( t \right) &= \mathtt{Sys1.c_2}\left( t \right) \\ \mathtt{Sys1.Sys2\_ddt\_c_{1}ˍt}\left( t \right) &= \mathtt{Sys2.Sys2\_ddt\_c_{1}ˍt}\left( t \right) \\ \mathtt{Sys1.Sys2\_ddt\_c_{2}ˍt}\left( t \right) &= \mathtt{Sys2.Sys2\_ddt\_c_{2}ˍt}\left( t \right) \end{align} \]

We can also run simulations using this system:

odeprob = ODEProblem(simplified_sys, [], (0.0,10.0), [])
 odesol = solve(odeprob)
-plot(odesol)
Example block output

Once we've confirmed that our model works in a 0D "box model" setting, we can expand it to 1, 2, or 3 dimensions using by adding in initial and boundary conditions. We will also add in advection using constant-velocity wind fields add the same time.

x_min = y_min = t_min = 0.0
+plot(odesol)
Example block output

Once we've confirmed that our model works in a 0D "box model" setting, we can expand it to 1, 2, or 3 dimensions using by adding in initial and boundary conditions. We will also add in advection using constant-velocity wind fields add the same time.

x_min = y_min = t_min = 0.0
 x_max = y_max = t_max = 1.0
 domain = DomainInfo(
     constIC(4.0, t ∈ Interval(t_min, t_max)),
@@ -78,14 +78,14 @@
 sys_pde = couple(sys, domain, ConstantWind(t, 1.0, 1.0), Advection())
 
 sys_pde_mtk = convert(PDESystem, sys_pde)

\[ \begin{align} -\frac{\mathrm{d}}{\mathrm{d}t} \mathtt{Sys1.c_1}\left( t, x, y \right) &= - 2 \mathtt{Sys1.c_1}\left( t, x, y \right) + \mathtt{Sys1.Sys2\_ddt\_c_{1}ˍt}\left( t, x, y \right) - \frac{\mathrm{d}}{\mathrm{d}y} \mathtt{Sys1.c_1}\left( t, x, y \right) \mathtt{MeanWind.v\_y}\left( t, x, y \right) - \frac{\mathrm{d}}{\mathrm{d}x} \mathtt{Sys1.c_1}\left( t, x, y \right) \mathtt{MeanWind.v\_x}\left( t, x, y \right) \\ -\frac{\mathrm{d}}{\mathrm{d}t} \mathtt{Sys1.c_2}\left( t, x, y \right) &= 2 \mathtt{Sys1.c_1}\left( t, x, y \right) + \mathtt{Sys1.Sys2\_ddt\_c_{2}ˍt}\left( t, x, y \right) - \mathtt{MeanWind.v\_y}\left( t, x, y \right) \frac{\mathrm{d}}{\mathrm{d}y} \mathtt{Sys1.c_2}\left( t, x, y \right) - \mathtt{MeanWind.v\_x}\left( t, x, y \right) \frac{\mathrm{d}}{\mathrm{d}x} \mathtt{Sys1.c_2}\left( t, x, y \right) \\ \mathtt{Sys1.Sys2\_ddt\_c_{1}ˍt}\left( t, x, y \right) &= \mathtt{Sys2.Sys2\_ddt\_c_{1}ˍt}\left( t, x, y \right) \\ \mathtt{Sys1.c_1}\left( t, x, y \right) &= \mathtt{Sys2.c_1}\left( t, x, y \right) \\ \mathtt{Sys1.Sys2\_ddt\_c_{2}ˍt}\left( t, x, y \right) &= \mathtt{Sys2.Sys2\_ddt\_c_{2}ˍt}\left( t, x, y \right) \\ \mathtt{Sys1.c_2}\left( t, x, y \right) &= \mathtt{Sys2.c_2}\left( t, x, y \right) \\ \mathtt{MeanWind.v\_x}\left( t, x, y \right) &= \mathtt{ConstantWind.v\_1}\left( t, x, y \right) \\ \mathtt{MeanWind.v\_y}\left( t, x, y \right) &= \mathtt{ConstantWind.v\_2}\left( t, x, y \right) \\ +\frac{\mathrm{d}}{\mathrm{d}t} \mathtt{Sys1.c_1}\left( t, x, y \right) &= - 2 \mathtt{Sys1.c_1}\left( t, x, y \right) + \mathtt{Sys1.Sys2\_ddt\_c_{1}ˍt}\left( t, x, y \right) - \frac{\mathrm{d}}{\mathrm{d}y} \mathtt{Sys1.c_1}\left( t, x, y \right) \mathtt{MeanWind.v\_y}\left( t, x, y \right) - \frac{\mathrm{d}}{\mathrm{d}x} \mathtt{Sys1.c_1}\left( t, x, y \right) \mathtt{MeanWind.v\_x}\left( t, x, y \right) \\ +\frac{\mathrm{d}}{\mathrm{d}t} \mathtt{Sys1.c_2}\left( t, x, y \right) &= 2 \mathtt{Sys1.c_1}\left( t, x, y \right) + \mathtt{Sys1.Sys2\_ddt\_c_{2}ˍt}\left( t, x, y \right) - \mathtt{MeanWind.v\_y}\left( t, x, y \right) \frac{\mathrm{d}}{\mathrm{d}y} \mathtt{Sys1.c_2}\left( t, x, y \right) - \mathtt{MeanWind.v\_x}\left( t, x, y \right) \frac{\mathrm{d}}{\mathrm{d}x} \mathtt{Sys1.c_2}\left( t, x, y \right) \\ \mathtt{Sys2.Sys2\_ddt\_c_{1}ˍt}\left( t, x, y \right) &= - \mathtt{Sys2.p{_1}} \\ \mathtt{Sys2.Sys2\_ddt\_c_{2}ˍt}\left( t, x, y \right) &= \mathtt{Sys2.p{_2}} \\ \mathtt{ConstantWind.v\_1}\left( t, x, y \right) &= \mathtt{ConstantWind.c\_v1} \\ @@ -94,13 +94,13 @@ \]

Now we can inspect this new system that we've created:

sys_pde_mtk.dvs

\[ \begin{equation} \left[ \begin{array}{c} -\mathtt{Sys1.c_1}\left( t, x, y \right) \\ -\mathtt{Sys1.c_2}\left( t, x, y \right) \\ \mathtt{Sys1.Sys2\_ddt\_c_{1}ˍt}\left( t, x, y \right) \\ \mathtt{Sys2.Sys2\_ddt\_c_{1}ˍt}\left( t, x, y \right) \\ +\mathtt{Sys1.c_1}\left( t, x, y \right) \\ \mathtt{Sys2.c_1}\left( t, x, y \right) \\ \mathtt{Sys1.Sys2\_ddt\_c_{2}ˍt}\left( t, x, y \right) \\ \mathtt{Sys2.Sys2\_ddt\_c_{2}ˍt}\left( t, x, y \right) \\ +\mathtt{Sys1.c_2}\left( t, x, y \right) \\ \mathtt{Sys2.c_2}\left( t, x, y \right) \\ \mathtt{MeanWind.v\_x}\left( t, x, y \right) \\ \mathtt{ConstantWind.v\_1}\left( t, x, y \right) \\ @@ -132,4 +132,4 @@ p2 = heatmap(solc2[k, 1:end-1, 1:end-1], title="c₂ t=\$(discrete_t[k])", clim=(0,7.0), lab=:none) plot(p1, p2, layout=(1,2), size=(800,400)) end -gif(anim, fps = 8)Example block output

Because our system is a system of ordinary differential equations rather than partial differential equations, all of the grid cells in the animation above have the same value. Refer to the advection example for an example of a system of partial differential equations.

+gif(anim, fps = 8)Example block output

Because our system is a system of ordinary differential equations rather than partial differential equations, all of the grid cells in the animation above have the same value. Refer to the advection example for an example of a system of partial differential equations.

diff --git a/dev/icbc/index.html b/dev/icbc/index.html index 829b9627..587128dc 100644 --- a/dev/icbc/index.html +++ b/dev/icbc/index.html @@ -50,4 +50,4 @@ \mathtt{Docs.Example.v}\left( t, x, 0 \right) &= 16 \\ \mathtt{Docs.Example.v}\left( t, x, 1 \right) &= 16 \end{align} - \]

+ \]

diff --git a/dev/index.html b/dev/index.html index 49befd5b..c2fe3d42 100644 --- a/dev/index.html +++ b/dev/index.html @@ -6,14 +6,14 @@ [0c46a032] DifferentialEquations v7.14.0 [e30172f5] Documenter v1.7.0 [5b8099bc] DomainSets v0.7.14 - [06fc5a27] DynamicQuantities v1.1.0 + [06fc5a27] DynamicQuantities v1.2.0 [e53f1632] EarthSciMLBase v0.19.2 `~/work/EarthSciMLBase.jl/EarthSciMLBase.jl` [1ecd5474] GraphMakie v0.5.12 [fa8bd995] MetaGraphsNext v0.7.1 [94925ecb] MethodOfLines v0.11.6 [961ee093] ModelingToolkit v9.49.0 [91a5bcdd] Plots v1.40.8 - [c0aeaf25] SciMLOperators v0.3.11 + [c0aeaf25] SciMLOperators v0.3.12 [0c5d862f] Symbolics v6.17.0
and using this machine and Julia version.
Julia Version 1.11.1
 Commit 8f5b7ca12ad (2024-10-16 10:53 UTC)
 Build Info:
@@ -100,7 +100,7 @@
   [e30172f5] Documenter v1.7.0
   [5b8099bc] DomainSets v0.7.14
   [7c1d4256] DynamicPolynomials v0.6.0
-  [06fc5a27] DynamicQuantities v1.1.0
+  [06fc5a27] DynamicQuantities v1.2.0
   [e53f1632] EarthSciMLBase v0.19.2 `~/work/EarthSciMLBase.jl/EarthSciMLBase.jl`
   [4e289a0a] EnumX v1.0.4
   [f151be2c] EnzymeCore v0.8.5
@@ -125,7 +125,7 @@
   [6a86dc24] FiniteDiff v2.26.0
   [53c48c17] FixedPointNumbers v0.8.5
   [1fa38f19] Format v1.3.7
-  [f6369f11] ForwardDiff v0.10.36
+  [f6369f11] ForwardDiff v0.10.37
   [b38be410] FreeType v4.1.1
   [663a7486] FreeTypeAbstraction v0.10.4
   [069b7b12] FunctionWrappers v1.1.3
@@ -179,7 +179,7 @@
   [b964fa9f] LaTeXStrings v1.4.0
   [23fbe1c1] Latexify v0.16.5
   [10f19ff3] LayoutPointers v0.1.17
-  [0e77f7df] LazilyInitializedFields v1.2.2
+  [0e77f7df] LazilyInitializedFields v1.3.0
   [5078a376] LazyArrays v2.2.1
   [8cdb02fc] LazyModules v0.3.1
   [2d8b4e74] LevyArea v1.0.0
@@ -214,7 +214,7 @@
   [77ba4419] NaNMath v1.0.2
   [f09324ee] Netpbm v1.1.1
   [46757867] NetworkLayout v0.4.7
-  [8913a72c] NonlinearSolve v3.15.1
+⌅ [8913a72c] NonlinearSolve v3.15.1
   [510215fc] Observables v0.5.5
   [6fe1bfb0] OffsetArrays v1.14.1
   [52e1d378] OpenEXR v0.3.2
@@ -263,7 +263,7 @@
   [b98c9c47] Pipe v1.3.0
   [eebad327] PkgVersion v0.3.3
   [ccf2f8ad] PlotThemes v3.3.0
-  [995b91a9] PlotUtils v1.4.2
+  [995b91a9] PlotUtils v1.4.3
   [91a5bcdd] Plots v1.40.8
   [e409e4f3] PoissonRandom v0.4.4
   [f517fe37] Polyester v0.7.16
@@ -302,7 +302,7 @@
   [476501e8] SLEEFPirates v0.6.43
   [0bca4576] SciMLBase v2.58.0
   [19f34311] SciMLJacobianOperators v0.1.1
-  [c0aeaf25] SciMLOperators v0.3.11
+  [c0aeaf25] SciMLOperators v0.3.12
   [53ae85a6] SciMLStructures v1.5.0
   [6c6a2e73] Scratch v1.2.1
   [efcf1570] Setfield v1.1.1
@@ -310,7 +310,7 @@
   [992d4aef] Showoff v1.0.3
   [73760f76] SignedDistanceFields v0.4.0
   [777ac1f9] SimpleBufferStream v1.2.0
-  [727e6d20] SimpleNonlinearSolve v1.12.3
+⌅ [727e6d20] SimpleNonlinearSolve v1.12.3
   [699a6c99] SimpleTraits v0.9.4
   [ce78b400] SimpleUnPack v1.1.0
   [45858cf5] Sixel v0.1.3
@@ -514,4 +514,4 @@
   [3f19e933] p7zip_jll v17.4.0+2
 Info Packages marked with ⌅ have new versions available but compatibility constraints restrict them from upgrading. To see why use `status --outdated -m`
You can also download the manifest file and the -project file. +project file. diff --git a/dev/operator/index.html b/dev/operator/index.html index 710cad41..b136dfa2 100644 --- a/dev/operator/index.html +++ b/dev/operator/index.html @@ -76,4 +76,4 @@ heatmap(u[1, :, :, 1]), ) end -gif(anim, fps = 15)Example block output +gif(anim, fps = 15)Example block output diff --git a/dev/operator_compose/index.html b/dev/operator_compose/index.html index d474774a..417def88 100644 --- a/dev/operator_compose/index.html +++ b/dev/operator_compose/index.html @@ -39,9 +39,9 @@ combined = couple(sys1, sys2) combined_mtk = convert(ODESystem, combined)

\[ \begin{align} -\frac{\mathrm{d} \mathtt{ExampleSys.x}\left( t \right)}{\mathrm{d}t} &= \mathtt{ExampleSys.p} + \mathtt{ExampleSys.ExampleSys2\_ddt\_xˍt}\left( t \right) \\ \mathtt{ExampleSys.ExampleSys2\_ddt\_xˍt}\left( t \right) &= \mathtt{ExampleSys2.ExampleSys2\_ddt\_xˍt}\left( t \right) \\ \mathtt{ExampleSys.x}\left( t \right) &= \mathtt{ExampleSys2.x}\left( t \right) \\ +\frac{\mathrm{d} \mathtt{ExampleSys.x}\left( t \right)}{\mathrm{d}t} &= \mathtt{ExampleSys.p} + \mathtt{ExampleSys.ExampleSys2\_ddt\_xˍt}\left( t \right) \\ \mathtt{ExampleSys2.ExampleSys2\_ddt\_xˍt}\left( t \right) &= 2 \mathtt{ExampleSys2.p} \end{align} \]

The simplified equation should be D(x) = p + sys2_xˍt:

combined_simplified = structural_simplify(combined_mtk)

\[ \begin{align} @@ -125,4 +125,4 @@ \mathtt{Docs.ExampleSysNonODE2.y}\left( t \right) &= \mathtt{Docs.ExampleSysNonODE2.p} \\ \mathtt{ExampleSys.Docs.ExampleSysNonODE2\_y}\left( t \right) &= 6 \mathtt{Docs.ExampleSysNonODE2.y}\left( t \right) \end{align} - \]

+ \]

diff --git a/dev/param_to_var/index.html b/dev/param_to_var/index.html index ec1d81d5..b1510d55 100644 --- a/dev/param_to_var/index.html +++ b/dev/param_to_var/index.html @@ -36,8 +36,8 @@ variable_loss = couple(l, temp) convert(ODESystem, variable_loss)

\[ \begin{align} +\mathtt{Loss.T}\left( t \right) &= \mathtt{Temperature.T}\left( t \right) \\ \frac{\mathrm{d} \mathtt{Loss.A}\left( t \right)}{\mathrm{d}t} &= - \mathtt{Loss.k} e^{\frac{\mathtt{Loss.T}\left( t \right)}{\mathtt{Loss.T{_0}}}} \mathtt{Loss.A}\left( t \right) \\ -\frac{\mathrm{d} \mathtt{Temperature.T}\left( t \right)}{\mathrm{d}t} &= \mathtt{Temperature.Tc} \sin\left( \frac{t}{\mathtt{Temperature.tc}} \right) \\ -\mathtt{Loss.T}\left( t \right) &= \mathtt{Temperature.T}\left( t \right) +\frac{\mathrm{d} \mathtt{Temperature.T}\left( t \right)}{\mathrm{d}t} &= \mathtt{Temperature.Tc} \sin\left( \frac{t}{\mathtt{Temperature.tc}} \right) \end{align} - \]

If we wanted to, we could then run a simulation with the composed system.

+ \]

If we wanted to, we could then run a simulation with the composed system.

diff --git a/dev/search_index.js b/dev/search_index.js index f3d8d60a..bc3fb321 100644 --- a/dev/search_index.js +++ b/dev/search_index.js @@ -1,3 +1,3 @@ var documenterSearchIndex = {"docs": -[{"location":"api/#API-Index","page":"API Reference","title":"API Index","text":"","category":"section"},{"location":"api/","page":"API Reference","title":"API Reference","text":"","category":"page"},{"location":"api/#API-Documentation","page":"API Reference","title":"API Documentation","text":"","category":"section"},{"location":"api/","page":"API Reference","title":"API Reference","text":"Modules = [EarthSciMLBase]","category":"page"},{"location":"api/#EarthSciMLBase.Advection","page":"API Reference","title":"EarthSciMLBase.Advection","text":"Apply advection to a model.\n\n\n\n\n\n","category":"type"},{"location":"api/#EarthSciMLBase.ConnectorSystem","page":"API Reference","title":"EarthSciMLBase.ConnectorSystem","text":"A connector for two systems.\n\neqs\nfrom\nto\n\n\n\n\n\n","category":"type"},{"location":"api/#EarthSciMLBase.CoupledSystem","page":"API Reference","title":"EarthSciMLBase.CoupledSystem","text":"A system for composing together other systems using the couple function.\n\nsystems: Model components to be composed together\ndomaininfo: Initial and boundary conditions and other domain information\npdefunctions: A vector of functions where each function takes as an argument the resulting PDESystem after DomainInfo is added to this system, and returns a transformed PDESystem.\n\nops: A vector of operators to run during simulations.\n\ncallbacks: A vector of callbacks to run during simulations.\ninit_callbacks: Objects x with an init_callback(x, Simulator)::DECallback method.\n\nThings that can be added to a CoupledSystem: * ModelingToolkit.ODESystems. If the ODESystem has a field in the metadata called :coupletype (e.g. ModelingToolkit.get_metadata(sys)[:coupletype] returns a struct type with a single field called sys) then that type will be used to check for methods of EarthSciMLBase.couple that use that type. * Operators * DomainInfos * Callbacks * Types X that implement a EarthSciMLBase.init_callback(::X, sys::CoupledSystem, sys_mtk, obs_eqs, domain::DomainInfo)::DECallback method * Other CoupledSystems * Types X that implement a EarthSciMLBase.couple2(::X, ::CoupledSystem) or EarthSciMLBase.couple2(::CoupledSystem, ::X) method. * Tuples or AbstractVectors of any of the things above.\n\n\n\n\n\n","category":"type"},{"location":"api/#EarthSciMLBase.DomainInfo","page":"API Reference","title":"EarthSciMLBase.DomainInfo","text":"Domain information for a ModelingToolkit.jl PDESystem. It can be used with the + operator to add initial and boundary conditions and coordinate transforms to a ModelingToolkit.jl ODESystem or Catalyst.jl ReactionSystem.\n\nNOTE: The independent variable (usually time) must be first in the list of initial and boundary conditions.\n\npartial_derivative_funcs: Function that returns spatial derivatives of the partially-independent variables, optionally performing a coordinate transformation first.\nCurrent function options in this package are:\npartialderivatives_δxyδlonlat: Returns partial derivatives after transforming any variables named lat and lon\nfrom degrees to cartesian meters, assuming a spherical Earth.\nOther packages may implement additional functions. They are encouraged to use function names starting with partialderivatives_.\n\ngrid_spacing: The discretization intervals for the partial independent variables.\nicbc: The sets of initial and/or boundary conditions.\nspatial_ref: The spatial reference system for the domain.\ntime_offset: The time offset for the domain.\n\n\n\n\n\n","category":"type"},{"location":"api/#EarthSciMLBase.ICBCcomponent","page":"API Reference","title":"EarthSciMLBase.ICBCcomponent","text":"Initial and boundary condition components that can be combined to create an DomainInfo object.\n\n\n\n\n\n","category":"type"},{"location":"api/#EarthSciMLBase.Operator","page":"API Reference","title":"EarthSciMLBase.Operator","text":"Operators are objects that modify the current state of a Simulator system. Each operator should be define a function with the signature:\n\n`EarthSciMLBase.get_scimlop(::Operator, csys::CoupledSystem, mtk_sys, domain::DomainInfo, obs_functions, coordinate_transform_functions, u0, p)::AbstractSciMLOperator`\n\nwhich should return a SciMLOperators.AbstractSciMLOperator. Refer to the SciMLOperators.jl documentation for more information on how to define operators.\n\n\n\n\n\n","category":"type"},{"location":"api/#EarthSciMLBase.SolverIMEX","page":"API Reference","title":"EarthSciMLBase.SolverIMEX","text":"A solver strategy based on implicit-explicit (IMEX) time integration. See here for additional information.\n\nkwargs for ODEProblem constructor:\n\nstiff_scimlop: Whether the stiff ODE function should be implemented as a SciMLOperator.\nstiff_sparse: Whether the stiff ODE function should use a sparse Jacobian.\nstiff_jac: Whether the stiff ODE function should use an analytical Jacobian.\nstiffjacscimlop: Whether the stiff ODE function Jacobian should be implemented as a SciMLOperator. (Ignored if stiff_jac==false.)\nstiff_tgrad: Whether the stiff ODE function should use an analytical time gradient.\nu0: initial condtions; if \"nothing\", default values will be used.\np: parameters; if \"nothing\", default values will be used.\n\n\n\n\n\n","category":"type"},{"location":"api/#EarthSciMLBase.SolverStrang","page":"API Reference","title":"EarthSciMLBase.SolverStrang","text":"A simulator strategy based on Strang splitting. Choose either SimulatorStrangThreads or SimulatorStrangSerial to run the simulation.\n\nkwargs for ODEProblem constructor:\n\nu0: initial condtions; if \"nothing\", default values will be used.\np: parameters; if \"nothing\", default values will be used.\nnonstiff_params: parameters for the non-stiff ODE system.\nname: name of the system.\n\n\n\n\n\n","category":"type"},{"location":"api/#EarthSciMLBase.SolverStrangSerial","page":"API Reference","title":"EarthSciMLBase.SolverStrangSerial","text":"# Specify the stiff ODE solver algorithm.\n# `timestep` is the length of time for each splitting step.\nSimulatorStrangSerial(stiffalg, timestep; kwargs...)\n\nPerform a simulation using Strang splitting, where the MTK system is assumed to be stiff and the operators are assumed to be non-stiff. The solution will be calculated in serial.\n\nstiffalg: Stiff solver algorithm to use (see https://docs.sciml.ai/DiffEqDocs/stable/solvers/ode_solve/)\ntimestep: Length of each splitting time step\nstiff_kwargs: Keyword arguments for the stiff ODEProblem constructor and solver.\n\n\n\n\n\n","category":"type"},{"location":"api/#EarthSciMLBase.SolverStrangThreads","page":"API Reference","title":"EarthSciMLBase.SolverStrangThreads","text":"# Specify the number of threads and the stiff ODE solver algorithm.\n# `timestep` is the length of time for each splitting step.\nSimulatorStrangThreads(threads, stiffalg, timestep; kwargs...)\n# Use the default number of threads.\nSimulatorStrangThreads(stiffalg, timestep; kwargs...)\n\nPerform a simulation using Strang splitting, where the MTK system is assumed to be stiff and the operators are assumed to be non-stiff. The solution of the stiff ODE system is parallelized across grid cells using the specified number of threads.\n\nthreads: Number of threads to use\nstiffalg: Stiff solver algorithm to use (see https://docs.sciml.ai/DiffEqDocs/stable/solvers/ode_solve/)\ntimestep: Length of each splitting time step\nstiff_kwargs: Keyword arguments for the stiff ODEProblem constructor and solver.\n\n\n\n\n\n","category":"type"},{"location":"api/#EarthSciMLBase.SolverStrategy","page":"API Reference","title":"EarthSciMLBase.SolverStrategy","text":"SolverStrategy is an abstract type that defines the strategy for running a simulation.\n\n\n\n\n\n","category":"type"},{"location":"api/#EarthSciMLBase.constBC","page":"API Reference","title":"EarthSciMLBase.constBC","text":"Construct constant boundary conditions equal to the value specified by val.\n\nval: The value of the constant boundary conditions.\npartialdomains: The partial domains, e.g. [x ∈ Interval(x_min, x_max), y ∈ Interval(y_min, y_max)].\n\n\n\n\n\n","category":"type"},{"location":"api/#EarthSciMLBase.constIC","page":"API Reference","title":"EarthSciMLBase.constIC","text":"Construct constant initial conditions equal to the value specified by val.\n\nval: The value of the constant initial conditions.\nindepdomain: The independent domain, e.g. t ∈ Interval(t_min, t_max).\n\n\n\n\n\n","category":"type"},{"location":"api/#EarthSciMLBase.periodicBC","page":"API Reference","title":"EarthSciMLBase.periodicBC","text":"Construct periodic boundary conditions for the given partialdomains. Periodic boundary conditions are defined as when the value at one side of the domain is set equal to the value at the other side, so that the domain \"wraps around\" from one side to the other.\n\npartialdomains: The partial domains, e.g. [x ∈ Interval(x_min, x_max), y ∈ Interval(y_min, y_max)].\n\n\n\n\n\n","category":"type"},{"location":"api/#EarthSciMLBase.zerogradBC","page":"API Reference","title":"EarthSciMLBase.zerogradBC","text":"Construct zero-gradient boundary conditions for the given partialdomains.\n\npartialdomains: The partial domains, e.g. [x ∈ Interval(x_min, x_max), y ∈ Interval(y_min, y_max)].\n\n\n\n\n\n","category":"type"},{"location":"api/#Base.convert-Tuple{Type{<:ModelingToolkit.ODESystem}, CoupledSystem}","page":"API Reference","title":"Base.convert","text":"convert(, sys; name, simplify, kwargs...)\n\n\nGet the ODE ModelingToolkit ODESystem representation of a CoupledSystem.\n\nkwargs:\n\nname: The desired name for the resulting ODESystem\nsimplify: if true, the observed variables that are not needed to specify the state variables will be pruned and returned as a second return value after the ODESystem, which will be structurally simplified.\n\n\n\n\n\n","category":"method"},{"location":"api/#Base.convert-Tuple{Type{<:ModelingToolkit.PDESystem}, CoupledSystem}","page":"API Reference","title":"Base.convert","text":"convert(, sys; name, kwargs...)\n\n\nGet the ModelingToolkit PDESystem representation of a CoupledSystem.\n\n\n\n\n\n","category":"method"},{"location":"api/#EarthSciMLBase.ConstantWind-Tuple{Any, Vararg{Any}}","page":"API Reference","title":"EarthSciMLBase.ConstantWind","text":"ConstantWind(t, vals; name)\n\n\nConstruct a constant wind velocity model component with the given wind speed(s), which should include units. For example, ConstantWind(t, 1u\"m/s\", 2u\"m/s\").\n\n\n\n\n\n","category":"method"},{"location":"api/#EarthSciMLBase.MeanWind-Tuple{Any, DomainInfo}","page":"API Reference","title":"EarthSciMLBase.MeanWind","text":"MeanWind(t, domain)\n\n\nA model component that represents the mean wind velocity, where pvars is the partial dependent variables for the domain.\n\n\n\n\n\n","category":"method"},{"location":"api/#EarthSciMLBase.add_dims-Tuple{Any, AbstractVector, AbstractVector}","page":"API Reference","title":"EarthSciMLBase.add_dims","text":"add_dims(expression, vars, dims)\nadd_dims(equation, vars, dims)\n\nAdd the given dimensions to each variable in vars in the given expression or equation. Each variable in vars must be unidimensional, i.e. defined like @variables u(t) rather than @variables u(..).\n\nExample:\n\nusing EarthSciMLBase, ModelingToolkit\n\n@parameters x y k t\n@variables u(t) q(t)\nexp = 2u + 3k*q + 1\nEarthSciMLBase.add_dims(exp, [u, q], [x, y, t])\n\n# output\n1 + 2u(x, y, t) + 3k*q(x, y, t)\n\n\n\n\n\n","category":"method"},{"location":"api/#EarthSciMLBase.add_metadata-Tuple{Any, Any}","page":"API Reference","title":"EarthSciMLBase.add_metadata","text":"Add the units and description in the variable from to the variable to.\n\n\n\n\n\n","category":"method"},{"location":"api/#EarthSciMLBase.add_scope-Tuple{Any, Any, Any}","page":"API Reference","title":"EarthSciMLBase.add_scope","text":"add_scope(sys, v, iv)\n\n\nAdd a system scope to a variable name, for example so that x in system sys1 becomes sys1₊x. iv is the independent variable.\n\n\n\n\n\n","category":"method"},{"location":"api/#EarthSciMLBase.copy_with_change-Tuple{ModelingToolkit.ODESystem}","page":"API Reference","title":"EarthSciMLBase.copy_with_change","text":"copy_with_change(\n sys;\n eqs,\n name,\n metadata,\n continuous_events,\n discrete_events\n)\n\n\nCreate a copy of an ODESystem with the given changes.\n\n\n\n\n\n","category":"method"},{"location":"api/#EarthSciMLBase.couple-Tuple","page":"API Reference","title":"EarthSciMLBase.couple","text":"couple(systems...) -> CoupledSystem\n\n\nCouple multiple ModelingToolkit systems together.\n\nThe systems that are arguments to this system can be of type ModelingToolkit.AbstractSystem, CoupledSystem, DomainInfo, or any type T that has a method couple(::CoupledSystem, ::T)::CoupledSystem or a method couple(::T, ::CoupledSystem)::CoupledSystem defined for it.\n\n\n\n\n\n","category":"method"},{"location":"api/#EarthSciMLBase.couple2-Tuple{}","page":"API Reference","title":"EarthSciMLBase.couple2","text":"couple2()\n\n\nPerform bi-directional coupling for two equation systems.\n\nTo specify couplings for system pairs, create methods for this function with the signature:\n\nEarthSciMLBase.couple2(a::ACoupler, b::BCoupler)::ConnectorSystem\n\nwhere ACoupler and BCoupler are :coupletypes defined like this:\n\nstruct ACoupler sys end\n@named asys = ODESystem([], t, metadata=Dict(:coupletype=>ACoupler))\n\n\n\n\n\n","category":"method"},{"location":"api/#EarthSciMLBase.dims-Tuple{EarthSciMLBase.ICcomponent}","page":"API Reference","title":"EarthSciMLBase.dims","text":"dims(\n icbc::EarthSciMLBase.ICcomponent\n) -> Vector{Symbolics.Num}\n\n\nReturns the dimensions of the independent and partial domains associated with these initial or boundary conditions.\n\n\n\n\n\n","category":"method"},{"location":"api/#EarthSciMLBase.domains-Tuple{EarthSciMLBase.ICcomponent}","page":"API Reference","title":"EarthSciMLBase.domains","text":"domains(icbc::EarthSciMLBase.ICcomponent) -> Vector\n\n\nReturns the domains associated with these initial or boundary conditions.\n\n\n\n\n\n","category":"method"},{"location":"api/#EarthSciMLBase.dtype-Union{Tuple{DomainInfo{T}}, Tuple{T}} where T","page":"API Reference","title":"EarthSciMLBase.dtype","text":"dtype(_)\n\n\nReturn the data type of the state variables for this domain, based on the data types of the boundary conditions domain intervals.\n\n\n\n\n\n","category":"method"},{"location":"api/#EarthSciMLBase.endpoints-Union{Tuple{DomainInfo{T}}, Tuple{T}} where T","page":"API Reference","title":"EarthSciMLBase.endpoints","text":"endpoints(d)\n\n\nReturn the endpoints of the partial independent variables for this domain.\n\n\n\n\n\n","category":"method"},{"location":"api/#EarthSciMLBase.get_coupletype-Tuple{ModelingToolkit.AbstractSystem}","page":"API Reference","title":"EarthSciMLBase.get_coupletype","text":"Return the coupling type associated with the given system.\n\n\n\n\n\n","category":"method"},{"location":"api/#EarthSciMLBase.get_dv-Tuple{Any, Any}","page":"API Reference","title":"EarthSciMLBase.get_dv","text":"Return the dependent variable, which is the first argument of the term, unless the term is a time derivative, in which case the dependent variable is the argument of the time derivative.\n\n\n\n\n\n","category":"method"},{"location":"api/#EarthSciMLBase.get_needed_vars-Tuple{ModelingToolkit.ODESystem}","page":"API Reference","title":"EarthSciMLBase.get_needed_vars","text":"get_needed_vars(sys)\n\n\nReturn the indexes of the system variables that the state variables of the final simplified system depend on. This should be done before running structural_simplify on the system.\n\n\n\n\n\n","category":"method"},{"location":"api/#EarthSciMLBase.graph-Tuple{CoupledSystem}","page":"API Reference","title":"EarthSciMLBase.graph","text":"Create a graph from a CoupledSystem using the MetaGraphsNext package.\n\n\n\n\n\n","category":"method"},{"location":"api/#EarthSciMLBase.grid-Union{Tuple{DomainInfo{T}}, Tuple{T}} where T","page":"API Reference","title":"EarthSciMLBase.grid","text":"grid(d)\n\n\nReturn the ranges representing the discretization of the partial independent variables for this domain, based on the discretization intervals given in Δs.\n\n\n\n\n\n","category":"method"},{"location":"api/#EarthSciMLBase.icbc-Tuple{DomainInfo, AbstractVector}","page":"API Reference","title":"EarthSciMLBase.icbc","text":"icbc(di, states)\n\n\nReturn a vector of equations that define the initial and boundary conditions for the given state variables.\n\n\n\n\n\n","category":"method"},{"location":"api/#EarthSciMLBase.init_callback-Tuple{}","page":"API Reference","title":"EarthSciMLBase.init_callback","text":"Types that implement an:\n\ninit_callback(x, sys::CoupledSystem, obs_eqs, domain::DomainInfo)::DECallback\n\nmethod can also be coupled into a CoupledSystem. The init_callback function will be run before the simulator is run to get the callback.\n\n\n\n\n\n","category":"method"},{"location":"api/#EarthSciMLBase.init_u-Tuple{ModelingToolkit.ODESystem, DomainInfo}","page":"API Reference","title":"EarthSciMLBase.init_u","text":"Initialize the state variables.\n\n\n\n\n\n","category":"method"},{"location":"api/#EarthSciMLBase.ivar-Tuple{DomainInfo}","page":"API Reference","title":"EarthSciMLBase.ivar","text":"ivar(di::DomainInfo) -> Any\n\n\nReturn the independent variable associated with these initial and boundary conditions.\n\n\n\n\n\n","category":"method"},{"location":"api/#EarthSciMLBase.observed_expression-Tuple{Any, Any}","page":"API Reference","title":"EarthSciMLBase.observed_expression","text":"observed_expression(eqs, x)\n\n\nReturn an expression for the observed value of a variable x after substituting in the constants observed values of other variables. extra_eqs is a list of additional equations to use in the substitution.\n\n\n\n\n\n","category":"method"},{"location":"api/#EarthSciMLBase.observed_function-Tuple{Any, Any, Any}","page":"API Reference","title":"EarthSciMLBase.observed_function","text":"observed_function(eqs, x, coords)\n\n\nReturn a function to for the observed value of a variable x based on the input arguments in coords. extra_eqs is a list of additional equations to use to determine the value of x.\n\n\n\n\n\n","category":"method"},{"location":"api/#EarthSciMLBase.ode_step!-Tuple{Any, SolverStrangThreads, Vararg{Any, 5}}","page":"API Reference","title":"EarthSciMLBase.ode_step!","text":"Take a step using the ODE solver.\n\n\n\n\n\n","category":"method"},{"location":"api/#EarthSciMLBase.operator_compose","page":"API Reference","title":"EarthSciMLBase.operator_compose","text":"operator_compose(a, b)\noperator_compose(a, b, translate)\n\n\nCompose to systems of equations together by adding the right-hand side terms together of equations that have matching left-hand sides. The left hand sides of two equations will be considered matching if:\n\nThey are both time derivatives of the same variable.\nThe first one is a time derivative of a variable and the second one is the variable itself.\nThere is an entry in the optional translate dictionary that maps the dependent variable in the first system to the dependent variable in the second system, e.g. Dict(sys1.sys.x => sys2.sys.y).\nThere is an entry in the optional translate dictionary that maps the dependent variable in the first system to the dependent variable in the second system, with a conversion factor, e.g. Dict(sys1.sys.x => sys2.sys.y => 6).\n\n\n\n\n\n","category":"function"},{"location":"api/#EarthSciMLBase.param_to_var-Tuple{ModelingToolkit.AbstractSystem, Vararg{Symbol}}","page":"API Reference","title":"EarthSciMLBase.param_to_var","text":"Replace the parameter p in the system sys with a new variable that has the same name, units, and description as p.\n\nparam_to_var(sys, ps)\n\n\nThis can be useful to replace a parameter that does not change in time in a model component with one specified by another system that does change in time (or space). For example, the code below specifies a first-order loss equation, and then changes the temperature (which determines the loss rate) with a temperature value that varies in time. ```\n\n\n\n\n\n","category":"method"},{"location":"api/#EarthSciMLBase.partialderivative_transforms-Tuple{DomainInfo}","page":"API Reference","title":"EarthSciMLBase.partialderivative_transforms","text":"partialderivative_transforms(di::DomainInfo) -> Vector{Any}\n\n\nReturn transform factor to multiply each partial derivative operator by, for example to convert from degrees to meters.\n\n\n\n\n\n","category":"method"},{"location":"api/#EarthSciMLBase.partialderivatives-Tuple{DomainInfo}","page":"API Reference","title":"EarthSciMLBase.partialderivatives","text":"partialderivatives(di::DomainInfo) -> Any\n\n\nReturn the partial derivative operators for the given domain.\n\n\n\n\n\n","category":"method"},{"location":"api/#EarthSciMLBase.partialderivatives_δxyδlonlat-Tuple{AbstractVector}","page":"API Reference","title":"EarthSciMLBase.partialderivatives_δxyδlonlat","text":"partialderivatives_δxyδlonlat(pvars; default_lat)\n\n\nReturn partial derivative operator transform factors corresponding for the given partial-independent variables after converting variables named lon and lat from degrees to x and y meters, assuming they represent longitude and latitude on a spherical Earth.\n\n\n\n\n\n","category":"method"},{"location":"api/#EarthSciMLBase.prune_observed-Tuple{ModelingToolkit.ODESystem}","page":"API Reference","title":"EarthSciMLBase.prune_observed","text":"prune_observed(sys)\n\n\nRemove equations from an ODESystem where the variable in the LHS is not present in any of the equations for the state variables. This can be used to remove computationally intensive equations that are not used in the final model.\n\n\n\n\n\n","category":"method"},{"location":"api/#EarthSciMLBase.pvars-Tuple{DomainInfo}","page":"API Reference","title":"EarthSciMLBase.pvars","text":"pvars(di::DomainInfo) -> Any\n\n\nReturn the partial independent variables associated with these initial and boundary conditions.\n\n\n\n\n\n","category":"method"},{"location":"api/#EarthSciMLBase.single_ode_step!-NTuple{6, Any}","page":"API Reference","title":"EarthSciMLBase.single_ode_step!","text":"Take a step using the ODE solver with the given IIchunk (grid cell interator) and integrator.\n\n\n\n\n\n","category":"method"},{"location":"api/#EarthSciMLBase.steplength-Tuple{Any}","page":"API Reference","title":"EarthSciMLBase.steplength","text":"steplength(timesteps)\n\n\nReturn the time step length common to all of the given timesteps. Throw an error if not all timesteps are the same length.\n\n\n\n\n\n","category":"method"},{"location":"api/#EarthSciMLBase.stiff_callback-Union{Tuple{T}, Tuple{Any, AbstractArray{T}, EarthSciMLBase.SolverStrang, Any, Any}} where T","page":"API Reference","title":"EarthSciMLBase.stiff_callback","text":"A callback to periodically run the stiff solver.\n\n\n\n\n\n","category":"method"},{"location":"api/#EarthSciMLBase.threaded_ode_step!-NTuple{6, Any}","page":"API Reference","title":"EarthSciMLBase.threaded_ode_step!","text":"Take a step using the ODE solver.\n\n\n\n\n\n","category":"method"},{"location":"api/#EarthSciMLBase.timesteps-Union{Tuple{Vararg{AbstractVector{T}}}, Tuple{T}} where T<:AbstractFloat","page":"API Reference","title":"EarthSciMLBase.timesteps","text":"timesteps(tsteps)\n\n\nReturn the time points during which integration should be stopped to run the operators.\n\n\n\n\n\n","category":"method"},{"location":"api/#EarthSciMLBase.tspan-Union{Tuple{DomainInfo{T}}, Tuple{T}} where T<:AbstractFloat","page":"API Reference","title":"EarthSciMLBase.tspan","text":"tspan(d)\n\n\nReturn the time range associated with this domain.\n\n\n\n\n\n","category":"method"},{"location":"icbc/#Initial-and-Boundary-conditions","page":"Initial and Boundary Conditions","title":"Initial and Boundary conditions","text":"","category":"section"},{"location":"icbc/","page":"Initial and Boundary Conditions","title":"Initial and Boundary Conditions","text":"Oftentimes we will want to do a 1, 2, or 3-dimensional simulation, rather than the 0-dimensional simulation we get by default with a system of ordinary differential equations. In these cases, we will need to specify initial and boundary conditions for the system.","category":"page"},{"location":"icbc/","page":"Initial and Boundary Conditions","title":"Initial and Boundary Conditions","text":"To demonstrate how to do this, we will use the following simple system of ordinary differential equations:","category":"page"},{"location":"icbc/","page":"Initial and Boundary Conditions","title":"Initial and Boundary Conditions","text":"using EarthSciMLBase\nusing ModelingToolkit\nusing ModelingToolkit: t_nounits, D_nounits\nt = t_nounits\nD = D_nounits\n\n@parameters x y\n\nfunction ExampleSys()\n @variables u(t) v(t)\n eqs = [\n D(u) ~ √abs(v),\n D(v) ~ √abs(u),\n ]\n ODESystem(eqs, t; name=:Docs₊Example)\nend\n\nExampleSys()","category":"page"},{"location":"icbc/","page":"Initial and Boundary Conditions","title":"Initial and Boundary Conditions","text":"Next, we specify our initial and boundary conditions using the DomainInfo type. We initialize DomainInfo with sets of initial and boundary conditions. In the example below, we set constant initial conditions using constIC and constant boundary conditions using constBC.","category":"page"},{"location":"icbc/","page":"Initial and Boundary Conditions","title":"Initial and Boundary Conditions","text":"using DomainSets\n\nx_min = y_min = t_min = 0.0\nx_max = y_max = t_max = 1.0\n\n# Create constant initial conditions = 16.0 and boundary conditions = 4.0.\nicbc = DomainInfo(\n constIC(4.0, t ∈ Interval(t_min, t_max)),\n constBC(16.0,\n x ∈ Interval(x_min, x_max),\n y ∈ Interval(y_min, y_max),\n ),\n)\nnothing # hide","category":"page"},{"location":"icbc/","page":"Initial and Boundary Conditions","title":"Initial and Boundary Conditions","text":"It is also possible to use periodic boundary conditions with periodicBC and zero-gradient boundary conditions with zerogradBC.","category":"page"},{"location":"icbc/","page":"Initial and Boundary Conditions","title":"Initial and Boundary Conditions","text":"Finally, we combine our initial and boundary conditions with our system of equations using the couple function.","category":"page"},{"location":"icbc/","page":"Initial and Boundary Conditions","title":"Initial and Boundary Conditions","text":"model = couple(ExampleSys(), icbc)\n\neq_sys = convert(PDESystem, model)","category":"page"},{"location":"icbc/","page":"Initial and Boundary Conditions","title":"Initial and Boundary Conditions","text":"We can also look at the expanded boundary conditions of the resulting equation system:","category":"page"},{"location":"icbc/","page":"Initial and Boundary Conditions","title":"Initial and Boundary Conditions","text":"eq_sys.bcs","category":"page"},{"location":"operator_compose/#Operator-Composition","page":"Operator Composition","title":"Operator Composition","text":"","category":"section"},{"location":"operator_compose/","page":"Operator Composition","title":"Operator Composition","text":"There are a lot of cases where there are two different \"processes\" or \"operators\" that change the same variable. For example, CO2 in the atmosphere can be emitted by human activity, and it can also be absorbed by the ocean. In models, typically the emission and removal are considered separate processes which are represented by separate model components. However, when we want to combine these two components into a single model, we need to be able to compose them together.","category":"page"},{"location":"operator_compose/","page":"Operator Composition","title":"Operator Composition","text":"We can use the operator_compose function for this. It composes to systems of equations together by adding the right-hand side terms together of equations that have matching left-hand sides. The left hand sides of two equations will be considered matching if:","category":"page"},{"location":"operator_compose/","page":"Operator Composition","title":"Operator Composition","text":"They are both time derivatives of the same variable.\nThe first one is a time derivative of a variable and the second one is the variable itself.\nThere is an entry in the optional translate dictionary that maps the dependent variable in the first system to the dependent variable in the second system, e.g. Dict(sys1.sys.x => sys2.sys.y).\nThere is an entry in the optional translate dictionary that maps the dependent variable in the first system to the dependent variable in the second system, with a conversion factor, e.g. Dict(sys1.sys.x => sys2.sys.y => 6).","category":"page"},{"location":"operator_compose/","page":"Operator Composition","title":"Operator Composition","text":"Perhaps we can make this somewhat clearer with some examples.","category":"page"},{"location":"operator_compose/#Examples","page":"Operator Composition","title":"Examples","text":"","category":"section"},{"location":"operator_compose/#Example-with-matching-variable-time-derivatives","page":"Operator Composition","title":"Example with matching variable time derivatives","text":"","category":"section"},{"location":"operator_compose/","page":"Operator Composition","title":"Operator Composition","text":"The example below shows that when we operator_compose two systems together that are both equal to D(x) = p, the resulting system is equal to D(x) = 2p.","category":"page"},{"location":"operator_compose/","page":"Operator Composition","title":"Operator Composition","text":"using EarthSciMLBase\nusing ModelingToolkit\nusing ModelingToolkit: t_nounits, D_nounits\nt = t_nounits\nD = D_nounits\n\n\nstruct ExampleSysCoupler sys end\nfunction ExampleSys()\n @variables x(t)\n @parameters p\n ODESystem([D(x) ~ p], t; name=:ExampleSys,\n metadata=Dict(:coupletype=>ExampleSysCoupler))\nend\n\nExampleSys()","category":"page"},{"location":"operator_compose/","page":"Operator Composition","title":"Operator Composition","text":"struct ExampleSys2Coupler sys end\nfunction ExampleSys2()\n @variables x(t)\n @parameters p\n ODESystem([D(x) ~ 2p], t; name=:ExampleSys2,\n metadata=Dict(:coupletype=>ExampleSys2Coupler))\nend\n\nExampleSys2()","category":"page"},{"location":"operator_compose/","page":"Operator Composition","title":"Operator Composition","text":"sys1 = ExampleSys()\nsys2 = ExampleSys2()\n\nfunction EarthSciMLBase.couple2(sys1::ExampleSysCoupler, sys2::ExampleSys2Coupler)\n sys1, sys2 = sys1.sys, sys2.sys\n operator_compose(sys1, sys2)\nend\n\ncombined = couple(sys1, sys2)\n\ncombined_mtk = convert(ODESystem, combined)","category":"page"},{"location":"operator_compose/","page":"Operator Composition","title":"Operator Composition","text":"The simplified equation should be D(x) = p + sys2_xˍt:","category":"page"},{"location":"operator_compose/","page":"Operator Composition","title":"Operator Composition","text":"combined_simplified = structural_simplify(combined_mtk)","category":"page"},{"location":"operator_compose/","page":"Operator Composition","title":"Operator Composition","text":"where sys2_xˍt is also equal to p:","category":"page"},{"location":"operator_compose/","page":"Operator Composition","title":"Operator Composition","text":"observed(combined_simplified)","category":"page"},{"location":"operator_compose/#Example-with-non-matching-variables","page":"Operator Composition","title":"Example with non-matching variables","text":"","category":"section"},{"location":"operator_compose/","page":"Operator Composition","title":"Operator Composition","text":"This example demonstrates a case where one variable in the first system is equal to another variable in the second system:","category":"page"},{"location":"operator_compose/","page":"Operator Composition","title":"Operator Composition","text":"struct ExampleSys3Coupler sys end\nfunction ExampleSys3()\n @variables y(t)\n @parameters p\n ODESystem([D(y) ~ p], t; name=:ExampleSys3,\n metadata=Dict(:coupletype=>ExampleSys3Coupler))\nend\n\nsys1 = ExampleSys()\nsys2 = ExampleSys3()\n\nfunction EarthSciMLBase.couple2(sys1::ExampleSysCoupler, sys2::ExampleSys3Coupler)\n sys1, sys2 = sys1.sys, sys2.sys\n operator_compose(sys1, sys2, Dict(sys1.x => sys2.y))\nend\n\ncombined = couple(sys1, sys2)\ncombined_simplified = structural_simplify(convert(ODESystem, combined))","category":"page"},{"location":"operator_compose/","page":"Operator Composition","title":"Operator Composition","text":"observed(combined_simplified)","category":"page"},{"location":"operator_compose/#Example-with-a-non-ODE-system","page":"Operator Composition","title":"Example with a non-ODE system","text":"","category":"section"},{"location":"operator_compose/","page":"Operator Composition","title":"Operator Composition","text":"In the second case above, we might have a variable in the second system that is equal to a rate, but it is not a time derivative. This could happen if we are extracting emissions from a file, and those emissions are already in units of kg/s, or something similar. The example below demonstrates this case. (Note that this case can also be used with the conversion factors shown in the last example.)","category":"page"},{"location":"operator_compose/","page":"Operator Composition","title":"Operator Composition","text":"struct ExampleSysNonODECoupler sys end\nfunction ExampleSysNonODE()\n @variables y(t)\n @parameters p\n ODESystem([y ~ p], t; name=:ExampleSysNonODE,\n metadata=Dict(:coupletype=>ExampleSysNonODECoupler))\nend\n\nsys1 = ExampleSys()\nsys2 = ExampleSysNonODE()\n\nfunction EarthSciMLBase.couple2(sys1::ExampleSysCoupler, sys2::ExampleSysNonODECoupler)\n sys1, sys2 = sys1.sys, sys2.sys\n operator_compose(sys1, sys2, Dict(sys1.x => sys2.y))\nend\n\ncombined = couple(sys1, sys2)\nsys_combined = structural_simplify(convert(ODESystem, combined))","category":"page"},{"location":"operator_compose/","page":"Operator Composition","title":"Operator Composition","text":"observed(sys_combined)","category":"page"},{"location":"operator_compose/#Example-with-non-matching-variables-and-a-conversion-factor","page":"Operator Composition","title":"Example with non-matching variables and a conversion factor","text":"","category":"section"},{"location":"operator_compose/","page":"Operator Composition","title":"Operator Composition","text":"Finally, this last example shows the fourth case, where a conversion factor is included in the translation dictionary.","category":"page"},{"location":"operator_compose/","page":"Operator Composition","title":"Operator Composition","text":"struct ExampleSysNonODE2Coupler sys end\nfunction ExampleSysNonODE2()\n @variables y(t)\n @parameters p\n ODESystem([y ~ p], t; name=:Docs₊ExampleSysNonODE2,\n metadata=Dict(:coupletype=>ExampleSysNonODE2Coupler))\nend\n\nsys1 = ExampleSys()\nsys2 = ExampleSysNonODE2()\n\nfunction EarthSciMLBase.couple2(sys1::ExampleSysCoupler, sys2::ExampleSysNonODE2Coupler)\n sys1, sys2 = sys1.sys, sys2.sys\n operator_compose(sys1, sys2, Dict(sys1.x => sys2.y => 6.0))\nend\n\ncombined = couple(sys1, sys2)\ncombined_simplified = structural_simplify(convert(ODESystem, combined))","category":"page"},{"location":"operator_compose/","page":"Operator Composition","title":"Operator Composition","text":"observed(combined_simplified)","category":"page"},{"location":"param_to_var/#Converting-parameters-to-variables","page":"Parameter Replacement","title":"Converting parameters to variables","text":"","category":"section"},{"location":"param_to_var/","page":"Parameter Replacement","title":"Parameter Replacement","text":"This can be useful to replace a parameter that does not change in time in a model component with one specified by another system that does change in time (or space). For example, the code below specifies a first-order loss equation, and then changes the temperature (which determines the loss rate) with a temperature value that varies in time.","category":"page"},{"location":"param_to_var/","page":"Parameter Replacement","title":"Parameter Replacement","text":"As an example, we will create a loss equation that depends on the temperature, starting with a constant temperature. We will then create a temperature equation that varies in time, and use the param_to_var function to replace the constant temperature in the loss equation with the time-varying temperature.","category":"page"},{"location":"param_to_var/","page":"Parameter Replacement","title":"Parameter Replacement","text":"So first, let's specify the original system with constant temperature.","category":"page"},{"location":"param_to_var/","page":"Parameter Replacement","title":"Parameter Replacement","text":"using ModelingToolkit, EarthSciMLBase, DynamicQuantities\nusing ModelingToolkit: t, D\n\nstruct LossCoupler sys end\nfunction Loss()\n @variables A(t)=1 [unit=u\"kg\"]\n @parameters k=1 [unit=u\"s^-1\"]\n @parameters T=300 [unit=u\"K\"]\n @constants T₀=300 [unit=u\"K\"]\n eq = D(A) ~ -k*exp(T/T₀) * A\n ODESystem([eq], t; name=:Loss, metadata=Dict(:coupletype=>LossCoupler))\nend\n\nLoss()","category":"page"},{"location":"param_to_var/","page":"Parameter Replacement","title":"Parameter Replacement","text":"Next, we specify the temperature that varies in time.","category":"page"},{"location":"param_to_var/","page":"Parameter Replacement","title":"Parameter Replacement","text":"struct TemperatureCoupler sys end\nfunction Temperature()\n @variables T(t)=300 [unit=u\"K\"]\n @constants Tc=1.0 [unit=u\"K/s\"]\n @constants tc=1.0 [unit=u\"s\"]\n eq = D(T) ~ sin(t/tc)*Tc\n ODESystem([eq], t; name=:Temperature, metadata=Dict(:coupletype=>TemperatureCoupler))\nend\n\nTemperature()","category":"page"},{"location":"param_to_var/","page":"Parameter Replacement","title":"Parameter Replacement","text":"Now, we specify how to compose the two systems using param_to_var.","category":"page"},{"location":"param_to_var/","page":"Parameter Replacement","title":"Parameter Replacement","text":"function EarthSciMLBase.couple2(loss::LossCoupler, temp::TemperatureCoupler)\n loss, temp = loss.sys, temp.sys\n loss = param_to_var(loss, :T)\n ConnectorSystem([loss.T ~ temp.T], loss, temp)\nend","category":"page"},{"location":"param_to_var/","page":"Parameter Replacement","title":"Parameter Replacement","text":"Finally, we create the system components and the composed system.","category":"page"},{"location":"param_to_var/","page":"Parameter Replacement","title":"Parameter Replacement","text":"l = Loss()\ntemp = Temperature()\nvariable_loss = couple(l, temp)\n\nconvert(ODESystem, variable_loss)","category":"page"},{"location":"param_to_var/","page":"Parameter Replacement","title":"Parameter Replacement","text":"If we wanted to, we could then run a simulation with the composed system.","category":"page"},{"location":"operator/#Operators-for-Large-scale-3D-simulations","page":"Operators","title":"Operators for Large-scale 3D simulations","text":"","category":"section"},{"location":"operator/","page":"Operators","title":"Operators","text":"In this documentation so far, we have talked about creating systems of ordinary differential equations in ModelingToolkit and then converting them to systems of partial differential equations to perform 1-, 2-, or 3-dimensional simulations. However, currently this does not work for large scale simulations.","category":"page"},{"location":"operator/","page":"Operators","title":"Operators","text":"While this ModelingToolkit functionality is being built, we have a different solution based on the Operator type in this package. Using this system, we still define systems of ODEs to define behavior in a single grid cell, and we also have Operator processes that define behavior between grid cells.","category":"page"},{"location":"operator/#ODE-System","page":"Operators","title":"ODE System","text":"","category":"section"},{"location":"operator/","page":"Operators","title":"Operators","text":"As an example, let's first define a system of ODEs:","category":"page"},{"location":"operator/","page":"Operators","title":"Operators","text":"using EarthSciMLBase\nusing ModelingToolkit, DifferentialEquations\nusing SciMLOperators, Plots\nusing DomainSets\nusing ModelingToolkit: t_nounits, D_nounits\nt = t_nounits\nD = D_nounits\n\n@parameters y lon = 0.0 lat = 0.0 lev = 1.0 α = 10.0\n@constants p = 1.0\n@variables(\n u(t) = 1.0, v(t) = 1.0, x(t) = 1.0, y(t) = 1.0, windspeed(t) = 1.0\n)\n\neqs = [D(u) ~ -α * √abs(v) + lon,\n D(v) ~ -α * √abs(u) + lat + 1e-14 * lev,\n windspeed ~ lat + lon + lev,\n]\nsys = ODESystem(eqs, t; name=:sys)","category":"page"},{"location":"operator/","page":"Operators","title":"Operators","text":"The equations above don't really have any physical meaning, but they include two state variables, some parameters, and a constant. There is also a variable windspeed which is \"observed\" based on the parameters, rather than being a state variable, which will be important later.","category":"page"},{"location":"operator/#Operator","page":"Operators","title":"Operator","text":"","category":"section"},{"location":"operator/","page":"Operators","title":"Operators","text":"Next, we define an operator. To do so, first we create a new type that is a subtype of Operator:","category":"page"},{"location":"operator/","page":"Operators","title":"Operators","text":"mutable struct ExampleOp <: Operator\n α::Num # Multiplier from ODESystem\nend","category":"page"},{"location":"operator/","page":"Operators","title":"Operators","text":"In the case above, we're setting up our operator so that it can hold a parameter from our ODE system.","category":"page"},{"location":"operator/","page":"Operators","title":"Operators","text":"Next, we need to define a method of EarthSciMLBase.get_scimlop for our operator. This method will be called to get a SciMLOperators.AbstractSciMLOperator that will be used conjunction with the ModelingToolkit system above to integrate the simulation forward in time.","category":"page"},{"location":"operator/","page":"Operators","title":"Operators","text":"function EarthSciMLBase.get_scimlop(op::ExampleOp, csys::CoupledSystem, mtk_sys, domain::DomainInfo, obs_functions, coordinate_transform_functions, u0, p)\n obs_f = obs_functions(op.α)\n grd = EarthSciMLBase.grid(domain)\n function run(du, u, p, t)\n u = reshape(u, size(u0)...)\n du = reshape(du, size(u0)...)\n for ix ∈ 1:size(u, 1)\n for (i, c1) ∈ enumerate(grd[1])\n for (j, c2) ∈ enumerate(grd[2])\n for (k, c3) ∈ enumerate(grd[3])\n # Demonstrate coordinate transforms\n t1 = coordinate_transform_functions[1](t, c1, c2, c3)\n t2 = coordinate_transform_functions[2](t, c1, c2, c3)\n t3 = coordinate_transform_functions[3](t, c1, c2, c3)\n # Demonstrate calculating observed value.\n fv = obs_f(t, c1, c2, c3)\n # Set derivative value.\n du[ix, i, j, k] = (t1 + t2 + t3) * fv\n end\n end\n end\n end\n nothing\n end\n FunctionOperator(run, u0[:], p=p)\nend","category":"page"},{"location":"operator/","page":"Operators","title":"Operators","text":"The function above also doesn't have any physical meaning, but it demonstrates some functionality of the Operator \"s\". First, it retrieves a function to get the current value of an observed variable in our ODE system using the obs_functions argement, and it demonstrates how to call the resulting function to get that value. It also demonstrates how to get coordinate transforms using the coordinate_transform_functions argument. Coordinate transforms are discussed in more detail in the documentation for the DomainInfo type.","category":"page"},{"location":"operator/#Domain","page":"Operators","title":"Domain","text":"","category":"section"},{"location":"operator/","page":"Operators","title":"Operators","text":"Once we have an ODE system and an operator, the final component we need is a domain to run the simulation on. Defining a domain is covered in more depth in the documentation for the DomainInfo type, but for now we'll just define a simple domain:","category":"page"},{"location":"operator/","page":"Operators","title":"Operators","text":"t_min = 0.0\nlon_min, lon_max = -π, π\nlat_min, lat_max = -0.45π, 0.45π\nt_max = 11.5\n\nindepdomain = t ∈ Interval(t_min, t_max)\n\npartialdomains = [lon ∈ Interval(lon_min, lon_max),\n lat ∈ Interval(lat_min, lat_max),\n lev ∈ Interval(1, 3)]\n\ndomain = DomainInfo(\n partialderivatives_δxyδlonlat,\n constIC(16.0, indepdomain), constBC(16.0, partialdomains...);\n grid_spacing = [0.1π, 0.1π, 1])\nnothing #hide","category":"page"},{"location":"operator/","page":"Operators","title":"Operators","text":"Note that our domain includes a coordinate transform to convert from degrees latitude and longitude to meters. Our domain specification also includes grid spacing the the lon, lat, and lev coordinates, which we set as 0.1π, 0.1π, and 1, respectively.","category":"page"},{"location":"operator/","page":"Operators","title":"Operators","text":"warning: Warning\nInitial and boundary conditions are not fully implemented for this case, so regardless of the conditions you specify, the initial conditions will be the default values of the variables in the ODE system, and the boundary conditions will be zero.","category":"page"},{"location":"operator/#Coupling-and-Running-the-Simulation","page":"Operators","title":"Coupling and Running the Simulation","text":"","category":"section"},{"location":"operator/","page":"Operators","title":"Operators","text":"Next, initialize our operator, giving the the windspeed observed variable, and we can couple our ODESystem, Operator, and Domain together into a single model:","category":"page"},{"location":"operator/","page":"Operators","title":"Operators","text":"op = ExampleOp(sys.windspeed)\n\ncsys = couple(sys, op, domain)","category":"page"},{"location":"operator/","page":"Operators","title":"Operators","text":"Finally, we can choose a EarthSciMLBase.SolverStrategy and run the simulation. We choose the SolverStrangThreads strategy, which needs us to specify an ODE solver from the options available in DifferentialEquations.jl for both the MTK system. We choose the Tsit5 solver. Then we create an ODEProblem which can be used to run the simulation. Finally, we solve the problem using the solve function. At this point we need to choose a solver for the Operator part of the system, and we choose the Euler solver. We also choose a splitting time step of 1.0 seconds, which we pass both to our SolverStrangThreads strategy and to the solve function.","category":"page"},{"location":"operator/","page":"Operators","title":"Operators","text":"dt = 1.0 # Splitting time step\nst = SolverStrangThreads(Tsit5(), 1.0)\n\nprob = ODEProblem(csys, st)\nsol = solve(prob, Euler(); dt=1.0)\nnothing #hide","category":"page"},{"location":"operator/","page":"Operators","title":"Operators","text":"After the simulation finishes, we can plot the result:","category":"page"},{"location":"operator/","page":"Operators","title":"Operators","text":"anim = @animate for i ∈ 1:length(sol.u)\n u = reshape(sol.u[i], :, size(domain)...)\n plot(\n heatmap(u[1, :, :, 1]),\n heatmap(u[1, :, :, 1]),\n )\nend\ngif(anim, fps = 15)","category":"page"},{"location":"coord_transforms/#Coordinate-Transforms-for-Partial-Derivatives","page":"Coordinate Transforms","title":"Coordinate Transforms for Partial Derivatives","text":"","category":"section"},{"location":"coord_transforms/","page":"Coordinate Transforms","title":"Coordinate Transforms","text":"Often, the coordinates of a grid may be defined in a different coordinate system than the one in which the partial derivatives are desired. For example, grids are often defined in latitude and longitude, but partial derivatives may be required in units of meters to correspond with wind speeds in meters per second.","category":"page"},{"location":"coord_transforms/","page":"Coordinate Transforms","title":"Coordinate Transforms","text":"To handle this, the DomainInfo type can be used to define coordinate system transformations. To use it, a coordinate transform function first needs to be defined, for example partialderivatives_δxyδlonlat which transforms partial derivatives from longitude and latitude to meters:","category":"page"},{"location":"coord_transforms/","page":"Coordinate Transforms","title":"Coordinate Transforms","text":"using EarthSciMLBase\nusing ModelingToolkit\nusing ModelingToolkit: t, D\nusing DomainSets\nusing DynamicQuantities\n\n@parameters lon [unit = u\"rad\"]\n@parameters lat [unit = u\"rad\"]\n@parameters lev\n\npartialderivatives_δxyδlonlat([lev, lon, lat])","category":"page"},{"location":"coord_transforms/","page":"Coordinate Transforms","title":"Coordinate Transforms","text":"As you can see in the output of the code above, the function should take as arguments a list of the coordinates describing the grid (in the case above we have a 3-dimensional grid with vertical level, latitude, and longitude), and return a Dictionary relating the index of each coordinate with a factor to multiply the partial derivative by to convert it to the desired units. The function only needs to return factors for the coordinates that are being transformed.","category":"page"},{"location":"coord_transforms/","page":"Coordinate Transforms","title":"Coordinate Transforms","text":"To include a coordinate transform in our domain, we include the function in the DomainInfo constructor:","category":"page"},{"location":"coord_transforms/","page":"Coordinate Transforms","title":"Coordinate Transforms","text":"deg2rad(x) = x * π / 180.0\ndomain = DomainInfo(\n partialderivatives_δxyδlonlat,\n constIC(0.0, t ∈ Interval(0.0f0, 3600.0f0)),\n periodicBC(lat ∈ Interval(deg2rad(-90.0f0), deg2rad(90.0f0))),\n periodicBC(lon ∈ Interval(deg2rad(-180.0f0), deg2rad(180.0f0))),\n zerogradBC(lev ∈ Interval(1.0f0, 10.0f0)),\n);","category":"page"},{"location":"coord_transforms/","page":"Coordinate Transforms","title":"Coordinate Transforms","text":"Multiple functions can be included in the DomainInfo constructor, just by including them as a vector, e.g.:","category":"page"},{"location":"coord_transforms/","page":"Coordinate Transforms","title":"Coordinate Transforms","text":"domain = DomainInfo(\n [transform1, transform2, ...],\n constIC(0.0, t ∈ Interval(0.0f0, 3600.0f0)),\n ...\n)","category":"page"},{"location":"coord_transforms/","page":"Coordinate Transforms","title":"Coordinate Transforms","text":"After we include the coordinate transforms in our domain, in general everything should be handled automatically. The coordinate transforms may also be automatically added when different model components are coupled together, so you may not need to worry about them at all in many cases. However, if you would like to use the transformed partial derivatives, for example to create a new PDE equation system, you can get them using the EarthSciMLBase.partialderivatives function:","category":"page"},{"location":"coord_transforms/","page":"Coordinate Transforms","title":"Coordinate Transforms","text":"δs = partialderivatives(domain)","category":"page"},{"location":"coord_transforms/","page":"Coordinate Transforms","title":"Coordinate Transforms","text":"This returns a list of functions, one corresponding to each coordinate in our domain. Then we can calculate the symbolic partial derivative of a variable by just calling each function:","category":"page"},{"location":"coord_transforms/","page":"Coordinate Transforms","title":"Coordinate Transforms","text":"@variables u(t)\n\n[δs[i](u) for i ∈ eachindex(δs)]","category":"page"},{"location":"coord_transforms/","page":"Coordinate Transforms","title":"Coordinate Transforms","text":"You can see an example of how this is implemented in the source code for the Advection model component.","category":"page"},{"location":"coord_transforms/","page":"Coordinate Transforms","title":"Coordinate Transforms","text":"Additional transformation functions may be defined in other packages. We recommend that the names of these functions start with partialderivatives_ to make it clear that they are intended to be used in this context.","category":"page"},{"location":"example_all_together/","page":"All Together","title":"All Together","text":"CurrentModule = EarthSciMLBase","category":"page"},{"location":"example_all_together/#Example-using-different-components-of-EarthSciMLBase-together","page":"All Together","title":"Example using different components of EarthSciMLBase together","text":"","category":"section"},{"location":"example_all_together/","page":"All Together","title":"All Together","text":"This example shows how to define and couple different components of EarthSciMLBase together to create a more complex model. First, we create several model components.","category":"page"},{"location":"example_all_together/","page":"All Together","title":"All Together","text":"Our first example system is a simple reaction system:","category":"page"},{"location":"example_all_together/","page":"All Together","title":"All Together","text":"using EarthSciMLBase\nusing ModelingToolkit, Catalyst, DomainSets, MethodOfLines, DifferentialEquations\nusing ModelingToolkit: t_nounits, D_nounits\nt = t_nounits\nD = D_nounits\nusing Plots\n\n# Create our independent variable `t` and our partially-independent variables `x` and `y`.\n@parameters x y\n\nstruct ExampleSys1Coupler sys end\nfunction ExampleSys1()\n @species c₁(t)=5.0 c₂(t)=5.0\n rs = ReactionSystem(\n [Reaction(2.0, [c₁], [c₂])],\n t; name=:Sys1, combinatoric_ratelaws=false)\n convert(ODESystem, complete(rs), metadata=Dict(:coupletype=>ExampleSys1Coupler))\nend\n\nExampleSys1()","category":"page"},{"location":"example_all_together/","page":"All Together","title":"All Together","text":"Our second example system is a simple ODE system, with the same two variables.","category":"page"},{"location":"example_all_together/","page":"All Together","title":"All Together","text":"struct ExampleSys2Coupler sys end\nfunction ExampleSys2()\n @variables c₁(t)=5.0 c₂(t)=5.0\n @parameters p₁=1.0 p₂=0.5\n ODESystem(\n [D(c₁) ~ -p₁, D(c₂) ~ p₂],\n t; name=:Sys2, metadata=Dict(:coupletype=>ExampleSys2Coupler))\nend\n\nExampleSys2()","category":"page"},{"location":"example_all_together/","page":"All Together","title":"All Together","text":"Now, we specify what should happen when we couple the two systems together. In this case, we want the the derivative of the composed system to be equal to the sum of the derivatives of the two systems. We can do that using the operator_compose function from this package.","category":"page"},{"location":"example_all_together/","page":"All Together","title":"All Together","text":"function EarthSciMLBase.couple2(sys1::ExampleSys1Coupler, sys2::ExampleSys2Coupler)\n sys1, sys2 = sys1.sys, sys2.sys\n sys1 = convert(ODESystem, sys1)\n operator_compose(sys1, sys2)\nend\nnothing # hide","category":"page"},{"location":"example_all_together/","page":"All Together","title":"All Together","text":"Once we specify all of the above, it is simple to create our two individual systems and then couple them together. ","category":"page"},{"location":"example_all_together/","page":"All Together","title":"All Together","text":"sys1 = ExampleSys1()\nsys2 = ExampleSys2()\nsys = couple(sys1, sys2)\n\nconvert(ODESystem, sys)","category":"page"},{"location":"example_all_together/","page":"All Together","title":"All Together","text":"At this point we have an ODE system that is composed of two other ODE systems. We can inspect its observed variables using the observed function.","category":"page"},{"location":"example_all_together/","page":"All Together","title":"All Together","text":"simplified_sys = structural_simplify(convert(ODESystem, sys))","category":"page"},{"location":"example_all_together/","page":"All Together","title":"All Together","text":"observed(simplified_sys)","category":"page"},{"location":"example_all_together/","page":"All Together","title":"All Together","text":"We can also run simulations using this system:","category":"page"},{"location":"example_all_together/","page":"All Together","title":"All Together","text":"odeprob = ODEProblem(simplified_sys, [], (0.0,10.0), [])\nodesol = solve(odeprob)\nplot(odesol)","category":"page"},{"location":"example_all_together/","page":"All Together","title":"All Together","text":"Once we've confirmed that our model works in a 0D \"box model\" setting, we can expand it to 1, 2, or 3 dimensions using by adding in initial and boundary conditions. We will also add in advection using constant-velocity wind fields add the same time.","category":"page"},{"location":"example_all_together/","page":"All Together","title":"All Together","text":"x_min = y_min = t_min = 0.0\nx_max = y_max = t_max = 1.0\ndomain = DomainInfo(\n constIC(4.0, t ∈ Interval(t_min, t_max)),\n periodicBC(x ∈ Interval(x_min, x_max)),\n zerogradBC(y ∈ Interval(y_min, y_max)),\n)\n\nsys_pde = couple(sys, domain, ConstantWind(t, 1.0, 1.0), Advection())\n\nsys_pde_mtk = convert(PDESystem, sys_pde)","category":"page"},{"location":"example_all_together/","page":"All Together","title":"All Together","text":"Now we can inspect this new system that we've created:","category":"page"},{"location":"example_all_together/","page":"All Together","title":"All Together","text":"sys_pde_mtk.dvs","category":"page"},{"location":"example_all_together/","page":"All Together","title":"All Together","text":"sys_pde_mtk.bcs","category":"page"},{"location":"example_all_together/","page":"All Together","title":"All Together","text":"Finally, we can run a simulation using this system:","category":"page"},{"location":"example_all_together/","page":"All Together","title":"All Together","text":"discretization = MOLFiniteDifference([x=>10, y=>10], t, approx_order=2)\n@time pdeprob = discretize(sys_pde_mtk, discretization)\n@time pdesol = solve(pdeprob, Tsit5(), saveat=0.1)\n\n# Plot the solution.\ndiscrete_x, discrete_y, discrete_t = pdesol[x], pdesol[y], pdesol[t]\n@variables Sys1₊c₁(..) Sys1₊c₂(..)\nsolc1, solc2 = pdesol[Sys1₊c₁(t, x, y)], pdesol[Sys1₊c₂(t, x, y)]\nanim = @animate for k in 1:length(discrete_t)\n p1 = heatmap(solc1[k, 1:end-1, 1:end-1], title=\"c₁ t=\\$(discrete_t[k])\", clim=(0,4.0), lab=:none)\n p2 = heatmap(solc2[k, 1:end-1, 1:end-1], title=\"c₂ t=\\$(discrete_t[k])\", clim=(0,7.0), lab=:none)\n plot(p1, p2, layout=(1,2), size=(800,400))\nend\ngif(anim, fps = 8)","category":"page"},{"location":"example_all_together/","page":"All Together","title":"All Together","text":"Because our system is a system of ordinary differential equations rather than partial differential equations, all of the grid cells in the animation above have the same value. Refer to the advection example for an example of a system of partial differential equations.","category":"page"},{"location":"composition/#Model-Composition","page":"Composition","title":"Model Composition","text":"","category":"section"},{"location":"composition/","page":"Composition","title":"Composition","text":"A main goal of the EarthSciMLBase package is to allow model components to be created independently and composed together. We achieve this by creating coupletypes that allow us to use multiple dispatch on the EarthSciMLBase.couple2 function to specify how particular systems should be coupled together.","category":"page"},{"location":"composition/","page":"Composition","title":"Composition","text":"To demonstrate how this works, below we define three model components, Photolysis, Chemistry, and Emissions, which represent different processes in the atmosphere.","category":"page"},{"location":"composition/","page":"Composition","title":"Composition","text":"using ModelingToolkit, Catalyst, EarthSciMLBase\nusing ModelingToolkit: t_nounits\nt = t_nounits\n\nstruct PhotolysisCoupler sys end\nfunction Photolysis(; name=:Photolysis)\n @variables j_NO2(t)\n eqs = [\n j_NO2 ~ max(sin(t/86400),0)\n ]\n ODESystem(eqs, t, [j_NO2], [], name=name,\n metadata=Dict(:coupletype=>PhotolysisCoupler))\nend\n\nPhotolysis()","category":"page"},{"location":"composition/","page":"Composition","title":"Composition","text":"You can see that the system defined above is mostly a standard ModelingToolkit ODE system, except for two things:","category":"page"},{"location":"composition/","page":"Composition","title":"Composition","text":"The first unique part is that we define a PhotolysisCoupler type:","category":"page"},{"location":"composition/","page":"Composition","title":"Composition","text":"struct PhotolysisCoupler sys end","category":"page"},{"location":"composition/","page":"Composition","title":"Composition","text":"It is important that this type is a struct, and that the struct has a single field named sys. When defining your own coupled systems, you can just copy the line of code above but change the name of the type (i.e., change it to something besides PhotolysisCoupler).","category":"page"},{"location":"composition/","page":"Composition","title":"Composition","text":"The second unique part is that we define some metadata for our ODESystem to tell it what coupling type to use:","category":"page"},{"location":"composition/","page":"Composition","title":"Composition","text":"metadata=Dict(:coupletype=>PhotolysisCoupler)","category":"page"},{"location":"composition/","page":"Composition","title":"Composition","text":"Again, when making your own components, just copy the code above but change PhotolysisCoupler to something else.","category":"page"},{"location":"composition/","page":"Composition","title":"Composition","text":"Let's follow the same process for some additional components:","category":"page"},{"location":"composition/","page":"Composition","title":"Composition","text":"struct ChemistryCoupler sys end\nfunction Chemistry(; name=:Chemistry)\n @parameters jNO2\n @species NO2(t)\n rxs = [\n Reaction(jNO2, [NO2], [], [1], [])\n ]\n rsys = ReactionSystem(rxs, t, [NO2], [jNO2]; \n combinatoric_ratelaws=false, name=name)\n convert(ODESystem, complete(rsys), metadata=Dict(:coupletype=>ChemistryCoupler))\nend\n\nChemistry()","category":"page"},{"location":"composition/","page":"Composition","title":"Composition","text":"For our chemistry component above, because it's is originally a ReactionSystem instead of an ODESystem, we convert it to an ODESystem before adding the metadata.","category":"page"},{"location":"composition/","page":"Composition","title":"Composition","text":"struct EmissionsCoupler sys end\nfunction Emissions(; name=:Emissions)\n @parameters emis = 1\n @variables NO2(t)\n eqs = [NO2 ~ emis]\n ODESystem(eqs, t; name=name,\n metadata=Dict(:coupletype=>EmissionsCoupler))\nend\n\nEmissions()","category":"page"},{"location":"composition/","page":"Composition","title":"Composition","text":"Now, we need to define ways to couple the model components together. We can do this by defining a coupling function (as a method of EarthSciMLBase.couple2) for each pair of model components that we want to couple. Each coupling function should have the signature EarthSciMLBase.couple2(a::ACoupler, b::BCoupler)::ConnectorSystem, and should assume that the two ODESystems are inside their respective couplers in the sys field. It should make any edits to the components as needed and return a ConnectorSystem which defines the relationship between the two components.","category":"page"},{"location":"composition/","page":"Composition","title":"Composition","text":"The code below defines a method for coupling the Chemistry and Photolysis components. First, it uses the param_to_var function to convert the photolysis rate parameter jNO2 from the Chemistry component to a variable, then it creates a new Chemistry component with the updated photolysis rate, and finally, it creates a ConnectorSystem object that sets the j_NO2 variable from the Photolysis component equal to the jNO2 variable from the Chemistry component. The next effect is that the photolysis rate in the Chemistry component is now controlled by the Photolysis component.","category":"page"},{"location":"composition/","page":"Composition","title":"Composition","text":"function EarthSciMLBase.couple2(c::ChemistryCoupler, p::PhotolysisCoupler)\n c, p = c.sys, p.sys\n c = param_to_var(convert(ODESystem, c), :jNO2)\n ConnectorSystem([c.jNO2 ~ p.j_NO2], c, p)\nend\nnothing # hide","category":"page"},{"location":"composition/","page":"Composition","title":"Composition","text":"Next, we define a method for coupling the Chemistry and Emissions components. To do this we use the operator_compose function to add the NO2 species from the Emissions component to the time derivative of NO2 in the Chemistry component.","category":"page"},{"location":"composition/","page":"Composition","title":"Composition","text":"function EarthSciMLBase.couple2(c::ChemistryCoupler, emis::EmissionsCoupler)\n c, emis = c.sys, emis.sys\n operator_compose(convert(ODESystem, c), emis, Dict(\n c.NO2 => emis.NO2,\n ))\nend\nnothing # hide","category":"page"},{"location":"composition/","page":"Composition","title":"Composition","text":"Finally, we can compose the model components together to create our complete model. To do so, we just initialize each model component and add the components together using the couple function. We can use the convert function to convert the composed model to a ModelingToolkit ODESystem so we can see what the combined equations look like.","category":"page"},{"location":"composition/","page":"Composition","title":"Composition","text":"model = couple(Photolysis(), Chemistry(), Emissions())\nconvert(ODESystem, model)","category":"page"},{"location":"composition/","page":"Composition","title":"Composition","text":"Finally, we can use the graph function to visualize the model components and their connections.","category":"page"},{"location":"composition/","page":"Composition","title":"Composition","text":"using MetaGraphsNext\nusing CairoMakie, GraphMakie\n\ng = graph(model)\n\nf, ax, p = graphplot(g; ilabels=labels(g))\nhidedecorations!(ax); hidespines!(ax); ax.aspect = DataAspect()\n\nf","category":"page"},{"location":"","page":"Home","title":"Home","text":"CurrentModule = EarthSciMLBase","category":"page"},{"location":"#EarthSciMLBase:-Utilities-for-Symbolic-Earth-Science-Modeling-and-Machine-Learning","page":"Home","title":"EarthSciMLBase: Utilities for Symbolic Earth Science Modeling and Machine Learning","text":"","category":"section"},{"location":"","page":"Home","title":"Home","text":"Documentation for EarthSciMLBase.","category":"page"},{"location":"","page":"Home","title":"Home","text":"This package contains utilities for constructing Earth Science models in Julia using ModelingToolkit.jl.","category":"page"},{"location":"#Installation","page":"Home","title":"Installation","text":"","category":"section"},{"location":"","page":"Home","title":"Home","text":"using Pkg\nPkg.add(\"EarthSciMLBase\")","category":"page"},{"location":"#Feature-Summary","page":"Home","title":"Feature Summary","text":"","category":"section"},{"location":"","page":"Home","title":"Home","text":"This package contains types and functions designed to simplify the process of constructing and composing symbolically-defined Earth Science model components together.","category":"page"},{"location":"#Feature-List","page":"Home","title":"Feature List","text":"","category":"section"},{"location":"","page":"Home","title":"Home","text":"Operations to compose ModelingToolkit.jl equation systems together.\nOperations to add intitial and boundary conditions to systems and to turn ODE systems into PDE systems, and to provide coordinate transformations.\nOperations to add Advection terms to systems.\nA Simulator type for running large-scale simulations.","category":"page"},{"location":"#Contributing","page":"Home","title":"Contributing","text":"","category":"section"},{"location":"","page":"Home","title":"Home","text":"Please refer to the SciML ColPrac: Contributor's Guide on Collaborative Practices for Community Packages for guidance on PRs, issues, and other matters relating to contributing.","category":"page"},{"location":"#Reproducibility","page":"Home","title":"Reproducibility","text":"","category":"section"},{"location":"","page":"Home","title":"Home","text":"
The documentation of this EarthSciML package was built using these direct dependencies,","category":"page"},{"location":"","page":"Home","title":"Home","text":"using Pkg # hide\nPkg.status() # hide","category":"page"},{"location":"","page":"Home","title":"Home","text":"
","category":"page"},{"location":"","page":"Home","title":"Home","text":"
and using this machine and Julia version.","category":"page"},{"location":"","page":"Home","title":"Home","text":"using InteractiveUtils # hide\nversioninfo() # hide","category":"page"},{"location":"","page":"Home","title":"Home","text":"
","category":"page"},{"location":"","page":"Home","title":"Home","text":"
A more complete overview of all dependencies and their versions is also provided.","category":"page"},{"location":"","page":"Home","title":"Home","text":"using Pkg # hide\nPkg.status(;mode = PKGMODE_MANIFEST) # hide","category":"page"},{"location":"","page":"Home","title":"Home","text":"
","category":"page"},{"location":"","page":"Home","title":"Home","text":"You can also download the \nmanifest file and the\nproject file.","category":"page"},{"location":"advection/#Advection","page":"Advection","title":"Advection","text":"","category":"section"},{"location":"advection/","page":"Advection","title":"Advection","text":"The Advection function adds advection to a system of equations. This is useful for modeling the transport of a substance by a fluid. Advection is implemented with the Advection type.","category":"page"},{"location":"advection/","page":"Advection","title":"Advection","text":"!warning Fully symbolic partial differential equations like those shown here don't currently work on domains that have a large number of grid cells. See here for additional information.","category":"page"},{"location":"advection/","page":"Advection","title":"Advection","text":"To demonstrate how this can work, we will start with a simple system of equations:","category":"page"},{"location":"advection/","page":"Advection","title":"Advection","text":"using EarthSciMLBase, ModelingToolkit\nusing ModelingToolkit: t_nounits, D_nounits\nt = t_nounits\nD = D_nounits\n\nfunction ExampleSys()\n @variables y(t)\n @parameters p=2.0\n ODESystem([D(y) ~ p], t; name=:ExampleSys)\nend\n\nExampleSys()","category":"page"},{"location":"advection/","page":"Advection","title":"Advection","text":"We also need to create our initial and boundary conditions.","category":"page"},{"location":"advection/","page":"Advection","title":"Advection","text":"using DomainSets\n@parameters x\ndomain = DomainInfo(constIC(0.0, t ∈ Interval(0, 1.0)), constBC(1.0, x ∈ Interval(0, 1.0)))\nnothing # hide","category":"page"},{"location":"advection/","page":"Advection","title":"Advection","text":"Now we convert add advection to each of the state variables. We're also adding a constant wind (ConstantWind) in the x-direction, with a speed of 1.0.","category":"page"},{"location":"advection/","page":"Advection","title":"Advection","text":"sys_advection = couple(ExampleSys(), domain, ConstantWind(t, 1.0), Advection())\nsys_mtk = convert(PDESystem, sys_advection)","category":"page"},{"location":"advection/","page":"Advection","title":"Advection","text":"Finally, we can discretize the system and solve it:","category":"page"},{"location":"advection/","page":"Advection","title":"Advection","text":"using MethodOfLines, DifferentialEquations, Plots\ndiscretization = MOLFiniteDifference([x=>10], t, approx_order=2)\n@time prob = discretize(sys_mtk, discretization)\n@time sol = solve(prob, Tsit5(), saveat=0.1)\n\n\n# Plot the solution.\ndiscrete_x = sol[x]\ndiscrete_t = sol[t]\nsoly = sol[sys_mtk.dvs[3]]\nanim = @animate for k in 1:length(discrete_t)\n plot(discrete_x, soly[k, 1:end], title=\"t=\\$(discrete_t[k])\", ylim=(0,2.5), lab=:none)\nend\ngif(anim, fps = 8)","category":"page"}] +[{"location":"api/#API-Index","page":"API Reference","title":"API Index","text":"","category":"section"},{"location":"api/","page":"API Reference","title":"API Reference","text":"","category":"page"},{"location":"api/#API-Documentation","page":"API Reference","title":"API Documentation","text":"","category":"section"},{"location":"api/","page":"API Reference","title":"API Reference","text":"Modules = [EarthSciMLBase]","category":"page"},{"location":"api/#EarthSciMLBase.Advection","page":"API Reference","title":"EarthSciMLBase.Advection","text":"Apply advection to a model.\n\n\n\n\n\n","category":"type"},{"location":"api/#EarthSciMLBase.ConnectorSystem","page":"API Reference","title":"EarthSciMLBase.ConnectorSystem","text":"A connector for two systems.\n\neqs\nfrom\nto\n\n\n\n\n\n","category":"type"},{"location":"api/#EarthSciMLBase.CoupledSystem","page":"API Reference","title":"EarthSciMLBase.CoupledSystem","text":"A system for composing together other systems using the couple function.\n\nsystems: Model components to be composed together\ndomaininfo: Initial and boundary conditions and other domain information\npdefunctions: A vector of functions where each function takes as an argument the resulting PDESystem after DomainInfo is added to this system, and returns a transformed PDESystem.\n\nops: A vector of operators to run during simulations.\n\ncallbacks: A vector of callbacks to run during simulations.\ninit_callbacks: Objects x with an init_callback(x, Simulator)::DECallback method.\n\nThings that can be added to a CoupledSystem: * ModelingToolkit.ODESystems. If the ODESystem has a field in the metadata called :coupletype (e.g. ModelingToolkit.get_metadata(sys)[:coupletype] returns a struct type with a single field called sys) then that type will be used to check for methods of EarthSciMLBase.couple that use that type. * Operators * DomainInfos * Callbacks * Types X that implement a EarthSciMLBase.init_callback(::X, sys::CoupledSystem, sys_mtk, obs_eqs, domain::DomainInfo)::DECallback method * Other CoupledSystems * Types X that implement a EarthSciMLBase.couple2(::X, ::CoupledSystem) or EarthSciMLBase.couple2(::CoupledSystem, ::X) method. * Tuples or AbstractVectors of any of the things above.\n\n\n\n\n\n","category":"type"},{"location":"api/#EarthSciMLBase.DomainInfo","page":"API Reference","title":"EarthSciMLBase.DomainInfo","text":"Domain information for a ModelingToolkit.jl PDESystem. It can be used with the + operator to add initial and boundary conditions and coordinate transforms to a ModelingToolkit.jl ODESystem or Catalyst.jl ReactionSystem.\n\nNOTE: The independent variable (usually time) must be first in the list of initial and boundary conditions.\n\npartial_derivative_funcs: Function that returns spatial derivatives of the partially-independent variables, optionally performing a coordinate transformation first.\nCurrent function options in this package are:\npartialderivatives_δxyδlonlat: Returns partial derivatives after transforming any variables named lat and lon\nfrom degrees to cartesian meters, assuming a spherical Earth.\nOther packages may implement additional functions. They are encouraged to use function names starting with partialderivatives_.\n\ngrid_spacing: The discretization intervals for the partial independent variables.\nicbc: The sets of initial and/or boundary conditions.\nspatial_ref: The spatial reference system for the domain.\ntime_offset: The time offset for the domain.\n\n\n\n\n\n","category":"type"},{"location":"api/#EarthSciMLBase.ICBCcomponent","page":"API Reference","title":"EarthSciMLBase.ICBCcomponent","text":"Initial and boundary condition components that can be combined to create an DomainInfo object.\n\n\n\n\n\n","category":"type"},{"location":"api/#EarthSciMLBase.Operator","page":"API Reference","title":"EarthSciMLBase.Operator","text":"Operators are objects that modify the current state of a Simulator system. Each operator should be define a function with the signature:\n\n`EarthSciMLBase.get_scimlop(::Operator, csys::CoupledSystem, mtk_sys, domain::DomainInfo, obs_functions, coordinate_transform_functions, u0, p)::AbstractSciMLOperator`\n\nwhich should return a SciMLOperators.AbstractSciMLOperator. Refer to the SciMLOperators.jl documentation for more information on how to define operators.\n\n\n\n\n\n","category":"type"},{"location":"api/#EarthSciMLBase.SolverIMEX","page":"API Reference","title":"EarthSciMLBase.SolverIMEX","text":"A solver strategy based on implicit-explicit (IMEX) time integration. See here for additional information.\n\nkwargs:\n\nstiff_scimlop: Whether the stiff ODE function should be implemented as a SciMLOperator.\nstiff_sparse: Whether the stiff ODE function should use a sparse Jacobian.\nstiff_jac: Whether the stiff ODE function should use an analytical Jacobian.\nstiffjacscimlop: Whether the stiff ODE function Jacobian should be implemented as a SciMLOperator. (Ignored if stiff_jac==false.)\nstiff_tgrad: Whether the stiff ODE function should use an analytical time gradient.\n\nAdditional kwargs for ODEProblem constructor:\n\nu0: initial condtions; if \"nothing\", default values will be used.\np: parameters; if \"nothing\", default values will be used.\nname: name of the model.\n\n\n\n\n\n","category":"type"},{"location":"api/#EarthSciMLBase.SolverStrang","page":"API Reference","title":"EarthSciMLBase.SolverStrang","text":"A simulator strategy based on Strang splitting. Choose either SimulatorStrangThreads or SimulatorStrangSerial to run the simulation.\n\nkwargs for ODEProblem constructor:\n\nu0: initial condtions; if \"nothing\", default values will be used.\np: parameters; if \"nothing\", default values will be used.\nnonstiff_params: parameters for the non-stiff ODE system.\nname: name of the system.\n\n\n\n\n\n","category":"type"},{"location":"api/#EarthSciMLBase.SolverStrangSerial","page":"API Reference","title":"EarthSciMLBase.SolverStrangSerial","text":"# Specify the stiff ODE solver algorithm.\n# `timestep` is the length of time for each splitting step.\nSimulatorStrangSerial(stiffalg, timestep; kwargs...)\n\nPerform a simulation using Strang splitting, where the MTK system is assumed to be stiff and the operators are assumed to be non-stiff. The solution will be calculated in serial.\n\nAdditional kwargs for ODEProblem constructor:\n\nu0: initial condtions; if \"nothing\", default values will be used.\np: parameters; if \"nothing\", default values will be used.\nnonstiff_params: parameters for the Operators.\n\nstiffalg: Stiff solver algorithm to use (see https://docs.sciml.ai/DiffEqDocs/stable/solvers/ode_solve/)\ntimestep: Length of each splitting time step\nstiff_kwargs: Keyword arguments for the stiff ODEProblem constructor and solver.\n\n\n\n\n\n","category":"type"},{"location":"api/#EarthSciMLBase.SolverStrangThreads","page":"API Reference","title":"EarthSciMLBase.SolverStrangThreads","text":"# Specify the number of threads and the stiff ODE solver algorithm.\n# `timestep` is the length of time for each splitting step.\nSimulatorStrangThreads(threads, stiffalg, timestep; kwargs...)\n# Use the default number of threads.\nSimulatorStrangThreads(stiffalg, timestep; kwargs...)\n\nPerform a simulation using Strang splitting, where the MTK system is assumed to be stiff and the operators are assumed to be non-stiff. The solution of the stiff ODE system is parallelized across grid cells using the specified number of threads.\n\nthreads: Number of threads to use\nstiffalg: Stiff solver algorithm to use (see https://docs.sciml.ai/DiffEqDocs/stable/solvers/ode_solve/)\ntimestep: Length of each splitting time step\nstiff_kwargs: Keyword arguments for the stiff ODEProblem constructor and solver.\n\n\n\n\n\n","category":"type"},{"location":"api/#EarthSciMLBase.SolverStrategy","page":"API Reference","title":"EarthSciMLBase.SolverStrategy","text":"SolverStrategy is an abstract type that defines the strategy for running a simulation.\n\n\n\n\n\n","category":"type"},{"location":"api/#EarthSciMLBase.constBC","page":"API Reference","title":"EarthSciMLBase.constBC","text":"Construct constant boundary conditions equal to the value specified by val.\n\nval: The value of the constant boundary conditions.\npartialdomains: The partial domains, e.g. [x ∈ Interval(x_min, x_max), y ∈ Interval(y_min, y_max)].\n\n\n\n\n\n","category":"type"},{"location":"api/#EarthSciMLBase.constIC","page":"API Reference","title":"EarthSciMLBase.constIC","text":"Construct constant initial conditions equal to the value specified by val.\n\nval: The value of the constant initial conditions.\nindepdomain: The independent domain, e.g. t ∈ Interval(t_min, t_max).\n\n\n\n\n\n","category":"type"},{"location":"api/#EarthSciMLBase.periodicBC","page":"API Reference","title":"EarthSciMLBase.periodicBC","text":"Construct periodic boundary conditions for the given partialdomains. Periodic boundary conditions are defined as when the value at one side of the domain is set equal to the value at the other side, so that the domain \"wraps around\" from one side to the other.\n\npartialdomains: The partial domains, e.g. [x ∈ Interval(x_min, x_max), y ∈ Interval(y_min, y_max)].\n\n\n\n\n\n","category":"type"},{"location":"api/#EarthSciMLBase.zerogradBC","page":"API Reference","title":"EarthSciMLBase.zerogradBC","text":"Construct zero-gradient boundary conditions for the given partialdomains.\n\npartialdomains: The partial domains, e.g. [x ∈ Interval(x_min, x_max), y ∈ Interval(y_min, y_max)].\n\n\n\n\n\n","category":"type"},{"location":"api/#Base.convert-Tuple{Type{<:ModelingToolkit.ODESystem}, CoupledSystem}","page":"API Reference","title":"Base.convert","text":"convert(, sys; name, simplify, kwargs...)\n\n\nGet the ODE ModelingToolkit ODESystem representation of a CoupledSystem.\n\nkwargs:\n\nname: The desired name for the resulting ODESystem\nsimplify: if true, the observed variables that are not needed to specify the state variables will be pruned and returned as a second return value after the ODESystem, which will be structurally simplified.\n\n\n\n\n\n","category":"method"},{"location":"api/#Base.convert-Tuple{Type{<:ModelingToolkit.PDESystem}, CoupledSystem}","page":"API Reference","title":"Base.convert","text":"convert(, sys; name, kwargs...)\n\n\nGet the ModelingToolkit PDESystem representation of a CoupledSystem.\n\n\n\n\n\n","category":"method"},{"location":"api/#EarthSciMLBase.ConstantWind-Tuple{Any, Vararg{Any}}","page":"API Reference","title":"EarthSciMLBase.ConstantWind","text":"ConstantWind(t, vals; name)\n\n\nConstruct a constant wind velocity model component with the given wind speed(s), which should include units. For example, ConstantWind(t, 1u\"m/s\", 2u\"m/s\").\n\n\n\n\n\n","category":"method"},{"location":"api/#EarthSciMLBase.MeanWind-Tuple{Any, DomainInfo}","page":"API Reference","title":"EarthSciMLBase.MeanWind","text":"MeanWind(t, domain)\n\n\nA model component that represents the mean wind velocity, where pvars is the partial dependent variables for the domain.\n\n\n\n\n\n","category":"method"},{"location":"api/#EarthSciMLBase.add_dims-Tuple{Any, AbstractVector, AbstractVector}","page":"API Reference","title":"EarthSciMLBase.add_dims","text":"add_dims(expression, vars, dims)\nadd_dims(equation, vars, dims)\n\nAdd the given dimensions to each variable in vars in the given expression or equation. Each variable in vars must be unidimensional, i.e. defined like @variables u(t) rather than @variables u(..).\n\nExample:\n\nusing EarthSciMLBase, ModelingToolkit\n\n@parameters x y k t\n@variables u(t) q(t)\nexp = 2u + 3k*q + 1\nEarthSciMLBase.add_dims(exp, [u, q], [x, y, t])\n\n# output\n1 + 2u(x, y, t) + 3k*q(x, y, t)\n\n\n\n\n\n","category":"method"},{"location":"api/#EarthSciMLBase.add_metadata-Tuple{Any, Any}","page":"API Reference","title":"EarthSciMLBase.add_metadata","text":"Add the units and description in the variable from to the variable to.\n\n\n\n\n\n","category":"method"},{"location":"api/#EarthSciMLBase.add_scope-Tuple{Any, Any, Any}","page":"API Reference","title":"EarthSciMLBase.add_scope","text":"add_scope(sys, v, iv)\n\n\nAdd a system scope to a variable name, for example so that x in system sys1 becomes sys1₊x. iv is the independent variable.\n\n\n\n\n\n","category":"method"},{"location":"api/#EarthSciMLBase.copy_with_change-Tuple{ModelingToolkit.ODESystem}","page":"API Reference","title":"EarthSciMLBase.copy_with_change","text":"copy_with_change(\n sys;\n eqs,\n name,\n unknowns,\n parameters,\n metadata,\n continuous_events,\n discrete_events\n)\n\n\nCreate a copy of an ODESystem with the given changes.\n\n\n\n\n\n","category":"method"},{"location":"api/#EarthSciMLBase.couple-Tuple","page":"API Reference","title":"EarthSciMLBase.couple","text":"couple(systems...) -> CoupledSystem\n\n\nCouple multiple ModelingToolkit systems together.\n\nThe systems that are arguments to this system can be of type ModelingToolkit.AbstractSystem, CoupledSystem, DomainInfo, or any type T that has a method couple(::CoupledSystem, ::T)::CoupledSystem or a method couple(::T, ::CoupledSystem)::CoupledSystem defined for it.\n\n\n\n\n\n","category":"method"},{"location":"api/#EarthSciMLBase.couple2-Tuple{}","page":"API Reference","title":"EarthSciMLBase.couple2","text":"couple2()\n\n\nPerform bi-directional coupling for two equation systems.\n\nTo specify couplings for system pairs, create methods for this function with the signature:\n\nEarthSciMLBase.couple2(a::ACoupler, b::BCoupler)::ConnectorSystem\n\nwhere ACoupler and BCoupler are :coupletypes defined like this:\n\nstruct ACoupler sys end\n@named asys = ODESystem([], t, metadata=Dict(:coupletype=>ACoupler))\n\n\n\n\n\n","category":"method"},{"location":"api/#EarthSciMLBase.dims-Tuple{EarthSciMLBase.ICcomponent}","page":"API Reference","title":"EarthSciMLBase.dims","text":"dims(\n icbc::EarthSciMLBase.ICcomponent\n) -> Vector{Symbolics.Num}\n\n\nReturns the dimensions of the independent and partial domains associated with these initial or boundary conditions.\n\n\n\n\n\n","category":"method"},{"location":"api/#EarthSciMLBase.domains-Tuple{EarthSciMLBase.ICcomponent}","page":"API Reference","title":"EarthSciMLBase.domains","text":"domains(icbc::EarthSciMLBase.ICcomponent) -> Vector\n\n\nReturns the domains associated with these initial or boundary conditions.\n\n\n\n\n\n","category":"method"},{"location":"api/#EarthSciMLBase.dtype-Union{Tuple{DomainInfo{T}}, Tuple{T}} where T","page":"API Reference","title":"EarthSciMLBase.dtype","text":"dtype(_)\n\n\nReturn the data type of the state variables for this domain, based on the data types of the boundary conditions domain intervals.\n\n\n\n\n\n","category":"method"},{"location":"api/#EarthSciMLBase.endpoints-Union{Tuple{DomainInfo{T}}, Tuple{T}} where T","page":"API Reference","title":"EarthSciMLBase.endpoints","text":"endpoints(d)\n\n\nReturn the endpoints of the partial independent variables for this domain.\n\n\n\n\n\n","category":"method"},{"location":"api/#EarthSciMLBase.get_coupletype-Tuple{ModelingToolkit.AbstractSystem}","page":"API Reference","title":"EarthSciMLBase.get_coupletype","text":"Return the coupling type associated with the given system.\n\n\n\n\n\n","category":"method"},{"location":"api/#EarthSciMLBase.get_dv-Tuple{Any, Any}","page":"API Reference","title":"EarthSciMLBase.get_dv","text":"Return the dependent variable, which is the first argument of the term, unless the term is a time derivative, in which case the dependent variable is the argument of the time derivative.\n\n\n\n\n\n","category":"method"},{"location":"api/#EarthSciMLBase.get_needed_vars-Tuple{ModelingToolkit.ODESystem}","page":"API Reference","title":"EarthSciMLBase.get_needed_vars","text":"get_needed_vars(sys)\n\n\nReturn the indexes of the system variables that the state variables of the final simplified system depend on. This should be done before running structural_simplify on the system.\n\n\n\n\n\n","category":"method"},{"location":"api/#EarthSciMLBase.graph-Tuple{CoupledSystem}","page":"API Reference","title":"EarthSciMLBase.graph","text":"Create a graph from a CoupledSystem using the MetaGraphsNext package.\n\n\n\n\n\n","category":"method"},{"location":"api/#EarthSciMLBase.grid-Union{Tuple{DomainInfo{T}}, Tuple{T}} where T","page":"API Reference","title":"EarthSciMLBase.grid","text":"grid(d)\n\n\nReturn the ranges representing the discretization of the partial independent variables for this domain, based on the discretization intervals given in Δs.\n\n\n\n\n\n","category":"method"},{"location":"api/#EarthSciMLBase.icbc-Tuple{DomainInfo, AbstractVector}","page":"API Reference","title":"EarthSciMLBase.icbc","text":"icbc(di, states)\n\n\nReturn a vector of equations that define the initial and boundary conditions for the given state variables.\n\n\n\n\n\n","category":"method"},{"location":"api/#EarthSciMLBase.init_callback-Tuple{}","page":"API Reference","title":"EarthSciMLBase.init_callback","text":"Types that implement an:\n\ninit_callback(x, sys::CoupledSystem, obs_eqs, domain::DomainInfo)::DECallback\n\nmethod can also be coupled into a CoupledSystem. The init_callback function will be run before the simulator is run to get the callback.\n\n\n\n\n\n","category":"method"},{"location":"api/#EarthSciMLBase.init_u-Tuple{ModelingToolkit.ODESystem, DomainInfo}","page":"API Reference","title":"EarthSciMLBase.init_u","text":"Initialize the state variables.\n\n\n\n\n\n","category":"method"},{"location":"api/#EarthSciMLBase.ivar-Tuple{DomainInfo}","page":"API Reference","title":"EarthSciMLBase.ivar","text":"ivar(di::DomainInfo) -> Any\n\n\nReturn the independent variable associated with these initial and boundary conditions.\n\n\n\n\n\n","category":"method"},{"location":"api/#EarthSciMLBase.observed_expression-Tuple{Any, Any}","page":"API Reference","title":"EarthSciMLBase.observed_expression","text":"observed_expression(eqs, x)\n\n\nReturn an expression for the observed value of a variable x after substituting in the constants observed values of other variables. extra_eqs is a list of additional equations to use in the substitution.\n\n\n\n\n\n","category":"method"},{"location":"api/#EarthSciMLBase.observed_function-Tuple{Any, Any, Any}","page":"API Reference","title":"EarthSciMLBase.observed_function","text":"observed_function(eqs, x, coords)\n\n\nReturn a function to for the observed value of a variable x based on the input arguments in coords. extra_eqs is a list of additional equations to use to determine the value of x.\n\n\n\n\n\n","category":"method"},{"location":"api/#EarthSciMLBase.ode_step!-Tuple{Any, SolverStrangThreads, Vararg{Any, 5}}","page":"API Reference","title":"EarthSciMLBase.ode_step!","text":"Take a step using the ODE solver.\n\n\n\n\n\n","category":"method"},{"location":"api/#EarthSciMLBase.operator_compose","page":"API Reference","title":"EarthSciMLBase.operator_compose","text":"operator_compose(a, b)\noperator_compose(a, b, translate)\n\n\nCompose to systems of equations together by adding the right-hand side terms together of equations that have matching left-hand sides. The left hand sides of two equations will be considered matching if:\n\nThey are both time derivatives of the same variable.\nThe first one is a time derivative of a variable and the second one is the variable itself.\nThere is an entry in the optional translate dictionary that maps the dependent variable in the first system to the dependent variable in the second system, e.g. Dict(sys1.sys.x => sys2.sys.y).\nThere is an entry in the optional translate dictionary that maps the dependent variable in the first system to the dependent variable in the second system, with a conversion factor, e.g. Dict(sys1.sys.x => sys2.sys.y => 6).\n\n\n\n\n\n","category":"function"},{"location":"api/#EarthSciMLBase.param_to_var-Tuple{ModelingToolkit.AbstractSystem, Vararg{Symbol}}","page":"API Reference","title":"EarthSciMLBase.param_to_var","text":"Replace the parameter p in the system sys with a new variable that has the same name, units, and description as p.\n\nparam_to_var(sys, ps)\n\n\nThis can be useful to replace a parameter that does not change in time in a model component with one specified by another system that does change in time (or space). For example, the code below specifies a first-order loss equation, and then changes the temperature (which determines the loss rate) with a temperature value that varies in time. ```\n\n\n\n\n\n","category":"method"},{"location":"api/#EarthSciMLBase.partialderivative_transforms-Tuple{DomainInfo}","page":"API Reference","title":"EarthSciMLBase.partialderivative_transforms","text":"partialderivative_transforms(di::DomainInfo) -> Vector{Any}\n\n\nReturn transform factor to multiply each partial derivative operator by, for example to convert from degrees to meters.\n\n\n\n\n\n","category":"method"},{"location":"api/#EarthSciMLBase.partialderivatives-Tuple{DomainInfo}","page":"API Reference","title":"EarthSciMLBase.partialderivatives","text":"partialderivatives(di::DomainInfo) -> Any\n\n\nReturn the partial derivative operators for the given domain.\n\n\n\n\n\n","category":"method"},{"location":"api/#EarthSciMLBase.partialderivatives_δxyδlonlat-Tuple{AbstractVector}","page":"API Reference","title":"EarthSciMLBase.partialderivatives_δxyδlonlat","text":"partialderivatives_δxyδlonlat(pvars; default_lat)\n\n\nReturn partial derivative operator transform factors corresponding for the given partial-independent variables after converting variables named lon and lat from degrees to x and y meters, assuming they represent longitude and latitude on a spherical Earth.\n\n\n\n\n\n","category":"method"},{"location":"api/#EarthSciMLBase.prune_observed-Tuple{ModelingToolkit.ODESystem}","page":"API Reference","title":"EarthSciMLBase.prune_observed","text":"prune_observed(sys)\n\n\nRemove equations from an ODESystem where the variable in the LHS is not present in any of the equations for the state variables. This can be used to remove computationally intensive equations that are not used in the final model.\n\n\n\n\n\n","category":"method"},{"location":"api/#EarthSciMLBase.pvars-Tuple{DomainInfo}","page":"API Reference","title":"EarthSciMLBase.pvars","text":"pvars(di::DomainInfo) -> Any\n\n\nReturn the partial independent variables associated with these initial and boundary conditions.\n\n\n\n\n\n","category":"method"},{"location":"api/#EarthSciMLBase.single_ode_step!-NTuple{6, Any}","page":"API Reference","title":"EarthSciMLBase.single_ode_step!","text":"Take a step using the ODE solver with the given IIchunk (grid cell interator) and integrator.\n\n\n\n\n\n","category":"method"},{"location":"api/#EarthSciMLBase.steplength-Tuple{Any}","page":"API Reference","title":"EarthSciMLBase.steplength","text":"steplength(timesteps)\n\n\nReturn the time step length common to all of the given timesteps. Throw an error if not all timesteps are the same length.\n\n\n\n\n\n","category":"method"},{"location":"api/#EarthSciMLBase.stiff_callback-Union{Tuple{T}, Tuple{Any, AbstractArray{T}, EarthSciMLBase.SolverStrang, Any, Any}} where T","page":"API Reference","title":"EarthSciMLBase.stiff_callback","text":"A callback to periodically run the stiff solver.\n\n\n\n\n\n","category":"method"},{"location":"api/#EarthSciMLBase.threaded_ode_step!-NTuple{6, Any}","page":"API Reference","title":"EarthSciMLBase.threaded_ode_step!","text":"Take a step using the ODE solver.\n\n\n\n\n\n","category":"method"},{"location":"api/#EarthSciMLBase.timesteps-Union{Tuple{Vararg{AbstractVector{T}}}, Tuple{T}} where T<:AbstractFloat","page":"API Reference","title":"EarthSciMLBase.timesteps","text":"timesteps(tsteps)\n\n\nReturn the time points during which integration should be stopped to run the operators.\n\n\n\n\n\n","category":"method"},{"location":"api/#EarthSciMLBase.tspan-Union{Tuple{DomainInfo{T}}, Tuple{T}} where T<:AbstractFloat","page":"API Reference","title":"EarthSciMLBase.tspan","text":"tspan(d)\n\n\nReturn the time range associated with this domain.\n\n\n\n\n\n","category":"method"},{"location":"icbc/#Initial-and-Boundary-conditions","page":"Initial and Boundary Conditions","title":"Initial and Boundary conditions","text":"","category":"section"},{"location":"icbc/","page":"Initial and Boundary Conditions","title":"Initial and Boundary Conditions","text":"Oftentimes we will want to do a 1, 2, or 3-dimensional simulation, rather than the 0-dimensional simulation we get by default with a system of ordinary differential equations. In these cases, we will need to specify initial and boundary conditions for the system.","category":"page"},{"location":"icbc/","page":"Initial and Boundary Conditions","title":"Initial and Boundary Conditions","text":"To demonstrate how to do this, we will use the following simple system of ordinary differential equations:","category":"page"},{"location":"icbc/","page":"Initial and Boundary Conditions","title":"Initial and Boundary Conditions","text":"using EarthSciMLBase\nusing ModelingToolkit\nusing ModelingToolkit: t_nounits, D_nounits\nt = t_nounits\nD = D_nounits\n\n@parameters x y\n\nfunction ExampleSys()\n @variables u(t) v(t)\n eqs = [\n D(u) ~ √abs(v),\n D(v) ~ √abs(u),\n ]\n ODESystem(eqs, t; name=:Docs₊Example)\nend\n\nExampleSys()","category":"page"},{"location":"icbc/","page":"Initial and Boundary Conditions","title":"Initial and Boundary Conditions","text":"Next, we specify our initial and boundary conditions using the DomainInfo type. We initialize DomainInfo with sets of initial and boundary conditions. In the example below, we set constant initial conditions using constIC and constant boundary conditions using constBC.","category":"page"},{"location":"icbc/","page":"Initial and Boundary Conditions","title":"Initial and Boundary Conditions","text":"using DomainSets\n\nx_min = y_min = t_min = 0.0\nx_max = y_max = t_max = 1.0\n\n# Create constant initial conditions = 16.0 and boundary conditions = 4.0.\nicbc = DomainInfo(\n constIC(4.0, t ∈ Interval(t_min, t_max)),\n constBC(16.0,\n x ∈ Interval(x_min, x_max),\n y ∈ Interval(y_min, y_max),\n ),\n)\nnothing # hide","category":"page"},{"location":"icbc/","page":"Initial and Boundary Conditions","title":"Initial and Boundary Conditions","text":"It is also possible to use periodic boundary conditions with periodicBC and zero-gradient boundary conditions with zerogradBC.","category":"page"},{"location":"icbc/","page":"Initial and Boundary Conditions","title":"Initial and Boundary Conditions","text":"Finally, we combine our initial and boundary conditions with our system of equations using the couple function.","category":"page"},{"location":"icbc/","page":"Initial and Boundary Conditions","title":"Initial and Boundary Conditions","text":"model = couple(ExampleSys(), icbc)\n\neq_sys = convert(PDESystem, model)","category":"page"},{"location":"icbc/","page":"Initial and Boundary Conditions","title":"Initial and Boundary Conditions","text":"We can also look at the expanded boundary conditions of the resulting equation system:","category":"page"},{"location":"icbc/","page":"Initial and Boundary Conditions","title":"Initial and Boundary Conditions","text":"eq_sys.bcs","category":"page"},{"location":"operator_compose/#Operator-Composition","page":"Operator Composition","title":"Operator Composition","text":"","category":"section"},{"location":"operator_compose/","page":"Operator Composition","title":"Operator Composition","text":"There are a lot of cases where there are two different \"processes\" or \"operators\" that change the same variable. For example, CO2 in the atmosphere can be emitted by human activity, and it can also be absorbed by the ocean. In models, typically the emission and removal are considered separate processes which are represented by separate model components. However, when we want to combine these two components into a single model, we need to be able to compose them together.","category":"page"},{"location":"operator_compose/","page":"Operator Composition","title":"Operator Composition","text":"We can use the operator_compose function for this. It composes to systems of equations together by adding the right-hand side terms together of equations that have matching left-hand sides. The left hand sides of two equations will be considered matching if:","category":"page"},{"location":"operator_compose/","page":"Operator Composition","title":"Operator Composition","text":"They are both time derivatives of the same variable.\nThe first one is a time derivative of a variable and the second one is the variable itself.\nThere is an entry in the optional translate dictionary that maps the dependent variable in the first system to the dependent variable in the second system, e.g. Dict(sys1.sys.x => sys2.sys.y).\nThere is an entry in the optional translate dictionary that maps the dependent variable in the first system to the dependent variable in the second system, with a conversion factor, e.g. Dict(sys1.sys.x => sys2.sys.y => 6).","category":"page"},{"location":"operator_compose/","page":"Operator Composition","title":"Operator Composition","text":"Perhaps we can make this somewhat clearer with some examples.","category":"page"},{"location":"operator_compose/#Examples","page":"Operator Composition","title":"Examples","text":"","category":"section"},{"location":"operator_compose/#Example-with-matching-variable-time-derivatives","page":"Operator Composition","title":"Example with matching variable time derivatives","text":"","category":"section"},{"location":"operator_compose/","page":"Operator Composition","title":"Operator Composition","text":"The example below shows that when we operator_compose two systems together that are both equal to D(x) = p, the resulting system is equal to D(x) = 2p.","category":"page"},{"location":"operator_compose/","page":"Operator Composition","title":"Operator Composition","text":"using EarthSciMLBase\nusing ModelingToolkit\nusing ModelingToolkit: t_nounits, D_nounits\nt = t_nounits\nD = D_nounits\n\n\nstruct ExampleSysCoupler sys end\nfunction ExampleSys()\n @variables x(t)\n @parameters p\n ODESystem([D(x) ~ p], t; name=:ExampleSys,\n metadata=Dict(:coupletype=>ExampleSysCoupler))\nend\n\nExampleSys()","category":"page"},{"location":"operator_compose/","page":"Operator Composition","title":"Operator Composition","text":"struct ExampleSys2Coupler sys end\nfunction ExampleSys2()\n @variables x(t)\n @parameters p\n ODESystem([D(x) ~ 2p], t; name=:ExampleSys2,\n metadata=Dict(:coupletype=>ExampleSys2Coupler))\nend\n\nExampleSys2()","category":"page"},{"location":"operator_compose/","page":"Operator Composition","title":"Operator Composition","text":"sys1 = ExampleSys()\nsys2 = ExampleSys2()\n\nfunction EarthSciMLBase.couple2(sys1::ExampleSysCoupler, sys2::ExampleSys2Coupler)\n sys1, sys2 = sys1.sys, sys2.sys\n operator_compose(sys1, sys2)\nend\n\ncombined = couple(sys1, sys2)\n\ncombined_mtk = convert(ODESystem, combined)","category":"page"},{"location":"operator_compose/","page":"Operator Composition","title":"Operator Composition","text":"The simplified equation should be D(x) = p + sys2_xˍt:","category":"page"},{"location":"operator_compose/","page":"Operator Composition","title":"Operator Composition","text":"combined_simplified = structural_simplify(combined_mtk)","category":"page"},{"location":"operator_compose/","page":"Operator Composition","title":"Operator Composition","text":"where sys2_xˍt is also equal to p:","category":"page"},{"location":"operator_compose/","page":"Operator Composition","title":"Operator Composition","text":"observed(combined_simplified)","category":"page"},{"location":"operator_compose/#Example-with-non-matching-variables","page":"Operator Composition","title":"Example with non-matching variables","text":"","category":"section"},{"location":"operator_compose/","page":"Operator Composition","title":"Operator Composition","text":"This example demonstrates a case where one variable in the first system is equal to another variable in the second system:","category":"page"},{"location":"operator_compose/","page":"Operator Composition","title":"Operator Composition","text":"struct ExampleSys3Coupler sys end\nfunction ExampleSys3()\n @variables y(t)\n @parameters p\n ODESystem([D(y) ~ p], t; name=:ExampleSys3,\n metadata=Dict(:coupletype=>ExampleSys3Coupler))\nend\n\nsys1 = ExampleSys()\nsys2 = ExampleSys3()\n\nfunction EarthSciMLBase.couple2(sys1::ExampleSysCoupler, sys2::ExampleSys3Coupler)\n sys1, sys2 = sys1.sys, sys2.sys\n operator_compose(sys1, sys2, Dict(sys1.x => sys2.y))\nend\n\ncombined = couple(sys1, sys2)\ncombined_simplified = structural_simplify(convert(ODESystem, combined))","category":"page"},{"location":"operator_compose/","page":"Operator Composition","title":"Operator Composition","text":"observed(combined_simplified)","category":"page"},{"location":"operator_compose/#Example-with-a-non-ODE-system","page":"Operator Composition","title":"Example with a non-ODE system","text":"","category":"section"},{"location":"operator_compose/","page":"Operator Composition","title":"Operator Composition","text":"In the second case above, we might have a variable in the second system that is equal to a rate, but it is not a time derivative. This could happen if we are extracting emissions from a file, and those emissions are already in units of kg/s, or something similar. The example below demonstrates this case. (Note that this case can also be used with the conversion factors shown in the last example.)","category":"page"},{"location":"operator_compose/","page":"Operator Composition","title":"Operator Composition","text":"struct ExampleSysNonODECoupler sys end\nfunction ExampleSysNonODE()\n @variables y(t)\n @parameters p\n ODESystem([y ~ p], t; name=:ExampleSysNonODE,\n metadata=Dict(:coupletype=>ExampleSysNonODECoupler))\nend\n\nsys1 = ExampleSys()\nsys2 = ExampleSysNonODE()\n\nfunction EarthSciMLBase.couple2(sys1::ExampleSysCoupler, sys2::ExampleSysNonODECoupler)\n sys1, sys2 = sys1.sys, sys2.sys\n operator_compose(sys1, sys2, Dict(sys1.x => sys2.y))\nend\n\ncombined = couple(sys1, sys2)\nsys_combined = structural_simplify(convert(ODESystem, combined))","category":"page"},{"location":"operator_compose/","page":"Operator Composition","title":"Operator Composition","text":"observed(sys_combined)","category":"page"},{"location":"operator_compose/#Example-with-non-matching-variables-and-a-conversion-factor","page":"Operator Composition","title":"Example with non-matching variables and a conversion factor","text":"","category":"section"},{"location":"operator_compose/","page":"Operator Composition","title":"Operator Composition","text":"Finally, this last example shows the fourth case, where a conversion factor is included in the translation dictionary.","category":"page"},{"location":"operator_compose/","page":"Operator Composition","title":"Operator Composition","text":"struct ExampleSysNonODE2Coupler sys end\nfunction ExampleSysNonODE2()\n @variables y(t)\n @parameters p\n ODESystem([y ~ p], t; name=:Docs₊ExampleSysNonODE2,\n metadata=Dict(:coupletype=>ExampleSysNonODE2Coupler))\nend\n\nsys1 = ExampleSys()\nsys2 = ExampleSysNonODE2()\n\nfunction EarthSciMLBase.couple2(sys1::ExampleSysCoupler, sys2::ExampleSysNonODE2Coupler)\n sys1, sys2 = sys1.sys, sys2.sys\n operator_compose(sys1, sys2, Dict(sys1.x => sys2.y => 6.0))\nend\n\ncombined = couple(sys1, sys2)\ncombined_simplified = structural_simplify(convert(ODESystem, combined))","category":"page"},{"location":"operator_compose/","page":"Operator Composition","title":"Operator Composition","text":"observed(combined_simplified)","category":"page"},{"location":"param_to_var/#Converting-parameters-to-variables","page":"Parameter Replacement","title":"Converting parameters to variables","text":"","category":"section"},{"location":"param_to_var/","page":"Parameter Replacement","title":"Parameter Replacement","text":"This can be useful to replace a parameter that does not change in time in a model component with one specified by another system that does change in time (or space). For example, the code below specifies a first-order loss equation, and then changes the temperature (which determines the loss rate) with a temperature value that varies in time.","category":"page"},{"location":"param_to_var/","page":"Parameter Replacement","title":"Parameter Replacement","text":"As an example, we will create a loss equation that depends on the temperature, starting with a constant temperature. We will then create a temperature equation that varies in time, and use the param_to_var function to replace the constant temperature in the loss equation with the time-varying temperature.","category":"page"},{"location":"param_to_var/","page":"Parameter Replacement","title":"Parameter Replacement","text":"So first, let's specify the original system with constant temperature.","category":"page"},{"location":"param_to_var/","page":"Parameter Replacement","title":"Parameter Replacement","text":"using ModelingToolkit, EarthSciMLBase, DynamicQuantities\nusing ModelingToolkit: t, D\n\nstruct LossCoupler sys end\nfunction Loss()\n @variables A(t)=1 [unit=u\"kg\"]\n @parameters k=1 [unit=u\"s^-1\"]\n @parameters T=300 [unit=u\"K\"]\n @constants T₀=300 [unit=u\"K\"]\n eq = D(A) ~ -k*exp(T/T₀) * A\n ODESystem([eq], t; name=:Loss, metadata=Dict(:coupletype=>LossCoupler))\nend\n\nLoss()","category":"page"},{"location":"param_to_var/","page":"Parameter Replacement","title":"Parameter Replacement","text":"Next, we specify the temperature that varies in time.","category":"page"},{"location":"param_to_var/","page":"Parameter Replacement","title":"Parameter Replacement","text":"struct TemperatureCoupler sys end\nfunction Temperature()\n @variables T(t)=300 [unit=u\"K\"]\n @constants Tc=1.0 [unit=u\"K/s\"]\n @constants tc=1.0 [unit=u\"s\"]\n eq = D(T) ~ sin(t/tc)*Tc\n ODESystem([eq], t; name=:Temperature, metadata=Dict(:coupletype=>TemperatureCoupler))\nend\n\nTemperature()","category":"page"},{"location":"param_to_var/","page":"Parameter Replacement","title":"Parameter Replacement","text":"Now, we specify how to compose the two systems using param_to_var.","category":"page"},{"location":"param_to_var/","page":"Parameter Replacement","title":"Parameter Replacement","text":"function EarthSciMLBase.couple2(loss::LossCoupler, temp::TemperatureCoupler)\n loss, temp = loss.sys, temp.sys\n loss = param_to_var(loss, :T)\n ConnectorSystem([loss.T ~ temp.T], loss, temp)\nend","category":"page"},{"location":"param_to_var/","page":"Parameter Replacement","title":"Parameter Replacement","text":"Finally, we create the system components and the composed system.","category":"page"},{"location":"param_to_var/","page":"Parameter Replacement","title":"Parameter Replacement","text":"l = Loss()\ntemp = Temperature()\nvariable_loss = couple(l, temp)\n\nconvert(ODESystem, variable_loss)","category":"page"},{"location":"param_to_var/","page":"Parameter Replacement","title":"Parameter Replacement","text":"If we wanted to, we could then run a simulation with the composed system.","category":"page"},{"location":"operator/#Operators-for-Large-scale-3D-simulations","page":"Operators","title":"Operators for Large-scale 3D simulations","text":"","category":"section"},{"location":"operator/","page":"Operators","title":"Operators","text":"In this documentation so far, we have talked about creating systems of ordinary differential equations in ModelingToolkit and then converting them to systems of partial differential equations to perform 1-, 2-, or 3-dimensional simulations. However, currently this does not work for large scale simulations.","category":"page"},{"location":"operator/","page":"Operators","title":"Operators","text":"While this ModelingToolkit functionality is being built, we have a different solution based on the Operator type in this package. Using this system, we still define systems of ODEs to define behavior in a single grid cell, and we also have Operator processes that define behavior between grid cells.","category":"page"},{"location":"operator/#ODE-System","page":"Operators","title":"ODE System","text":"","category":"section"},{"location":"operator/","page":"Operators","title":"Operators","text":"As an example, let's first define a system of ODEs:","category":"page"},{"location":"operator/","page":"Operators","title":"Operators","text":"using EarthSciMLBase\nusing ModelingToolkit, DifferentialEquations\nusing SciMLOperators, Plots\nusing DomainSets\nusing ModelingToolkit: t_nounits, D_nounits\nt = t_nounits\nD = D_nounits\n\n@parameters y lon = 0.0 lat = 0.0 lev = 1.0 α = 10.0\n@constants p = 1.0\n@variables(\n u(t) = 1.0, v(t) = 1.0, x(t) = 1.0, y(t) = 1.0, windspeed(t) = 1.0\n)\n\neqs = [D(u) ~ -α * √abs(v) + lon,\n D(v) ~ -α * √abs(u) + lat + 1e-14 * lev,\n windspeed ~ lat + lon + lev,\n]\nsys = ODESystem(eqs, t; name=:sys)","category":"page"},{"location":"operator/","page":"Operators","title":"Operators","text":"The equations above don't really have any physical meaning, but they include two state variables, some parameters, and a constant. There is also a variable windspeed which is \"observed\" based on the parameters, rather than being a state variable, which will be important later.","category":"page"},{"location":"operator/#Operator","page":"Operators","title":"Operator","text":"","category":"section"},{"location":"operator/","page":"Operators","title":"Operators","text":"Next, we define an operator. To do so, first we create a new type that is a subtype of Operator:","category":"page"},{"location":"operator/","page":"Operators","title":"Operators","text":"mutable struct ExampleOp <: Operator\n α::Num # Multiplier from ODESystem\nend","category":"page"},{"location":"operator/","page":"Operators","title":"Operators","text":"In the case above, we're setting up our operator so that it can hold a parameter from our ODE system.","category":"page"},{"location":"operator/","page":"Operators","title":"Operators","text":"Next, we need to define a method of EarthSciMLBase.get_scimlop for our operator. This method will be called to get a SciMLOperators.AbstractSciMLOperator that will be used conjunction with the ModelingToolkit system above to integrate the simulation forward in time.","category":"page"},{"location":"operator/","page":"Operators","title":"Operators","text":"function EarthSciMLBase.get_scimlop(op::ExampleOp, csys::CoupledSystem, mtk_sys, domain::DomainInfo, obs_functions, coordinate_transform_functions, u0, p)\n obs_f = obs_functions(op.α)\n grd = EarthSciMLBase.grid(domain)\n function run(du, u, p, t)\n u = reshape(u, size(u0)...)\n du = reshape(du, size(u0)...)\n for ix ∈ 1:size(u, 1)\n for (i, c1) ∈ enumerate(grd[1])\n for (j, c2) ∈ enumerate(grd[2])\n for (k, c3) ∈ enumerate(grd[3])\n # Demonstrate coordinate transforms\n t1 = coordinate_transform_functions[1](t, c1, c2, c3)\n t2 = coordinate_transform_functions[2](t, c1, c2, c3)\n t3 = coordinate_transform_functions[3](t, c1, c2, c3)\n # Demonstrate calculating observed value.\n fv = obs_f(t, c1, c2, c3)\n # Set derivative value.\n du[ix, i, j, k] = (t1 + t2 + t3) * fv\n end\n end\n end\n end\n nothing\n end\n FunctionOperator(run, u0[:], p=p)\nend","category":"page"},{"location":"operator/","page":"Operators","title":"Operators","text":"The function above also doesn't have any physical meaning, but it demonstrates some functionality of the Operator \"s\". First, it retrieves a function to get the current value of an observed variable in our ODE system using the obs_functions argement, and it demonstrates how to call the resulting function to get that value. It also demonstrates how to get coordinate transforms using the coordinate_transform_functions argument. Coordinate transforms are discussed in more detail in the documentation for the DomainInfo type.","category":"page"},{"location":"operator/#Domain","page":"Operators","title":"Domain","text":"","category":"section"},{"location":"operator/","page":"Operators","title":"Operators","text":"Once we have an ODE system and an operator, the final component we need is a domain to run the simulation on. Defining a domain is covered in more depth in the documentation for the DomainInfo type, but for now we'll just define a simple domain:","category":"page"},{"location":"operator/","page":"Operators","title":"Operators","text":"t_min = 0.0\nlon_min, lon_max = -π, π\nlat_min, lat_max = -0.45π, 0.45π\nt_max = 11.5\n\nindepdomain = t ∈ Interval(t_min, t_max)\n\npartialdomains = [lon ∈ Interval(lon_min, lon_max),\n lat ∈ Interval(lat_min, lat_max),\n lev ∈ Interval(1, 3)]\n\ndomain = DomainInfo(\n partialderivatives_δxyδlonlat,\n constIC(16.0, indepdomain), constBC(16.0, partialdomains...);\n grid_spacing = [0.1π, 0.1π, 1])\nnothing #hide","category":"page"},{"location":"operator/","page":"Operators","title":"Operators","text":"Note that our domain includes a coordinate transform to convert from degrees latitude and longitude to meters. Our domain specification also includes grid spacing the the lon, lat, and lev coordinates, which we set as 0.1π, 0.1π, and 1, respectively.","category":"page"},{"location":"operator/","page":"Operators","title":"Operators","text":"warning: Warning\nInitial and boundary conditions are not fully implemented for this case, so regardless of the conditions you specify, the initial conditions will be the default values of the variables in the ODE system, and the boundary conditions will be zero.","category":"page"},{"location":"operator/#Coupling-and-Running-the-Simulation","page":"Operators","title":"Coupling and Running the Simulation","text":"","category":"section"},{"location":"operator/","page":"Operators","title":"Operators","text":"Next, initialize our operator, giving the the windspeed observed variable, and we can couple our ODESystem, Operator, and Domain together into a single model:","category":"page"},{"location":"operator/","page":"Operators","title":"Operators","text":"op = ExampleOp(sys.windspeed)\n\ncsys = couple(sys, op, domain)","category":"page"},{"location":"operator/","page":"Operators","title":"Operators","text":"Finally, we can choose a EarthSciMLBase.SolverStrategy and run the simulation. We choose the SolverStrangThreads strategy, which needs us to specify an ODE solver from the options available in DifferentialEquations.jl for both the MTK system. We choose the Tsit5 solver. Then we create an ODEProblem which can be used to run the simulation. Finally, we solve the problem using the solve function. At this point we need to choose a solver for the Operator part of the system, and we choose the Euler solver. We also choose a splitting time step of 1.0 seconds, which we pass both to our SolverStrangThreads strategy and to the solve function.","category":"page"},{"location":"operator/","page":"Operators","title":"Operators","text":"dt = 1.0 # Splitting time step\nst = SolverStrangThreads(Tsit5(), 1.0)\n\nprob = ODEProblem(csys, st)\nsol = solve(prob, Euler(); dt=1.0)\nnothing #hide","category":"page"},{"location":"operator/","page":"Operators","title":"Operators","text":"After the simulation finishes, we can plot the result:","category":"page"},{"location":"operator/","page":"Operators","title":"Operators","text":"anim = @animate for i ∈ 1:length(sol.u)\n u = reshape(sol.u[i], :, size(domain)...)\n plot(\n heatmap(u[1, :, :, 1]),\n heatmap(u[1, :, :, 1]),\n )\nend\ngif(anim, fps = 15)","category":"page"},{"location":"coord_transforms/#Coordinate-Transforms-for-Partial-Derivatives","page":"Coordinate Transforms","title":"Coordinate Transforms for Partial Derivatives","text":"","category":"section"},{"location":"coord_transforms/","page":"Coordinate Transforms","title":"Coordinate Transforms","text":"Often, the coordinates of a grid may be defined in a different coordinate system than the one in which the partial derivatives are desired. For example, grids are often defined in latitude and longitude, but partial derivatives may be required in units of meters to correspond with wind speeds in meters per second.","category":"page"},{"location":"coord_transforms/","page":"Coordinate Transforms","title":"Coordinate Transforms","text":"To handle this, the DomainInfo type can be used to define coordinate system transformations. To use it, a coordinate transform function first needs to be defined, for example partialderivatives_δxyδlonlat which transforms partial derivatives from longitude and latitude to meters:","category":"page"},{"location":"coord_transforms/","page":"Coordinate Transforms","title":"Coordinate Transforms","text":"using EarthSciMLBase\nusing ModelingToolkit\nusing ModelingToolkit: t, D\nusing DomainSets\nusing DynamicQuantities\n\n@parameters lon [unit = u\"rad\"]\n@parameters lat [unit = u\"rad\"]\n@parameters lev\n\npartialderivatives_δxyδlonlat([lev, lon, lat])","category":"page"},{"location":"coord_transforms/","page":"Coordinate Transforms","title":"Coordinate Transforms","text":"As you can see in the output of the code above, the function should take as arguments a list of the coordinates describing the grid (in the case above we have a 3-dimensional grid with vertical level, latitude, and longitude), and return a Dictionary relating the index of each coordinate with a factor to multiply the partial derivative by to convert it to the desired units. The function only needs to return factors for the coordinates that are being transformed.","category":"page"},{"location":"coord_transforms/","page":"Coordinate Transforms","title":"Coordinate Transforms","text":"To include a coordinate transform in our domain, we include the function in the DomainInfo constructor:","category":"page"},{"location":"coord_transforms/","page":"Coordinate Transforms","title":"Coordinate Transforms","text":"deg2rad(x) = x * π / 180.0\ndomain = DomainInfo(\n partialderivatives_δxyδlonlat,\n constIC(0.0, t ∈ Interval(0.0f0, 3600.0f0)),\n periodicBC(lat ∈ Interval(deg2rad(-90.0f0), deg2rad(90.0f0))),\n periodicBC(lon ∈ Interval(deg2rad(-180.0f0), deg2rad(180.0f0))),\n zerogradBC(lev ∈ Interval(1.0f0, 10.0f0)),\n);","category":"page"},{"location":"coord_transforms/","page":"Coordinate Transforms","title":"Coordinate Transforms","text":"Multiple functions can be included in the DomainInfo constructor, just by including them as a vector, e.g.:","category":"page"},{"location":"coord_transforms/","page":"Coordinate Transforms","title":"Coordinate Transforms","text":"domain = DomainInfo(\n [transform1, transform2, ...],\n constIC(0.0, t ∈ Interval(0.0f0, 3600.0f0)),\n ...\n)","category":"page"},{"location":"coord_transforms/","page":"Coordinate Transforms","title":"Coordinate Transforms","text":"After we include the coordinate transforms in our domain, in general everything should be handled automatically. The coordinate transforms may also be automatically added when different model components are coupled together, so you may not need to worry about them at all in many cases. However, if you would like to use the transformed partial derivatives, for example to create a new PDE equation system, you can get them using the EarthSciMLBase.partialderivatives function:","category":"page"},{"location":"coord_transforms/","page":"Coordinate Transforms","title":"Coordinate Transforms","text":"δs = partialderivatives(domain)","category":"page"},{"location":"coord_transforms/","page":"Coordinate Transforms","title":"Coordinate Transforms","text":"This returns a list of functions, one corresponding to each coordinate in our domain. Then we can calculate the symbolic partial derivative of a variable by just calling each function:","category":"page"},{"location":"coord_transforms/","page":"Coordinate Transforms","title":"Coordinate Transforms","text":"@variables u(t)\n\n[δs[i](u) for i ∈ eachindex(δs)]","category":"page"},{"location":"coord_transforms/","page":"Coordinate Transforms","title":"Coordinate Transforms","text":"You can see an example of how this is implemented in the source code for the Advection model component.","category":"page"},{"location":"coord_transforms/","page":"Coordinate Transforms","title":"Coordinate Transforms","text":"Additional transformation functions may be defined in other packages. We recommend that the names of these functions start with partialderivatives_ to make it clear that they are intended to be used in this context.","category":"page"},{"location":"example_all_together/","page":"All Together","title":"All Together","text":"CurrentModule = EarthSciMLBase","category":"page"},{"location":"example_all_together/#Example-using-different-components-of-EarthSciMLBase-together","page":"All Together","title":"Example using different components of EarthSciMLBase together","text":"","category":"section"},{"location":"example_all_together/","page":"All Together","title":"All Together","text":"This example shows how to define and couple different components of EarthSciMLBase together to create a more complex model. First, we create several model components.","category":"page"},{"location":"example_all_together/","page":"All Together","title":"All Together","text":"Our first example system is a simple reaction system:","category":"page"},{"location":"example_all_together/","page":"All Together","title":"All Together","text":"using EarthSciMLBase\nusing ModelingToolkit, Catalyst, DomainSets, MethodOfLines, DifferentialEquations\nusing ModelingToolkit: t_nounits, D_nounits\nt = t_nounits\nD = D_nounits\nusing Plots\n\n# Create our independent variable `t` and our partially-independent variables `x` and `y`.\n@parameters x y\n\nstruct ExampleSys1Coupler sys end\nfunction ExampleSys1()\n @species c₁(t)=5.0 c₂(t)=5.0\n rs = ReactionSystem(\n [Reaction(2.0, [c₁], [c₂])],\n t; name=:Sys1, combinatoric_ratelaws=false)\n convert(ODESystem, complete(rs), metadata=Dict(:coupletype=>ExampleSys1Coupler))\nend\n\nExampleSys1()","category":"page"},{"location":"example_all_together/","page":"All Together","title":"All Together","text":"Our second example system is a simple ODE system, with the same two variables.","category":"page"},{"location":"example_all_together/","page":"All Together","title":"All Together","text":"struct ExampleSys2Coupler sys end\nfunction ExampleSys2()\n @variables c₁(t)=5.0 c₂(t)=5.0\n @parameters p₁=1.0 p₂=0.5\n ODESystem(\n [D(c₁) ~ -p₁, D(c₂) ~ p₂],\n t; name=:Sys2, metadata=Dict(:coupletype=>ExampleSys2Coupler))\nend\n\nExampleSys2()","category":"page"},{"location":"example_all_together/","page":"All Together","title":"All Together","text":"Now, we specify what should happen when we couple the two systems together. In this case, we want the the derivative of the composed system to be equal to the sum of the derivatives of the two systems. We can do that using the operator_compose function from this package.","category":"page"},{"location":"example_all_together/","page":"All Together","title":"All Together","text":"function EarthSciMLBase.couple2(sys1::ExampleSys1Coupler, sys2::ExampleSys2Coupler)\n sys1, sys2 = sys1.sys, sys2.sys\n sys1 = convert(ODESystem, sys1)\n operator_compose(sys1, sys2)\nend\nnothing # hide","category":"page"},{"location":"example_all_together/","page":"All Together","title":"All Together","text":"Once we specify all of the above, it is simple to create our two individual systems and then couple them together. ","category":"page"},{"location":"example_all_together/","page":"All Together","title":"All Together","text":"sys1 = ExampleSys1()\nsys2 = ExampleSys2()\nsys = couple(sys1, sys2)\n\nconvert(ODESystem, sys)","category":"page"},{"location":"example_all_together/","page":"All Together","title":"All Together","text":"At this point we have an ODE system that is composed of two other ODE systems. We can inspect its observed variables using the observed function.","category":"page"},{"location":"example_all_together/","page":"All Together","title":"All Together","text":"simplified_sys = structural_simplify(convert(ODESystem, sys))","category":"page"},{"location":"example_all_together/","page":"All Together","title":"All Together","text":"observed(simplified_sys)","category":"page"},{"location":"example_all_together/","page":"All Together","title":"All Together","text":"We can also run simulations using this system:","category":"page"},{"location":"example_all_together/","page":"All Together","title":"All Together","text":"odeprob = ODEProblem(simplified_sys, [], (0.0,10.0), [])\nodesol = solve(odeprob)\nplot(odesol)","category":"page"},{"location":"example_all_together/","page":"All Together","title":"All Together","text":"Once we've confirmed that our model works in a 0D \"box model\" setting, we can expand it to 1, 2, or 3 dimensions using by adding in initial and boundary conditions. We will also add in advection using constant-velocity wind fields add the same time.","category":"page"},{"location":"example_all_together/","page":"All Together","title":"All Together","text":"x_min = y_min = t_min = 0.0\nx_max = y_max = t_max = 1.0\ndomain = DomainInfo(\n constIC(4.0, t ∈ Interval(t_min, t_max)),\n periodicBC(x ∈ Interval(x_min, x_max)),\n zerogradBC(y ∈ Interval(y_min, y_max)),\n)\n\nsys_pde = couple(sys, domain, ConstantWind(t, 1.0, 1.0), Advection())\n\nsys_pde_mtk = convert(PDESystem, sys_pde)","category":"page"},{"location":"example_all_together/","page":"All Together","title":"All Together","text":"Now we can inspect this new system that we've created:","category":"page"},{"location":"example_all_together/","page":"All Together","title":"All Together","text":"sys_pde_mtk.dvs","category":"page"},{"location":"example_all_together/","page":"All Together","title":"All Together","text":"sys_pde_mtk.bcs","category":"page"},{"location":"example_all_together/","page":"All Together","title":"All Together","text":"Finally, we can run a simulation using this system:","category":"page"},{"location":"example_all_together/","page":"All Together","title":"All Together","text":"discretization = MOLFiniteDifference([x=>10, y=>10], t, approx_order=2)\n@time pdeprob = discretize(sys_pde_mtk, discretization)\n@time pdesol = solve(pdeprob, Tsit5(), saveat=0.1)\n\n# Plot the solution.\ndiscrete_x, discrete_y, discrete_t = pdesol[x], pdesol[y], pdesol[t]\n@variables Sys1₊c₁(..) Sys1₊c₂(..)\nsolc1, solc2 = pdesol[Sys1₊c₁(t, x, y)], pdesol[Sys1₊c₂(t, x, y)]\nanim = @animate for k in 1:length(discrete_t)\n p1 = heatmap(solc1[k, 1:end-1, 1:end-1], title=\"c₁ t=\\$(discrete_t[k])\", clim=(0,4.0), lab=:none)\n p2 = heatmap(solc2[k, 1:end-1, 1:end-1], title=\"c₂ t=\\$(discrete_t[k])\", clim=(0,7.0), lab=:none)\n plot(p1, p2, layout=(1,2), size=(800,400))\nend\ngif(anim, fps = 8)","category":"page"},{"location":"example_all_together/","page":"All Together","title":"All Together","text":"Because our system is a system of ordinary differential equations rather than partial differential equations, all of the grid cells in the animation above have the same value. Refer to the advection example for an example of a system of partial differential equations.","category":"page"},{"location":"composition/#Model-Composition","page":"Composition","title":"Model Composition","text":"","category":"section"},{"location":"composition/","page":"Composition","title":"Composition","text":"A main goal of the EarthSciMLBase package is to allow model components to be created independently and composed together. We achieve this by creating coupletypes that allow us to use multiple dispatch on the EarthSciMLBase.couple2 function to specify how particular systems should be coupled together.","category":"page"},{"location":"composition/","page":"Composition","title":"Composition","text":"To demonstrate how this works, below we define three model components, Photolysis, Chemistry, and Emissions, which represent different processes in the atmosphere.","category":"page"},{"location":"composition/","page":"Composition","title":"Composition","text":"using ModelingToolkit, Catalyst, EarthSciMLBase\nusing ModelingToolkit: t_nounits\nt = t_nounits\n\nstruct PhotolysisCoupler sys end\nfunction Photolysis(; name=:Photolysis)\n @variables j_NO2(t)\n eqs = [\n j_NO2 ~ max(sin(t/86400),0)\n ]\n ODESystem(eqs, t, [j_NO2], [], name=name,\n metadata=Dict(:coupletype=>PhotolysisCoupler))\nend\n\nPhotolysis()","category":"page"},{"location":"composition/","page":"Composition","title":"Composition","text":"You can see that the system defined above is mostly a standard ModelingToolkit ODE system, except for two things:","category":"page"},{"location":"composition/","page":"Composition","title":"Composition","text":"The first unique part is that we define a PhotolysisCoupler type:","category":"page"},{"location":"composition/","page":"Composition","title":"Composition","text":"struct PhotolysisCoupler sys end","category":"page"},{"location":"composition/","page":"Composition","title":"Composition","text":"It is important that this type is a struct, and that the struct has a single field named sys. When defining your own coupled systems, you can just copy the line of code above but change the name of the type (i.e., change it to something besides PhotolysisCoupler).","category":"page"},{"location":"composition/","page":"Composition","title":"Composition","text":"The second unique part is that we define some metadata for our ODESystem to tell it what coupling type to use:","category":"page"},{"location":"composition/","page":"Composition","title":"Composition","text":"metadata=Dict(:coupletype=>PhotolysisCoupler)","category":"page"},{"location":"composition/","page":"Composition","title":"Composition","text":"Again, when making your own components, just copy the code above but change PhotolysisCoupler to something else.","category":"page"},{"location":"composition/","page":"Composition","title":"Composition","text":"Let's follow the same process for some additional components:","category":"page"},{"location":"composition/","page":"Composition","title":"Composition","text":"struct ChemistryCoupler sys end\nfunction Chemistry(; name=:Chemistry)\n @parameters jNO2\n @species NO2(t)\n rxs = [\n Reaction(jNO2, [NO2], [], [1], [])\n ]\n rsys = ReactionSystem(rxs, t, [NO2], [jNO2]; \n combinatoric_ratelaws=false, name=name)\n convert(ODESystem, complete(rsys), metadata=Dict(:coupletype=>ChemistryCoupler))\nend\n\nChemistry()","category":"page"},{"location":"composition/","page":"Composition","title":"Composition","text":"For our chemistry component above, because it's is originally a ReactionSystem instead of an ODESystem, we convert it to an ODESystem before adding the metadata.","category":"page"},{"location":"composition/","page":"Composition","title":"Composition","text":"struct EmissionsCoupler sys end\nfunction Emissions(; name=:Emissions)\n @parameters emis = 1\n @variables NO2(t)\n eqs = [NO2 ~ emis]\n ODESystem(eqs, t; name=name,\n metadata=Dict(:coupletype=>EmissionsCoupler))\nend\n\nEmissions()","category":"page"},{"location":"composition/","page":"Composition","title":"Composition","text":"Now, we need to define ways to couple the model components together. We can do this by defining a coupling function (as a method of EarthSciMLBase.couple2) for each pair of model components that we want to couple. Each coupling function should have the signature EarthSciMLBase.couple2(a::ACoupler, b::BCoupler)::ConnectorSystem, and should assume that the two ODESystems are inside their respective couplers in the sys field. It should make any edits to the components as needed and return a ConnectorSystem which defines the relationship between the two components.","category":"page"},{"location":"composition/","page":"Composition","title":"Composition","text":"The code below defines a method for coupling the Chemistry and Photolysis components. First, it uses the param_to_var function to convert the photolysis rate parameter jNO2 from the Chemistry component to a variable, then it creates a new Chemistry component with the updated photolysis rate, and finally, it creates a ConnectorSystem object that sets the j_NO2 variable from the Photolysis component equal to the jNO2 variable from the Chemistry component. The next effect is that the photolysis rate in the Chemistry component is now controlled by the Photolysis component.","category":"page"},{"location":"composition/","page":"Composition","title":"Composition","text":"function EarthSciMLBase.couple2(c::ChemistryCoupler, p::PhotolysisCoupler)\n c, p = c.sys, p.sys\n c = param_to_var(convert(ODESystem, c), :jNO2)\n ConnectorSystem([c.jNO2 ~ p.j_NO2], c, p)\nend\nnothing # hide","category":"page"},{"location":"composition/","page":"Composition","title":"Composition","text":"Next, we define a method for coupling the Chemistry and Emissions components. To do this we use the operator_compose function to add the NO2 species from the Emissions component to the time derivative of NO2 in the Chemistry component.","category":"page"},{"location":"composition/","page":"Composition","title":"Composition","text":"function EarthSciMLBase.couple2(c::ChemistryCoupler, emis::EmissionsCoupler)\n c, emis = c.sys, emis.sys\n operator_compose(convert(ODESystem, c), emis, Dict(\n c.NO2 => emis.NO2,\n ))\nend\nnothing # hide","category":"page"},{"location":"composition/","page":"Composition","title":"Composition","text":"Finally, we can compose the model components together to create our complete model. To do so, we just initialize each model component and add the components together using the couple function. We can use the convert function to convert the composed model to a ModelingToolkit ODESystem so we can see what the combined equations look like.","category":"page"},{"location":"composition/","page":"Composition","title":"Composition","text":"model = couple(Photolysis(), Chemistry(), Emissions())\nconvert(ODESystem, model)","category":"page"},{"location":"composition/","page":"Composition","title":"Composition","text":"Finally, we can use the graph function to visualize the model components and their connections.","category":"page"},{"location":"composition/","page":"Composition","title":"Composition","text":"using MetaGraphsNext\nusing CairoMakie, GraphMakie\n\ng = graph(model)\n\nf, ax, p = graphplot(g; ilabels=labels(g))\nhidedecorations!(ax); hidespines!(ax); ax.aspect = DataAspect()\n\nf","category":"page"},{"location":"","page":"Home","title":"Home","text":"CurrentModule = EarthSciMLBase","category":"page"},{"location":"#EarthSciMLBase:-Utilities-for-Symbolic-Earth-Science-Modeling-and-Machine-Learning","page":"Home","title":"EarthSciMLBase: Utilities for Symbolic Earth Science Modeling and Machine Learning","text":"","category":"section"},{"location":"","page":"Home","title":"Home","text":"Documentation for EarthSciMLBase.","category":"page"},{"location":"","page":"Home","title":"Home","text":"This package contains utilities for constructing Earth Science models in Julia using ModelingToolkit.jl.","category":"page"},{"location":"#Installation","page":"Home","title":"Installation","text":"","category":"section"},{"location":"","page":"Home","title":"Home","text":"using Pkg\nPkg.add(\"EarthSciMLBase\")","category":"page"},{"location":"#Feature-Summary","page":"Home","title":"Feature Summary","text":"","category":"section"},{"location":"","page":"Home","title":"Home","text":"This package contains types and functions designed to simplify the process of constructing and composing symbolically-defined Earth Science model components together.","category":"page"},{"location":"#Feature-List","page":"Home","title":"Feature List","text":"","category":"section"},{"location":"","page":"Home","title":"Home","text":"Operations to compose ModelingToolkit.jl equation systems together.\nOperations to add intitial and boundary conditions to systems and to turn ODE systems into PDE systems, and to provide coordinate transformations.\nOperations to add Advection terms to systems.\nA Simulator type for running large-scale simulations.","category":"page"},{"location":"#Contributing","page":"Home","title":"Contributing","text":"","category":"section"},{"location":"","page":"Home","title":"Home","text":"Please refer to the SciML ColPrac: Contributor's Guide on Collaborative Practices for Community Packages for guidance on PRs, issues, and other matters relating to contributing.","category":"page"},{"location":"#Reproducibility","page":"Home","title":"Reproducibility","text":"","category":"section"},{"location":"","page":"Home","title":"Home","text":"
The documentation of this EarthSciML package was built using these direct dependencies,","category":"page"},{"location":"","page":"Home","title":"Home","text":"using Pkg # hide\nPkg.status() # hide","category":"page"},{"location":"","page":"Home","title":"Home","text":"
","category":"page"},{"location":"","page":"Home","title":"Home","text":"
and using this machine and Julia version.","category":"page"},{"location":"","page":"Home","title":"Home","text":"using InteractiveUtils # hide\nversioninfo() # hide","category":"page"},{"location":"","page":"Home","title":"Home","text":"
","category":"page"},{"location":"","page":"Home","title":"Home","text":"
A more complete overview of all dependencies and their versions is also provided.","category":"page"},{"location":"","page":"Home","title":"Home","text":"using Pkg # hide\nPkg.status(;mode = PKGMODE_MANIFEST) # hide","category":"page"},{"location":"","page":"Home","title":"Home","text":"
","category":"page"},{"location":"","page":"Home","title":"Home","text":"You can also download the \nmanifest file and the\nproject file.","category":"page"},{"location":"advection/#Advection","page":"Advection","title":"Advection","text":"","category":"section"},{"location":"advection/","page":"Advection","title":"Advection","text":"The Advection function adds advection to a system of equations. This is useful for modeling the transport of a substance by a fluid. Advection is implemented with the Advection type.","category":"page"},{"location":"advection/","page":"Advection","title":"Advection","text":"!warning Fully symbolic partial differential equations like those shown here don't currently work on domains that have a large number of grid cells. See here for additional information.","category":"page"},{"location":"advection/","page":"Advection","title":"Advection","text":"To demonstrate how this can work, we will start with a simple system of equations:","category":"page"},{"location":"advection/","page":"Advection","title":"Advection","text":"using EarthSciMLBase, ModelingToolkit\nusing ModelingToolkit: t_nounits, D_nounits\nt = t_nounits\nD = D_nounits\n\nfunction ExampleSys()\n @variables y(t)\n @parameters p=2.0\n ODESystem([D(y) ~ p], t; name=:ExampleSys)\nend\n\nExampleSys()","category":"page"},{"location":"advection/","page":"Advection","title":"Advection","text":"We also need to create our initial and boundary conditions.","category":"page"},{"location":"advection/","page":"Advection","title":"Advection","text":"using DomainSets\n@parameters x\ndomain = DomainInfo(constIC(0.0, t ∈ Interval(0, 1.0)), constBC(1.0, x ∈ Interval(0, 1.0)))\nnothing # hide","category":"page"},{"location":"advection/","page":"Advection","title":"Advection","text":"Now we convert add advection to each of the state variables. We're also adding a constant wind (ConstantWind) in the x-direction, with a speed of 1.0.","category":"page"},{"location":"advection/","page":"Advection","title":"Advection","text":"sys_advection = couple(ExampleSys(), domain, ConstantWind(t, 1.0), Advection())\nsys_mtk = convert(PDESystem, sys_advection)","category":"page"},{"location":"advection/","page":"Advection","title":"Advection","text":"Finally, we can discretize the system and solve it:","category":"page"},{"location":"advection/","page":"Advection","title":"Advection","text":"using MethodOfLines, DifferentialEquations, Plots\ndiscretization = MOLFiniteDifference([x=>10], t, approx_order=2)\n@time prob = discretize(sys_mtk, discretization)\n@time sol = solve(prob, Tsit5(), saveat=0.1)\n\n\n# Plot the solution.\ndiscrete_x = sol[x]\ndiscrete_t = sol[t]\nsoly = sol[sys_mtk.dvs[3]]\nanim = @animate for k in 1:length(discrete_t)\n plot(discrete_x, soly[k, 1:end], title=\"t=\\$(discrete_t[k])\", ylim=(0,2.5), lab=:none)\nend\ngif(anim, fps = 8)","category":"page"}] }