Fix convoluted/messy output from F_and_∇ₓF and load (#85)

* JOSS submission WIP

* Move from F_and_∇ₓF to AIBECSFunction API

* Remove default He output from OCIM2.load

.gitignore

@@ -1,6 +1,7 @@
 Manifest.toml
 docs/build/
 *html
+*pdf
 
 docs/src/tutorials/
 docs/src/howtos/

Project.toml

@@ -1,7 +1,7 @@
 name = "AIBECS"
 uuid = "ace601d6-714c-11e9-04e5-89b7fad23838"
 authors = ["Benoit Pasquier <[email protected]>"]
-version = "0.10.12"
+version = "0.11.0"
 
 [deps]
 Bijectors = "76274a88-744f-5084-9051-94815aaf08c4"

docs/lit/tutorials/1_ideal_age.jl

@@ -82,13 +82,11 @@ p = IdealAgeParameters(1.0, 30.0)
 
 # We now use the AIBECS to generate the state function $\boldsymbol{F}$ (and its Jacobian) via
 
-F, ∇ₓF = F_and_∇ₓF(TOCIM2, G)
+F = AIBECSFunction(TOCIM2, G)
 
-# (`∇ₓF` is the **Jacobian** of the state function $\nabla_{\boldsymbol{x}}\boldsymbol{F}$, calculated automatically using dual numbers.)
-
-# Now that `F(x,p)`, and `p` are defined, we are going to solve for the steady-state.
-# But first, we must create a `SteadyStateProblem` object that contains `F`, `∇ₓF`, `p`, and an initial guess `x_init` for the age.
-# (`SteadyStateProblem` is specialized from [DiffEqBase](https://github.com/JuliaDiffEq/DiffEqBase.jl) for AIBECS models.)
+# Now that `F` and `p` are defined, we are going to solve for the steady-state.
+# But first, we must create a `SteadyStateProblem` object that contains `F`, `p`, and an initial guess `x_init` for the age.
+# (`SteadyStateProblem` comes from [DiffEqBase](https://github.com/JuliaDiffEq/DiffEqBase.jl).)
 
 # Let's make a vector of 0's for our initial guess.
 
@@ -97,7 +95,7 @@ x_init = zeros(nb)    # Start with age = 0 everywhere
 
 # Now we can create our `SteadyStateProblem` instance
 
-prob = SteadyStateProblem(F, ∇ₓF, x_init, p)
+prob = SteadyStateProblem(F, x_init, p)
 
 # And finally, we can `solve` this problem, using the AIBECS `CTKAlg()` algorithm,
 

docs/lit/tutorials/2_radiocarbon.jl

@@ -89,11 +89,11 @@ p = RadiocarbonParameters(λ = 50u"m"/10u"yr",
 #nb # > The parameters are converted to SI units when unpacked.
 #nb # >  When you specify units for your parameters, you must supply their values in that unit.
 
-# We generate the state function and its Jacobian, generate the corresponding steady-state problem, and solve it, via
+# We build the state function `F` and the corresponding steady-state problem (and solve it) via
 
-F, ∇ₓF = F_and_∇ₓF(T_OCCA, G)
+F = AIBECSFunction(T_OCCA, G)
 x = zeros(length(z)) # an initial guess
-prob = SteadyStateProblem(F, ∇ₓF, x, p)
+prob = SteadyStateProblem(F, x, p)
 R = solve(prob, CTKAlg()).u
 
 # This should take a few seconds on a laptop.

docs/lit/tutorials/3_Pmodel.jl

@@ -33,8 +33,6 @@
 #
 # $$R(x_\mathsf{POP}) = \frac{x_\mathsf{POP}}{\tau_\mathsf{POP}}.$$
 
-
-
 # We start by telling Julia we want to use the AIBECS and the OCIM0.1 circulation for DIP.
 
 using AIBECS
@@ -116,16 +114,16 @@ end
 
 p = PmodelParameters()
 
-# We generate the state function `F` and its Jacobian `∇ₓF`,
+# We generate the state function `F`,
 
 nb = sum(iswet(grd))
-F, ∇ₓF = F_and_∇ₓF((T_DIP, T_POP), (G_DIP, G_POP), nb)
+F = AIBECSFunction((T_DIP, T_POP), (G_DIP, G_POP), nb)
 
 # generate the steady-state problem,
 
 @unpack DIP_geo = p
 x = DIP_geo * ones(2nb) # initial guess
-prob = SteadyStateProblem(F, ∇ₓF, x, p)
+prob = SteadyStateProblem(F, x, p)
 
 # and solve it
 
@@ -165,8 +163,8 @@ T_POP2(p) = transportoperator(grd, z -> w(z,p); frac_seafloor=f_topo)
 
 # With this new vertical transport for POP, we can recreate our problem, solve it again
 
-F2, ∇ₓF2 = F_and_∇ₓF((T_DIP, T_POP2), (G_DIP, G_POP), nb)
-prob2 = SteadyStateProblem(F2, ∇ₓF2, x, p)
+F2 = AIBECSFunction((T_DIP, T_POP2), (G_DIP, G_POP), nb)
+prob2 = SteadyStateProblem(F2, x, p)
 sol2 = solve(prob2, CTKAlg()).u
 DIP2, POP2 = state_to_tracers(sol2, grd) # unpack tracers
 

docs/lit/tutorials/4_dustmodel.jl

@@ -104,15 +104,15 @@ end
 
 p = DustModelParameters()
 
-# We generate the state function `F` and its Jacobian `∇ₓF`,
+# We build the state function `F`,
 
 nb = count(iswet(grd))
-F, ∇ₓF = F_and_∇ₓF(T_dust, G_dust, nb)
+F = AIBECSFunction(T_dust, G_dust, nb)
 
-# generate the steady-state problem,
+# the steady-state problem,
 
 x = ones(nb) # initial guess
-prob = SteadyStateProblem(F, ∇ₓF, x, p)
+prob = SteadyStateProblem(F, x, p)
 
 # and solve it
 
@@ -149,7 +149,7 @@ profile_plot = plothorizontalmean(sol * u"g/m^3" .|> u"mg/m^3", grd)
 # For example, impose a larger fraction that stays in the bottom or change other parameters via
 
 p2 = DustModelParameters(w₀=0.1u"km/yr", w′=0.1u"km/yr/km", fsedremin=95.0u"percent")
-prob2 = SteadyStateProblem(F, ∇ₓF, x, p2)
+prob2 = SteadyStateProblem(F, x, p2)
 sol2 = solve(prob2, CTKAlg()).u
 plotverticalmean(sol2 * u"g/m^3", grd, color=cgrad(:turbid, rev=true, scale=:exp))
 

docs/lit/tutorials/5_river_discharge.jl

@@ -93,42 +93,42 @@ end
 
 p = RadioRiversParameters()
 
-# We generate the state function `F` and its Jacobian `∇ₓF`,
+# We generate the state function `F` that gives the tendencies,
 
-F, ∇ₓF = F_and_∇ₓF(T_OCIM2, G_radiorivers)
+F = AIBECSFunction(T_OCIM2, G_radiorivers)
 
 # generate the steady-state problem `prob`,
 
 nb = sum(iswet(grd))
 x = ones(nb) # initial guess
-prob = SteadyStateProblem(F, ∇ₓF, x, p)
+prob = SteadyStateProblem(F, x, p)
 
 # and solve it
 
-s = solve(prob, CTKAlg()).u * u"mol/m^3"
+sol = solve(prob, CTKAlg()).u * u"mol/m^3"
 
 # Let's now run some visualizations using the plot recipes.
 # Taking a horizontal slice of the 3D field at 200m gives
 
 cmap = :viridis
-plothorizontalslice(s, grd, zunit=u"μmol/m^3", depth=200, color=cmap)
+plothorizontalslice(sol, grd, zunit=u"μmol/m^3", depth=200, color=cmap)
 
 # and at 500m:
 
-plothorizontalslice(s, grd, zunit=u"μmol/m^3", depth=500, color=cmap, clim=(0,0.05))
+plothorizontalslice(sol, grd, zunit=u"μmol/m^3", depth=500, color=cmap, clim=(0,0.05))
 
 # Or we can change the timescale and watch the tracer fill the oceans:
 
 p = RadioRiversParameters(τ = 50.0u"yr")
-prob = SteadyStateProblem(F, ∇ₓF, x, p)
-s_τ50 = solve(prob, CTKAlg()).u * u"mol/m^3"
-plothorizontalslice(s_τ50, grd, zunit=u"μmol/m^3", depth=500, color=cmap, clim=(0,1))
+prob = SteadyStateProblem(F, x, p)
+sol_τ50 = solve(prob, CTKAlg()).u * u"mol/m^3"
+plothorizontalslice(sol_τ50, grd, zunit=u"μmol/m^3", depth=500, color=cmap, clim=(0,1))
 
 # Point-wise sources like the rivers used here can be problematic numerically because these
 # can generate large gradients and numerical noise with the given spatial discretization:
 
 cmap = :RdYlBu_4
-plothorizontalslice(s_τ50, grd, zunit=u"μmol/m^3", depth=0, color=cmap, clim=(-10,10))
+plothorizontalslice(sol_τ50, grd, zunit=u"μmol/m^3", depth=0, color=cmap, clim=(-10,10))
 
 # Note, for an example of numerical noise, the negative values appearing in red.
 # Hence, it might be wise to smooth the 3D field of the river sources.
@@ -139,16 +139,17 @@ plothorizontalslice(s_τ50, grd, zunit=u"μmol/m^3", depth=0, color=cmap, clim=(
 
 S = smooth_operator(grd, T_OCIM2)
 s_0 = S * (S * s_0) # smooth river sources twice
-s_smooth = solve(prob, CTKAlg()).u * u"mol/m^3"
-plothorizontalslice(s_smooth, grd, zunit=u"μmol/m^3", depth=0, color=cmap, clim=(-10,10))
+sol_smooth = solve(prob, CTKAlg()).u * u"mol/m^3"
+plothorizontalslice(sol_smooth, grd, zunit=u"μmol/m^3", depth=0, color=cmap, clim=(-10,10))
 
 # Note that such numerical noise generally dampens out with distance and depth
 
 cmap = :viridis
 plot(
-    plothorizontalslice(s_τ50, grd, zunit=u"μmol/m^3", depth=200, color=cmap, clim=(0,2)),
-    plothorizontalslice(s_smooth, grd, zunit=u"μmol/m^3", depth=200, color=cmap, clim=(0,2)),
-    layout=(2,1)
+    plothorizontalslice(sol_τ50, grd, zunit=u"μmol/m^3", depth=200, color=cmap, clim=(0,2)),
+    plothorizontalslice(sol_smooth, grd, zunit=u"μmol/m^3", depth=200, color=cmap, clim=(0,2)),
+    layout=(2,1),
+    size=(600,600)
 )
 
-# Nevertheless, it is always good to reduce noise as much as possible!
+# Nevertheless, it can be good to reduce noise as much as possible!

docs/lit/tutorials/6_groundwater_discharge.jl

@@ -91,33 +91,33 @@ end
 
 p = GroundWatersParameters()
 
-# We generate the state function `F` and its Jacobian `∇ₓF`,
+# We build the state function `F`,
 
-F, ∇ₓF = F_and_∇ₓF(T_OCIM2, G_gw)
+F = AIBECSFunction(T_OCIM2, G_gw)
 
-# generate the steady-state problem `prob`,
+# the steady-state problem `prob`,
 
 nb = sum(iswet(grd))
 x = ones(nb) # initial guess
-prob = SteadyStateProblem(F, ∇ₓF, x, p)
+prob = SteadyStateProblem(F, x, p)
 
 # and solve it
 
-s = solve(prob, CTKAlg()).u * u"mol/m^3"
+sol = solve(prob, CTKAlg()).u * u"mol/m^3"
 
 # Let's now run some visualizations using the plot recipes.
 # Taking a horizontal slice of the 3D field at 200m gives
 
 cmap = :viridis
-plothorizontalslice(s, grd, zunit=u"μmol/m^3", depth=200, color=cmap, clim=(0,100))
+plothorizontalslice(sol, grd, zunit=u"μmol/m^3", depth=200, color=cmap, clim=(0,100))
 
 # and at 500m:
 
-plothorizontalslice(s, grd, zunit=u"μmol/m^3", depth=500, color=cmap, clim=(0,25))
+plothorizontalslice(sol, grd, zunit=u"μmol/m^3", depth=500, color=cmap, clim=(0,25))
 
 # Or we can increase the decay timescale (×10) and decrease the groundwater concentration (÷10) to get a different (more well-mixed) tracer distribution:
 
 p = GroundWatersParameters(τ = 200.0u"yr", C_gw = 0.1u"mol/m^3")
-prob = SteadyStateProblem(F, ∇ₓF, x, p)
-s_τ50 = solve(prob, CTKAlg()).u * u"mol/m^3"
-plothorizontalslice(s_τ50, grd, zunit=u"μmol/m^3", depth=500, color=cmap, clim=(0,25))
+prob = SteadyStateProblem(F, x, p)
+sol_τ50 = solve(prob, CTKAlg()).u * u"mol/m^3"
+plothorizontalslice(sol_τ50, grd, zunit=u"μmol/m^3", depth=500, color=cmap, clim=(0,25))