diff --git a/Project.toml b/Project.toml index e5a0e6be..0fc0c5f5 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "AIBECS" uuid = "ace601d6-714c-11e9-04e5-89b7fad23838" authors = ["Benoit Pasquier "] -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" diff --git a/ReadMe.md b/ReadMe.md index 422827a0..e18a4feb 100644 --- a/ReadMe.md +++ b/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! diff --git a/docs/Project.toml b/docs/Project.toml index 048cd950..6d811c4a 100644 --- a/docs/Project.toml +++ b/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" diff --git a/docs/lit/tutorials/5_river_discharge.jl b/docs/lit/tutorials/5_river_discharge.jl index a8605e8e..d83e245f 100644 --- a/docs/lit/tutorials/5_river_discharge.jl +++ b/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! diff --git a/src/AIBECS.jl b/src/AIBECS.jl index 6e9ef5ca..8f2c9017 100644 --- a/src/AIBECS.jl +++ b/src/AIBECS.jl @@ -25,6 +25,7 @@ using NearestNeighbors @reexport using MetadataArrays using Shapefile using Bijectors +using ImageFiltering diff --git a/src/Rivers.jl b/src/Rivers.jl index 56de7aa7..fbac7f71 100644 --- a/src/Rivers.jl +++ b/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)) diff --git a/src/diagnostics.jl b/src/diagnostics.jl index ba929597..1b28d3eb 100644 --- a/src/diagnostics.jl +++ b/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