diff --git a/.gitignore b/.gitignore index a14808df..88bc7020 100644 --- a/.gitignore +++ b/.gitignore @@ -4,3 +4,4 @@ *.jl.mem /docs/build/ *~ +*nc \ No newline at end of file diff --git a/Manifest.toml b/Manifest.toml index a6c0344a..2197c88c 100644 --- a/Manifest.toml +++ b/Manifest.toml @@ -2,7 +2,7 @@ julia_version = "1.8.5" manifest_format = "2.0" -project_hash = "df12c94795f35fec1a6c70b73daa8a083a2a1473" +project_hash = "c02643bd8035ebad2cc399fe86fb6c70dd21fbf9" [[deps.ANSIColoredPrinters]] git-tree-sha1 = "574baf8110975760d391c710b6341da1afa48d8c" @@ -77,12 +77,42 @@ git-tree-sha1 = "844b061c104c408b24537482469400af6075aae4" uuid = "9e997f8a-9a97-42d5-a9f1-ce6bfc15e2c0" version = "0.1.5" +[[deps.CodeTracking]] +deps = ["InteractiveUtils", "UUIDs"] +git-tree-sha1 = "d57c99cc7e637165c81b30eb268eabe156a45c49" +uuid = "da1fd8a2-8d9e-5ec2-8556-3022fb5608a2" +version = "1.2.2" + [[deps.CodecZlib]] deps = ["TranscodingStreams", "Zlib_jll"] git-tree-sha1 = "9c209fb7536406834aa938fb149964b985de6c83" uuid = "944b1d66-785c-5afd-91f1-9de20f533193" version = "0.7.1" +[[deps.ColorSchemes]] +deps = ["ColorTypes", "ColorVectorSpace", "Colors", "FixedPointNumbers", "Random", "SnoopPrecompile"] +git-tree-sha1 = "aa3edc8f8dea6cbfa176ee12f7c2fc82f0608ed3" +uuid = "35d6a980-a343-548e-a6ea-1d62b119f2f4" +version = "3.20.0" + +[[deps.ColorTypes]] +deps = ["FixedPointNumbers", "Random"] +git-tree-sha1 = "eb7f0f8307f71fac7c606984ea5fb2817275d6e4" +uuid = "3da002f7-5984-5a60-b8a6-cbb66c0b333f" +version = "0.11.4" + +[[deps.ColorVectorSpace]] +deps = ["ColorTypes", "FixedPointNumbers", "LinearAlgebra", "SpecialFunctions", "Statistics", "TensorCore"] +git-tree-sha1 = "600cc5508d66b78aae350f7accdb58763ac18589" +uuid = "c3611d14-8923-5661-9e6a-0046d554d3a4" +version = "0.9.10" + +[[deps.Colors]] +deps = ["ColorTypes", "FixedPointNumbers", "Reexport"] +git-tree-sha1 = "fc08e5930ee9a4e03f84bfb5211cb54e7769758a" +uuid = "5ae59095-9a9b-59fe-a467-6f913c188581" +version = "0.12.10" + [[deps.CommonSubexpressions]] deps = ["MacroTools", "Test"] git-tree-sha1 = "7b8a93dba8af7e3b42fecabf646260105ac373f7" @@ -112,6 +142,11 @@ git-tree-sha1 = "fb21ddd70a051d882a1686a5a550990bbe371a95" uuid = "187b0558-2788-49d3-abe0-74a17ed4e7c9" version = "1.4.1" +[[deps.Contour]] +git-tree-sha1 = "d05d9e7b7aedff4e5b51a029dced05cfb6125781" +uuid = "d38c429a-6771-53c6-b99e-75d170b6e991" +version = "0.6.2" + [[deps.Crayons]] git-tree-sha1 = "249fe38abf76d48563e2f4556bebd215aa317e15" uuid = "a8cc5b0e-0ffa-5ad4-8c14-923d3ee1735f" @@ -237,6 +272,12 @@ git-tree-sha1 = "04ed1f0029b6b3af88343e439b995141cb0d0b8d" uuid = "6a86dc24-6348-571c-b903-95158fe2bd41" version = "2.17.0" +[[deps.FixedPointNumbers]] +deps = ["Statistics"] +git-tree-sha1 = "335bfdceacc84c5cdf16aadc768aa5ddfc5383cc" +uuid = "53c48c17-4a7d-5ca2-90c5-79b7896eea93" +version = "0.8.4" + [[deps.Formatting]] deps = ["Printf"] git-tree-sha1 = "8339d61043228fdd3eb658d86c926cb282ae72a8" @@ -357,6 +398,12 @@ git-tree-sha1 = "3c837543ddb02250ef42f4738347454f95079d4e" uuid = "682c06a0-de6a-54ab-a142-c8b1cf79cde6" version = "0.21.3" +[[deps.JuliaInterpreter]] +deps = ["CodeTracking", "InteractiveUtils", "Random", "UUIDs"] +git-tree-sha1 = "d9ae7a9081d9b1a3b2a5c1d3dac5e2fdaafbd538" +uuid = "aa1ae85d-cabe-5617-a682-6adf51b2e16a" +version = "0.9.22" + [[deps.LaTeXStrings]] git-tree-sha1 = "f2355693d6778a178ade15952b7ac47a4ff97996" uuid = "b964fa9f-0449-5b57-a5c2-d3ea65f4040f" @@ -415,6 +462,12 @@ git-tree-sha1 = "cedb76b37bc5a6c702ade66be44f831fa23c681e" uuid = "e6f89c97-d47a-5376-807f-9c37f3926c36" version = "1.0.0" +[[deps.LoweredCodeUtils]] +deps = ["JuliaInterpreter"] +git-tree-sha1 = "60168780555f3e663c536500aa790b6368adc02a" +uuid = "6f1432cf-f94c-5a45-995e-cdbf5db27b0b" +version = "2.3.0" + [[deps.MAT]] deps = ["BufferedStreams", "CodecZlib", "HDF5", "SparseArrays"] git-tree-sha1 = "971be550166fe3f604d28715302b58a3f7293160" @@ -427,6 +480,12 @@ git-tree-sha1 = "42324d08725e200c23d4dfb549e0d5d89dede2d2" uuid = "1914dd2f-81c6-5fcd-8719-6d5c9610ff09" version = "0.5.10" +[[deps.MarchingCubes]] +deps = ["SnoopPrecompile", "StaticArrays"] +git-tree-sha1 = "b198463d1a631e8771709bc8e011ba329da9ad38" +uuid = "299715c1-40a9-479a-aaf9-4a633d36f717" +version = "0.1.7" + [[deps.Markdown]] deps = ["Base64"] uuid = "d6f4376e-aef5-505a-96c1-9c027394607a" @@ -616,6 +675,12 @@ git-tree-sha1 = "838a3a4188e2ded87a4f9f184b4b0d78a1e91cb7" uuid = "ae029012-a4dd-5104-9daa-d747884805df" version = "1.3.0" +[[deps.Revise]] +deps = ["CodeTracking", "Distributed", "FileWatching", "JuliaInterpreter", "LibGit2", "LoweredCodeUtils", "OrderedCollections", "Pkg", "REPL", "Requires", "UUIDs", "Unicode"] +git-tree-sha1 = "90cb983381a9dc7d3dff5fb2d1ee52cd59877412" +uuid = "295af30f-e4ad-537b-8983-00126c2a3abe" +version = "3.5.1" + [[deps.Rmath]] deps = ["Random", "Rmath_jll"] git-tree-sha1 = "f65dcb5fa46aee0cf9ed6274ccbd597adc49aa7b" @@ -740,6 +805,12 @@ deps = ["ArgTools", "SHA"] uuid = "a4e569a6-e804-4fa4-b0f3-eef7a1d5b13e" version = "1.10.1" +[[deps.TensorCore]] +deps = ["LinearAlgebra"] +git-tree-sha1 = "1feb45f88d133a655e001435632f019a9a1bcdb6" +uuid = "62fd8b95-f654-4bbd-a8a5-9c27f68ccd50" +version = "0.1.1" + [[deps.Test]] deps = ["InteractiveUtils", "Logging", "Random", "Serialization"] uuid = "8dfed614-e22c-5e08-85e1-65c5234f0b40" @@ -767,6 +838,12 @@ version = "1.0.2" [[deps.Unicode]] uuid = "4ec0a83e-493e-50e2-b9ac-8f72acf5a8f5" +[[deps.UnicodePlots]] +deps = ["ColorSchemes", "ColorTypes", "Contour", "Crayons", "Dates", "LinearAlgebra", "MarchingCubes", "NaNMath", "Printf", "Requires", "SnoopPrecompile", "SparseArrays", "StaticArrays", "StatsBase"] +git-tree-sha1 = "516462456b9e54f2b94351f51186b3b8aedfb7fd" +uuid = "b8865327-cd53-5732-bb35-84acbb429228" +version = "3.4.2" + [[deps.VersionParsing]] git-tree-sha1 = "58d6e80b4ee071f5efd07fda82cb9fbe17200868" uuid = "81def892-9a0e-5fdd-b105-ffc91e053289" diff --git a/Project.toml b/Project.toml index 7b4dc12b..dc05664d 100644 --- a/Project.toml +++ b/Project.toml @@ -23,9 +23,11 @@ NetCDF = "30363a11-5582-574a-97bb-aa9a979735b9" NetCDF_jll = "7243133f-43d8-5620-bbf4-c2c921802cf3" Optim = "429524aa-4258-5aef-a3af-852621145aeb" Pkg = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f" +Revise = "295af30f-e4ad-537b-8983-00126c2a3abe" SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2" StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91" +UnicodePlots = "b8865327-cd53-5732-bb35-84acbb429228" [compat] CSV = "0.9, 0.10" diff --git a/scripts/mat2netcdf.jl b/scripts/mat2netcdf.jl index 78685196..d6a2e678 100644 --- a/scripts/mat2netcdf.jl +++ b/scripts/mat2netcdf.jl @@ -12,6 +12,7 @@ TMIversion = "LGM_90x45x33_G14A" TMIversion = "LGM_90x45x33_GPLS1" TMIversion = "LGM_90x45x33_GPLS2" TMIversion = "LGM_90x45x33_OG18" +TMIversion = "nordic_201x115x46_B23" # original NetCDF version #A, Alu, γ, TMIfile = config(TMIversion) diff --git a/src/TMI.jl b/src/TMI.jl index 931dd741..cb6c0b95 100644 --- a/src/TMI.jl +++ b/src/TMI.jl @@ -14,6 +14,7 @@ using LineSearches using MAT using NCDatasets using DataStructures +using UnicodePlots export config, config_from_mat, config_from_nc, download_ncfile, download_matfile, @@ -29,6 +30,7 @@ export config, config_from_mat, config_from_nc, linearindex, nearestneighbor, updatelinearindex, nearestneighbormask, horizontaldistance, readtracer, readfield, writefield, + readsource, cartesianindex, Γ, costfunction_gridded_obs, costfunction_gridded_obs!, costfunction_point_obs, costfunction_point_obs!, @@ -45,10 +47,15 @@ export config, config_from_mat, config_from_nc, surface_oxygensaturation, oxygen, location_obs, getsurfaceboundary, zerosurfaceboundary, onesurfaceboundary, setboundarycondition!, - adjustboundarycondition, adjustboundarycondition!, + adjustboundarycondition!, gsetboundarycondition, setsource!, - zeros, one, oneunit, ones, maximum, minimum, (+), (-), (*), dot, - Grid, Field, BoundaryCondition, vec, unvec!, unvec + zeros, one, oneunit, ones, + maximum, minimum, (+), (-), (*), dot, + zerosource, onesource, + Grid, Field, BoundaryCondition, vec, unvec!, unvec, + adjustsource, adjustsource!, + zerowestboundary, zeronorthboundary, zeroeastboundary, zerosouthboundary, +onewestboundary, onenorthboundary, oneeastboundary, onesouthboundary import Base: zeros, one, oneunit, ones, maximum, minimum, (\) import Base: (+), (-), (*), (/), vec @@ -63,10 +70,28 @@ struct Grid lon::Vector{Float64} lat::Vector{Float64} depth::Vector{Float64} - I::Vector{CartesianIndex{3}} # index - R::Array{Int,3} -# R::LinearIndices{3, Tuple{UnitRange{Int64}, UnitRange{Int64}, UnitRange{Int64}}} wet::BitArray{3} + interior::BitArray{3} + # I::Vector{CartesianIndex{3}} # index + # R::Array{Int,3} + # R::LinearIndices{3, Tuple{UnitRange{Int64}, UnitRange{Int64}, UnitRange{Int64}}} +end + +""" + Base.propertynames(γ::Grid) = (I,R,fieldnames(typeof(x))...) + + Do not store Cartesian and linear indices. + Compute them on demand. +""" +Base.propertynames(γ::Grid) = (:I,:R,fieldnames(typeof(γ))...) +function Base.getproperty(γ::Grid, d::Symbol) + if d === :I + return cartesianindex(γ.wet) + elseif d === :R + return linearindex(γ.wet) + else + return getfield(γ, d) + end end """ @@ -108,6 +133,19 @@ end # return Field(tracer,γ,:none,"unknown","unknown") #end +function Base.show(io::IO, mime::MIME{Symbol("text/plain")}, x::Field) + summary(io, x); println(io) + println(io, "Field size") + println(io, size(x.tracer)) + println(io, "Surface view") + show(io,mime,heatmap(transpose(x.tracer[:,:,1]))) + #show(io, mime, x.x) + # println(io, "\nsingular values:") + # show(io, mime, F.S) + # println(io, "\nV (right singular vectors):") + # show(io, mime, F.V) +end + """ struct BoundaryCondition @@ -157,6 +195,22 @@ function BoundaryCondition(tracer::Array{Float64,2},i::Vector{Float64},j::Vector return BoundaryCondition(tracer,i,j,k,dim,dimval,wet,:none,"unknown","unknown") end +""" + struct Source + + This structure describes Sources, which are + similar to Fields, but they may be + 1) non-negative + 2) have only interior mask +""" +struct Source{T} + tracer::Array{T,3} + γ::Grid + name::Symbol + longname::String + units::String + logscale::Bool +end # Credit to DrWatson.jl for these functions # Didn't want to add dependency for these small functions @@ -172,7 +226,7 @@ pkgdatadir(args...) = joinpath(pkgdatadir(), args...) pkgsrcdir() = joinpath(pkgdir(),"src") pkgsrcdir(args...) = joinpath(pkgsrcdir(), args...) -include(pkgsrcdir("TMIconfig.jl")) +include(pkgsrcdir("config.jl")) """ function trackpathways(TMIversion,latbox,lonbox) @@ -417,9 +471,10 @@ end - `fg!`: compute cost function and gradient in place - `γ`: grid """ -function sparsedatamap(u₀::Vector{T},Alu,b::Union{BoundaryCondition{T},NamedTuple{<:Any, NTuple{N1,BoundaryCondition{T}}}},u::Union{BoundaryCondition{T},NamedTuple{<:Any, NTuple{N2,BoundaryCondition{T}}}},y::Vector{T},W⁻,wis::Vector{Tuple{Interpolations.WeightedAdjIndex{2,T}, Interpolations.WeightedAdjIndex{2,T}, Interpolations.WeightedAdjIndex{2,T}}},locs,Q⁻,γ::Grid,iterations=10) where {N1, N2, T <: Real} +function sparsedatamap(u₀::Vector,Alu,b::Union{BoundaryCondition,NamedTuple},u::Union{BoundaryCondition,NamedTuple},y::Vector,W⁻,wis::Vector,locs,Q⁻,γ::Grid;q = nothing, r = 1.0,iterations=10) #where {N1, N2, T <: Real} +#function sparsedatamap(u₀::Vector{T},Alu,b::Union{BoundaryCondition{T},NamedTuple{<:Any, NTuple{N1,BoundaryCondition{T}}}},u::Union{BoundaryCondition{T},NamedTuple{<:Any, NTuple{N2,BoundaryCondition{T}}}},y::Vector{T},W⁻,wis::Vector{Tuple{Interpolations.WeightedAdjIndex{2,T}, Interpolations.WeightedAdjIndex{2,T}, Interpolations.WeightedAdjIndex{2,T}}},locs,Q⁻,γ::Grid,iterations=10) where {N1, N2, T <: Real} - fg!(F,G,x) = costfunction_point_obs!(F,G,x,Alu,b,u,y,W⁻,wis,locs,Q⁻,γ) + fg!(F,G,x) = costfunction_point_obs!(F,G,x,Alu,b,u,y,W⁻,wis,locs,Q⁻,γ,q₀=q,r=r) # a first guess: observed surface boundary conditions are perfect. # set surface boundary condition to the observations. @@ -487,8 +542,8 @@ end - `c`::Field """ function readfield(file,tracername,γ::Grid) - # The mode "r" stands for read-only. The mode "r" is the default mode and the parameter can be omitted. + # The mode "r" stands for read-only. The mode "r" is the default mode and the parameter can be omitted. ds = Dataset(file,"r") v = ds[tracername] @@ -504,13 +559,16 @@ function readfield(file,tracername,γ::Grid) # perform a check of file compatibility # with grid if sum(isnan.(tracer[γ.wet])) > 0 - println("readfield warning: NaN on grid") + error("readfield warning: NaN on grid") end + # check for non NaN or nonzero off grid # Need to rethink how to do this. - # if sum( !iszero(tracer[.!γ.wet])) > 0 - # println("readfield warning: nonzero value off grid") - # end + if sum(isnan.(tracer[.!(γ.wet)])) < length(isnan.(tracer[.!(γ.wet)])) + println("readfield warning: non-NaN value off grid") + println("resetting to NaN") + tracer[.!(γ.wet)] .= NaN + end #c = Field(tracer,γ) c = Field(tracer,γ,Symbol(tracername),longname,units) @@ -583,6 +641,17 @@ function writefield(file,field::Field{T}) where T <: Real return nothing end +function readsource(file,tracername,γ::Grid;logscale=false) + c = readfield(file,tracername,γ) + if logscale + ct = log.(c.tracer) + else + ct = c.tracer + end + q = Source(ct,c.γ,c.name,c.longname,c.units,logscale) + return q +end + """ function depthindex(I) Get the k-index (depth level) from the Cartesian index @@ -958,6 +1027,22 @@ function zeros(γ::Grid,name=:none,longname="unknown",units="unknown")::Field return d end +function zerosource(γ::Grid,name=:none,longname="unknown",units="unknown";logscale=false)::Source + T = eltype(γ.depth) + tracer = Array{T}(undef,size(γ.interior)) + tracer[γ.interior] .= zero(T) + tracer[.!γ.interior] .= zero(T)/zero(T) + return Source(tracer,γ,name,longname,units,logscale) +end + +function onesource(γ::Grid,name=:none,longname="unknown",units="unknown";logscale=false)::Source + T = eltype(γ.depth) + tracer = Array{T}(undef,size(γ.interior)) + tracer[γ.interior] .= one(T) + tracer[.!γ.interior] .= zero(T)/zero(T) + return Source(tracer,γ,name,longname,units,logscale) +end + """ function ones(γ::Grid,name=:none,longname="unknown",units="unknown")::Field @@ -979,7 +1064,7 @@ function ones(γ::Grid,name=:none,longname="unknown",units="unknown")::Field # set ocean to zero, land to NaN # consider whether land should be nothing or missing - println("calling one with ",T) + #println("calling one with ",T) tracer[γ.wet] .= Base.one(T) # add Base: error "should import Base" tracer[.!γ.wet] .= zero(T)/zero(T) # NaNs with right type @@ -1003,7 +1088,7 @@ function one(field::Field{T})::Field{T} where T <: Real # set ocean to zero, land to NaN # consider whether land should be nothing or missing - println("calling one with ",T) + #println("calling one with ",T) tracer[field.γ.wet] .= Base.one(T) # add Base: error "should import Base" tracer[.!field.γ.wet] .= zero(T)/zero(T) # NaNs with right type @@ -1036,7 +1121,7 @@ function one(field::Field{Float64}) # set ocean to zero, land to NaN # consider whether land should be nothing or missing - println("calling one with ",T) + #println("calling one with ",T) tracer[field.γ.wet] .= one(T) # add Base: error "should import Base" tracer[.!field.γ.wet] .= zero(T)/zero(T) # NaNs with right type @@ -1058,7 +1143,7 @@ function Field(field::Field{Float64}) # set ocean to zero, land to NaN # consider whether land should be nothing or missing - println("calling one with ",T) + # println("calling one with ",T) tracer[field.γ.wet] .= one(T) # add Base: error "should import Base" tracer[.!field.γ.wet] .= zero(T)/zero(T) # NaNs with right type @@ -1075,7 +1160,7 @@ function /(field::Field{Float64},scalar::Float64) # set ocean to zero, land to NaN # consider whether land should be nothing or missing - println("calling one with ",T) + # println("calling one with ",T) tracer[field.γ.wet] ./= scalar # add Base: error "should import Base" tracer[.!field.γ.wet] .= zero(T)/zero(T) # NaNs with right type @@ -1166,7 +1251,7 @@ end Initialize boundary condition with ones """ -function ones(dim::Int64,dimval::Int64,γ::Grid)::BoundaryCondition +function ones(dim::Int64,dimval::Int64,γ::Grid,name::Symbol,longname::String,units::String)::BoundaryCondition i,j,k,wet = boundaryconditionatts(dim,dimval,γ) @@ -1253,16 +1338,12 @@ function getboundarycondition(field::Field,dim,dimval)::BoundaryCondition end # Define maximum for Field to not include NaNs -maximum(c::Field) = maximum(c.tracer[c.γ.wet]) -minimum(c::Field) = minimum(c.tracer[c.γ.wet]) +Base.maximum(c::Union{Field,Source,BoundaryCondition}) = maximum(c.tracer[wet(c)]) +Base.minimum(c::Union{Field,Source,BoundaryCondition}) = minimum(c.tracer[wet(c)]) # uncomment this #mean(x::Field) = mean(x.tracer[x.γ.wet]) -# Define max/min for BoundaryCondition -maximum(b::BoundaryCondition) = maximum(b.tracer[b.wet]) -minimum(b::BoundaryCondition) = minimum(b.tracer[b.wet]) - """ `function \\(A,d::Field)::Field` Define left division for Fields @@ -1279,112 +1360,96 @@ end `function +(c::BoundaryCondition,d::BoundaryCondition)::BoundaryCondition` Define addition for Fields """ -function +(c::BoundaryCondition{T},d::BoundaryCondition{T})::BoundaryCondition{T} where T <: Real - - if c.wet != d.wet # check conformability - error("BoundaryCondition's not conformable for addition") +function add!(c::T,d::T) where T <: Union{Source,Field,BoundaryCondition} + if wet(c) != wet(d) # check conformability + error("TMI type not conformable for addition") end - array = zeros(c.wet) - e = BoundaryCondition(array,c.i,c.j,c.k,c.dim,c.dimval,c.wet) - - # a strange formulation to get - # return e to be correct - e.tracer[e.wet] += c.tracer[c.wet] - e.tracer[e.wet] += d.tracer[d.wet] - return e + # a strange formulation to do in-place addition + c.tracer[wet(c)] += d.tracer[wet(d)] end -""" - `function +(c::Field,d::Field)::Field` - Define addition for Fields -""" -function +(c::Field{T},d::Field{T})::Field{T} where T <: Real - # initialize output - if c.γ.wet != d.γ.wet # check conformability - error("Fields not conformable for addition") - end - - if !isequal(d.units,c.units) - error("Units not consistent:",d.units," vs ",c.units) - end - - e = zeros(d.γ,d.name,d.longname,d.units) - - # a strange formulation to get - # return e to be correct - e.tracer[e.γ.wet] += c.tracer[c.γ.wet] - e.tracer[e.γ.wet] += d.tracer[d.γ.wet] +function Base.:+(c::T,d::T) where T <: Union{Source,Field,BoundaryCondition} + e = deepcopy(c) + add!(e,d) return e end -""" - `function -(c::Field,d::Field)::Field` - Define addition for Fields -""" -function -(c::Field{T},d::Field{T})::Field{T} where T <: Real - # initialize output - if c.γ.wet !== d.γ.wet # check conformability - error("Fields not conformable for addition") +function subtract!(c::T,d::T) where T <: Union{Source,Field,BoundaryCondition} + if wet(c) != wet(d) # check conformability + error("TMI type not conformable for addition") end - e = zeros(d.γ) - e.tracer[e.γ.wet] += c.tracer[c.γ.wet] - e.tracer[e.γ.wet] -= d.tracer[d.γ.wet] + # a strange formulation to do in-place addition + c.tracer[wet(c)] -= d.tracer[wet(d)] +end + +function Base.:-(c::T,d::T) where T <: Union{Source,Field,BoundaryCondition} + e = deepcopy(c) + subtract!(e,d) return e end """ `function *(C,d::Field)::Field` Define scalar or matrix multiplication for fields + + one argument can be a number, but not both (type piracy?) """ -function *(C,d::Field{T})::Field{T} where T <: Real +# the order doesn't matter when multiplying by a scalar - e = zeros(d.γ) - e.tracer[e.γ.wet] += C*d.tracer[d.γ.wet] +Base.:*(c::Number,d::Union{Field,BoundaryCondition,Source}) = d*c +# right matrix multiply not handled +function Base.:*(c::AbstractArray,d::Union{Field,BoundaryCondition,Source}) + e = deepcopy(d) + mul!(c,e) return e end - -""" - `function *(C,d::BoundaryCondition)::BoundaryCondition` - Define scalar or matrix multiplication for BoundaryCondition`s -""" -function *(C,d::BoundaryCondition{T})::BoundaryCondition{T} where T <: Real - array = zeros(d.wet) - e = BoundaryCondition(array,d.i,d.j,d.k,d.dim,d.dimval,d.wet) - e.tracer[e.wet] += C*d.tracer[d.wet] +function Base.:*(d::T,c::Union{Number,T}) where T <: Union{Field,BoundaryCondition,Source} + e = deepcopy(d) + mul!(e,c) return e end -""" - `function *(c::Field,d::Field)::Field` - Field by field multiplication is element-by-element. -""" -function *(c::Field{T},d::Field{T})::Field{T} where T <: Real +function mul!(d::Union{Field,BoundaryCondition,Source},C::Number) + d.tracer[wet(d)] *= C #*d.tracer[wet(d)] +end +function mul!(c::T,d::T) where T <: Union{Field,BoundaryCondition,Source} # initialize output - if c.γ.wet != d.γ.wet # check conformability - error("Fields not conformable for addition") + if wet(c) != wet(d) # check conformability + error("Fields not conformable for multiplication") end if !isequal(d.units,c.units) error("Units not consistent:",d.units," vs ",c.units) end - - e = zeros(d.γ,d.name,d.longname,d.units) - e.tracer[e.γ.wet] += c.tracer[c.γ.wet] .* d.tracer[d.γ.wet] - return e + c.tracer[wet(c)] .*= d.tracer[wet(d)] end +# matrix multiplication with arrays/matrices: order matters +# d is mutated +function mul!(c::AbstractArray,d::Union{Field,BoundaryCondition,Source}) + # initialize output + if size(c,2) != sum(wet(d)) # check conformability + error("Number/Matrix-Field/BC/Source multiplication: not right size") + end + d.tracer[wet(d)] = c*d.tracer[wet(d)] +end +function mul!(d::Union{Field,BoundaryCondition,Source},c::AbstractArray) + # initialize output + if size(c,1) != sum(wet(d)) # check conformability + error("Field/BC/Source-Number/Matrix multiplication: not right size") + end + d.tracer[wet(d)] = d.tracer[wet(d)]*c +end """ `function *(c::Field,d::Field)::Field` Field by field multiplication is element-by-element. """ function dot(c::Field{T},d::Field{T})::T where T <: Real - e = c.tracer[c.γ.wet]' * d.tracer[d.γ.wet] return e end - """ function setsource!(d::Field,q::Field,r::Number) apply interior source q to the equation constraints d @@ -1394,22 +1459,27 @@ end - `r`::Number, default = 1.0, stoichiometric ratio """ function setsource!(d::Field{T},q::Field{T},r=1.0) where T<: Real - - # warning d.γ.wet and q.γ.wet must match d.tracer[d.γ.wet] -= r * q.tracer[q.γ.wet] - +end +function setsource!(d::Field{T},q::Source{T},r=1.0) where T<: Real + d.tracer[d.γ.interior] -= r * q.tracer[q.γ.interior] end """ - function vec(u) + function gsetsource!(gq::Field{T},gd::Field{T},r=1.0) - Turn a collection of controls into a vector - for use with Optim.jl. - An in-place version of this function would be handy. + Adjoint to `setsource!` """ -function vec(u::BoundaryCondition{T}) where T <: Real - uvec = u.tracer[u.wet] - return uvec +function gsetsource!(gq::Field{T},gd::Field{T},r) where T<: Real + gq.tracer[gq.γ.wet] -= r * gd.tracer[gd.γ.wet] +end +function gsetsource!(gq::Source{T},gd::Field{T},r) where T<: Real + gq.tracer[gq.γ.interior] -= r * gd.tracer[gd.γ.interior] +end +function gsetsource(gd::Field{T},q::Union{Field,Source},r) where T<: Real + gq = 0.0 * q + gsetsource!(gq,gd,r) + return gq end """ @@ -1419,49 +1489,32 @@ end for use with Optim.jl. An in-place version of this function would be handy. """ -function vec(u::NamedTuple{<:Any,NTuple{N,BoundaryCondition{T}}}) where {N, T <: Real} +vec(u::Field) = u.tracer[u.γ.wet] +vec(u::BoundaryCondition) = u.tracer[u.wet] +vec(u::Source) = u.tracer[u.γ.interior] +function vec(u::NamedTuple) + T = eltype(values(u)[1].tracer) + #T = eltype(u) uvec = Vector{T}(undef,0) for v in u - append!(uvec,v.tracer[v.wet]) + #append!(uvec,v.tracer[v.wet]) + append!(uvec,vec(v)) end return uvec end """ - function unvec!(u,uvec) - - Undo the operations by vec(u) - Needs to update u because attributes of - u need to be known at runtime. -""" -function unvec!(u::NamedTuple{<:Any,NTuple{N,BoundaryCondition{T}}},uvec::Vector{T}) where {N, T <: Real} - - counter = 0 - for v in u - n = sum(v.wet) - v.tracer[v.wet] = uvec[counter+1:counter+n] - counter += n - end -end - -""" - function unvec(utemplate,uvec) + function unvec(u,uvec) + Replace u with new u Undo the operations by vec(u) Needs to update u because attributes of u need to be known at runtime. """ -function unvec(utemplate::NamedTuple{<:Any,NTuple{N,BoundaryCondition{T}}},uvec::Vector{T}) where {N, T <: Real} - - counter = 0 - vals = Vector{BoundaryCondition}(undef,N) - for (ii,vv) in enumerate(utemplate) - n = sum(vv.wet) - vals[ii] = unvec(vv, uvec[counter+1:counter+n]) - counter += n - end - u = (;zip(keys(utemplate), vals)...) +function unvec(u₀::Union{NamedTuple,Field,BoundaryCondition},uvec::Vector) #where T <: Real + u = deepcopy(u₀) + unvec!(u,uvec) return u end @@ -1472,33 +1525,20 @@ end Needs to update u because attributes of u need to be known at runtime. """ -function unvec!(u::BoundaryCondition{T},uvec::Vector{T}) where T <: Real - - I = findall(u.wet) - counter = 0 - for i in I - counter +=1 - # doesn't work to pass u back - #u.tracer[u.wet] = uvec - u.tracer[i] = uvec[counter] +function unvec!(u::Union{BoundaryCondition{T},Field{T},Source{T}},uvec::Vector{T}) where T <: Real + I = findall(wet(u)) # findall seems slow + for (ii,vv) in enumerate(I) + u.tracer[vv] = uvec[ii] end - end - -""" - function unvec(u,uvec) - - Replace u with new u - Undo the operations by vec(u) - Needs to update u because attributes of - u need to be known at runtime. -""" -function unvec(utemplate::BoundaryCondition{T},uvec::Vector{T}) where T <: Real - - tracer = zeros(utemplate.wet) - tracer[utemplate.wet] = uvec - u = BoundaryCondition(tracer,utemplate.i,utemplate.j,utemplate.k,utemplate.dim,utemplate.dimval,utemplate.wet) - return u +function unvec!(u::NamedTuple,uvec::Vector) #where {N, T <: Real} + nlo = 1 + nhi = 0 + for v in u + nhi += sum(wet(v)) + unvec!(v,uvec[nlo:nhi]) + nlo = nhi + 1 + end end """ @@ -1584,13 +1624,8 @@ function synthetic_observations(TMIversion,variable,γ) #ntrue += rand(Normal(),length(σθ[γ.wet])) .* σθ[γ.wet] ntrue = Field(rand(Normal(),size(γ.wet)),γ,θtrue.name,θtrue.longname,θtrue.units) - println(θtrue.units) - println(ntrue.units) - ntrue *= σθ - y = θtrue + ntrue - #y = θtrue .+ ntrue # get cost function (J) based on model misfit @@ -1799,11 +1834,11 @@ function costfunction_gridded_obs!(J,guvec,uvec::Vector{T},Alu,b₀::Union{Bound end """ - function costfunction_point_obs(uvec::Vector{T},Alu,b₀::BoundaryCondition{T},u₀::BoundaryCondition{T},y::Vector{T},Wⁱ::Diagonal{T, Vector{T}},wis,locs,Q⁻,γ::Grid) where T <: Real + function costfunction_point_obs(uvec::Vector{T},Alu,b₀::BoundaryCondition{T},u₀::BoundaryCondition{T},y::Vector{T},Wⁱ::Diagonal{T, Vector{T}},wis,locs,Q⁻,γ::Grid;q=nothing,r=1.0) where T <: Real - squared model-data misfit for pointwise data - controls are a vector input for Optim.jl - Issue #1: couldn't figure out how to nest with costfunction_obs! + Squared model-data misfit for pointwise data. + Controls are a vector input for Optim.jl. + Core numerics handled by `costfunction_point_obs`. # Arguments - `uvec`: controls, vector format @@ -1815,36 +1850,20 @@ end - `locs`: data locations (lon,lat,depth) - `Q⁻`: weights for control vector - `γ`: grid +# Optional +- `q::Field`: interior source +- `r::Number`: scalar factor for source # Output - `J`: cost function of sum of squared misfits - `gJ`: derivative of cost function wrt to controls """ -function costfunction_point_obs(uvec::Vector{T},Alu,b₀::Union{BoundaryCondition{T},NamedTuple{<:Any, NTuple{N1,BoundaryCondition{T}}}},u₀::Union{BoundaryCondition{T},NamedTuple{<:Any, NTuple{N2,BoundaryCondition{T}}}},y::Vector{T},Wⁱ::Diagonal{T, Vector{T}},wis,locs,Q⁻,γ::Grid) where {N1, N2, T <: Real} - - # control penalty and gradient - Jcontrol = uvec'*(Q⁻*uvec) - guvec = 2*(Q⁻*uvec) - - u = unvec(u₀,uvec) - b = adjustboundarycondition(b₀,u) # combine b₀, u +function costfunction_point_obs(uvec::Vector,Alu,b::Union{BoundaryCondition,NamedTuple},u::Union{BoundaryCondition,NamedTuple},y::Vector,Wⁱ::Diagonal,wis,locs,Q⁻,γ::Grid;q=nothing,r=1.0) - c = steadyinversion(Alu,b,γ) # gives the misfit - - # observe at right spots - ỹ = observe(c,wis,γ) - n = ỹ - y - - Jdata = n ⋅ (Wⁱ * n) # dot product - J = Jdata + Jcontrol - - gn = 2Wⁱ * n - gỹ = gn - - gc = gobserve(gỹ,c,locs) - gb = gsteadyinversion(gc, Alu, b, γ) - gu = gadjustboundarycondition(gb,u) - guvec += vec(gu) - + J = 0.0 + guvec = 0.0.*uvec # same size + #gu = unvec(u,0 .* uvec) + J = costfunction_point_obs!(J,guvec,uvec,Alu,b,u,y,Wⁱ,wis,locs,Q⁻,γ,q₀=q,r=r) + return J, guvec end @@ -1868,24 +1887,42 @@ end - `Q⁻`: weights for control vector - `γ`: grid """ -function costfunction_point_obs!(J,guvec,uvec::Vector{T},Alu,b₀::Union{BoundaryCondition{T},NamedTuple{<:Any, NTuple{N1,BoundaryCondition{T}}}},u₀::Union{BoundaryCondition{T},NamedTuple{<:Any, NTuple{N2,BoundaryCondition{T}}}},y::Vector{T},Wⁱ::Diagonal{T, Vector{T}},wis,locs,Q⁻,γ::Grid) where {N1, N2, T <: Real} +function costfunction_point_obs!(J,guvec::Union{Nothing,Vector},uvec::Vector,Alu,b₀::Union{BoundaryCondition,NamedTuple},u₀::Union{BoundaryCondition,NamedTuple},y::Vector,Wⁱ::Diagonal,wis,locs,Q⁻,γ::Grid;q₀=nothing,r=1.0) u = unvec(u₀,uvec) b = adjustboundarycondition(b₀,u) # combine b₀, u - c = steadyinversion(Alu,b,γ) # gives the misfit + + if !isnothing(q₀) + # careful with scope of c + q = adjustsource(q₀,u) + c = steadyinversion(Alu,b,γ,q=q,r=r) + else + c = steadyinversion(Alu,b,γ) + end # observe at right spots ỹ = observe(c,wis,γ) n = ỹ - y if guvec != nothing - # adjoint equations - gtmp = 2*(Q⁻*uvec) + ## start adjoint model + + # initialize to zero + gu = unvec(u,0 .* uvec) + + gtmp = 2*(Q⁻*uvec) # control penalty gradient gn = 2Wⁱ * n gỹ = gn gc = gobserve(gỹ,c,locs) - gb = gsteadyinversion(gc, Alu, b, γ) - gu = gadjustboundarycondition(gb,u) + + if !isnothing(q₀) + gb,gq = gsteadyinversion(gc,Alu,b,γ,q=q,r=r) + gadjustsource!(gu,gq,q) # pass q to linearize logscale version + else + gb = gsteadyinversion(gc, Alu, b, γ) + end + + gadjustboundarycondition!(gu,gb) gtmp += vec(gu) for (ii,vv) in enumerate(gtmp) guvec[ii] = vv @@ -1896,6 +1933,8 @@ function costfunction_point_obs!(J,guvec,uvec::Vector{T},Alu,b₀::Union{Boundar # control penalty and gradient Jcontrol = uvec'*(Q⁻*uvec) Jdata = n ⋅ (Wⁱ * n) # dot product + #println("Jcontrol:",Jcontrol) + #println("Jdata:",Jdata) return Jdata + Jcontrol end end @@ -1905,7 +1944,7 @@ end invert for a steady-state tracer distribution # Arguments - `Alu`: LU decomposition of water-mass matrix -- `b`: boundary condition +- `b`: boundary condition, assumed to be surface boundary condition - `γ`::Grid # Optional Arguments - `q`: interior sources/sinks of phosphate @@ -1915,8 +1954,6 @@ end """ function steadyinversion(Alu,b::BoundaryCondition{T},γ::Grid;q=nothing,r=1.0)::Field{T} where T <: Real - #println("running steady inversion") - # preallocate Field for equation constraints d = zeros(γ) @@ -1929,17 +1966,12 @@ function steadyinversion(Alu,b::BoundaryCondition{T},γ::Grid;q=nothing,r=1.0):: setsource!(d,q,r) end - # define ldiv with fields - # println(b.units) - #c = zeros(d.γ,b.name,b.longname,b.units) - #println(c.units) c = Alu \ d - return c end """ - function gsteadyinversion(Alu,b;q=nothing,r=1.0) + function gsteadyinversion(gc,Alu,b;q=nothing,r=1.0) ADJOINT invert for a steady-state tracer distribution @@ -1953,20 +1985,18 @@ end # Output - `c`::Field, steady-state tracer distribution """ -function gsteadyinversion(gc,Alu,b::BoundaryCondition{T},γ::Grid;q=nothing,r=1.0)::BoundaryCondition{T} where T <: Real - +function gsteadyinversion(gc::Field,Alu,b::Union{BoundaryCondition,NamedTuple},γ::Grid;q=nothing,r=1.0) #where T <: Real #println("running adjoint steady inversion") - gd = Alu' \ gc + gb = gsetboundarycondition(gd,b) - # still need to develop this - # and to find a way to return the value if !isnothing(q) - gq = gsetsource!(gd,d,q,r) + gq = gsetsource(gd,q,r) + return gb,gq + else + return gb end - gb = gsetboundarycondition(gd,b) - return gb end """ @@ -1989,7 +2019,6 @@ function steadyinversion(Alu,b::NamedTuple{<:Any, NTuple{N,BoundaryCondition{T}} end # Warning: doesn't this need to loop over the NamedTuple? - # define ldiv with fields #c = zeros(d.γ,b.name,b.longname,b.units) @@ -1998,29 +2027,29 @@ function steadyinversion(Alu,b::NamedTuple{<:Any, NTuple{N,BoundaryCondition{T}} return c end -""" - function gsteadyinversion(gc::Field{T},Alu,b::NamedTuple{<:Any, NTuple{N,BoundaryCondition{T}}},γ::Grid;q=nothing,r=1.0)::Field{T} where {N, T <: Real} +# """ +# function gsteadyinversion(gc::Field{T},Alu,b::NamedTuple{<:Any, NTuple{N,BoundaryCondition{T}}},γ::Grid;q=nothing,r=1.0)::Field{T} where {N, T <: Real} - ADDJOINT steady inversion for b::NamedTuple -""" -function gsteadyinversion(gc::Field{T},Alu,b::NamedTuple{<:Any, NTuple{N,BoundaryCondition{T}}},γ::Grid;q=nothing,r=1.0) where {N, T <: Real} +# ADDJOINT steady inversion for b::NamedTuple +# """ +# function gsteadyinversion(gc::Field{T},Alu,b::NamedTuple{<:Any, NTuple{N,BoundaryCondition{T}}},γ::Grid;q=nothing,r=1.0) where {N, T <: Real} - #println("running adjoint steady inversion") +# #println("running adjoint steady inversion") - # sensitivity of Field of equation constraints - gd = Alu' \ gc +# # sensitivity of Field of equation constraints +# gd = Alu' \ gc - # still need to develop this - # and to find a way to return the value - if !isnothing(q) - gq = gsetsource!(gd,d,q,r) - end +# # still need to develop this +# # and to find a way to return the value +# if !isnothing(q) +# gq = gsetsource!(gd,d,q,r) +# end - # update d with the boundary condition b - gb = gsetboundarycondition(gd,b) +# # update d with the boundary condition b +# gb = gsetboundarycondition(gd,b) - return gb -end +# return gb +# end """ function wetlocation(γ) @@ -2095,8 +2124,15 @@ zerosouthboundary(γ,name=:none,longname="unknown",units="unknown") = zeros(2,1, zerowestboundary(γ,name=:none,longname="unknown",units="unknown") = zeros(1,1,γ,name,longname,units)::BoundaryCondition -## Warning: need to update this next function -onesurfaceboundary(γ) = ones(3,1,γ)::BoundaryCondition +onesurfaceboundary(γ,name=:none,longname="unknown",units="unknown") = ones(3,1,γ,name,longname,units)::BoundaryCondition + +onenorthboundary(γ,name=:none,longname="unknown",units="unknown") = ones(2,maximum(latindex(γ.I)),γ,name,longname,units)::BoundaryCondition + +oneeastboundary(γ,name=:none,longname="unknown",units="unknown") = ones(1,maximum(lonindex(γ.I)),γ,name,longname,units)::BoundaryCondition + +onesouthboundary(γ,name=:none,longname="unknown",units="unknown") = ones(2,1,γ,name,longname,units)::BoundaryCondition + +onewestboundary(γ,name=:none,longname="unknown",units="unknown") = ones(1,1,γ,name,longname,units)::BoundaryCondition getsurfaceboundary(c::Field) = getboundarycondition(c,3,1)::BoundaryCondition getnorthboundary(c::Field) = getboundarycondition(c,2,maximum(latindex(c.γ.I)))::BoundaryCondition @@ -2183,25 +2219,46 @@ function gsetboundarycondition(gd::Field{T},b::NamedTuple{<:Any, NTuple{N,Bounda end """ - function adjustboundarycondition!(b::BoundaryCondition{T},u::BoundaryCondition{T}) where T <: Real + function adjustboundarycondition!(b::Union{BoundaryCondition,NamedTuple},u::Union{BoundaryCondition,NamedTuple}) - adjust the (one) boundary condition - problem: passes back a mutated b + adjust all boundary conditions b that are described in u """ -function adjustboundarycondition!(b::BoundaryCondition{T},u::BoundaryCondition{T}) where T <: Real - b.tracer[b.wet] += u.tracer[u.wet] # write it out so b changes when returned +adjustboundarycondition(b::BoundaryCondition,u::BoundaryCondition) = b + u +function adjustboundarycondition(b::Union{BoundaryCondition,NamedTuple},u::Union{BoundaryCondition,NamedTuple}) + bnew = deepcopy(b) + adjustboundarycondition!(bnew,u) + return bnew end """ - function adjustboundarycondition(b::BoundaryCondition{T},u::BoundaryCondition{T}) where T <: Real + function adjustboundarycondition!(b::Union{BoundaryCondition,NamedTuple},u::Union{BoundaryCondition,NamedTuple}) - adjust the (one) boundary condition -""" -function adjustboundarycondition(b::BoundaryCondition{T},u::BoundaryCondition{T}) where T <: Real + adjust all boundary conditions b that are described in u - bnew = b + u - return bnew + warning: if u doesn't contain any boundary condition adjustments, + nothing will change. +""" +function adjustboundarycondition!(b::BoundaryCondition,u::BoundaryCondition) + # write it out so b changes when returned + b.tracer[b.wet] += u.tracer[u.wet] +end +function adjustboundarycondition!(b::BoundaryCondition,u::NamedTuple) #where {N1, N2, T <: Real} + # if not explicit, assume a surface boundary condition + bkey = :surface # + if haskey(u,bkey) + adjustboundarycondition!(b,u[bkey]) + end +end +function adjustboundarycondition!(b::NamedTuple,u::NamedTuple) + for bkey in keys(b) + adjustboundarycondition!(b[bkey],u) + end +end +# Seems not to be general because gu overwritten. +function gadjustboundarycondition(gb::BoundaryCondition{T},u::BoundaryCondition{T}) where T <: Real + gu = gb + return gu end """ @@ -2212,53 +2269,115 @@ end Keep this function so that calling functions can look alike. Could probably combine with lower function, use Union type """ -function gadjustboundarycondition(gb::BoundaryCondition{T},u::BoundaryCondition{T}) where T <: Real - gu = gb - return gu +function gadjustboundarycondition!(gu::BoundaryCondition,gb::BoundaryCondition) + gu.tracer[gu.wet] += gb.tracer[gb.wet] end - -""" - function adjustboundarycondition!(b::NamedTuple{<:Any, NTuple{N1,BoundaryCondition{T}}},u::NamedTuple{<:Any, NTuple{N2,BoundaryCondition{T}}}) where N1, N2, T <: Real - - adjust all boundary conditions b that are described in u -""" -function adjustboundarycondition!(b::NamedTuple{<:Any, NTuple{N1,BoundaryCondition{T}}},u::NamedTuple{<:Any, NTuple{N2,BoundaryCondition{T}}}) where {N1, N2, T <: Real} - - ukeys = keys(u) - for ukey in keys(u) - b[ukey].tracer[b[ukey].wet] += u[ukey].tracer[b[ukey].wet] +function gadjustboundarycondition!(gu::NamedTuple,gb::BoundaryCondition) + bkey = :surface + if haskey(gu,bkey) + gadjustboundarycondition!(gu[bkey],gb) end - end - -""" - function adjustboundarycondition(b::NamedTuple{<:Any, NTuple{N1,BoundaryCondition{T}}},u::NamedTuple{<:Any, NTuple{N2,BoundaryCondition{T}}}) where N1, N2, T <: Real - - adjust all boundary conditions b that are described in u -""" -function adjustboundarycondition(b::NamedTuple{<:Any, NTuple{N1,BoundaryCondition{T}}},u::NamedTuple{<:Any, NTuple{N2,BoundaryCondition{T}}}) where {N1, N2, T <: Real} - - bnew = deepcopy(b) - ukeys = keys(u) - for ukey in keys(u) - bnew[ukey].tracer[bnew[ukey].wet] += u[ukey].tracer[bnew[ukey].wet] +function gadjustboundarycondition!(gu::NamedTuple,gb::NamedTuple) + for bkey in keys(gb) + gadjustboundarycondition!(gu,gb[bkey]) end - return bnew end -""" - function gadjustboundarycondition!(b::BoundaryCondition{T},u::BoundaryCondition{T}) where T <: Real - - ADJOINT CODE - adjust the (one) boundary condition - Just copy the variable. - Keep this function so that calling functions can look alike. -""" -function gadjustboundarycondition(gb::NamedTuple{<:Any, NTuple{N1,BoundaryCondition{T}}},u::NamedTuple{<:Any, NTuple{N2,BoundaryCondition{T}}}) where {N1, N2, T <: Real} +#function gadjustboundarycondition(gb::NamedTuple{<:Any, NTuple{N1,BoundaryCondition{T}}},u::NamedTuple{<:Any, NTuple{N2,BoundaryCondition{T}}}) where {N1, N2, T <: Real} +function gadjustboundarycondition(gb::Union{NamedTuple,BoundaryCondition},u::NamedTuple) #where {N1, N2, T <: Real} gu = gb[keys(u)] # grab the parts of the named tuple corresponding to u return gu end +function adjustsource(q₀::Union{Source,Field,NamedTuple},u::Union{Source,Field,NamedTuple}) + q = deepcopy(q₀) + adjustsource!(q,u) + return q +end + +function adjustsource!(q::Field,u::Field) + # write it out so b changes when returned + q.tracer[q.γ.wet] += u.tracer[u.γ.wet] +end +function adjustsource!(q::Source,u::Source) + if q.logscale && u.logscale + q.tracer[q.γ.interior] += u.tracer[u.γ.interior] + q.tracer[q.γ.interior] = exp.(q.tracer[q.γ.interior]) + q.logscale = false + elseif ~q.logscale && u.logscale + q.tracer[q.γ.interior] = log.(q.tracer[q.γ.interior]) + q.tracer[q.γ.interior] += u.tracer[u.γ.interior] + q.tracer[q.γ.interior] = exp.(q.tracer[q.γ.interior]) + elseif q.logscale && ~u.logscale + error("not implemented: logscale would not lead to non-negative source") + # u.tracer[u.γ.interior] = log.(u.tracer[u.γ.interior]) + # q.tracer[q.γ.interior] += u.tracer[u.γ.interior] + # q.tracer[q.γ.interior] = exp.(q.tracer[q.γ.interior]) + else + q.tracer[q.γ.interior] += u.tracer[u.γ.interior] + end +end +function adjustsource!(q::NamedTuple,u::NamedTuple) #where {N1, N2, T <: Real} + for qkey in keys(q) + if haskey(u,qkey) + adjustsource!(q[qkey],u[qkey]) + else + error("adjustsource!: u doesn't have qkey ",qkey) + end + end +end +function adjustsource!(q::Union{Field,Source},u::NamedTuple) #where {N1, N2, T <: Real} + qkey = :source + if haskey(u,qkey) + adjustsource!(q,u[qkey]) + else + error("adjustsource!: u doesn't have source info") + end +end + +function gadjustsource!(gu::Field,gq::Field) + # write it out so b changes when returned + gu.tracer[gu.γ.wet] += gq.tracer[gq.γ.wet] +end +function gadjustsource!(gu::Source,gq::Source,q₀::Source) + q = deepcopy(q₀) + if gq.logscale && gu.logscale + error("not implemented") + # gu.tracer[gu.γ.interior] += gq.tracer[gq.γ.interior] + # q.tracer[q.γ.interior] = exp.(q.tracer[q.γ.interior]) + # q.logscale = false + gq.logscale = false + + elseif gu.logscale + q.tracer[q.γ.interior] = log.(q.tracer[q.γ.interior]) + gq.tracer[gq.γ.interior] = gq.tracer[gq.γ.interior] .* exp.(q.tracer[q.γ.interior]) + gu.tracer[gu.γ.interior] += gq.tracer[gq.γ.interior] + # no need to adjoint this next line? + elseif gq.logscale + error("not implemented: logscale would not lead to non-negative source") + # u.tracer[u.γ.interior] = log.(u.tracer[u.γ.interior]) + # q.tracer[q.γ.interior] += u.tracer[u.γ.interior] + # q.tracer[q.γ.interior] = exp.(q.tracer[q.γ.interior]) + else + gu.tracer[gu.γ.interior] += gq.tracer[gq.γ.interior] + end +end + +function gadjustsource!(gu::NamedTuple,gq::T,q::T) where T <: Union{Source,Field} + qkey = :source + if haskey(gu,qkey) + gadjustsource!(gu[qkey],gq,q) + end +end +function gadjustsource!(gu::NamedTuple,gq::T,q::T) where T <: NamedTuple + for qkey in keys(gq) + if haskey(gu,qkey) + gadjustsource!(gu[qkey],gq[qkey],q[qkey]) + end + end +end + """ function section View latitude-depth slice of field @@ -2293,183 +2412,10 @@ function planview(c::Field{T},depth)::Array{T,2} where T <: Real return cplan end -""" - function control2state(tracer2D,γ) - turn 2D surface field into 3D field with zeroes below surface -# Arguments -- `tracer2D`:: 2D surface tracer field -- `wet`::BitArray mask of ocean points -# Output -- `tracer3D`:: 3d tracer field with NaN on dry points -""" -function control2state(tracer2D::Matrix{T},wet) where T<: Real - # preallocate - tracer3D = Array{T}(undef,size(wet)) - - # set ocean to zero, land to NaN - # consider whether land should be nothing or missing - tracer3D[wet] .= zero(T) - tracer3D[.!wet] .= zero(T)/zero(T) - tracer3D[:,:,1] = tracer2D - return tracer3D -end - -""" - function surfacecontrol2field(usfc,γ.wet) - turn surface control vector into 3D field with zeroes below surface -# Arguments -- `usfc`:: surface control vector -- `wet`::BitArray mask of ocean points -# Output -- `tracer3D`:: 3d tracer field with NaN on dry points -""" -function surfacecontrol2field(usfc::Vector{T},wet) where T<: Real - # preallocate - tracer3D = Array{T}(undef,size(wet)) - - # set ocean to zero, land to NaN - # consider whether land should be nothing or missing - tracer3D[wet] .= zero(T) - tracer3D[.!wet] .= zero(T)/zero(T) - tracer3D[:,:,1][wet[:,:,1]] = usfc - return tracer3D -end - -""" - function Γsfc - Γsfc anonymously calls surfacecontrol2field -""" -Γsfc = surfacecontrol2field +wet(a::BoundaryCondition) = a.wet +wet(a::Field) = a.γ.wet +wet(a::Source) = a.γ.interior -""" - function surfacecontrol2field!(c,u,γ) - Add surface control vector to existing 3D field -# Arguments -- `c`:: state field, 3d tracer field with NaN on dry points, modified by function -- `usfc`:: surface control vector -- `wet`::BitArray mask of ocean points -""" -function surfacecontrol2field!(c::Array{T,3},usfc::Vector{T},γ) where T<: Real - #c[:,:,1][wet[:,:,1]] .+= u # doesn't work -# [c[γ.I[ii][1],γ.I[ii][2],γ.I[ii][3]] += u[ii] for ii ∈ eachindex(γ.I) if γ.I[ii][3] == 1] - list = surfaceindex(γ.I) - [c[γ.I[ii]] += usfc[list[ii]] for ii ∈ eachindex(γ.I) if γ.I[ii][3] == 1] -end - -""" - function surfacecontrol2field!(c,u,γ) - Add surface control vector to tracer vector -# Arguments -- `c`:: state field, 3d tracer field with NaN on dry points, modified by function -- `u`:: surface control vector -- `wet`::BitArray mask of ocean points -""" -function surfacecontrol2field!(c::Vector{T},u::Vector{T},γ) where T<: Real - list = surfaceindex(γ.I) - [c[ii] += u[list[ii]] for ii ∈ eachindex(γ.I) if γ.I[ii][3] == 1] -end - -""" - function Γsfc! - Γsfc! anonymously calls surfacecontrol2field! -""" -Γsfc! = surfacecontrol2field! - -function field2obs(cvec,wis,γ) - # interpolate onto data points - N = length(wis) - sumwis = Vector{Float64}(undef,N) - list = vcat(1:length(γ.lon),1) - - # perhaps the most clever line in TMI.jl? - wetwrap = view(γ.wet,list,:,:) - - # some interpolation weights on land, oh no - # sum up all weights in ocean - [sumwis[i] = Interpolations.InterpGetindex(wetwrap)[wis[i]...] for i in eachindex(wis)] - - # reconstruct the observations - ỹ = Vector{Float64}(undef,N) - c̃ = tracerinit(γ.wet) - c̃[γ.wet] = cvec - replace!(c̃,NaN=>0.0) - cwrap = view(c̃,list,:,:) - - # divide by sum of all ocean weights so that this is still a true average - [ỹ[i] = Interpolations.InterpGetindex(cwrap)[wis[i]...]/sumwis[i] for i in eachindex(wis)] - return ỹ -end - -""" - function E - E anonymously calls field2obs -""" -E = field2obs - -""" - function tracerinit(wet,vec,I) - initialize tracer field on TMI grid - perhaps better to have a tracer struct and constructor -# Arguments -- `wet`:: BitArray mask of ocean points -- `vec`:: vector of values at wet points -- `I`:: Cartesian Index for vector -# Output -- `field`:: 3d tracer field with NaN on dry points -""" -function tracerinit(vec,I,wet) - - # preallocate - T = eltype(vec) - field = Array{T}(undef,size(wet)) - fill!(field,zero(T)/zero(T)) - - #- a comprehension - [field[I[n]]=vec[n] for n ∈ eachindex(I)] - return field -end - -""" - function vec2fld - Transfer a vector to a 3D field with accounting for ocean bathymetry -# Arguments -- `vector`: field in vector form (no land points) -- `I`: cartesian indices of ocean points -# Output -- `field`: field in 3d form including land points (NaN) -""" -function vec2fld(vector::Vector{T},I::Vector{CartesianIndex{3}}) where T<:Real - # choose NaN for now, zero better? nothing better? - fillvalue = zero(T)/zero(T) - - nx = maximum(I)[1] - ny = maximum(I)[2] - nz = maximum(I)[3] - - # faster instead to allocate as undef and then fill! ? - field = (NaN .* zero(T)) .* zeros(nx,ny,nz) - - # a comprehension - [field[I[n]]=vector[n] for n ∈ eachindex(I)] - return field -end - -""" - function fld2vec - Transfer 3D field with accounting for ocean bathymetry to a vector without land points. - This is done more easily with a BitArray mask, i.e., vector = field[mask]. - This function may be removed in the future. -# Arguments -- `field`: field in 3d form including land points (NaN) -- `I`: cartesian indices of ocean points -# Output -- `vector`: field in vector form (no land points) -""" -function fld2vec(field::Array{T,3},I::Vector{CartesianIndex{3}}) where T<:Real - vector = Vector{T}(undef,length(I)) - #- a comprehension - [vector[n] = field[I[n]] for n ∈ eachindex(I)]; - return vector - end +include("deprecated.jl") end diff --git a/src/TMIconfig.jl b/src/config.jl similarity index 90% rename from src/TMIconfig.jl rename to src/config.jl index 085986e1..8eeba8a4 100644 --- a/src/TMIconfig.jl +++ b/src/config.jl @@ -14,8 +14,6 @@ function config_from_nc(TMIversion) TMIfile = download_ncfile(TMIversion) - γ = Grid(TMIfile) - println("A") @time A = watermassmatrix(TMIfile) @@ -23,6 +21,8 @@ function config_from_nc(TMIversion) println("Alu") @time Alu = lu(A) + γ = Grid(TMIfile) + # would be good to make this optional println("L=") @time L = circulationmatrix(TMIfile,A,γ) @@ -63,6 +63,12 @@ function download_ncfile(TMIversion::String) shellscript = pkgsrcdir("read_nc_modern_180x90x33_GH11_GH12.sh") run(`sh $shellscript`) mv(joinpath(pwd(),"TMI_"*TMIversion*".nc"),TMIfile) + elseif TMIversion == "nordic_201x115x46_B23" + println("workaround for regional Nordic Seas file") + shellscript = pkgsrcdir("read_nc_nordic_201x115x46_B23.sh") + run(`sh $shellscript`) + TMIfile = TMIversion*".nc" + mv(joinpath(pwd(),"TMI_"*TMIversion*".nc"),TMIfile) else println("read via GoogleDrive.jl") #- `url`: Google Drive URL for data @@ -85,31 +91,48 @@ function gridinit(TMIfile) # Output - `γ::Grid`: TMI grid struct """ -function Grid(TMIfile::String) +function Grid(TMIfile::String; A = watermassmatrix(TMIfile)) # get properties of grid lon,lat,depth = gridprops(TMIfile) + + # make ocean mask + wet = wetmask(TMIfile,length(lon),length(lat),length(depth)) + + # make interior mask + interior = interiormask(A,wet,length(lon),length(lat),length(depth)) + + # do not store: compute on demand + #R = linearindex(wet) + #γ = Grid(lon,lat,depth,I,R,wet) + + return Grid(lon,lat,depth,wet,interior) +end + +function wetmask(TMIfile,nx,ny,nz) # read Cartesian Index from file. I = cartesianindex(TMIfile) - # make a mask # first choice: smaller but inconsistent with input grid #wet = falses(maximum(I)[1],maximum(I)[2],maximum(I)[3]) - wet = falses(length(lon),length(lat),length(depth)) + wet = falses(nx,ny,nz) wet[I] .= 1 + return wet +end - R = linearindex(wet) - - γ = Grid(lon,lat,depth,I,R,wet) - - return γ +function interiormask(A,wet,nx,ny,nz) + interior = falses(nx,ny,nz) + I = cartesianindex(wet) + list = findall(.!(isone.(sum(abs.(A),dims=2)))) + interior[I[list]] .= true + return interior end """ Configure TMI environment from original MATLAB output """ function config_from_mat(TMIversion) - + # repetitive and long TMIfile = download_matfile(TMIversion) # # make a sample field from zyx cartesian indices @@ -120,19 +143,21 @@ function config_from_mat(TMIversion) wet[Izyx] .= 1 I = cartesianindex(wet) - R = linearindex(wet) Azyx = watermassmatrix(TMIfile) + A = matrix_zyx2xyz(TMIfile,Azyx,R) # # LU factorization for efficient matrix inversions Alu = lu(A) - # get properties of grid - lat,lon,depth = gridprops(TMIfile) + γ = Grid(TMIfile,A=A) - γ = Grid(lon,lat,depth,I,R,wet) + # # get properties of grid + # lat,lon,depth = gridprops(TMIfile) + + # γ = Grid(lon,lat,depth,I,R,wet) # need to make this optional L = circulationmatrix(TMIfile,γ) @@ -175,6 +200,11 @@ function download_matfile(TMIversion::String) shellscript = pkgsrcdir("read_mat_modern_180x90x33_GH11_GH12.sh") run(`sh $shellscript`) mv(joinpath(pwd(),"TMI_"*TMIversion*".mat.gz"),TMIfilegz,force=true) + elseif TMIversion == "nordic_201x115x46_B23" + println("workaround for regional Nordic Seas file") + shellscript = pkgsrcdir("read_mat_nordic_201x115x46_B23.sh") + run(`sh $shellscript`) + mv(joinpath(pwd(),"TMI_"*TMIversion*".mat.gz"),TMIfilegz,force=true) else !isfile(TMIfilegz) & !isfile(TMIfile) && google_download(url,pkgdatadir()) end @@ -206,9 +236,25 @@ function cartesianindex(file::String) elseif file[end-2:end] == "mat" matobj = matopen(file) - haskey(matobj,"it") ? it=convert(Vector{Integer},vec(read(matobj,"it"))) : it = convert(Vector{Integer},vec(read(matobj,"i"))) - haskey(matobj,"jt") ? jt=convert(Vector{Integer},vec(read(matobj,"jt"))) : jt=convert(Vector{Integer},vec(read(matobj,"j"))) - haskey(matobj,"kt") ? kt=convert(Vector{Integer},vec(read(matobj,"kt"))) : kt=convert(Vector{Integer},vec(read(matobj,"k"))) + if haskey(matobj,"it") + it=convert(Vector{Integer},vec(read(matobj,"it"))) + jt=convert(Vector{Integer},vec(read(matobj,"jt"))) + kt=convert(Vector{Integer},vec(read(matobj,"kt"))) + elseif haskey(matobj,"i") + it = convert(Vector{Integer},vec(read(matobj,"i"))) + jt = convert(Vector{Integer},vec(read(matobj,"j"))) + kt = convert(Vector{Integer},vec(read(matobj,"k"))) + elseif haskey(matobj,"grd") + #grd = read(matobj,"grd") + it = convert(Vector{Integer},vec(read(matobj,"grd")["it"])) + jt = convert(Vector{Integer},vec(read(matobj,"grd")["jt"])) + kt = convert(Vector{Integer},vec(read(matobj,"grd")["kt"])) + else + error("grid index key not found") + end + + #haskey(matobj,"jt") ? jt=convert(Vector{Integer},vec(read(matobj,"jt"))) : jt=convert(Vector{Integer},vec(read(matobj,"j"))) + #haskey(matobj,"kt") ? kt=convert(Vector{Integer},vec(read(matobj,"kt"))) : kt=convert(Vector{Integer},vec(read(matobj,"k"))) close(matobj) I = CartesianIndex.(it,jt,kt) end @@ -231,15 +277,19 @@ function gridprops(file) depth = convert(Vector{Float64},ncread(file,"depth")) elseif file[end-2:end] == "mat" - + matobj = matopen(file) - lon=convert(Vector{Float64},vec(read(matobj,"LON"))) - lat=convert(Vector{Float64},vec(read(matobj,"LAT"))) - depth=convert(Vector{Float64},vec(read(matobj,"DEPTH"))) + if haskey(matobj,"grd") + lon = convert(Vector{Float64},vec(read(matobj,"grd")["LON"])) + lat = convert(Vector{Float64},vec(read(matobj,"grd")["LAT"])) + depth = convert(Vector{Float64},vec(read(matobj,"grd")["DEPTH"])) + else + lon=convert(Vector{Float64},vec(read(matobj,"LON"))) + lat=convert(Vector{Float64},vec(read(matobj,"LAT"))) + depth=convert(Vector{Float64},vec(read(matobj,"DEPTH"))) + end close(matobj) - end - return lon,lat,depth end @@ -444,6 +494,9 @@ end """ function updatelinearindex(izyx,Izyx,R) Linear index translated from z,y,x to x,y,z accounting + + get Izyx Cartesian index stored from legacy MATLAB code + # Arguments - `izyx`: index of interest in z,y,x accounting - `Izyx`: wet Cartesian Index for z,y,x @@ -451,12 +504,7 @@ end # Output - `ixyz`: index of interest in x,y,z accounting """ -function updatelinearindex(izyx,Izyx,R) - # get Izyx Cartesian index stored from legacy MATLAB code - ixyz = R[Izyx[izyx]] - return ixyz -end - +updatelinearindex(izyx,Izyx,R) = R[Izyx[izyx]] """ function surfaceregion(TMIversion::String,region::String,γ::Grid)::BoundaryCondition diff --git a/src/deprecated.jl b/src/deprecated.jl new file mode 100644 index 00000000..4d4449c9 --- /dev/null +++ b/src/deprecated.jl @@ -0,0 +1,339 @@ + +# """ +# function adjustboundarycondition(b::NamedTuple{<:Any, NTuple{N1,BoundaryCondition{T}}},u::NamedTuple{<:Any, NTuple{N2,BoundaryCondition{T}}}) where N1, N2, T <: Real + +# adjust all boundary conditions b that are described in u +# """ +# function adjustboundarycondition(b::NamedTuple{<:Any, NTuple{N1,BoundaryCondition{T}}},u::NamedTuple{<:Any, NTuple{N2,BoundaryCondition{T}}}) where {N1, N2, T <: Real} + +# bnew = deepcopy(b) +# ukeys = keys(u) +# for ukey in keys(u) +# bnew[ukey].tracer[bnew[ukey].wet] += u[ukey].tracer[bnew[ukey].wet] +# end +# return bnew +# end + +# function unvec!(u::BoundaryCondition{T},uvec::Vector{T}) where T <: Real +# I = findall(u.wet) +# counter = 0 +# for i in I +# counter +=1 +# # doesn't work to pass u back +# #u.tracer[u.wet] = uvec +# u.tracer[i] = uvec[counter] +# end +# end +# function unvec!(u::NamedTuple,uvec::Vector) #where {N, T <: Real} + +# #counter = 0 +# for v in u +# # if v isa BoundaryCondition +# # wet = v.wet +# # elseif v isa Field +# # wet = v.γ.wet +# # end +# n = sum(wet(v)) +# v.tracer[wet(v)] = uvec[counter+1:counter+n] +# counter += n +# end +# end + + +# function zerotemplate(utemplate,uvec) +# tracer = zeros(wet(utemplate)) +# tracer[wet(utemplate)] = uvec +# return tracer +# end + +# """ +# function unvec(u,uvec) + +# Replace u with new u +# Undo the operations by vec(u) +# Needs to update u because attributes of +# u need to be known at runtime. +# """ +# function unvec(utemplate::Union{Field{T},BoundaryCondition{T}},uvec::Vector{T}) where T <: Real +# tracer = zerotemplate(utemplate,uvec) +# u = BoundaryCondition(tracer,utemplate.i,utemplate.j,utemplate.k,utemplate.dim,utemplate.dimval,wet(utemplate)) +# return u +# end +# function unvec(utemplate::BoundaryCondition{T},uvec::Vector{T}) where T <: Real +# tracer = zeros(utemplate.wet) +# tracer[utemplate.wet] = uvec +# u = BoundaryCondition(tracer,utemplate.i,utemplate.j,utemplate.k,utemplate.dim,utemplate.dimval,utemplate.wet) +# return u +# end +# function unvec(utemplate::Field{T},uvec::Vector{T}) where T <: Real +# tracer = zeros(utemplate.γ) +# tracer.tracer[utemplate.γ.wet] .= uvec +# #u = BoundaryCondition(tracer,utemplate.i,utemplate.j,utemplate.k,utemplate.dim,utemplate.dimval,utemplate.wet) +# return tracer +# end +# function unvec(utemplate::NamedTuple{<:Any,NTuple{N,BoundaryCondition{T}}},uvec::Vector{T}) where {N, T <: Real} +# counter = 0 +# vals = Vector{BoundaryCondition}(undef,N) +# for (ii,vv) in enumerate(utemplate) +# n = sum(vv.wet) +# vals[ii] = unvec(vv, uvec[counter+1:counter+n]) +# counter += n +# end +# u = (;zip(keys(utemplate), vals)...) +# return u +# end +# function unvec(utemplate::NamedTuple{<:Any,NTuple{N,Field{T}}},uvec::Vector{T}) where {N, T <: Real} +# counter = 0 +# vals = Vector{Field}(undef,N) +# for (ii,vv) in enumerate(utemplate) +# n = sum(vv.γ.wet) +# vals[ii] = unvec(vv, uvec[counter+1:counter+n]) +# counter += n +# end +# u = (;zip(keys(utemplate), vals)...) +# return u +# end +# function unvec(utemplate,uvec::Vector{T}) where {T <: Real} +# counter = 0 +# vals = Vector{Any}(undef,length(utemplate)) +# for (ii,vv) in enumerate(utemplate) +# if vv isa Field +# n = sum(vv.γ.wet) +# elseif vv isa BoundaryCondition +# n = sum(vv.wet) +# end +# vals[ii] = unvec(vv, uvec[counter+1:counter+n]) +# counter += n +# end +# u = (;zip(keys(utemplate), vals)...) +# return u +# end + + +""" + function vec2fld + Transfer a vector to a 3D field with accounting for ocean bathymetry +# Arguments +- `vector`: field in vector form (no land points) +- `I`: cartesian indices of ocean points +# Output +- `field`: field in 3d form including land points (NaN) +""" +function vec2fld(vector::Vector{T},I::Vector{CartesianIndex{3}}) where T<:Real + # choose NaN for now, zero better? nothing better? + fillvalue = zero(T)/zero(T) + + nx = maximum(I)[1] + ny = maximum(I)[2] + nz = maximum(I)[3] + + # faster instead to allocate as undef and then fill! ? + field = (NaN .* zero(T)) .* zeros(nx,ny,nz) + + # a comprehension + [field[I[n]]=vector[n] for n ∈ eachindex(I)] + return field +end + + +""" + function tracerinit(wet,vec,I) + initialize tracer field on TMI grid + perhaps better to have a tracer struct and constructor +# Arguments +- `wet`:: BitArray mask of ocean points +- `vec`:: vector of values at wet points +- `I`:: Cartesian Index for vector +# Output +- `field`:: 3d tracer field with NaN on dry points +""" +function tracerinit(vec,I,wet) + + # preallocate + T = eltype(vec) + field = Array{T}(undef,size(wet)) + fill!(field,zero(T)/zero(T)) + + #- a comprehension + [field[I[n]]=vec[n] for n ∈ eachindex(I)] + return field +end + +""" + function Γsfc + Γsfc anonymously calls surfacecontrol2field +""" +#Γsfc = surfacecontrol2field + +""" + function surfacecontrol2field!(c,u,γ) + Add surface control vector to existing 3D field +# Arguments +- `c`:: state field, 3d tracer field with NaN on dry points, modified by function +- `usfc`:: surface control vector +- `wet`::BitArray mask of ocean points +""" +function surfacecontrol2field!(c::Array{T,3},usfc::Vector{T},γ) where T<: Real + #c[:,:,1][wet[:,:,1]] .+= u # doesn't work +# [c[γ.I[ii][1],γ.I[ii][2],γ.I[ii][3]] += u[ii] for ii ∈ eachindex(γ.I) if γ.I[ii][3] == 1] + list = surfaceindex(γ.I) + [c[γ.I[ii]] += usfc[list[ii]] for ii ∈ eachindex(γ.I) if γ.I[ii][3] == 1] +end + +""" + function surfacecontrol2field!(c,u,γ) + Add surface control vector to tracer vector +# Arguments +- `c`:: state field, 3d tracer field with NaN on dry points, modified by function +- `u`:: surface control vector +- `wet`::BitArray mask of ocean points +""" +function surfacecontrol2field!(c::Vector{T},u::Vector{T},γ) where T<: Real + list = surfaceindex(γ.I) + [c[ii] += u[list[ii]] for ii ∈ eachindex(γ.I) if γ.I[ii][3] == 1] +end + +""" + function Γsfc! + Γsfc! anonymously calls surfacecontrol2field! +""" +Γsfc! = surfacecontrol2field! + +function field2obs(cvec,wis,γ) + # interpolate onto data points + N = length(wis) + sumwis = Vector{Float64}(undef,N) + list = vcat(1:length(γ.lon),1) + + # perhaps the most clever line in TMI.jl? + wetwrap = view(γ.wet,list,:,:) + + # some interpolation weights on land, oh no + # sum up all weights in ocean + [sumwis[i] = Interpolations.InterpGetindex(wetwrap)[wis[i]...] for i in eachindex(wis)] + + # reconstruct the observations + ỹ = Vector{Float64}(undef,N) + c̃ = tracerinit(γ.wet) + c̃[γ.wet] = cvec + replace!(c̃,NaN=>0.0) + cwrap = view(c̃,list,:,:) + + # divide by sum of all ocean weights so that this is still a true average + [ỹ[i] = Interpolations.InterpGetindex(cwrap)[wis[i]...]/sumwis[i] for i in eachindex(wis)] + return ỹ +end + +""" + function control2state(tracer2D,γ) + turn 2D surface field into 3D field with zeroes below surface +# Arguments +- `tracer2D`:: 2D surface tracer field +- `wet`::BitArray mask of ocean points +# Output +- `tracer3D`:: 3d tracer field with NaN on dry points +""" +function control2state(tracer2D::Matrix{T},wet) where T<: Real + # preallocate + tracer3D = Array{T}(undef,size(wet)) + + # set ocean to zero, land to NaN + # consider whether land should be nothing or missing + tracer3D[wet] .= zero(T) + tracer3D[.!wet] .= zero(T)/zero(T) + tracer3D[:,:,1] = tracer2D + return tracer3D +end + +""" + function surfacecontrol2field(usfc,γ.wet) + turn surface control vector into 3D field with zeroes below surface +# Arguments +- `usfc`:: surface control vector +- `wet`::BitArray mask of ocean points +# Output +- `tracer3D`:: 3d tracer field with NaN on dry points +""" +function surfacecontrol2field(usfc::Vector{T},wet) where T<: Real + # preallocate + tracer3D = Array{T}(undef,size(wet)) + + # set ocean to zero, land to NaN + # consider whether land should be nothing or missing + tracer3D[wet] .= zero(T) + tracer3D[.!wet] .= zero(T)/zero(T) + tracer3D[:,:,1][wet[:,:,1]] = usfc + return tracer3D +end + +# """ +# `function +(c::BoundaryCondition,d::BoundaryCondition)::BoundaryCondition` +# Define addition for Fields +# """ +# function Base.:+(c::T,d::T)::T where T <: Union{Source,Field,BoundaryCondition} + +# if c.wet != d.wet # check conformability +# error("BoundaryCondition's not conformable for addition") +# end +# array = zeros(c.wet) +# e = BoundaryCondition(array,c.i,c.j,c.k,c.dim,c.dimval,c.wet) + +# # a strange formulation to get +# # return e to be correct +# e.tracer[e.wet] += c.tracer[c.wet] +# e.tracer[e.wet] += d.tracer[d.wet] +# return e +# end + +# """ +# `function +(c::Field,d::Field)::Field` +# Define addition for Fields +# """ +# function +(c::Field{T},d::Field{T})::Field{T} where T <: Real +# # initialize output +# if c.γ.wet != d.γ.wet # check conformability +# error("Fields not conformable for addition") +# end + +# if !isequal(d.units,c.units) +# error("Units not consistent:",d.units," vs ",c.units) +# end + +# e = zeros(d.γ,d.name,d.longname,d.units) + +# # a strange formulation to get +# # return e to be correct +# e.tracer[e.γ.wet] += c.tracer[c.γ.wet] +# e.tracer[e.γ.wet] += d.tracer[d.γ.wet] +# return e +# end + +# """ +# `function -(c::Field,d::Field)::Field` +# Define addition for Fields +# """ +# function -(c::Field{T},d::Field{T})::Field{T} where T <: Real +# # initialize output +# if c.γ.wet !== d.γ.wet # check conformability +# error("Fields not conformable for addition") +# end +# e = zeros(d.γ) +# e.tracer[e.γ.wet] += c.tracer[c.γ.wet] +# e.tracer[e.γ.wet] -= d.tracer[d.γ.wet] +# return e +# end + +# """ +# `function *(C,d::BoundaryCondition)::BoundaryCondition` +# Define scalar or matrix multiplication for BoundaryCondition`s +# """ +# function Base.:*(C,d::BoundaryCondition{T})::BoundaryCondition{T} where T <: Real +# array = zeros(d.wet) +# e = BoundaryCondition(array,d.i,d.j,d.k,d.dim,d.dimval,d.wet) +# e.tracer[e.wet] += C*d.tracer[d.wet] +# return e +# end + +# """ +# `function *(c::Field,d::Field)::Field` +# Field by field multiplication is element-by-element. diff --git a/src/read_mat_nordic_201x115x46_B23.sh b/src/read_mat_nordic_201x115x46_B23.sh new file mode 100644 index 00000000..70dd066e --- /dev/null +++ b/src/read_mat_nordic_201x115x46_B23.sh @@ -0,0 +1,6 @@ +#!/bin/sh +# Extract data from Google Drive using your favorite method. MATLAB's webread may be an alternative method. Google Drive may ask for spam confirmation for files bigger than 40 MB. Sometimes throws ERROR: cannot verify docs.google.com's certificate, but still works. + +# Old Google ID = 1VUkucxO7aKVsVD8XpHUT9GCN-_CHz9aU +# Google ID = 1LddFEsGVCor1IOr4sYZ1S8SAVR9nDQ4R +wget --load-cookies /tmp/cookies.txt "https://docs.google.com/uc?export=download&confirm=$(wget --quiet --save-cookies /tmp/cookies.txt --keep-session-cookies --no-check-certificate 'https://docs.google.com/uc?export=download&id=1LddFEsGVCor1IOr4sYZ1S8SAVR9nDQ4R' -O- | sed -rn 's/.*confirm=([0-9A-Za-z_]+).*/\1\n/p')&id=1LddFEsGVCor1IOr4sYZ1S8SAVR9nDQ4R" -O TMI_nordic_201x115x46_B23.mat.gz && rm -rf /tmp/cookies.txt diff --git a/src/read_nc_nordic_201x115x46_B23.sh b/src/read_nc_nordic_201x115x46_B23.sh new file mode 100644 index 00000000..0df817dd --- /dev/null +++ b/src/read_nc_nordic_201x115x46_B23.sh @@ -0,0 +1,7 @@ +#!/bin/sh +# Extract data from Google Drive using your favorite method. MATLAB's webread may be an alternative method. Google Drive may ask for spam confirmation for files bigger than 40 MB. Sometimes throws ERROR: cannot verify docs.google.com's certificate, but still works. + +# Old Google ID = 1mAmJKS3xeNFjKLLTMBpgSS-aWZ1WHNm- +# Google ID = 1aLCoNkAMSujC-ImX6Xw-mMK5qQeWwZHd + +wget --load-cookies /tmp/cookies.txt "https://docs.google.com/uc?export=download&confirm=$(wget --quiet --save-cookies /tmp/cookies.txt --keep-session-cookies --no-check-certificate 'https://docs.google.com/uc?export=download&id=1aLCoNkAMSujC-ImX6Xw-mMK5qQeWwZHd' -O- | sed -rn 's/.*confirm=([0-9A-Za-z_]+).*/\1\n/p')&id=1aLCoNkAMSujC-ImX6Xw-mMK5qQeWwZHd" -O TMI_nordic_201x115x46_B23.nc && rm -rf /tmp/cookies.txt diff --git a/test/runtests.jl b/test/runtests.jl index 3fbc50c7..eb8eef91 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -31,7 +31,32 @@ using TMI O₂ = steadyinversion(Alu,bO₂,γ,q=qPO₄,r=-170.0) end - + + @testset "steadyinversion with Source type" begin + + yPO₄ = readfield(TMIfile,"PO₄",γ) + bPO₄ = getsurfaceboundary(yPO₄) + PO₄pre = steadyinversion(Alu,bPO₄,γ) + + # this line changes + qPO₄ = readsource(TMIfile,"qPO₄",γ) + b₀ = zerosurfaceboundary(γ) + PO₄ᴿ = steadyinversion(Alu,b₀,γ,q=qPO₄) + PO₄total = PO₄ᴿ + PO₄pre + PO₄direct = steadyinversion(Alu,bPO₄,γ,q=qPO₄) + + ## how big is the maximum difference? + # could replace with abs function + @test maximum(PO₄direct - PO₄total) < 0.1 + @test minimum(PO₄direct - PO₄total) > -0.1 + + ## oxygen distribution, just be sure it runs + yO₂ = readfield(TMIfile,"O₂",γ) + bO₂ = getsurfaceboundary(yO₂) + O₂ = steadyinversion(Alu,bO₂,γ,q=qPO₄,r=-170.0) + + end + ############################ ## trackpathways @testset "trackpathways" begin @@ -198,7 +223,7 @@ using TMI fg(x) = costfunction_point_obs(x,Alu,b,u,y,W⁻,wis,locs,Q⁻,γ) f(x) = fg(x)[1] J0 = f(uvec) - J̃₀,gJ₀ = fg(uvec) + J̃₀,∂J₀∂u = fg(uvec) fg!(F,G,x) = costfunction_point_obs!(F,G,x,Alu,b,u,y,W⁻,wis,locs,Q⁻,γ) ϵ = 1e-3 # size of finite perturbation @@ -208,21 +233,137 @@ using TMI ∇f_finite = (f(δu) - f(uvec))/ϵ println("∇f_finite=",∇f_finite) - fg!(J̃₀,gJ₀,(uvec+δu)./2) # J̃₀ is not overwritten - ∇f = gJ₀[ii] + ∂J₀∂u = 0.0 .* uvec + fg!(J̃₀,∂J₀∂u,(uvec+δu)./2) # J̃₀ is not overwritten + ∇f = ∂J₀∂u[ii] println("∇f=",∇f) # error less than 10 percent? println("Percent error=",100*abs(∇f - ∇f_finite)/abs(∇f + ∇f_finite)) @test abs(∇f - ∇f_finite)/abs(∇f + ∇f_finite) < 0.1 - iterations = 5 - out = sparsedatamap(uvec,Alu,b,u,y,W⁻,wis,locs,Q⁻,γ,iterations) + iters = 5 + out = sparsedatamap(uvec,Alu,b,u,y,W⁻,wis,locs,Q⁻,γ,iterations=iters) + # was cost function decreased? @test out.minimum < J̃₀ # reconstruct by hand to double-check. ũ = out.minimizer - J̃,gJ̃ = fg(ũ) + J̃,∂J̃∂ũ = fg(ũ) @test J̃ < J̃₀ + + + end + end + @testset "sourcemap" begin + + using Statistics, Interpolations + + for lscale in (false,true) + for surfacetoo in (false,true) + + yPO₄ = readfield(TMIfile,"PO₄",γ) + bPO₄ = getsurfaceboundary(yPO₄) + qPO₄ = readsource(TMIfile,"qPO₄",γ) + q₀ = 1e-2*onesource(γ) + + N = 20 + y, W⁻, ctrue, ytrue, locs, wis = synthetic_observations(TMIversion,"PO₄",γ,N) + + #u = (; source = zerosource(γ)) + if surfacetoo + u = (; surface = zerosurfaceboundary(γ), source = zerosource(γ)) + else + u = (; source = zerosource(γ,logscale=lscale)) + end + + b = (; surface = bPO₄) # surface boundary condition + + PO₄true = steadyinversion(Alu,b,γ,q=qPO₄) + PO₄₀ = steadyinversion(Alu,b,γ,q=q₀) + uvec = vec(u) + + σq = 1.0 + Q⁻ = 1.0/(σq^2) # how well is q (source) known? + fg(x) = costfunction_point_obs(x,Alu,b,u,y,W⁻,wis,locs,Q⁻,γ,q=q₀) + f(x) = fg(x)[1] + J0 = f(uvec) + J̃₀,∂J₀∂u = fg(uvec) + fg!(F,G,x) = costfunction_point_obs!(F,G,x,Alu,b,u,y,W⁻,wis,locs,Q⁻,γ,q₀=q₀) + + ϵ = 1e-3 # size of finite perturbation + ii = rand(1:sum(γ.wet[:,:,1])) + println("gradient check location=",ii) + δu = copy(uvec); δu[ii] += ϵ + ∇f_finite = (f(δu) - f(uvec))/ϵ + println("∇f_finite=",∇f_finite) + + fg!(J̃₀,∂J₀∂u,(uvec+δu)./2) # J̃₀ is not overwritten + ∇f = ∂J₀∂u[ii] + println("∇f=",∇f) + + # error less than 10 percent? + println("Percent error=",100*abs(∇f - ∇f_finite)/abs(∇f + ∇f_finite)) + @test abs(∇f - ∇f_finite)/abs(∇f + ∇f_finite) < 0.1 + iters = 5 + out = sparsedatamap(uvec,Alu,b,u,y,W⁻,wis,locs,Q⁻,γ,q=qPO₄,r=1.0,iterations=iters) + # was cost function decreased? + @test out.minimum < J̃₀ + # reconstruct by hand to double-check. + ũ = out.minimizer + J̃,∂J̃∂ũ = fg(ũ) + @test J̃ < J̃₀ + #@test J̃ < 3N # too strict if first guess is bad + + if lscale + #b̃ = adjustboundarycondition(b₀,unvec(u,ũ)) # combine b₀, u + q̃ = TMI.adjustsource(q₀,unvec(u,ũ)) + @test minimum(q̃) ≥ 0 + end + end end end + + @testset "unvec" begin + u = (;surface = onesurfaceboundary(γ), + west = 2 * onewestboundary(γ), + east = 3 *oneeastboundary(γ)) + + vecu = vec(u) + + u2 = (;surface = zerosurfaceboundary(γ), + west = 2 * zerowestboundary(γ), + east = 3 * zeroeastboundary(γ)) + + u3 = unvec(u2,vecu) + @test u.surface.tracer[10,10] == u3.surface.tracer[10,10] + @test u.east.tracer[10,10] == u3.east.tracer[10,10] + @test u.west.tracer[10,10] == u3.west.tracer[10,10] + @test u.surface.tracer[10,10] == 1.0 + @test u.west.tracer[10,10] == 2.0 + @test u.east.tracer[10,10] == 3.0 + + unvec!(u2,vecu) + @test u.surface.tracer[10,10] == u2.surface.tracer[10,10] + @test u.east.tracer[10,10] == u2.east.tracer[10,10] + @test u.west.tracer[10,10] == u2.west.tracer[10,10] + @test u.surface.tracer[10,10] == 1.0 + @test u.west.tracer[10,10] == 2.0 + @test u.east.tracer[10,10] == 3.0 + + end + + @testset "regional" begin + TMIversion = "nordic_201x115x46_B23" + A, Alu, γ, TMIfile, L, B = config_from_mat(TMIversion) + +# shellscript = TMI.pkgsrcdir("read_nc_nordic_lowresolution.sh") +# run(`sh $shellscript`) + +# # NordicSeas_30-Dec-2020_low.mat +# TMIversion = "NordicSeas_30-Dec-2020_low"; +# TMIfile = TMIversion*".nc" +# !ispath(TMI.pkgdatadir()) && mkpath(TMI.pkgdatadir()) +# mv(joinpath(pwd(),TMIfile),TMI.pkgdatadir(TMIfile),force=true) + + end end