Add carbonate chemistry examples

docs/src/PALEOsediment_Reactions.md

@@ -1,6 +1,6 @@
 # PALEOsediment Reactions
 
-## Sediment
+## Sediment transport
 ```@meta
 CurrentModule = PALEOsediment.Sediment
 ```

examples/Project.toml

@@ -1,6 +1,7 @@
 [deps]
 BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"
 DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
+Infiltrator = "5903a43b-9cc3-4c30-8d17-598619ec4e9b"
 LineSearches = "d3d80556-e9d4-5f37-9878-2ab0fcc64255"
 Logging = "56ddb016-857b-54e1-b83d-db4d58db5568"
 NLsolve = "2774e3e8-f4cf-5e23-947b-6d7e65073b56"
@@ -10,6 +11,7 @@ PALEOmodel = "bf7b4fbe-ccb1-42c5-83c2-e6e9378b660c"
 PALEOsediment = "e0a37952-6f01-4236-91ff-62fdc855f67b"
 Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
 Revise = "295af30f-e4ad-537b-8983-00126c2a3abe"
+SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
 Sundials = "c3572dad-4567-51f8-b174-8c6c989267f4"
 
 [extras]

examples/boudreau1996/PALEO_examples_sediment_carb.jl

@@ -0,0 +1,79 @@
+using Logging
+using Sundials
+import LineSearches
+using BenchmarkTools
+
+using Plots
+
+import PALEOboxes as PB
+import PALEOmodel
+import PALEOaqchem
+import PALEOsediment
+
+global_logger(ConsoleLogger(stderr,Logging.Info))
+
+include("config_sediment_expts.jl")
+include("../plot_sediment.jl")
+
+
+model = PB.create_model_from_config(
+    joinpath(@__DIR__, "PALEO_examples_sediment_cfg.yaml"), 
+    "sediment_Corg_O2_carb",
+)
+
+#############################
+# Steady state solutions
+############################
+
+config_sediment_expts(model, 
+    [
+        "baseline",
+    ]
+) 
+
+tspan=(0.0, 10000.0)
+
+initial_state, modeldata = PALEOmodel.initialize!(model)
+
+run = PALEOmodel.Run(model=model, output = PALEOmodel.OutputWriters.OutputMemory())
+
+
+# PTC, Newton, no line search
+# Bounds and max step size for Newton solve. NB some tracers start at zero so set newton_max_ratio=Inf
+newton_min, newton_max, newton_min_ratio, newton_max_ratio = 1e-80, Inf, 0.1, Inf 
+PALEOmodel.SteadyState.steadystate_ptcForwardDiff(
+    run, initial_state, modeldata, tspan, 1e-3,
+    deltat_fac=2.0,
+    solvekwargs=(
+        ftol=1e-7,
+        iterations=20,
+        method=:newton,
+        linesearch=LineSearches.Static(),
+        apply_step! = PALEOmodel.SolverFunctions.StepClampMultAll!(newton_min, newton_max, newton_min_ratio, newton_max_ratio),
+        # store_trace=true,
+        # show_trace=true, 
+    ),
+    verbose=false,
+) 
+
+
+############################
+# Plot 
+############################
+
+# single plots
+# plotlyjs(size=(750, 565))
+# pager = PALEOmodel.DefaultPlotPager()
+
+# multiple plots per screen
+gr(size=(900, 600))
+pager = PALEOmodel.PlotPager((1,3), (legend_background_color=nothing, ))
+
+plot_Corg_O2(run.output; Corgs=["Corg1", "Corg2"], colT=[first(tspan), last(tspan)], pager)
+plot_solutes(run.output; solutes = ["P", "SO4", "H2S", "CH4", "DIC", "TAlk"], colT=[first(tspan), last(tspan)], pager)
+plot_rates(run.output; colT=[first(tspan), last(tspan)], pager=pager)
+plot_carbchem(run.output; include_constraint_error=true, colT=last(tspan), pager)
+
+pager(:newpage)
+
+