Merge pull request #158 from ggebbie/157-getboundarycondition-with-1-…

…dim-field

1-dimensional boundary conditions

Project.toml

@@ -1,7 +1,7 @@
 name = "TMI"
 uuid = "582500f6-28c8-4d8f-aabe-b197735ec1d4"
 authors = ["G Jake Gebbie <ggebbie@whoi.edu>"]
-version = "0.2.12"
+version = "0.2.13"
 
 [deps]
 COSMO = "1e616198-aa4e-51ec-90a2-23f7fbd31d8d"

README.md

@@ -18,7 +18,7 @@ Started by G Jake Gebbie, WHOI, ggebbie@whoi.edu
 
 # Requirements
 
-The built-in tests are automatically checked with Julia 1.8. 
+The built-in tests are automatically checked with Julia 1.10. 
 
 # Setting up project environment
 
@@ -57,28 +57,15 @@ Scripts can be run non-interactively like this:\
 `cd TMI.jl`\
 `julia --project=scripts scripts/ex1.trackpathways.jl`
 
-To show graphical results, TMI.jl uses `GGplot.jl` for plotting routines. In particular, matplotlib, cartopy and cmocean packages are handled by `GGplot` so that they are not dependencies in `TMI.jl`. Thus the `scripts` directory has its own environment distinct from the TMI project. If you are working interactively, try the following commands to activate the scripts environment: 
+To show graphical results, TMI.jl uses `GeoPythonPlot.jl` for plotting routines. In particular, matplotlib, cartopy and cmocean packages are handled by `GeoPythonPlot` so that they are not dependencies in `TMI.jl`. Thus the `scripts` directory has its own environment distinct from the TMI project. If you are working interactively, try the following commands to activate the scripts environment: 
 
 `cd("scripts")`\
 `import Pkg; Pkg.activate(".")`
 
-If you use the examples and `GGplot`, you will need a python environment in Julia. Direct the python environment to an existing system-wide version of python with these already installed:
-`ENV["PYTHON"]="python/directory/on/your/machine"`
-
-GGplot will use a package-specific python environment built from scratch using the `CondaPkg.jl` package. Check out the `GGplot/deps/build.jl` file to see how this Python environment is set up. In particular, it executes:
+GeoPythonPlot will use a package-specific python environment built from scratch using the `CondaPkg.jl` package. Check out the `GeoPythonPlot/deps/build.jl` file to see how this Python environment is set up. In particular, it executes:
 `ENV["PYTHON"]=""` 
 
-Rather than using the `PyCall.jl` package, `GGplot.jl` uses `PythonCall.jl` in order to minimize errors that occur due to incompatible Python environments. 
-`import Pkg; Pkg.add("PythonCall")`\
-`import Pkg; Pkg.build("PythonCall") # probably not necessary`
-
-In particular, GGplot uses the matplotlib, cartopy, and cmocean packages. The channel is automatically assumed to be conda-forge. \
-`using CondaPkg`\
-`] conda add matplotlib # from the package manager`\
-`] Conda add cartopy`\
-`] Conda add cmocean`
-
-This version of cartopy (v0.20.0+) can download coastlines from Amazon Web Services.
+Rather than using the `PyCall.jl` package, `GeoPythonPlot.jl` uses `PythonCall.jl` in order to minimize errors that occur due to incompatible Python environments. 
 
 # Data files
 

src/boundary_condition.jl

@@ -30,8 +30,7 @@ struct BoundaryCondition{T <: Real,
     units::String
 end
 
-Base.propertynames(b::BoundaryCondition) = (:i,:j,:k,fieldnames(typeof(γ))...)
-
+#Base.propertynames(b::BoundaryCondition) = (:i,:j,:k,fieldnames(typeof(γ))...)
 # function Base.getproperty(b::BoundaryCondition,
 #     d::Symbol)
 #     if d === :i 
@@ -240,7 +239,7 @@ function ones(dim::I,dimval::I,γ::Grid,name::Symbol,longname::String,units::Str
 end
 
 """
-   Get boundary condition by extracting from 3D tracer
+   Get boundary condition by extracting from N-dimensional tracer and returning (N-1)-dimensional array
 """
 function getboundarycondition(field::Field{T,R,N},dim::Integer,dimval::Integer,γ::Grid)::BoundaryCondition where {T<:Real,R<:Real,N}
 
@@ -250,17 +249,17 @@ function getboundarycondition(field::Field{T,R,N},dim::Integer,dimval::Integer,
             ind = deleteat!(collect(1:N),n)
             # boundary axes
             baxes = γ.axes[ind]
-            wet2d = copy(selectdim(γ.wet,dim,dimval))
-            tracer2d = copy(selectdim(field.tracer,
+            wet_nminus1 = copy(selectdim(γ.wet,dim,dimval))
+            tracer_nminus1 = copy(selectdim(field.tracer,
                 dim,
                 dimval))
 
-            return BoundaryCondition(tracer2d,
+            return BoundaryCondition(tracer_nminus1,
                 baxes,
                 k,
                 dim,
                 dimval,
-                wet2d)
+                wet_nminus1)
         end
     end
 end
@@ -320,7 +319,7 @@ vec(u::BoundaryCondition) = u.tracer[u.wet]
 # Output
 - `d`: vector that describes surface patch
 """
-function surfacepatch(lonbox,latbox,γ::Grid)::BoundaryCondition
+function surfacepatch(lonbox,latbox,γ::Grid{R,3})::BoundaryCondition where R
 
     # ternary operator to handle longitudinal wraparound
     lonbox[1] ≤ 0 ? lonbox[1] += 360 : nothing

src/mass_fractions.jl

@@ -588,7 +588,7 @@ function local_watermass_matrix(c::NamedTuple,
 end
 
 """
-`function watermassmatrix(m::NamedTuple, γ::Grid)`
+`function watermassmatrix(m::Union{NamedTuple,Vector}, γ::Grid)`
 
 Produce water-mass matrix from mass fractions and grid.
 
@@ -599,7 +599,7 @@ Produce water-mass matrix from mass fractions and grid.
 # Output
 - `A`: sparse water-mass matrix
 """
-function watermassmatrix(m::NamedTuple, γ::Grid)
+function watermassmatrix(m::Union{NamedTuple,Vector}, γ::Grid)
 
     # get total size of mass fractions
 

test/test_1d.jl

@@ -1,7 +1,5 @@
-@testset "mass fractions 1d" begin
+@testset "1-dimensional domain" begin
 
-    #using Revise
-    #using TMI
     ngrid = (50) # number of grid cells
     xmax = 1000.0 # domain size 
     lon = collect(range(0.0,1000.0,length=ngrid[1]))
@@ -28,5 +26,36 @@
     w = (c = 0.01,)
 
     m = massfractions(y, w)
+    A = watermassmatrix(m, γ)
 
+    ## following ex0
+    # manually pick out boundary conditions
+    dim =1
+    b = (west = TMI.getboundarycondition(c, dim, 1, γ),
+        east = TMI.getboundarycondition(c, 1, ngrid[dim], γ))
+
+        # no real need to worry about LU decomposition in small problem.
+    c̃ = steadyinversion(A,b,γ)
+    @test abs(maximum(c-c̃)) < 1e-4
+
+    # Introduce a nonconservative variables with a source
+    qfield = 1.0e-2 * ones(ngrid)
+    # requires negative sign which is counterintutive (needs to be fixed)
+    q = TMI.Source(-qfield, γ, :q, "remineralized stuff", "μmol/kg", false)
+    c_noncons = steadyinversion(A,b,γ, q=q)
+    Δc = c_noncons - c
+    @test iszero(sum(Δc .< 0.0))
+
+    # ex1: trackpathways
+    # no `trackpathways` function for 1D has been implemented
+    # Here is a manual tracking of pathways from western boundary
+
+     # ones() and trues(): trick to get a 0-dim array
+    bwest = (west = BoundaryCondition(ones(),(),0.0,dim,1,trues()),
+        east = BoundaryCondition(zeros(),(),xmax,dim,ngrid[1],trues()))
+
+    gwest = steadyinversion(A,bwest,γ) # fraction from west boundary
+    @test maximum(gwest) ≤ 1.0
+    @test minimum(gwest) ≥ 0.0
+
 end