Merge branch 'master' into inplacejac

* master:
  SmoothOperator (#64)
  Update Project.toml
  CompatHelper: bump compat for "BSON" to "0.3"
  Fix ReadMe typos
  Update Project.toml
  CompatHelper: bump compat for "Shapefile" to "0.7"
  Update Project.toml
  bump version
  CompatHelper: bump compat for "Reexport" to "1.0"
  CompatHelper: bump compat for "NCDatasets" to "0.11"

Project.toml

@@ -1,7 +1,7 @@
 name = "AIBECS"
 uuid = "ace601d6-714c-11e9-04e5-89b7fad23838"
 authors = ["Benoit Pasquier <[email protected]>"]
-version = "0.8.2"
+version = "0.8.6"
 
 [deps]
 BSON = "fbb218c0-5317-5bc6-957e-2ee96dd4b1f0"
@@ -15,6 +15,7 @@ Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
 FieldMetadata = "bf96fef3-21d2-5d20-8afa-0e7d4c32a885"
 Flatten = "4c728ea3-d9ee-5c9a-9642-b6f7d7dc04fa"
 ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
+ImageFiltering = "6a3955dd-da59-5b1f-98d4-e7296123deb5"
 Interpolations = "a98d9a8b-a2ab-59e6-89dd-64a1c18fca59"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 MetadataArrays = "49441bc9-da82-574f-b07c-a0d10dd4ac13"
@@ -32,7 +33,7 @@ Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d"
 UnitfulRecipes = "42071c24-d89e-48dd-8a24-8a12d9b8861f"
 
 [compat]
-BSON = "0.2"
+BSON = "0.2, 0.3"
 Bijectors = "0.8"
 DataDeps = "0.7"
 DataFrames = "0.20, 0.21, 0.22"
@@ -42,14 +43,15 @@ Distributions = "0.23, 0.24"
 FieldMetadata = "0.3"
 Flatten = "0.3, 0.4"
 ForwardDiff = "0.10"
+ImageFiltering = "0.6"
 Interpolations = "0.12, 0.13"
 MetadataArrays = "0.1"
-NCDatasets = "0.10"
+NCDatasets = "0.10, 0.11"
 NearestNeighbors = "0.4"
 OceanGrids = "0.3"
 RecipesBase = "1"
-Reexport = "0.2"
-Shapefile = "0.6"
+Reexport = "0.2, 1.0"
+Shapefile = "0.6, 0.7"
 UnPack = "1"
 Unitful = "1"
 UnitfulRecipes = "0.2, 1.0"

ReadMe.md

@@ -42,9 +42,9 @@
 
 
 
-**AIBECS** (for **A**lgebraic **I**mplicit **B**iogeochemical **E**lemental **C**ycling **S**ystem, pronounced like the cool [ibex](https://en.wikipedia.org/wiki/Ibex)) is a Julia package that provides ocean biogeochemistry modelers with an easy-to-use interface for creating and running models of the ocean system.
+**AIBECS** (for **A**lgebraic **I**mplicit **B**iogeochemical **E**lemental **C**ycling **S**ystem, pronounced like the cool [ibex](https://en.wikipedia.org/wiki/Ibex)) is a Julia package that provides ocean biogeochemistry modellers with an easy-to-use interface for creating and running models of the ocean system.
 
-AIBECS is a system because it allows you to chose some biogeochemical tracers, define their interactions, select an ocean circulation and *Voilà!* — your model is ready to run.
+AIBECS is a system because it allows you to choose some biogeochemical tracers, define their interactions, select an ocean circulation and *Voilà!* — your model is ready to run.
 
 ## Getting started
 
@@ -63,7 +63,7 @@ For a single tracer, *x* can be interpreted as the 3D field of its concentration
 In AIBECS, *x* is represented as a column vector.
 
 The operator 𝓣 is a spatial differential operator that represents the transport of tracers.
-For example, for a single tracer transported by the ocean circulation,
+For example, for a single tracer transported by ocean circulation,
 
 𝓣 = ∇ ⋅ (***u*** + **K**∇)
 
@@ -81,8 +81,10 @@ That's pretty much the whole concept!
 ## References
 
 If you use this package, please cite it.
-And if you use data with these package (like the ocean circulation from the OCIM) please also cite them.
-The references under bibtex format are available in the [CITATION.bib](./CITATION.bib) file.
+
+If you use data provided by this package (like the ocean circulation from the OCIM), please cite them as well.
+
+For convenience, all the references are available in [BibTeX](https://en.wikipedia.org/wiki/BibTeX) format in the [CITATION.bib](./CITATION.bib) file.
 
 Also, if you want to do research using the AIBECS, and you think I could help, do not hesitate to contact me directly (contacts on my [website](www.bpasquier.com)), I would be happy to contribute and collaborate!
 

docs/Project.toml

@@ -7,4 +7,4 @@ Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
 Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d"
 
 [compat]
-Documenter = "0.25"
+Documenter = "0.26"

docs/lit/tutorials/5_river_discharge.jl

@@ -123,3 +123,32 @@ 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))
+
+# 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))
+
+# 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.
+# We can do this using the `smooth_operator` function
+# (`smooth_operator` creates an operator that can be applied to smooth the AIBECS vector
+# along the stencil of the given transport matrix and in a volume-conservative way).
+# And we can check that the negative values disappear after smoothing of the sources:
+
+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))
+
+# 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)
+)
+
+# Nevertheless, it is always good to reduce noise as much as possible!

src/AIBECS.jl

@@ -25,6 +25,7 @@ using NearestNeighbors
 @reexport using MetadataArrays
 using Shapefile
 using Bijectors
+using ImageFiltering
 
 
 

src/Rivers.jl

@@ -241,7 +241,7 @@ citation() = """
 - For the formal publication that describe the dataset: Dai, A., and K. E. Trenberth, 2002: Estimates of freshwater discharge from continents: Latitudinal and seasonal variations. J. Hydrometeorol., 3, 660–687.
 - For descriptions of updates to the dataset:
     - Dai, A., T. Qian, K. E. Trenberth, and J. D Milliman, 2009: Changes in continental freshwater discharge from 1948-2004. J. Climate, 22, 2773–2791.
-    - Dai, A., 2016: Historical and future changes in streamflow and continental runoff: A review, in Terrestrial Water Cycle and Climate Change: Natural and Human-Induced Impacts, Geophysical Monograph 221. Ed. Qiuhong Tang and Taikan Oki, AGU, John Wiley & Sons, 17–37 (DOI: 10.1002/9781118971772).
+    - Dai, A., 2016: Historical and future changes in streamflow and continental runoff: A review, in Terrestrial Water Cycle and Climate Change: Natural and Human-Induced Impacts, Geophysical Monograph 221. Ed. Qiuhong Tang and Taikan Oki, AGU, John Wiley & Sons, 17–37 (DOI: 10.1002/9781118971772).
 """
 
 uconvert(u, r::River) = River(r.name, r.lon, r.lat, uconvert(u, r.VFR))

src/diagnostics.jl

@@ -85,7 +85,25 @@ export directional_transport, directional_transports
 
 
 
+"""
+    smooth_operator(grd, T; σs=(1.0, 1.0, 0.25))
 
+return a matrix of the same size and sparsity as `T` that smoothes data using
+a Gaussian Kernel for values, but conserving mass.
+"""
+function smooth_operator(grd, T; σs=(1.0, 1.0, 0.25))
+    st = stencil(grd, T)
+    weight(dir) = Kernel.gaussian(σs)[dir]
+    Isym, Jsym, _, _ = symmetric_IJVVᵀ(T)
+    dirs = directions(Isym, Jsym, grd)
+    A = sparse(Isym, Jsym, weight.(dirs))
+    # We add mass conservation to the Gaussian smoothing matrix
+    # by modifying A column by column such that A conserves mass,
+    # i.e., such that vᵀ A x = vᵀ x
+    v = depthvec(grd)
+    A * sparse(Diagonal(v ./ (A' * v)))
+end
+export smooth_operator