From 9606e6ecdcfa48580c0b82d0a7ba07fe70cc6f20 Mon Sep 17 00:00:00 2001 From: G Jake Gebbie Date: Mon, 27 Mar 2023 16:43:49 -0400 Subject: [PATCH 01/39] lost: minor changes --- .gitignore | 1 + src/TMI.jl | 2 +- 2 files changed, 2 insertions(+), 1 deletion(-) 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/src/TMI.jl b/src/TMI.jl index 931dd741..5f55d0ad 100644 --- a/src/TMI.jl +++ b/src/TMI.jl @@ -2224,7 +2224,7 @@ end """ function adjustboundarycondition!(b::NamedTuple{<:Any, NTuple{N1,BoundaryCondition{T}}},u::NamedTuple{<:Any, NTuple{N2,BoundaryCondition{T}}}) where {N1, N2, T <: Real} - ukeys = keys(u) + #ukeys = keys(u) for ukey in keys(u) b[ukey].tracer[b[ukey].wet] += u[ukey].tracer[b[ukey].wet] end From 52269779af577d38328ebafad15e3d805c2e155f Mon Sep 17 00:00:00 2001 From: G Jake Gebbie Date: Tue, 28 Mar 2023 21:53:35 -0400 Subject: [PATCH 02/39] shorten `updatelinearindex` --- src/TMI.jl | 11 ++--------- src/TMIconfig.jl | 10 ++++------ test/runtests.jl | 51 ++++++++++++++++++++++++++++++++++++++++++++++++ 3 files changed, 57 insertions(+), 15 deletions(-) diff --git a/src/TMI.jl b/src/TMI.jl index 5f55d0ad..3eaa5531 100644 --- a/src/TMI.jl +++ b/src/TMI.jl @@ -487,8 +487,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,7 +504,7 @@ 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. @@ -1915,8 +1915,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,12 +1927,7 @@ 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 diff --git a/src/TMIconfig.jl b/src/TMIconfig.jl index 085986e1..4efc8f7d 100644 --- a/src/TMIconfig.jl +++ b/src/TMIconfig.jl @@ -444,6 +444,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 +454,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/test/runtests.jl b/test/runtests.jl index 3fbc50c7..4c4a9779 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -225,4 +225,55 @@ using TMI @test J̃ < J̃₀ end end + @testset "sourcemap" begin + + using Statistics, Interpolations + + yPO₄ = readfield(TMIfile,"PO₄",γ) + bPO₄ = getsurfaceboundary(yPO₄) + + #PO₄pre = steadyinversion(Alu,bPO₄,γ) + qPO₄ = readfield(TMIfile,"qPO₄",γ) + #b₀ = bPO₄ #zerosurfaceboundary(γ) + #PO₄ᴿ = steadyinversion(Alu,b₀,γ,q=qPO₄) + + N = 20 + y, W⁻, ctrue, ytrue, locs, wis = synthetic_observations(TMIversion,"PO₄",γ,N) + + #u = (;surface = zerosurfaceboundary(γ)) + u = (; source = zeros(γ)) + b = (;surface = bPO₄) + + uvec = vec(u) + σb = 5.0 + Q⁻ = 1.0/(σb^2) + 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) + fg!(F,G,x) = costfunction_point_obs!(F,G,x,Alu,b,u,y,W⁻,wis,locs,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̃₀,gJ₀,(uvec+δu)./2) # J̃₀ is not overwritten + ∇f = gJ₀[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) + # was cost function decreased? + @test out.minimum < J̃₀ + # reconstruct by hand to double-check. + ũ = out.minimizer + J̃,gJ̃ = fg(ũ) + @test J̃ < J̃₀ + end end From 896eaf87c5d789a9dfa42a529e713084581cef01 Mon Sep 17 00:00:00 2001 From: G Jake Gebbie Date: Fri, 31 Mar 2023 16:50:21 -0400 Subject: [PATCH 03/39] wip: adjust- boundary and sources --- Manifest.toml | 26 +++++++- Project.toml | 1 + src/TMI.jl | 156 ++++++++++++++++++++++++----------------------- test/runtests.jl | 9 ++- 4 files changed, 113 insertions(+), 79 deletions(-) diff --git a/Manifest.toml b/Manifest.toml index a6c0344a..838c0f4f 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 = "9a28a4abf9f655c03a8cb01bab4e21c8a88412f5" [[deps.ANSIColoredPrinters]] git-tree-sha1 = "574baf8110975760d391c710b6341da1afa48d8c" @@ -77,6 +77,12 @@ 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" @@ -357,6 +363,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 +427,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" @@ -616,6 +634,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" diff --git a/Project.toml b/Project.toml index 7b4dc12b..10a2c2e8 100644 --- a/Project.toml +++ b/Project.toml @@ -23,6 +23,7 @@ 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" diff --git a/src/TMI.jl b/src/TMI.jl index 3eaa5531..2644fb7b 100644 --- a/src/TMI.jl +++ b/src/TMI.jl @@ -45,10 +45,12 @@ export config, config_from_mat, config_from_nc, surface_oxygensaturation, oxygen, location_obs, getsurfaceboundary, zerosurfaceboundary, onesurfaceboundary, setboundarycondition!, - adjustboundarycondition, adjustboundarycondition!, + #adjustboundarycondition, + adjustboundarycondition!, gsetboundarycondition, setsource!, zeros, one, oneunit, ones, maximum, minimum, (+), (-), (*), dot, - Grid, Field, BoundaryCondition, vec, unvec!, unvec + Grid, Field, BoundaryCondition, vec, unvec!, unvec, + adjustsource! import Base: zeros, one, oneunit, ones, maximum, minimum, (\) import Base: (+), (-), (*), (/), vec @@ -979,7 +981,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 +1005,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 +1038,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 +1060,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 +1077,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 @@ -1419,15 +1421,22 @@ 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} +#function vec(u::NamedTuple{<:Any,NTuple{N,BoundaryCondition{T}}}) where {N, T <: Real} +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 +vec(u::Field) = u.tracer[u.γ.wet] +vec(u::BoundaryCondition) = u.tracer[u.wet] + """ function unvec!(u,uvec) @@ -1453,7 +1462,6 @@ end 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) @@ -1464,6 +1472,17 @@ function unvec(utemplate::NamedTuple{<:Any,NTuple{N,BoundaryCondition{T}}},uvec: 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!(u,uvec) @@ -1494,12 +1513,17 @@ end 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 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 """ Surface oxygen saturation value and fraction of saturation value in field @@ -1584,13 +1608,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,7 +1818,7 @@ 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 @@ -1815,11 +1834,15 @@ 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} +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) +#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) @@ -1827,9 +1850,7 @@ function costfunction_point_obs(uvec::Vector{T},Alu,b₀::Union{BoundaryConditio u = unvec(u₀,uvec) b = adjustboundarycondition(b₀,u) # combine b₀, u - - c = steadyinversion(Alu,b,γ) # gives the misfit - + c = steadyinversion(Alu,b,γ,q,r) # gives the misfit # observe at right spots ỹ = observe(c,wis,γ) n = ỹ - y @@ -1982,7 +2003,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) @@ -2181,20 +2201,14 @@ end adjust the (one) boundary condition problem: passes back a mutated b """ -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 +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{T},u::BoundaryCondition{T}) where T <: Real - - adjust the (one) boundary condition -""" -function adjustboundarycondition(b::BoundaryCondition{T},u::BoundaryCondition{T}) where T <: Real - - bnew = b + u - return bnew - +function adjustsource!(b::Field,u::Field) + # write it out so b changes when returned + b.tracer[b.γ.wet] += u.tracer[u.γ.wet] end """ @@ -2215,28 +2229,42 @@ end 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 adjustboundarycondition!(b::NamedTuple{<:Any, NTuple{N1,BoundaryCondition{T}}},u::NamedTuple{<:Any, NTuple{N2,BoundaryCondition{T}}}) where {N1, N2, T <: Real} +function adjustboundarycondition!(b::NamedTuple,u::NamedTuple) #where {N1, N2, T <: Real} + for bkey in keys(b) + if haskey(u,bkey) + adjustboundarycondition!(b[bkey],u[bkey]) + #b[ukey].tracer[b[ukey].wet] += u[ukey].tracer[b[ukey].wet] + else + error("adjustboundarycondition: u doesn't have bkey ",bkey) + end + end +end +function adjustboundarycondition!(b::BoundaryCondition,u::NamedTuple) #where {N1, N2, T <: Real} + bkey = :surface # if not explicit, assume a surface boundary condition + if haskey(u,bkey) + adjustboundarycondition!(b,u[bkey]) + else + error("adjustboundarycondition: u doesn't have bkey ",bkey) 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 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::Field,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 - return bnew end """ @@ -2393,12 +2421,6 @@ function field2obs(cvec,wis,γ) return ỹ end -""" - function E - E anonymously calls field2obs -""" -E = field2obs - """ function tracerinit(wet,vec,I) initialize tracer field on TMI grid @@ -2447,22 +2469,6 @@ function vec2fld(vector::Vector{T},I::Vector{CartesianIndex{3}}) where T<:Real 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/test/runtests.jl b/test/runtests.jl index 4c4a9779..07f8ae0e 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -225,7 +225,7 @@ using TMI @test J̃ < J̃₀ end end - @testset "sourcemap" begin + @testset "adjustsource" begin using Statistics, Interpolations @@ -241,10 +241,13 @@ using TMI y, W⁻, ctrue, ytrue, locs, wis = synthetic_observations(TMIversion,"PO₄",γ,N) #u = (;surface = zerosurfaceboundary(γ)) - u = (; source = zeros(γ)) - b = (;surface = bPO₄) + u = (; source = ones(γ)) + b = (; surface = bPO₄) + b = bPO₄ # assume a surface boundary condition + PO₄ = steadyinversion(Alu,b,γ,q=qPO₄) uvec = vec(u) + σb = 5.0 Q⁻ = 1.0/(σb^2) fg(x) = costfunction_point_obs(x,Alu,b,u,y,W⁻,wis,locs,Q⁻,γ) From 1625f6c066953ced09b2310d5cfe52d1b9b3dc8a Mon Sep 17 00:00:00 2001 From: G Jake Gebbie Date: Fri, 31 Mar 2023 17:08:16 -0400 Subject: [PATCH 04/39] sparsedatamap fails --- src/TMI.jl | 54 +++++++++++++++++++++++++++++++++++------------------- 1 file changed, 35 insertions(+), 19 deletions(-) diff --git a/src/TMI.jl b/src/TMI.jl index 2644fb7b..f2606e7c 100644 --- a/src/TMI.jl +++ b/src/TMI.jl @@ -1934,7 +1934,7 @@ end # Output - `c`::Field, steady-state tracer distribution """ -function steadyinversion(Alu,b::BoundaryCondition{T},γ::Grid;q=nothing,r=1.0)::Field{T} where T <: Real +function steadyinversion(Alu,b::Union{NamedTuple,BoundaryCondition{T}},γ::Grid;q=nothing,r=1.0)::Field{T} where T <: Real # preallocate Field for equation constraints d = zeros(γ) @@ -2196,34 +2196,31 @@ 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::BoundaryCondition{T},u::BoundaryCondition{T}) where T <: Real adjust the (one) boundary condition - problem: passes back a mutated b """ -function adjustboundarycondition!(b::BoundaryCondition,u::BoundaryCondition) - # write it out so b changes when returned - b.tracer[b.wet] += u.tracer[u.wet] -end - -function adjustsource!(b::Field,u::Field) - # write it out so b changes when returned - b.tracer[b.γ.wet] += u.tracer[u.γ.wet] +adjustboundarycondition(b::BoundaryCondition,u::BoundaryCondition) = b + u +function adjustboundarycondition(b::Union{BoundaryCondition,NamedTuple},u::NamedTuple) + bnew = deepcopy(b) + adjustboundarycondition!(bnew,u) + # bnew[ukey] += u[ukey] + # #bnew[ukey].tracer[bnew[ukey].wet] += u[ukey].tracer[bnew[ukey].wet] + # end + # end + return bnew end """ - function gadjustboundarycondition!(b::BoundaryCondition{T},u::BoundaryCondition{T}) where T <: Real + function adjustboundarycondition!(b::BoundaryCondition{T},u::BoundaryCondition{T}) where T <: Real adjust the (one) boundary condition - Just copy the variable. - Keep this function so that calling functions can look alike. - Could probably combine with lower function, use Union type + problem: passes back a mutated b """ -function gadjustboundarycondition(gb::BoundaryCondition{T},u::BoundaryCondition{T}) where T <: Real - gu = gb - return gu +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::NamedTuple{<:Any, NTuple{N1,BoundaryCondition{T}}},u::NamedTuple{<:Any, NTuple{N2,BoundaryCondition{T}}}) where N1, N2, T <: Real @@ -2249,6 +2246,25 @@ function adjustboundarycondition!(b::BoundaryCondition,u::NamedTuple) #where {N1 end end +function adjustsource!(b::Field,u::Field) + # write it out so b changes when returned + b.tracer[b.γ.wet] += u.tracer[u.γ.wet] +end + +""" + function gadjustboundarycondition!(b::BoundaryCondition{T},u::BoundaryCondition{T}) where T <: Real + + adjust the (one) boundary condition + Just copy the variable. + 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 +end + + function adjustsource!(q::NamedTuple,u::NamedTuple) #where {N1, N2, T <: Real} for qkey in keys(q) if haskey(u,qkey) From 6cff74802bc1edd94694ea058d6d40b12628bb9e Mon Sep 17 00:00:00 2001 From: G Jake Gebbie Date: Sat, 1 Apr 2023 10:51:51 -0400 Subject: [PATCH 05/39] `sparsedatamap` passes tests again --- src/TMI.jl | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/TMI.jl b/src/TMI.jl index f2606e7c..a847a185 100644 --- a/src/TMI.jl +++ b/src/TMI.jl @@ -1850,7 +1850,7 @@ function costfunction_point_obs(uvec::Vector,Alu,b₀::Union{BoundaryCondition,N u = unvec(u₀,uvec) b = adjustboundarycondition(b₀,u) # combine b₀, u - c = steadyinversion(Alu,b,γ,q,r) # gives the misfit + c = steadyinversion(Alu,b,γ,q=q,r=r) # gives the misfit # observe at right spots ỹ = observe(c,wis,γ) n = ỹ - y @@ -1926,7 +1926,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 @@ -1934,7 +1934,7 @@ end # Output - `c`::Field, steady-state tracer distribution """ -function steadyinversion(Alu,b::Union{NamedTuple,BoundaryCondition{T}},γ::Grid;q=nothing,r=1.0)::Field{T} where T <: Real +function steadyinversion(Alu,b::BoundaryCondition{T},γ::Grid;q=nothing,r=1.0)::Field{T} where T <: Real # preallocate Field for equation constraints d = zeros(γ) From 28f3f2541ec061950676888b3ae7953ffc84edce Mon Sep 17 00:00:00 2001 From: G Jake Gebbie Date: Sun, 2 Apr 2023 18:41:21 -0400 Subject: [PATCH 06/39] reorganize `adjustboundarycondition` --- src/TMI.jl | 135 +++++++++++++++++++++++------------------------ test/runtests.jl | 9 ++-- 2 files changed, 71 insertions(+), 73 deletions(-) diff --git a/src/TMI.jl b/src/TMI.jl index a847a185..f2879c0c 100644 --- a/src/TMI.jl +++ b/src/TMI.jl @@ -1409,18 +1409,9 @@ end for use with Optim.jl. An in-place version of this function would be handy. """ -function vec(u::BoundaryCondition{T}) where T <: Real - uvec = u.tracer[u.wet] - return uvec -end - -""" - function vec(u) +vec(u::Field) = u.tracer[u.γ.wet] +vec(u::BoundaryCondition) = u.tracer[u.wet] - Turn a collection of controls into a vector - 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} function vec(u::NamedTuple) @@ -1434,9 +1425,6 @@ function vec(u::NamedTuple) return uvec end -vec(u::Field) = u.tracer[u.γ.wet] -vec(u::BoundaryCondition) = u.tracer[u.wet] - """ function unvec!(u,uvec) @@ -1444,10 +1432,16 @@ vec(u::BoundaryCondition) = u.tracer[u.wet] 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} +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(v.wet) v.tracer[v.wet] = uvec[counter+1:counter+n] counter += n @@ -1842,11 +1836,9 @@ end - `gJ`: derivative of cost function wrt to controls """ 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) -#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 + # control penalty Jcontrol = uvec'*(Q⁻*uvec) - guvec = 2*(Q⁻*uvec) u = unvec(u₀,uvec) b = adjustboundarycondition(b₀,u) # combine b₀, u @@ -1858,12 +1850,22 @@ function costfunction_point_obs(uvec::Vector,Alu,b₀::Union{BoundaryCondition,N Jdata = n ⋅ (Wⁱ * n) # dot product J = Jdata + Jcontrol + # start adjoint model + + # initialize adjoint control variable with zero + gu = unvec(u₀,0 .* uvec) + gn = 2Wⁱ * n gỹ = gn gc = gobserve(gỹ,c,locs) gb = gsteadyinversion(gc, Alu, b, γ) - gu = gadjustboundarycondition(gb,u) + #gu = gadjustboundarycondition(gb,u) + gadjustboundarycondition!(gu,gb) + + # control penalty gradient + guvec = 2*(Q⁻*uvec) + guvec += vec(gu) return J, guvec @@ -1889,10 +1891,11 @@ 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,uvec::Vector,Alu,b::Union{BoundaryCondition,NamedTuple},u::Union{BoundaryCondition,NamedTuple},y::Vector,Wⁱ::Diagonal,wis,locs,Q⁻,γ::Grid;q=nothing,r=1.0) #where {N1, N2, T <: Real} +#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} - u = unvec(u₀,uvec) - b = adjustboundarycondition(b₀,u) # combine b₀, u + unvec!(u,uvec) + adjustboundarycondition!(b,u) # combine b₀, u c = steadyinversion(Alu,b,γ) # gives the misfit # observe at right spots @@ -2196,59 +2199,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 + adjust all boundary conditions b that are described in u """ adjustboundarycondition(b::BoundaryCondition,u::BoundaryCondition) = b + u -function adjustboundarycondition(b::Union{BoundaryCondition,NamedTuple},u::NamedTuple) +function adjustboundarycondition(b::Union{BoundaryCondition,NamedTuple},u::Union{BoundaryCondition,NamedTuple}) bnew = deepcopy(b) adjustboundarycondition!(bnew,u) - # bnew[ukey] += u[ukey] - # #bnew[ukey].tracer[bnew[ukey].wet] += u[ukey].tracer[bnew[ukey].wet] - # end - # end 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 - problem: passes back a mutated b + adjust all boundary conditions b that are described in u + + 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::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} -function adjustboundarycondition!(b::NamedTuple,u::NamedTuple) #where {N1, N2, T <: Real} - for bkey in keys(b) - if haskey(u,bkey) - adjustboundarycondition!(b[bkey],u[bkey]) - #b[ukey].tracer[b[ukey].wet] += u[ukey].tracer[b[ukey].wet] - else - error("adjustboundarycondition: u doesn't have bkey ",bkey) - end - end -end function adjustboundarycondition!(b::BoundaryCondition,u::NamedTuple) #where {N1, N2, T <: Real} - bkey = :surface # if not explicit, assume a surface boundary condition + # if not explicit, assume a surface boundary condition + bkey = :surface # if haskey(u,bkey) adjustboundarycondition!(b,u[bkey]) - else - error("adjustboundarycondition: u doesn't have bkey ",bkey) + end +end +function adjustboundarycondition!(b::NamedTuple,u::NamedTuple) + for bkey in keys(b) + adjustboundarycondition!(b[bkey],u) end end -function adjustsource!(b::Field,u::Field) - # write it out so b changes when returned - b.tracer[b.γ.wet] += u.tracer[u.γ.wet] +# Seems not to be general because gu overwritten. +function gadjustboundarycondition(gb::BoundaryCondition{T},u::BoundaryCondition{T}) where T <: Real + gu = gb + return gu end """ @@ -2259,11 +2249,31 @@ 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 +function gadjustboundarycondition!(gu::BoundaryCondition,gb::BoundaryCondition) + gu.tracer[gu.wet] += gb.tracer[gb.wet] +end +function gadjustboundarycondition!(gu::NamedTuple,gb::BoundaryCondition) + bkey = :surface + if haskey(gu,bkey) + gadjustboundarycondition!(gu[bkey],gb) + end +end +function gadjustboundarycondition!(gu::NamedTuple,gb::NamedTuple) + for bkey in keys(gb) + gadjustboundarycondition!(gu,gb[bkey]) + end +end + +#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{<:Any, NTuple{N2,BoundaryCondition{T}}}) where {N1, N2, T <: Real} + gu = gb[keys(u)] # grab the parts of the named tuple corresponding to u return gu end +function adjustsource!(b::Field,u::Field) + # write it out so b changes when returned + b.tracer[b.γ.wet] += u.tracer[u.γ.wet] +end function adjustsource!(q::NamedTuple,u::NamedTuple) #where {N1, N2, T <: Real} for qkey in keys(q) @@ -2283,19 +2293,6 @@ function adjustsource!(q::Field,u::NamedTuple) #where {N1, N2, T <: Real} end 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} - gu = gb[keys(u)] # grab the parts of the named tuple corresponding to u - return gu -end - """ function section View latitude-depth slice of field diff --git a/test/runtests.jl b/test/runtests.jl index 07f8ae0e..b2bdcd03 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -241,16 +241,17 @@ using TMI y, W⁻, ctrue, ytrue, locs, wis = synthetic_observations(TMIversion,"PO₄",γ,N) #u = (;surface = zerosurfaceboundary(γ)) - u = (; source = ones(γ)) - b = (; surface = bPO₄) + u = (; source = zeros(γ)) + #b = (; surface = bPO₄) b = bPO₄ # assume a surface boundary condition PO₄ = steadyinversion(Alu,b,γ,q=qPO₄) uvec = vec(u) σb = 5.0 - Q⁻ = 1.0/(σb^2) - fg(x) = costfunction_point_obs(x,Alu,b,u,y,W⁻,wis,locs,Q⁻,γ) + σq = 0.2 + 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=qPO₄) f(x) = fg(x)[1] J0 = f(uvec) J̃₀,gJ₀ = fg(uvec) From 4a5df032107202e08a61b9325519a676f2d67642 Mon Sep 17 00:00:00 2001 From: G Jake Gebbie Date: Sun, 2 Apr 2023 20:32:36 -0400 Subject: [PATCH 07/39] core numerics in `costfunction_point_obs!` --- src/TMI.jl | 147 ++++++++++++++++++++++------------------------- test/runtests.jl | 18 +++--- 2 files changed, 80 insertions(+), 85 deletions(-) diff --git a/src/TMI.jl b/src/TMI.jl index f2879c0c..6ca4d27f 100644 --- a/src/TMI.jl +++ b/src/TMI.jl @@ -419,7 +419,8 @@ 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,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⁻,γ) @@ -1426,35 +1427,25 @@ function vec(u::NamedTuple) end """ - function unvec!(u,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!(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(v.wet) - v.tracer[v.wet] = uvec[counter+1:counter+n] - counter += n - 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,uvec) - - 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) @@ -1485,8 +1476,22 @@ 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 +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.tracer[wet] = uvec[counter+1:counter+n] + counter += n + end +end +function unvec!(u::BoundaryCondition{T},uvec::Vector{T}) where T <: Real I = findall(u.wet) counter = 0 for i in I @@ -1495,28 +1500,6 @@ function unvec!(u::BoundaryCondition{T},uvec::Vector{T}) where T <: Real #u.tracer[u.wet] = uvec u.tracer[i] = uvec[counter] 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 -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 """ @@ -1835,40 +1818,44 @@ end - `J`: cost function of sum of squared misfits - `gJ`: derivative of cost function wrt to controls """ -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) +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) - # control penalty - Jcontrol = uvec'*(Q⁻*uvec) + 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 - u = unvec(u₀,uvec) - b = adjustboundarycondition(b₀,u) # combine b₀, u - c = steadyinversion(Alu,b,γ,q=q,r=r) # gives the misfit - # observe at right spots - ỹ = observe(c,wis,γ) - n = ỹ - y - - Jdata = n ⋅ (Wⁱ * n) # dot product - J = Jdata + Jcontrol + # # control penalty + # Jcontrol = uvec'*(Q⁻*uvec) - # start adjoint model + # u = unvec(u₀,uvec) + # b = adjustboundarycondition(b₀,u) # combine b₀, u + # c = steadyinversion(Alu,b,γ,q=q,r=r) # gives the misfit + # # observe at right spots + # ỹ = observe(c,wis,γ) + # n = ỹ - y + + # Jdata = n ⋅ (Wⁱ * n) # dot product + # J = Jdata + Jcontrol - # initialize adjoint control variable with zero - gu = unvec(u₀,0 .* uvec) + # ## start adjoint model + # # initialize adjoint control variable with zero + # gu = unvec(u₀,0 .* uvec) - gn = 2Wⁱ * n - gỹ = gn + # gn = 2Wⁱ * n + # gỹ = gn - gc = gobserve(gỹ,c,locs) - gb = gsteadyinversion(gc, Alu, b, γ) - #gu = gadjustboundarycondition(gb,u) - gadjustboundarycondition!(gu,gb) + # gc = gobserve(gỹ,c,locs) + # gb = gsteadyinversion(gc, Alu, b, γ) + # #gu = gadjustboundarycondition(gb,u) + # gadjustboundarycondition!(gu,gb) - # control penalty gradient - guvec = 2*(Q⁻*uvec) + # # control penalty gradient + # guvec = 2*(Q⁻*uvec) - guvec += vec(gu) - - return J, guvec + # guvec += vec(gu) end """ @@ -1891,12 +1878,13 @@ end - `Q⁻`: weights for control vector - `γ`: grid """ -function costfunction_point_obs!(J,guvec,uvec::Vector,Alu,b::Union{BoundaryCondition,NamedTuple},u::Union{BoundaryCondition,NamedTuple},y::Vector,Wⁱ::Diagonal,wis,locs,Q⁻,γ::Grid;q=nothing,r=1.0) #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) #where {N1, N2, T <: Real} #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} - unvec!(u,uvec) - adjustboundarycondition!(b,u) # combine b₀, u - c = steadyinversion(Alu,b,γ) # gives the misfit + #unvec!(u,uvec) # causes issues with Optim.jl + u = unvec(u₀,uvec) + b = adjustboundarycondition(b₀,u) # combine b₀, u + c = steadyinversion(Alu,b,γ,q=q,r=r) # observe at right spots ỹ = observe(c,wis,γ) @@ -1904,6 +1892,7 @@ function costfunction_point_obs!(J,guvec,uvec::Vector,Alu,b::Union{BoundaryCondi if guvec != nothing # adjoint equations + #guvec = 0.0.*uvec # re-initialize gtmp = 2*(Q⁻*uvec) gn = 2Wⁱ * n gỹ = gn @@ -1914,9 +1903,13 @@ function costfunction_point_obs!(J,guvec,uvec::Vector,Alu,b::Union{BoundaryCondi for (ii,vv) in enumerate(gtmp) guvec[ii] = vv end + println(maximum(guvec)) + println(minimum(guvec)) end if J !=nothing + # reinitialize for Optim.jl + #J = 0.0 # control penalty and gradient Jcontrol = uvec'*(Q⁻*uvec) Jdata = n ⋅ (Wⁱ * n) # dot product diff --git a/test/runtests.jl b/test/runtests.jl index b2bdcd03..4de68583 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -198,7 +198,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,8 +208,9 @@ 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? @@ -217,6 +218,7 @@ using TMI @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) + # was cost function decreased? @test out.minimum < J̃₀ # reconstruct by hand to double-check. @@ -254,8 +256,8 @@ using TMI fg(x) = costfunction_point_obs(x,Alu,b,u,y,W⁻,wis,locs,Q⁻,γ,q=qPO₄) f(x) = fg(x)[1] J0 = f(uvec) - J̃₀,gJ₀ = fg(uvec) - fg!(F,G,x) = costfunction_point_obs!(F,G,x,Alu,b,u,y,W⁻,wis,locs,Q⁻,γ) + J̃₀,∂J₀∂u = fg(uvec) + fg!(F,G,x) = costfunction_point_obs!(F,G,x,Alu,b,u,y,W⁻,wis,locs,Q⁻,γ,q=qPO₄) ϵ = 1e-3 # size of finite perturbation ii = rand(1:sum(γ.wet[:,:,1])) @@ -264,8 +266,8 @@ 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] + fg!(J̃₀,∂J₀∂u,(uvec+δu)./2) # J̃₀ is not overwritten + ∇f = ∂J₀∂u[ii] println("∇f=",∇f) # error less than 10 percent? @@ -277,7 +279,7 @@ using TMI @test out.minimum < J̃₀ # reconstruct by hand to double-check. ũ = out.minimizer - J̃,gJ̃ = fg(ũ) + J̃,∂J̃∂ũ = fg(ũ) @test J̃ < J̃₀ end end From fea74cbccec8095ec44dd19af501d7965ffa17ce Mon Sep 17 00:00:00 2001 From: G Jake Gebbie Date: Sun, 2 Apr 2023 20:36:24 -0400 Subject: [PATCH 08/39] `sparsedatamap` cleanup --- src/TMI.jl | 46 ++++++---------------------------------------- 1 file changed, 6 insertions(+), 40 deletions(-) diff --git a/src/TMI.jl b/src/TMI.jl index 6ca4d27f..2c4c5ef0 100644 --- a/src/TMI.jl +++ b/src/TMI.jl @@ -1797,9 +1797,9 @@ 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;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 @@ -1826,36 +1826,6 @@ function costfunction_point_obs(uvec::Vector,Alu,b::Union{BoundaryCondition,Name J = costfunction_point_obs!(J,guvec,uvec,Alu,b,u,y,Wⁱ,wis,locs,Q⁻,γ,q=q,r=r) return J, guvec - - # # control penalty - # Jcontrol = uvec'*(Q⁻*uvec) - - # u = unvec(u₀,uvec) - # b = adjustboundarycondition(b₀,u) # combine b₀, u - # c = steadyinversion(Alu,b,γ,q=q,r=r) # gives the misfit - # # observe at right spots - # ỹ = observe(c,wis,γ) - # n = ỹ - y - - # Jdata = n ⋅ (Wⁱ * n) # dot product - # J = Jdata + Jcontrol - - # ## start adjoint model - # # initialize adjoint control variable with zero - # gu = unvec(u₀,0 .* uvec) - - # gn = 2Wⁱ * n - # gỹ = gn - - # gc = gobserve(gỹ,c,locs) - # gb = gsteadyinversion(gc, Alu, b, γ) - # #gu = gadjustboundarycondition(gb,u) - # gadjustboundarycondition!(gu,gb) - - # # control penalty gradient - # guvec = 2*(Q⁻*uvec) - - # guvec += vec(gu) end """ @@ -1891,9 +1861,9 @@ function costfunction_point_obs!(J,guvec::Union{Nothing,Vector},uvec::Vector,Alu n = ỹ - y if guvec != nothing - # adjoint equations - #guvec = 0.0.*uvec # re-initialize - gtmp = 2*(Q⁻*uvec) + ## start adjoint model + + gtmp = 2*(Q⁻*uvec) # control penalty gradient gn = 2Wⁱ * n gỹ = gn gc = gobserve(gỹ,c,locs) @@ -1903,13 +1873,9 @@ function costfunction_point_obs!(J,guvec::Union{Nothing,Vector},uvec::Vector,Alu for (ii,vv) in enumerate(gtmp) guvec[ii] = vv end - println(maximum(guvec)) - println(minimum(guvec)) end if J !=nothing - # reinitialize for Optim.jl - #J = 0.0 # control penalty and gradient Jcontrol = uvec'*(Q⁻*uvec) Jdata = n ⋅ (Wⁱ * n) # dot product From 2c5bd464a0fe31ddfcf8111d807d969dd8ba5e03 Mon Sep 17 00:00:00 2001 From: G Jake Gebbie Date: Mon, 3 Apr 2023 08:57:20 -0400 Subject: [PATCH 09/39] wip: error in adjustsources test --- test/runtests.jl | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/test/runtests.jl b/test/runtests.jl index 4de68583..bf52532d 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -245,12 +245,11 @@ using TMI #u = (;surface = zerosurfaceboundary(γ)) u = (; source = zeros(γ)) #b = (; surface = bPO₄) - b = bPO₄ # assume a surface boundary condition + b = (; surface = bPO₄) # assume a surface boundary condition PO₄ = steadyinversion(Alu,b,γ,q=qPO₄) uvec = vec(u) - σb = 5.0 σq = 0.2 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=qPO₄) From b92fd0c21fed674719653f8095960ad4600dd502 Mon Sep 17 00:00:00 2001 From: G Jake Gebbie Date: Mon, 3 Apr 2023 09:32:51 -0400 Subject: [PATCH 10/39] use `gadjustboundarycondition!` --- src/TMI.jl | 8 ++++++-- 1 file changed, 6 insertions(+), 2 deletions(-) diff --git a/src/TMI.jl b/src/TMI.jl index 2c4c5ef0..a2b6c336 100644 --- a/src/TMI.jl +++ b/src/TMI.jl @@ -1862,13 +1862,17 @@ function costfunction_point_obs!(J,guvec::Union{Nothing,Vector},uvec::Vector,Alu if guvec != nothing ## 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) + #gu = gadjustboundarycondition(gb,u) + gadjustboundarycondition!(gu,gb) gtmp += vec(gu) for (ii,vv) in enumerate(gtmp) guvec[ii] = vv @@ -2224,7 +2228,7 @@ function gadjustboundarycondition!(gu::NamedTuple,gb::NamedTuple) end #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{<: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 From e79469a5ab655f966fc760de38944ff2fffefad2 Mon Sep 17 00:00:00 2001 From: G Jake Gebbie Date: Mon, 3 Apr 2023 10:40:59 -0400 Subject: [PATCH 11/39] `adjustsource` has nested `adjustsource!` --- src/TMI.jl | 27 +++++++++++++++++++-------- test/runtests.jl | 6 +++--- 2 files changed, 22 insertions(+), 11 deletions(-) diff --git a/src/TMI.jl b/src/TMI.jl index a2b6c336..c017497b 100644 --- a/src/TMI.jl +++ b/src/TMI.jl @@ -419,10 +419,10 @@ end - `fg!`: compute cost function and gradient in place - `γ`: grid """ -function sparsedatamap(u₀::Vector,Alu,b::Union{BoundaryCondition,NamedTuple},u::Union{BoundaryCondition,NamedTuple},y::Vector,W⁻,wis::Vector,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. @@ -1823,7 +1823,7 @@ function costfunction_point_obs(uvec::Vector,Alu,b::Union{BoundaryCondition,Name 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) + J = costfunction_point_obs!(J,guvec,uvec,Alu,b,u,y,Wⁱ,wis,locs,Q⁻,γ,q₀=q,r=r) return J, guvec end @@ -1848,12 +1848,17 @@ end - `Q⁻`: weights for control vector - `γ`: grid """ -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) #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) #where {N1, N2, T <: Real} #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} #unvec!(u,uvec) # causes issues with Optim.jl u = unvec(u₀,uvec) b = adjustboundarycondition(b₀,u) # combine b₀, u + + if !isnothing(q₀) + q = adjustsource(q₀,u) + end + c = steadyinversion(Alu,b,γ,q=q,r=r) # observe at right spots @@ -1871,7 +1876,6 @@ function costfunction_point_obs!(J,guvec::Union{Nothing,Vector},uvec::Vector,Alu gỹ = gn gc = gobserve(gỹ,c,locs) gb = gsteadyinversion(gc, Alu, b, γ) - #gu = gadjustboundarycondition(gb,u) gadjustboundarycondition!(gu,gb) gtmp += vec(gu) for (ii,vv) in enumerate(gtmp) @@ -1883,6 +1887,8 @@ function costfunction_point_obs!(J,guvec::Union{Nothing,Vector},uvec::Vector,Alu # 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 @@ -2233,11 +2239,16 @@ function gadjustboundarycondition(gb::Union{NamedTuple,BoundaryCondition},u::Nam return gu end -function adjustsource!(b::Field,u::Field) +function adjustsource(q₀::Union{Field,NamedTuple},u::Union{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 - b.tracer[b.γ.wet] += u.tracer[u.γ.wet] + q.tracer[q.γ.wet] += u.tracer[u.γ.wet] end - function adjustsource!(q::NamedTuple,u::NamedTuple) #where {N1, N2, T <: Real} for qkey in keys(q) if haskey(u,qkey) diff --git a/test/runtests.jl b/test/runtests.jl index bf52532d..56292afb 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -223,7 +223,7 @@ using TMI @test out.minimum < J̃₀ # reconstruct by hand to double-check. ũ = out.minimizer - J̃,gJ̃ = fg(ũ) + J̃,∂J̃∂ũ = fg(ũ) @test J̃ < J̃₀ end end @@ -250,7 +250,7 @@ using TMI PO₄ = steadyinversion(Alu,b,γ,q=qPO₄) uvec = vec(u) - σq = 0.2 + σq = 0.1 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=qPO₄) f(x) = fg(x)[1] @@ -273,7 +273,7 @@ using TMI 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) + out = sparsedatamap(uvec,Alu,b,u,y,W⁻,wis,locs,Q⁻,γ,q=qPO₄,iterations=iterations) # was cost function decreased? @test out.minimum < J̃₀ # reconstruct by hand to double-check. From b45a994934d44524ac5c186476dae4eedf3fe557 Mon Sep 17 00:00:00 2001 From: G Jake Gebbie Date: Mon, 3 Apr 2023 14:28:28 -0400 Subject: [PATCH 12/39] adjustsources test passes --- src/TMI.jl | 108 +++++++++++++++++++++++++++++++---------------- test/runtests.jl | 10 ++--- 2 files changed, 76 insertions(+), 42 deletions(-) diff --git a/src/TMI.jl b/src/TMI.jl index c017497b..c5ea1331 100644 --- a/src/TMI.jl +++ b/src/TMI.jl @@ -422,7 +422,7 @@ end 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⁻,γ,q=q,r=r) + 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. @@ -1387,7 +1387,6 @@ function dot(c::Field{T},d::Field{T})::T where T <: Real return e end - """ function setsource!(d::Field,q::Field,r::Number) apply interior source q to the equation constraints d @@ -1397,10 +1396,21 @@ 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 gsetsource!(gq::Field{T},gd::Field{T},r=1.0) + Adjoint to `setsource!` +""" +function gsetsource!(gq::Field{T},gd::Field{T},r) where T<: Real + gq.tracer[gq.γ.wet] -= r * gd.tracer[gd.γ.wet] +end +function gsetsource(gd::Field{T},q,r) where T<: Real + gq = 0.0 * q + gsetsource!(gq,gd,r) + return gq end """ @@ -1848,18 +1858,19 @@ end - `Q⁻`: weights for control vector - `γ`: grid """ -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) #where {N1, N2, T <: Real} -#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) #unvec!(u,uvec) # causes issues with Optim.jl u = unvec(u₀,uvec) b = adjustboundarycondition(b₀,u) # combine b₀, u 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 - - c = steadyinversion(Alu,b,γ,q=q,r=r) # observe at right spots ỹ = observe(c,wis,γ) @@ -1875,7 +1886,14 @@ function costfunction_point_obs!(J,guvec::Union{Nothing,Vector},uvec::Vector,Alu gn = 2Wⁱ * n gỹ = gn gc = gobserve(gỹ,c,locs) - gb = gsteadyinversion(gc, Alu, b, γ) + + if !isnothing(q₀) + gb,gq = gsteadyinversion(gc,Alu,b,γ,q=q,r=r) + gadjustsource!(gu,gq) + else + gb = gsteadyinversion(gc, Alu, b, γ) + end + gadjustboundarycondition!(gu,gb) gtmp += vec(gu) for (ii,vv) in enumerate(gtmp) @@ -1887,8 +1905,8 @@ function costfunction_point_obs!(J,guvec::Union{Nothing,Vector},uvec::Vector,Alu # control penalty and gradient Jcontrol = uvec'*(Q⁻*uvec) Jdata = n ⋅ (Wⁱ * n) # dot product - println("Jcontrol:",Jcontrol) - println("Jdata:",Jdata) + #println("Jcontrol:",Jcontrol) + #println("Jdata:",Jdata) return Jdata + Jcontrol end end @@ -1925,7 +1943,7 @@ function steadyinversion(Alu,b::BoundaryCondition{T},γ::Grid;q=nothing,r=1.0):: 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 @@ -1939,20 +1957,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 """ @@ -1983,29 +1999,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(γ) @@ -2267,6 +2283,24 @@ function adjustsource!(q::Field,u::NamedTuple) #where {N1, N2, T <: Real} 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::NamedTuple,gq::Field) #where {N1, N2, T <: Real} + qkey = :source + if haskey(gu,qkey) + gadjustsource!(gu[qkey],gq) + end +end +function gadjustsource!(gu::NamedTuple,gq::NamedTuple) + for qkey in keys(gq) + if haskey(gu,qkey) + gadjustsource!(gu[qkey],gq[qkey]) + end + end +end + """ function section View latitude-depth slice of field diff --git a/test/runtests.jl b/test/runtests.jl index 56292afb..d2965f76 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -216,8 +216,8 @@ using TMI # 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̃₀ @@ -256,7 +256,7 @@ using TMI 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=qPO₄) + fg!(F,G,x) = costfunction_point_obs!(F,G,x,Alu,b,u,y,W⁻,wis,locs,Q⁻,γ,q₀=qPO₄) ϵ = 1e-3 # size of finite perturbation ii = rand(1:sum(γ.wet[:,:,1])) @@ -272,8 +272,8 @@ using TMI # 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⁻,γ,q=qPO₄,iterations=iterations) + iters = 5 + out = sparsedatamap(uvec,Alu,b,u,y,W⁻,wis,locs,Q⁻,γ,q=qPO₄,r=1.0,iterations=iterations) # was cost function decreased? @test out.minimum < J̃₀ # reconstruct by hand to double-check. From 55d3b51d360952c1a621aff358bb1e86d1703440 Mon Sep 17 00:00:00 2001 From: G Jake Gebbie Date: Mon, 3 Apr 2023 14:28:53 -0400 Subject: [PATCH 13/39] testname: sourcemap --- test/runtests.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/runtests.jl b/test/runtests.jl index d2965f76..72302959 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -227,7 +227,7 @@ using TMI @test J̃ < J̃₀ end end - @testset "adjustsource" begin + @testset "sourcemap" begin using Statistics, Interpolations From 63f1c91306b748f639ab1de30697b0367b88f9d1 Mon Sep 17 00:00:00 2001 From: G Jake Gebbie Date: Mon, 3 Apr 2023 17:08:57 -0400 Subject: [PATCH 14/39] cleanup `unvec!` again next integrate with `unvec` proper --- src/TMI.jl | 94 +++++++++++++++++++++++++++++++----------------- test/runtests.jl | 78 ++++++++++++++++++++++++++++++++++++++++ 2 files changed, 140 insertions(+), 32 deletions(-) diff --git a/src/TMI.jl b/src/TMI.jl index c5ea1331..eba554fe 100644 --- a/src/TMI.jl +++ b/src/TMI.jl @@ -50,7 +50,9 @@ export config, config_from_mat, config_from_nc, gsetboundarycondition, setsource!, zeros, one, oneunit, ones, maximum, minimum, (+), (-), (*), dot, 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 @@ -1169,7 +1171,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,γ) @@ -1436,6 +1438,12 @@ function vec(u::NamedTuple) return uvec end +function zerotemplate(utemplate,uvec) + tracer = zeros(wet(utemplate)) + tracer[wet(utemplate)] = uvec + return tracer +end + """ function unvec(u,uvec) @@ -1444,18 +1452,23 @@ end 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) +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::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::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) @@ -1478,6 +1491,21 @@ function unvec(utemplate::NamedTuple{<:Any,NTuple{N,Field{T}}},uvec::Vector{T}) 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 unvec!(u,uvec) @@ -1488,27 +1516,19 @@ end """ function unvec!(u::NamedTuple,uvec::Vector) #where {N, T <: Real} - counter = 0 + nlo = 1 + nhi = 0 for v in u - if v isa BoundaryCondition - wet = v.wet - elseif v isa Field - wet = v.γ.wet - end - - n = sum(wet) - v.tracer[wet] = uvec[counter+1:counter+n] - counter += n + nhi += sum(wet(v)) + unvec!(v,uvec[nlo:nhi]) + #v.tracer[wet(v)] = uvec[nlo:nhi] + nlo = nhi + 1 end 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] +function unvec!(u::Union{BoundaryCondition{T},Field{T}},uvec::Vector{T}) where T <: Real + I = findall(wet(u)) + for (ii,vv) in enumerate(I) + u.tracer[vv] = uvec[ii] end end @@ -2096,8 +2116,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 @@ -2490,6 +2517,9 @@ function vec2fld(vector::Vector{T},I::Vector{CartesianIndex{3}}) where T<:Real return field end +wet(a::BoundaryCondition) = a.wet +wet(a::Field) = a.γ.wet + include("deprecated.jl") end diff --git a/test/runtests.jl b/test/runtests.jl index 72302959..44f4580b 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -281,4 +281,82 @@ using TMI J̃,∂J̃∂ũ = fg(ũ) @test J̃ < J̃₀ end + @testset "source_and_surfacemap" begin + + using Statistics, Interpolations + + yPO₄ = readfield(TMIfile,"PO₄",γ) + bPO₄ = getsurfaceboundary(yPO₄) + + qPO₄ = readfield(TMIfile,"qPO₄",γ) + + N = 20 + y, W⁻, ctrue, ytrue, locs, wis = synthetic_observations(TMIversion,"PO₄",γ,N) + + u = (; surface = zerosurfaceboundary(γ), source = zeros(γ)) + #b = (; surface = bPO₄) + b = (; surface = bPO₄) # assume a surface boundary condition + + PO₄ = steadyinversion(Alu,b,γ,q=qPO₄) + uvec = vec(u) + + σq = 0.1 + 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=qPO₄) + 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₀=qPO₄) + + ϵ = 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=iterations) + # was cost function decreased? + @test out.minimum < J̃₀ + # reconstruct by hand to double-check. + ũ = out.minimizer + J̃,∂J̃∂ũ = fg(ũ) + @test J̃ < J̃₀ + 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(γ)) + + 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 + end From 4c582a4f5c958a891d43e1672b2e3fc8d1670138 Mon Sep 17 00:00:00 2001 From: G Jake Gebbie Date: Mon, 3 Apr 2023 20:44:24 -0400 Subject: [PATCH 15/39] `source and surface map`: passes test --- src/TMI.jl | 73 ++++++------------------------------------------ test/runtests.jl | 19 +++++++------ 2 files changed, 19 insertions(+), 73 deletions(-) diff --git a/src/TMI.jl b/src/TMI.jl index eba554fe..089a76d7 100644 --- a/src/TMI.jl +++ b/src/TMI.jl @@ -1438,12 +1438,6 @@ function vec(u::NamedTuple) return uvec end -function zerotemplate(utemplate,uvec) - tracer = zeros(wet(utemplate)) - tracer[wet(utemplate)] = uvec - return tracer -end - """ function unvec(u,uvec) @@ -1452,58 +1446,9 @@ end 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)...) +function unvec(u₀::Union{NamedTuple,Field,BoundaryCondition},uvec::Vector) #where T <: Real + u = deepcopy(u₀) + unvec!(u,uvec) return u end @@ -1514,6 +1459,12 @@ end Needs to update u because attributes of u need to be known at runtime. """ +function unvec!(u::Union{BoundaryCondition{T},Field{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::NamedTuple,uvec::Vector) #where {N, T <: Real} nlo = 1 @@ -1525,12 +1476,6 @@ function unvec!(u::NamedTuple,uvec::Vector) #where {N, T <: Real} nlo = nhi + 1 end end -function unvec!(u::Union{BoundaryCondition{T},Field{T}},uvec::Vector{T}) where T <: Real - I = findall(wet(u)) - for (ii,vv) in enumerate(I) - u.tracer[vv] = uvec[ii] - end -end """ Surface oxygen saturation value and fraction of saturation value in field diff --git a/test/runtests.jl b/test/runtests.jl index 44f4580b..25b625d6 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -307,9 +307,9 @@ using TMI 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₀=qPO₄) - + ϵ = 1e-3 # size of finite perturbation - ii = rand(1:sum(γ.wet[:,:,1])) + ii = rand(1:length(uvec)) println("gradient check location=",ii) δu = copy(uvec); δu[ii] += ϵ ∇f_finite = (f(δu) - f(uvec))/ϵ @@ -343,20 +343,21 @@ using TMI 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 - end From c18ee86500d67ab7e2eabb22c814d850a3b7297b Mon Sep 17 00:00:00 2001 From: G Jake Gebbie Date: Mon, 3 Apr 2023 20:48:11 -0400 Subject: [PATCH 16/39] old code in deprecated.jl --- src/deprecated.jl | 111 ++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 111 insertions(+) create mode 100644 src/deprecated.jl diff --git a/src/deprecated.jl b/src/deprecated.jl new file mode 100644 index 00000000..d9a0cf75 --- /dev/null +++ b/src/deprecated.jl @@ -0,0 +1,111 @@ + +# """ +# 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 + From 11cec7a0edb1a862d271f058e47544e96fb7668b Mon Sep 17 00:00:00 2001 From: G Jake Gebbie Date: Mon, 3 Apr 2023 21:55:18 -0400 Subject: [PATCH 17/39] TMIconfig -> config --- src/TMI.jl | 2 +- src/{TMIconfig.jl => config.jl} | 0 2 files changed, 1 insertion(+), 1 deletion(-) rename src/{TMIconfig.jl => config.jl} (100%) diff --git a/src/TMI.jl b/src/TMI.jl index 089a76d7..2d612b82 100644 --- a/src/TMI.jl +++ b/src/TMI.jl @@ -176,7 +176,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) diff --git a/src/TMIconfig.jl b/src/config.jl similarity index 100% rename from src/TMIconfig.jl rename to src/config.jl From 377d88ddb1be5c88b01ba5906391208dcab1d3db Mon Sep 17 00:00:00 2001 From: G Jake Gebbie Date: Mon, 3 Apr 2023 22:09:47 -0400 Subject: [PATCH 18/39] read regional Nordic input --- test/runtests.jl | 13 +++++++++++++ 1 file changed, 13 insertions(+) diff --git a/test/runtests.jl b/test/runtests.jl index 25b625d6..8eb49125 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -361,3 +361,16 @@ using TMI end end + +@testset "regional" begin + + 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 From 9e40e6dc5a5d38e7167444786fdf08a70f7b351b Mon Sep 17 00:00:00 2001 From: G Jake Gebbie Date: Tue, 4 Apr 2023 14:12:39 -0400 Subject: [PATCH 19/39] shell script for reading regional matrix --- src/read_nc_nordic_lowresolution.sh | 6 ++++++ 1 file changed, 6 insertions(+) create mode 100644 src/read_nc_nordic_lowresolution.sh diff --git a/src/read_nc_nordic_lowresolution.sh b/src/read_nc_nordic_lowresolution.sh new file mode 100644 index 00000000..be4a59cf --- /dev/null +++ b/src/read_nc_nordic_lowresolution.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. + +# Google ID = 1mAmJKS3xeNFjKLLTMBpgSS-aWZ1WHNm- + +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=1mAmJKS3xeNFjKLLTMBpgSS-aWZ1WHNm-' -O- | sed -rn 's/.*confirm=([0-9A-Za-z_]+).*/\1\n/p')&id=1mAmJKS3xeNFjKLLTMBpgSS-aWZ1WHNm-" -O NordicSeas_30-Dec-2020_low.nc && rm -rf /tmp/cookies.txt From 24d02b03ad589e165fe884a133793e13979dd8ba Mon Sep 17 00:00:00 2001 From: G Jake Gebbie Date: Tue, 4 Apr 2023 14:13:28 -0400 Subject: [PATCH 20/39] minor change --- test/runtests.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/test/runtests.jl b/test/runtests.jl index 8eb49125..6d2a748d 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -273,7 +273,7 @@ using TMI 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=iterations) + 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. @@ -323,7 +323,7 @@ using TMI 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=iterations) + 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. From a29b0e72f843b1643e071339a30768289a9ad85a Mon Sep 17 00:00:00 2001 From: G Jake Gebbie Date: Tue, 4 Apr 2023 16:41:47 -0400 Subject: [PATCH 21/39] wip: UnicodePlots and read NordicSeas in src --- Manifest.toml | 55 +++++++++++++++++++++++++++++++++++++++++++- Project.toml | 1 + src/TMI.jl | 59 +++++++++++++++++++++++++++++++++++++++++++++++- src/config.jl | 7 ++++++ test/runtests.jl | 25 +++++++++++++++++++- 5 files changed, 144 insertions(+), 3 deletions(-) diff --git a/Manifest.toml b/Manifest.toml index 838c0f4f..2197c88c 100644 --- a/Manifest.toml +++ b/Manifest.toml @@ -2,7 +2,7 @@ julia_version = "1.8.5" manifest_format = "2.0" -project_hash = "9a28a4abf9f655c03a8cb01bab4e21c8a88412f5" +project_hash = "c02643bd8035ebad2cc399fe86fb6c70dd21fbf9" [[deps.ANSIColoredPrinters]] git-tree-sha1 = "574baf8110975760d391c710b6341da1afa48d8c" @@ -89,6 +89,30 @@ 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" @@ -118,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" @@ -243,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" @@ -445,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" @@ -764,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" @@ -791,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 10a2c2e8..dc05664d 100644 --- a/Project.toml +++ b/Project.toml @@ -27,6 +27,7 @@ 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/src/TMI.jl b/src/TMI.jl index 2d612b82..4a5942d1 100644 --- a/src/TMI.jl +++ b/src/TMI.jl @@ -71,6 +71,7 @@ struct Grid R::Array{Int,3} # R::LinearIndices{3, Tuple{UnitRange{Int64}, UnitRange{Int64}, UnitRange{Int64}}} wet::BitArray{3} + #interior::BitArray{3} end """ @@ -161,6 +162,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} + interior::BitArray{3} + name::Symbol + longname::String + units::String + nonnegative::Bool +end # Credit to DrWatson.jl for these functions # Didn't want to add dependency for these small functions @@ -588,6 +605,42 @@ function writefield(file,field::Field{T}) where T <: Real return nothing end +function readsource(file,tracername,γ::Grid) + + # 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] + + # load all data + tracer = v[:,:,:] + #println(tracer) + + # Make an interior mask + + + # load an attribute + units = v.attrib["units"] + longname = v.attrib["longname"] + + close(ds) + + # perform a check of file compatibility + # with grid + if sum(isnan.(tracer[γ.wet])) > 0 + 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 + + #c = Field(tracer,γ) + c = Field(tracer,γ,Symbol(tracername),longname,units) + + return c +end + """ function depthindex(I) Get the k-index (depth level) from the Cartesian index @@ -1460,7 +1513,11 @@ end u need to be known at runtime. """ function unvec!(u::Union{BoundaryCondition{T},Field{T}},uvec::Vector{T}) where T <: Real - I = findall(wet(u)) # findall seems slow + if u isa BoundaryCondition + I = findall(wet(u)) # findall seems slow + elseif u isa Field + I = u.γ.I + end for (ii,vv) in enumerate(I) u.tracer[vv] = uvec[ii] end diff --git a/src/config.jl b/src/config.jl index 4efc8f7d..94ec5226 100644 --- a/src/config.jl +++ b/src/config.jl @@ -175,6 +175,13 @@ 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 == "NordicSeas_30-Dec-2020_low" + # problems with file name inconsistency here + println("workaround for regional Nordic Seas file") + shellscript = pkgsrcdir("read_mat_nordic_lowresolution.sh") + run(`sh $shellscript`) + TMIfile = TMIversion*".nc" + mv(joinpath(pwd(),TMIfilegz),pkgdatadir(TMIfile),force=true) else !isfile(TMIfilegz) & !isfile(TMIfile) && google_download(url,pkgdatadir()) end diff --git a/test/runtests.jl b/test/runtests.jl index 6d2a748d..ab85aa47 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -31,7 +31,30 @@ 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₄,γ) + qPO₄ = readfield(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 From 9550f529b9dc02bbddd2dde7ed2863a2957984ae Mon Sep 17 00:00:00 2001 From: G Jake Gebbie Date: Tue, 4 Apr 2023 21:46:33 -0400 Subject: [PATCH 22/39] read netcdf nordic with consistent naming --- src/config.jl | 5 ++--- ...ion.sh => read_nc_nordic_201x115x46_B23.sh} | 5 +++-- test/runtests.jl | 18 +++++++++--------- 3 files changed, 14 insertions(+), 14 deletions(-) rename src/{read_nc_nordic_lowresolution.sh => read_nc_nordic_201x115x46_B23.sh} (59%) diff --git a/src/config.jl b/src/config.jl index 94ec5226..3069ca4f 100644 --- a/src/config.jl +++ b/src/config.jl @@ -175,10 +175,9 @@ 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 == "NordicSeas_30-Dec-2020_low" - # problems with file name inconsistency here + elseif TMIversion == "nordic_201x115x46_B23" println("workaround for regional Nordic Seas file") - shellscript = pkgsrcdir("read_mat_nordic_lowresolution.sh") + shellscript = pkgsrcdir("read_mat_nordic_201x115x46_B23.sh") run(`sh $shellscript`) TMIfile = TMIversion*".nc" mv(joinpath(pwd(),TMIfilegz),pkgdatadir(TMIfile),force=true) diff --git a/src/read_nc_nordic_lowresolution.sh b/src/read_nc_nordic_201x115x46_B23.sh similarity index 59% rename from src/read_nc_nordic_lowresolution.sh rename to src/read_nc_nordic_201x115x46_B23.sh index be4a59cf..0df817dd 100644 --- a/src/read_nc_nordic_lowresolution.sh +++ b/src/read_nc_nordic_201x115x46_B23.sh @@ -1,6 +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. -# Google ID = 1mAmJKS3xeNFjKLLTMBpgSS-aWZ1WHNm- +# 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=1mAmJKS3xeNFjKLLTMBpgSS-aWZ1WHNm-' -O- | sed -rn 's/.*confirm=([0-9A-Za-z_]+).*/\1\n/p')&id=1mAmJKS3xeNFjKLLTMBpgSS-aWZ1WHNm-" -O NordicSeas_30-Dec-2020_low.nc && rm -rf /tmp/cookies.txt +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 ab85aa47..5510ad90 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -385,15 +385,15 @@ using TMI end end -@testset "regional" begin +# @testset "regional" begin - shellscript = TMI.pkgsrcdir("read_nc_nordic_lowresolution.sh") - run(`sh $shellscript`) +# 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) +# # 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 From 2d6fa91c1d1dd0635e68694b111ee5e3de954ed2 Mon Sep 17 00:00:00 2001 From: G Jake Gebbie Date: Tue, 4 Apr 2023 21:49:32 -0400 Subject: [PATCH 23/39] read nordic MATLAB file --- src/read_mat_nordic_201x115x46_B23.sh | 6 ++++++ 1 file changed, 6 insertions(+) create mode 100644 src/read_mat_nordic_201x115x46_B23.sh 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 From d1dfb67a6f47085fa4b365009ae99476e7f704b0 Mon Sep 17 00:00:00 2001 From: G Jake Gebbie Date: Tue, 4 Apr 2023 21:59:46 -0400 Subject: [PATCH 24/39] bugfix `download_matfile` for nordic --- src/config.jl | 9 +++++++-- 1 file changed, 7 insertions(+), 2 deletions(-) diff --git a/src/config.jl b/src/config.jl index 3069ca4f..9b7792c5 100644 --- a/src/config.jl +++ b/src/config.jl @@ -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 @@ -179,8 +185,7 @@ function download_matfile(TMIversion::String) println("workaround for regional Nordic Seas file") shellscript = pkgsrcdir("read_mat_nordic_201x115x46_B23.sh") run(`sh $shellscript`) - TMIfile = TMIversion*".nc" - mv(joinpath(pwd(),TMIfilegz),pkgdatadir(TMIfile),force=true) + mv(joinpath(pwd(),TMIfilegz),TMIfile,force=true) else !isfile(TMIfilegz) & !isfile(TMIfile) && google_download(url,pkgdatadir()) end From 526714dd259aabbe634a5f3550a83a2cb2e178d0 Mon Sep 17 00:00:00 2001 From: G Jake Gebbie Date: Wed, 5 Apr 2023 15:00:40 -0400 Subject: [PATCH 25/39] grid elements computed on demand --- src/TMI.jl | 24 +++++++++++++++++++++--- src/config.jl | 23 ++++++++++++++--------- 2 files changed, 35 insertions(+), 12 deletions(-) diff --git a/src/TMI.jl b/src/TMI.jl index 4a5942d1..9b76dc99 100644 --- a/src/TMI.jl +++ b/src/TMI.jl @@ -67,13 +67,31 @@ 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 + + """ struct Field diff --git a/src/config.jl b/src/config.jl index 9b7792c5..36b6dc6e 100644 --- a/src/config.jl +++ b/src/config.jl @@ -95,22 +95,27 @@ function Grid(TMIfile::String) # get properties of grid lon,lat,depth = gridprops(TMIfile) + # make ocean mask + wet = wetmask(TMIfile,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) +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 - - R = linearindex(wet) - - γ = Grid(lon,lat,depth,I,R,wet) - - return γ + return wet end - + """ Configure TMI environment from original MATLAB output """ From a54262c7a1c9692b01622310280763ff6cf9da1e Mon Sep 17 00:00:00 2001 From: G Jake Gebbie Date: Wed, 5 Apr 2023 16:52:14 -0400 Subject: [PATCH 26/39] add `interior` mask to Grid --- src/TMI.jl | 5 ++--- src/config.jl | 22 +++++++++++++++++----- 2 files changed, 19 insertions(+), 8 deletions(-) diff --git a/src/TMI.jl b/src/TMI.jl index 9b76dc99..e6c96c6a 100644 --- a/src/TMI.jl +++ b/src/TMI.jl @@ -68,7 +68,7 @@ struct Grid lat::Vector{Float64} depth::Vector{Float64} wet::BitArray{3} - #interior::BitArray{3} + interior::BitArray{3} # I::Vector{CartesianIndex{3}} # index # R::Array{Int,3} # R::LinearIndices{3, Tuple{UnitRange{Int64}, UnitRange{Int64}, UnitRange{Int64}}} @@ -91,7 +91,6 @@ function Base.getproperty(γ::Grid, d::Symbol) end end - """ struct Field @@ -190,7 +189,7 @@ end """ struct Source{T} tracer::Array{T,3} - interior::BitArray{3} + γ::Grid name::Symbol longname::String units::String diff --git a/src/config.jl b/src/config.jl index 36b6dc6e..1d8a2a2a 100644 --- a/src/config.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,γ) @@ -91,18 +91,22 @@ 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) + return Grid(lon,lat,depth,wet,interior) end function wetmask(TMIfile,nx,ny,nz) @@ -115,7 +119,15 @@ function wetmask(TMIfile,nx,ny,nz) wet[I] .= 1 return wet end - + +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 """ From e05ade1b5dd396b0724c476068fe69f9c3d83d58 Mon Sep 17 00:00:00 2001 From: G Jake Gebbie Date: Wed, 5 Apr 2023 17:14:30 -0400 Subject: [PATCH 27/39] `readsource` plus a source type --- src/TMI.jl | 52 ++++++++++++++---------------------------------- test/runtests.jl | 2 +- 2 files changed, 16 insertions(+), 38 deletions(-) diff --git a/src/TMI.jl b/src/TMI.jl index e6c96c6a..cd7fad5c 100644 --- a/src/TMI.jl +++ b/src/TMI.jl @@ -193,7 +193,7 @@ struct Source{T} name::Symbol longname::String units::String - nonnegative::Bool + logscale::Bool end # Credit to DrWatson.jl for these functions @@ -545,11 +545,14 @@ function readfield(file,tracername,γ::Grid) if sum(isnan.(tracer[γ.wet])) > 0 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) @@ -622,40 +625,15 @@ function writefield(file,field::Field{T}) where T <: Real return nothing end -function readsource(file,tracername,γ::Grid) - - # 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] - - # load all data - tracer = v[:,:,:] - #println(tracer) - - # Make an interior mask - - - # load an attribute - units = v.attrib["units"] - longname = v.attrib["longname"] - - close(ds) - - # perform a check of file compatibility - # with grid - if sum(isnan.(tracer[γ.wet])) > 0 - error("readfield warning: NaN on grid") +function readsource(file,tracername,γ::Grid;logscale=false) + c = readfield(file,tracername,γ) + if logscale + ct = log.(c.tracer) + else + ct = c.tracer 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 - - #c = Field(tracer,γ) - c = Field(tracer,γ,Symbol(tracername),longname,units) - - return c + q = Source(ct,c.γ,c.name,c.longname,c.units,logscale) + return q end """ diff --git a/test/runtests.jl b/test/runtests.jl index 5510ad90..47d83d95 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -14,7 +14,7 @@ using TMI yPO₄ = readfield(TMIfile,"PO₄",γ) bPO₄ = getsurfaceboundary(yPO₄) PO₄pre = steadyinversion(Alu,bPO₄,γ) - qPO₄ = readfield(TMIfile,"qPO₄",γ) + qPO₄ = TMI.readsource(TMIfile,"qPO₄",γ) b₀ = zerosurfaceboundary(γ) PO₄ᴿ = steadyinversion(Alu,b₀,γ,q=qPO₄) PO₄total = PO₄ᴿ + PO₄pre From 58c7f22640ac7faccb6224016040edfd99f7e661 Mon Sep 17 00:00:00 2001 From: G Jake Gebbie Date: Wed, 5 Apr 2023 21:57:31 -0400 Subject: [PATCH 28/39] wip: logarithmic source adjustment --- src/TMI.jl | 223 +++++++++++----------------------------------- src/deprecated.jl | 156 ++++++++++++++++++++++++++++++++ test/runtests.jl | 22 ++--- 3 files changed, 222 insertions(+), 179 deletions(-) diff --git a/src/TMI.jl b/src/TMI.jl index cd7fad5c..07c0eed6 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,11 @@ export config, config_from_mat, config_from_nc, surface_oxygensaturation, oxygen, location_obs, getsurfaceboundary, zerosurfaceboundary, onesurfaceboundary, setboundarycondition!, - #adjustboundarycondition, adjustboundarycondition!, gsetboundarycondition, setsource!, - zeros, one, oneunit, ones, maximum, minimum, (+), (-), (*), dot, + zeros, one, oneunit, ones, + maximum, minimum, (+), (-), (*), dot, + zerosource, Grid, Field, BoundaryCondition, vec, unvec!, unvec, adjustsource!, zerowestboundary, zeronorthboundary, zeroeastboundary, zerosouthboundary, @@ -130,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 @@ -1011,6 +1027,16 @@ 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 ones(γ::Grid,name=:none,longname="unknown",units="unknown")::Field @@ -1448,6 +1474,9 @@ end function setsource!(d::Field{T},q::Field{T},r=1.0) where T<: Real 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 gsetsource!(gq::Field{T},gd::Field{T},r=1.0) @@ -1472,8 +1501,7 @@ end """ vec(u::Field) = u.tracer[u.γ.wet] vec(u::BoundaryCondition) = u.tracer[u.wet] - -#function vec(u::NamedTuple{<:Any,NTuple{N,BoundaryCondition{T}}}) where {N, T <: Real} +vec(u::Source) = u.tracer[u.γ.interior] function vec(u::NamedTuple) T = eltype(values(u)[1].tracer) @@ -1507,24 +1535,18 @@ end Needs to update u because attributes of u need to be known at runtime. """ -function unvec!(u::Union{BoundaryCondition{T},Field{T}},uvec::Vector{T}) where T <: Real - if u isa BoundaryCondition - I = findall(wet(u)) # findall seems slow - elseif u isa Field - I = u.γ.I - end +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::NamedTuple,uvec::Vector) #where {N, T <: Real} - nlo = 1 nhi = 0 for v in u nhi += sum(wet(v)) unvec!(v,uvec[nlo:nhi]) - #v.tracer[wet(v)] = uvec[nlo:nhi] nlo = nhi + 1 end end @@ -2279,7 +2301,7 @@ function gadjustboundarycondition(gb::Union{NamedTuple,BoundaryCondition},u::Nam return gu end -function adjustsource(q₀::Union{Field,NamedTuple},u::Union{Field,NamedTuple}) +function adjustsource(q₀::Union{Source,Field,NamedTuple},u::Union{Source,Field,NamedTuple}) q = deepcopy(q₀) adjustsource!(q,u) return q @@ -2289,6 +2311,23 @@ 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 + 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) @@ -2298,7 +2337,7 @@ function adjustsource!(q::NamedTuple,u::NamedTuple) #where {N1, N2, T <: Real} end end end -function adjustsource!(q::Field,u::NamedTuple) #where {N1, N2, T <: Real} +function adjustsource!(q::Union{Field,Source},u::NamedTuple) #where {N1, N2, T <: Real} qkey = :source if haskey(u,qkey) adjustsource!(q,u[qkey]) @@ -2359,163 +2398,9 @@ 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 - -""" - 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 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 - wet(a::BoundaryCondition) = a.wet wet(a::Field) = a.γ.wet +wet(a::Source) = a.γ.interior include("deprecated.jl") diff --git a/src/deprecated.jl b/src/deprecated.jl index d9a0cf75..dad35e89 100644 --- a/src/deprecated.jl +++ b/src/deprecated.jl @@ -109,3 +109,159 @@ # 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 diff --git a/test/runtests.jl b/test/runtests.jl index 47d83d95..8baf998b 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -14,7 +14,7 @@ using TMI yPO₄ = readfield(TMIfile,"PO₄",γ) bPO₄ = getsurfaceboundary(yPO₄) PO₄pre = steadyinversion(Alu,bPO₄,γ) - qPO₄ = TMI.readsource(TMIfile,"qPO₄",γ) + qPO₄ = readfield(TMIfile,"qPO₄",γ) b₀ = zerosurfaceboundary(γ) PO₄ᴿ = steadyinversion(Alu,b₀,γ,q=qPO₄) PO₄total = PO₄ᴿ + PO₄pre @@ -37,7 +37,9 @@ using TMI yPO₄ = readfield(TMIfile,"PO₄",γ) bPO₄ = getsurfaceboundary(yPO₄) PO₄pre = steadyinversion(Alu,bPO₄,γ) - qPO₄ = readfield(TMIfile,"qPO₄",γ) + + # this line changes + qPO₄ = readsource(TMIfile,"qPO₄",γ) b₀ = zerosurfaceboundary(γ) PO₄ᴿ = steadyinversion(Alu,b₀,γ,q=qPO₄) PO₄total = PO₄ᴿ + PO₄pre @@ -257,18 +259,18 @@ using TMI yPO₄ = readfield(TMIfile,"PO₄",γ) bPO₄ = getsurfaceboundary(yPO₄) - #PO₄pre = steadyinversion(Alu,bPO₄,γ) - qPO₄ = readfield(TMIfile,"qPO₄",γ) - #b₀ = bPO₄ #zerosurfaceboundary(γ) - #PO₄ᴿ = steadyinversion(Alu,b₀,γ,q=qPO₄) + if rand([1,2]) == 1 + qPO₄ = readfield(TMIfile,"qPO₄",γ) + else + qPO₄ = readsource(TMIfile,"qPO₄",γ) + end N = 20 y, W⁻, ctrue, ytrue, locs, wis = synthetic_observations(TMIversion,"PO₄",γ,N) - #u = (;surface = zerosurfaceboundary(γ)) - u = (; source = zeros(γ)) - #b = (; surface = bPO₄) - b = (; surface = bPO₄) # assume a surface boundary condition + #u = (; source = zeros(γ)) + u = (; source = zerosource(γ)) + b = (; surface = bPO₄) # surface boundary condition PO₄ = steadyinversion(Alu,b,γ,q=qPO₄) uvec = vec(u) From 1f2bd628140ecee06bcedff45128d95145a3f415 Mon Sep 17 00:00:00 2001 From: G Jake Gebbie Date: Thu, 6 Apr 2023 08:02:57 -0400 Subject: [PATCH 29/39] streamline addition/subtraction --- src/TMI.jl | 131 +++++++++++++++++++++++++++++----------------- src/deprecated.jl | 2 +- 2 files changed, 83 insertions(+), 50 deletions(-) diff --git a/src/TMI.jl b/src/TMI.jl index 07c0eed6..e605e28a 100644 --- a/src/TMI.jl +++ b/src/TMI.jl @@ -1358,64 +1358,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::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::Field)::Field` Define scalar or matrix multiplication for fields """ -function *(C,d::Field{T})::Field{T} where T <: Real +function Base.:*(C,d::Field{T})::Field{T} where T <: Real e = zeros(d.γ) e.tracer[e.γ.wet] += C*d.tracer[d.γ.wet] @@ -1426,7 +1458,7 @@ end `function *(C,d::BoundaryCondition)::BoundaryCondition` Define scalar or matrix multiplication for BoundaryCondition`s """ -function *(C,d::BoundaryCondition{T})::BoundaryCondition{T} where T <: Real +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] @@ -1437,7 +1469,7 @@ 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 Base.:*(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") @@ -2321,9 +2353,10 @@ function adjustsource!(q::Source,u::Source) q.tracer[q.γ.interior] += u.tracer[u.γ.interior] q.tracer[q.γ.interior] = exp.(q.tracer[q.γ.interior]) elseif q.logscale && ~u.logscale - 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]) + 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 diff --git a/src/deprecated.jl b/src/deprecated.jl index dad35e89..2f373625 100644 --- a/src/deprecated.jl +++ b/src/deprecated.jl @@ -163,7 +163,7 @@ end function Γsfc Γsfc anonymously calls surfacecontrol2field """ -Γsfc = surfacecontrol2field +#Γsfc = surfacecontrol2field """ function surfacecontrol2field!(c,u,γ) From bbd118a15e9005e381ecf165a4791ec909799442 Mon Sep 17 00:00:00 2001 From: G Jake Gebbie Date: Thu, 6 Apr 2023 08:05:51 -0400 Subject: [PATCH 30/39] move old +/- to deprecated --- src/TMI.jl | 56 ---------------------------------------------- src/deprecated.jl | 57 +++++++++++++++++++++++++++++++++++++++++++++++ test/runtests.jl | 7 +----- 3 files changed, 58 insertions(+), 62 deletions(-) diff --git a/src/TMI.jl b/src/TMI.jl index e605e28a..54d7ff13 100644 --- a/src/TMI.jl +++ b/src/TMI.jl @@ -1386,62 +1386,6 @@ function Base.:-(c::T,d::T) where T <: Union{Source,Field,BoundaryCondition} return e 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::Field)::Field` diff --git a/src/deprecated.jl b/src/deprecated.jl index 2f373625..f4c2290e 100644 --- a/src/deprecated.jl +++ b/src/deprecated.jl @@ -265,3 +265,60 @@ function surfacecontrol2field(usfc::Vector{T},wet) where T<: Real 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 diff --git a/test/runtests.jl b/test/runtests.jl index 8baf998b..08fe6bb3 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -258,12 +258,7 @@ using TMI yPO₄ = readfield(TMIfile,"PO₄",γ) bPO₄ = getsurfaceboundary(yPO₄) - - if rand([1,2]) == 1 - qPO₄ = readfield(TMIfile,"qPO₄",γ) - else - qPO₄ = readsource(TMIfile,"qPO₄",γ) - end + qPO₄ = readsource(TMIfile,"qPO₄",γ) N = 20 y, W⁻, ctrue, ytrue, locs, wis = synthetic_observations(TMIversion,"PO₄",γ,N) From 0b7e1ce62e1754d08e36da234674e200b9c75749 Mon Sep 17 00:00:00 2001 From: G Jake Gebbie Date: Thu, 6 Apr 2023 08:18:18 -0400 Subject: [PATCH 31/39] streamline multiplication for BC,Source,Field --- src/TMI.jl | 39 ++++++--------------------------------- src/deprecated.jl | 30 ++++++++++++++++++++++++++++++ 2 files changed, 36 insertions(+), 33 deletions(-) diff --git a/src/TMI.jl b/src/TMI.jl index 54d7ff13..2a43be3f 100644 --- a/src/TMI.jl +++ b/src/TMI.jl @@ -1386,49 +1386,22 @@ function Base.:-(c::T,d::T) where T <: Union{Source,Field,BoundaryCondition} return e end - """ `function *(C,d::Field)::Field` Define scalar or matrix multiplication for fields """ -function Base.:*(C,d::Field{T})::Field{T} where T <: Real +Base.:*(C,d::Union{Field,BoundaryCondition,Source}) = d*C - e = zeros(d.γ) - e.tracer[e.γ.wet] += C*d.tracer[d.γ.wet] +function Base.:*(d::Union{Field,BoundaryCondition,Source},C) + e = deepcopy(d) + mul!(e,C) 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 +function mul!(d::Union{Field,BoundaryCondition,Source},C) + d.tracer[wet(d)] *= C #*d.tracer[wet(d)] end -""" - `function *(c::Field,d::Field)::Field` - Field by field multiplication is element-by-element. -""" -function Base.:*(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) - e.tracer[e.γ.wet] += c.tracer[c.γ.wet] .* 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/deprecated.jl b/src/deprecated.jl index f4c2290e..aa76af10 100644 --- a/src/deprecated.jl +++ b/src/deprecated.jl @@ -322,3 +322,33 @@ end # 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. +# """ +# function Base.:*(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) +# e.tracer[e.γ.wet] += c.tracer[c.γ.wet] .* d.tracer[d.γ.wet] +# return e +# end From be56fa298f8285e5d383d4492899178fdf8e78ed Mon Sep 17 00:00:00 2001 From: G Jake Gebbie Date: Thu, 6 Apr 2023 09:35:35 -0400 Subject: [PATCH 32/39] refactor multiplication --- src/TMI.jl | 46 +++++++++++++++++++++++++++++++++++++++++----- src/deprecated.jl | 15 --------------- 2 files changed, 41 insertions(+), 20 deletions(-) diff --git a/src/TMI.jl b/src/TMI.jl index 2a43be3f..aeca143a 100644 --- a/src/TMI.jl +++ b/src/TMI.jl @@ -1389,25 +1389,61 @@ end """ `function *(C,d::Field)::Field` Define scalar or matrix multiplication for fields + + one argument can be a number, but not both (type piracy?) """ -Base.:*(C,d::Union{Field,BoundaryCondition,Source}) = d*C +# the order doesn't matter when multiplying by a scalar -function Base.:*(d::Union{Field,BoundaryCondition,Source},C) +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!(e,C) + mul!(c,e) + return e +end +function Base.:*(d::T,c::Union{Number,T}) where T <: Union{Field,BoundaryCondition,Source} + e = deepcopy(d) + mul!(e,c) return e end -function mul!(d::Union{Field,BoundaryCondition,Source},C) +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 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 + 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 diff --git a/src/deprecated.jl b/src/deprecated.jl index aa76af10..4d4449c9 100644 --- a/src/deprecated.jl +++ b/src/deprecated.jl @@ -337,18 +337,3 @@ end # """ # `function *(c::Field,d::Field)::Field` # Field by field multiplication is element-by-element. -# """ -# function Base.:*(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) -# e.tracer[e.γ.wet] += c.tracer[c.γ.wet] .* d.tracer[d.γ.wet] -# return e -# end From 2dcee19c2204385d91f55960e1e2ed7d75e12cbb Mon Sep 17 00:00:00 2001 From: G Jake Gebbie Date: Thu, 6 Apr 2023 11:08:00 -0400 Subject: [PATCH 33/39] gadjustsource updated for Source --- src/TMI.jl | 28 ++++++++++++++++++++++++++-- 1 file changed, 26 insertions(+), 2 deletions(-) diff --git a/src/TMI.jl b/src/TMI.jl index aeca143a..be2b0e7d 100644 --- a/src/TMI.jl +++ b/src/TMI.jl @@ -1471,7 +1471,10 @@ end function gsetsource!(gq::Field{T},gd::Field{T},r) where T<: Real gq.tracer[gq.γ.wet] -= r * gd.tracer[gd.γ.wet] end -function gsetsource(gd::Field{T},q,r) where T<: Real +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 @@ -2336,7 +2339,28 @@ 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::NamedTuple,gq::Field) #where {N1, N2, T <: Real} +function gadjustsource!(gu::Source,gq::Source) + 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 + elseif gu.logscale + error("not implemented") + # 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 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::Union{Source,Field}) #where {N1, N2, T <: Real} qkey = :source if haskey(gu,qkey) gadjustsource!(gu[qkey],gq) From 4d2a865e1379179054e43de26a955a106653c01e Mon Sep 17 00:00:00 2001 From: G Jake Gebbie Date: Thu, 6 Apr 2023 14:59:50 -0400 Subject: [PATCH 34/39] maximum, minimum update --- src/TMI.jl | 44 ++++++++++++++++++++++++-------------------- test/runtests.jl | 26 +++++++++++++++++++------- 2 files changed, 43 insertions(+), 27 deletions(-) diff --git a/src/TMI.jl b/src/TMI.jl index be2b0e7d..cb6c0b95 100644 --- a/src/TMI.jl +++ b/src/TMI.jl @@ -51,9 +51,9 @@ export config, config_from_mat, config_from_nc, gsetboundarycondition, setsource!, zeros, one, oneunit, ones, maximum, minimum, (+), (-), (*), dot, - zerosource, + zerosource, onesource, Grid, Field, BoundaryCondition, vec, unvec!, unvec, - adjustsource!, + adjustsource, adjustsource!, zerowestboundary, zeronorthboundary, zeroeastboundary, zerosouthboundary, onewestboundary, onenorthboundary, oneeastboundary, onesouthboundary @@ -1028,13 +1028,19 @@ function zeros(γ::Grid,name=:none,longname="unknown",units="unknown")::Field 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 """ @@ -1332,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 @@ -1887,7 +1889,6 @@ end """ 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) - #unvec!(u,uvec) # causes issues with Optim.jl u = unvec(u₀,uvec) b = adjustboundarycondition(b₀,u) # combine b₀, u @@ -1916,7 +1917,7 @@ function costfunction_point_obs!(J,guvec::Union{Nothing,Vector},uvec::Vector,Alu if !isnothing(q₀) gb,gq = gsteadyinversion(gc,Alu,b,γ,q=q,r=r) - gadjustsource!(gu,gq) + gadjustsource!(gu,gq,q) # pass q to linearize logscale version else gb = gsteadyinversion(gc, Alu, b, γ) end @@ -2339,17 +2340,20 @@ 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) +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 - error("not implemented") - # 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]) + 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]) @@ -2360,16 +2364,16 @@ function gadjustsource!(gu::Source,gq::Source) end end -function gadjustsource!(gu::NamedTuple,gq::Union{Source,Field}) #where {N1, N2, T <: Real} +function gadjustsource!(gu::NamedTuple,gq::T,q::T) where T <: Union{Source,Field} qkey = :source if haskey(gu,qkey) - gadjustsource!(gu[qkey],gq) + gadjustsource!(gu[qkey],gq,q) end end -function gadjustsource!(gu::NamedTuple,gq::NamedTuple) +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]) + gadjustsource!(gu[qkey],gq[qkey],q[qkey]) end end end diff --git a/test/runtests.jl b/test/runtests.jl index 08fe6bb3..e6ed92b2 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -250,6 +250,8 @@ using TMI ũ = out.minimizer J̃,∂J̃∂ũ = fg(ũ) @test J̃ < J̃₀ + + end end @testset "sourcemap" begin @@ -259,24 +261,26 @@ using TMI 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 = zeros(γ)) - u = (; source = zerosource(γ)) + #u = (; source = zerosource(γ)) + u = (; source = zerosource(γ,logscale=true)) b = (; surface = bPO₄) # surface boundary condition - PO₄ = steadyinversion(Alu,b,γ,q=qPO₄) + PO₄true = steadyinversion(Alu,b,γ,q=qPO₄) + PO₄₀ = steadyinversion(Alu,b,γ,q=q₀) uvec = vec(u) - σq = 0.1 + σ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=qPO₄) + 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₀=qPO₄) + 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])) @@ -300,6 +304,14 @@ using TMI ũ = out.minimizer J̃,∂J̃∂ũ = fg(ũ) @test J̃ < J̃₀ + #@test J̃ < 3N # too strict if first guess is bad + + #b̃ = adjustboundarycondition(b₀,unvec(u,ũ)) # combine b₀, u + q̃ = TMI.adjustsource(q₀,unvec(u,ũ)) + + @test + + end @testset "source_and_surfacemap" begin From 5bb716276b4eb248b6db6680350af1b1d156f7fb Mon Sep 17 00:00:00 2001 From: G Jake Gebbie Date: Thu, 6 Apr 2023 15:23:14 -0400 Subject: [PATCH 35/39] `surfacetoo`: test surface and source adjustment together --- test/runtests.jl | 207 ++++++++++++++++++++++++----------------------- 1 file changed, 108 insertions(+), 99 deletions(-) diff --git a/test/runtests.jl b/test/runtests.jl index e6ed92b2..4b11a120 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -258,111 +258,120 @@ using TMI using Statistics, Interpolations - 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(γ)) - u = (; source = zerosource(γ,logscale=true)) - 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 - - #b̃ = adjustboundarycondition(b₀,unvec(u,ũ)) # combine b₀, u - q̃ = TMI.adjustsource(q₀,unvec(u,ũ)) - - @test - - + 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 "source_and_surfacemap" begin + # @testset "source_and_surfacemap" begin - using Statistics, Interpolations + # using Statistics, Interpolations - yPO₄ = readfield(TMIfile,"PO₄",γ) - bPO₄ = getsurfaceboundary(yPO₄) + # yPO₄ = readfield(TMIfile,"PO₄",γ) + # bPO₄ = getsurfaceboundary(yPO₄) + # qPO₄ = readfield(TMIfile,"qPO₄",γ) + # q₀ = 1e-2*onesource(γ) - qPO₄ = readfield(TMIfile,"qPO₄",γ) + # N = 20 + # y, W⁻, ctrue, ytrue, locs, wis = synthetic_observations(TMIversion,"PO₄",γ,N) - N = 20 - y, W⁻, ctrue, ytrue, locs, wis = synthetic_observations(TMIversion,"PO₄",γ,N) - - u = (; surface = zerosurfaceboundary(γ), source = zeros(γ)) - #b = (; surface = bPO₄) - b = (; surface = bPO₄) # assume a surface boundary condition - - PO₄ = steadyinversion(Alu,b,γ,q=qPO₄) - uvec = vec(u) - - σq = 0.1 - 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=qPO₄) - 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₀=qPO₄) + # u = (; surface = zerosurfaceboundary(γ), source = zeros(γ)) + # #b = (; surface = bPO₄) + # b = (; surface = bPO₄) # assume a surface boundary condition + + # PO₄ = steadyinversion(Alu,b,γ,q=qPO₄) + # uvec = vec(u) + + # σq = 0.1 + # 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=qPO₄) + # 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₀=qPO₄) - ϵ = 1e-3 # size of finite perturbation - ii = rand(1:length(uvec)) - 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̃₀ - end + # ϵ = 1e-3 # size of finite perturbation + # ii = rand(1:length(uvec)) + # 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̃₀ + # end @testset "unvec" begin u = (;surface = onesurfaceboundary(γ), From f6c806099660c5ad7e75efd7ec05d675eb12eca5 Mon Sep 17 00:00:00 2001 From: G Jake Gebbie Date: Thu, 6 Apr 2023 15:29:38 -0400 Subject: [PATCH 36/39] remove old source and surface test --- test/runtests.jl | 58 +++++------------------------------------------- 1 file changed, 5 insertions(+), 53 deletions(-) diff --git a/test/runtests.jl b/test/runtests.jl index 4b11a120..83b438ad 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -322,56 +322,6 @@ using TMI end end end - # @testset "source_and_surfacemap" begin - - # using Statistics, Interpolations - - # yPO₄ = readfield(TMIfile,"PO₄",γ) - # bPO₄ = getsurfaceboundary(yPO₄) - # qPO₄ = readfield(TMIfile,"qPO₄",γ) - # q₀ = 1e-2*onesource(γ) - - # N = 20 - # y, W⁻, ctrue, ytrue, locs, wis = synthetic_observations(TMIversion,"PO₄",γ,N) - - # u = (; surface = zerosurfaceboundary(γ), source = zeros(γ)) - # #b = (; surface = bPO₄) - # b = (; surface = bPO₄) # assume a surface boundary condition - - # PO₄ = steadyinversion(Alu,b,γ,q=qPO₄) - # uvec = vec(u) - - # σq = 0.1 - # 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=qPO₄) - # 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₀=qPO₄) - - # ϵ = 1e-3 # size of finite perturbation - # ii = rand(1:length(uvec)) - # 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̃₀ - # end @testset "unvec" begin u = (;surface = onesurfaceboundary(γ), @@ -401,9 +351,10 @@ using TMI @test u.east.tracer[10,10] == 3.0 end -end -# @testset "regional" begin + @testset "regional" begin + TMIversion = "nordic_201x115x46_B23" + A, Alu, γ, TMIfile, L, B = config_from_nc(TMIversion) # shellscript = TMI.pkgsrcdir("read_nc_nordic_lowresolution.sh") # run(`sh $shellscript`) @@ -414,4 +365,5 @@ end # !ispath(TMI.pkgdatadir()) && mkpath(TMI.pkgdatadir()) # mv(joinpath(pwd(),TMIfile),TMI.pkgdatadir(TMIfile),force=true) -# end + end +end From 5eb58aabda73d0228d3c13e8b0a6e79b3089ddb2 Mon Sep 17 00:00:00 2001 From: G Jake Gebbie Date: Thu, 6 Apr 2023 16:06:24 -0400 Subject: [PATCH 37/39] cartesianindex gridprops for regional case --- scripts/mat2netcdf.jl | 1 + src/config.jl | 40 ++++++++++++++++++++++++++++++---------- 2 files changed, 31 insertions(+), 10 deletions(-) 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/config.jl b/src/config.jl index 1d8a2a2a..1426d3bd 100644 --- a/src/config.jl +++ b/src/config.jl @@ -202,7 +202,7 @@ function download_matfile(TMIversion::String) println("workaround for regional Nordic Seas file") shellscript = pkgsrcdir("read_mat_nordic_201x115x46_B23.sh") run(`sh $shellscript`) - mv(joinpath(pwd(),TMIfilegz),TMIfile,force=true) + mv(joinpath(pwd(),"TMI_"*TMIversion*".mat.gz"),TMIfilegz,force=true) else !isfile(TMIfilegz) & !isfile(TMIfile) && google_download(url,pkgdatadir()) end @@ -234,9 +234,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 @@ -259,15 +275,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 From 62c776a49a8042089cfb7771dbc6c7f4800d455d Mon Sep 17 00:00:00 2001 From: G Jake Gebbie Date: Thu, 6 Apr 2023 16:23:59 -0400 Subject: [PATCH 38/39] update `config_from_mat` still redundant --- src/config.jl | 12 +++++++----- 1 file changed, 7 insertions(+), 5 deletions(-) diff --git a/src/config.jl b/src/config.jl index 1426d3bd..8eeba8a4 100644 --- a/src/config.jl +++ b/src/config.jl @@ -132,7 +132,7 @@ 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 @@ -143,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) + + # # get properties of grid + # lat,lon,depth = gridprops(TMIfile) - γ = Grid(lon,lat,depth,I,R,wet) + # γ = Grid(lon,lat,depth,I,R,wet) # need to make this optional L = circulationmatrix(TMIfile,γ) From 093e654155928fbde7ae9b4982d4764207cc8fb0 Mon Sep 17 00:00:00 2001 From: G Jake Gebbie Date: Thu, 6 Apr 2023 16:24:34 -0400 Subject: [PATCH 39/39] test nordic with mat file later try nc file --- test/runtests.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/runtests.jl b/test/runtests.jl index 83b438ad..eb8eef91 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -354,7 +354,7 @@ using TMI @testset "regional" begin TMIversion = "nordic_201x115x46_B23" - A, Alu, γ, TMIfile, L, B = config_from_nc(TMIversion) + A, Alu, γ, TMIfile, L, B = config_from_mat(TMIversion) # shellscript = TMI.pkgsrcdir("read_nc_nordic_lowresolution.sh") # run(`sh $shellscript`)