From cb55693ffb60a229035a41cebb4380bde5e80e46 Mon Sep 17 00:00:00 2001 From: mmorligh Date: Sat, 14 Sep 2024 19:31:28 -0400 Subject: [PATCH] CHG: implemented most of transient inputs, need to work on the final scaling for general case --- src/core/analyses/levelsetanalysis.jl | 6 +- src/core/analyses/masstransportanalysis.jl | 12 +-- src/core/analyses/stressbalanceanalysis.jl | 9 +- src/core/elements.jl | 12 +-- src/core/inputs.jl | 95 ++++++++++++++++++---- src/core/modules.jl | 14 ++-- 6 files changed, 112 insertions(+), 36 deletions(-) diff --git a/src/core/analyses/levelsetanalysis.jl b/src/core/analyses/levelsetanalysis.jl index d04d24f..571f9dd 100644 --- a/src/core/analyses/levelsetanalysis.jl +++ b/src/core/analyses/levelsetanalysis.jl @@ -7,6 +7,8 @@ function CreateConstraints(analysis::LevelsetAnalysis,constraints::Vector{Constr #load constraints from model spclevelset = md.levelset.spclevelset + @assert size(spclevelset,1)==md.mesh.numberofvertices + @assert size(spclevelset,2)==1 count = 1 for i in 1:md.mesh.numberofvertices @@ -37,8 +39,8 @@ function UpdateElements(analysis::LevelsetAnalysis,elements::Vector{Tria}, input #Add necessary inputs to perform this analysis FetchDataToInput(md,inputs,elements,md.mask.ice_levelset,MaskIceLevelsetEnum) FetchDataToInput(md,inputs,elements,md.mask.ocean_levelset,MaskOceanLevelsetEnum) - FetchDataToInput(md,inputs,elements,md.initialization.vx./md.constants.yts,VxEnum) - FetchDataToInput(md,inputs,elements,md.initialization.vy./md.constants.yts,VyEnum) + FetchDataToInput(md,inputs,elements,md.initialization.vx,VxEnum, 1.0/md.constants.yts) + FetchDataToInput(md,inputs,elements,md.initialization.vy,VyEnum, 1.0/md.constants.yts) FetchDataToInput(md,inputs,elements,md.geometry.thickness,ThicknessEnum) FetchDataToInput(md,inputs,elements,md.geometry.surface,SurfaceEnum) diff --git a/src/core/analyses/masstransportanalysis.jl b/src/core/analyses/masstransportanalysis.jl index 42c6024..6039399 100644 --- a/src/core/analyses/masstransportanalysis.jl +++ b/src/core/analyses/masstransportanalysis.jl @@ -7,6 +7,8 @@ function CreateConstraints(analysis::MasstransportAnalysis,constraints::Vector{C #load constraints from model spcthickness = md.masstransport.spcthickness + @assert size(spcthickness,1)==md.mesh.numberofvertices + @assert size(spcthickness,2)==1 count = 1 for i in 1:md.mesh.numberofvertices @@ -39,13 +41,13 @@ function UpdateElements(analysis::MasstransportAnalysis,elements::Vector{Tria}, FetchDataToInput(md,inputs,elements,md.geometry.surface,SurfaceEnum) FetchDataToInput(md,inputs,elements,md.geometry.base,BaseEnum) FetchDataToInput(md,inputs,elements,md.geometry.bed,BedEnum) - FetchDataToInput(md,inputs,elements,md.basalforcings.groundedice_melting_rate./md.constants.yts,BasalforcingsGroundediceMeltingRateEnum) - FetchDataToInput(md,inputs,elements,md.basalforcings.floatingice_melting_rate./md.constants.yts,BasalforcingsFloatingiceMeltingRateEnum) - FetchDataToInput(md,inputs,elements,md.smb.mass_balance./md.constants.yts,SmbMassBalanceEnum) + FetchDataToInput(md,inputs,elements,md.basalforcings.groundedice_melting_rate,BasalforcingsGroundediceMeltingRateEnum, 1.0/md.constants.yts) + FetchDataToInput(md,inputs,elements,md.basalforcings.floatingice_melting_rate,BasalforcingsFloatingiceMeltingRateEnum, 1.0/md.constants.yts) + FetchDataToInput(md,inputs,elements,md.smb.mass_balance,SmbMassBalanceEnum, 1.0/md.constants.yts) FetchDataToInput(md,inputs,elements,md.mask.ice_levelset, MaskIceLevelsetEnum) FetchDataToInput(md,inputs,elements,md.mask.ocean_levelset, MaskOceanLevelsetEnum) - FetchDataToInput(md,inputs,elements,md.initialization.vx./md.constants.yts,VxEnum) - FetchDataToInput(md,inputs,elements,md.initialization.vy./md.constants.yts,VyEnum) + FetchDataToInput(md,inputs,elements,md.initialization.vx,VxEnum, 1.0/md.constants.yts) + FetchDataToInput(md,inputs,elements,md.initialization.vy,VyEnum, 1.0/md.constants.yts) return nothing end#}}} diff --git a/src/core/analyses/stressbalanceanalysis.jl b/src/core/analyses/stressbalanceanalysis.jl index d82dade..3b9a382 100644 --- a/src/core/analyses/stressbalanceanalysis.jl +++ b/src/core/analyses/stressbalanceanalysis.jl @@ -9,6 +9,11 @@ function CreateConstraints(analysis::StressbalanceAnalysis,constraints::Vector{C spcvx = md.stressbalance.spcvx spcvy = md.stressbalance.spcvy + @assert size(spcvx,1)==md.mesh.numberofvertices + @assert size(spcvx,2)==1 + @assert size(spcvy,1)==md.mesh.numberofvertices + @assert size(spcvy,2)==1 + count = 1 for i in 1:md.mesh.numberofvertices if ~isnan(spcvx[i]) @@ -44,8 +49,8 @@ function UpdateElements(analysis::StressbalanceAnalysis,elements::Vector{Tria}, FetchDataToInput(md,inputs,elements,md.geometry.thickness,ThicknessEnum) FetchDataToInput(md,inputs,elements,md.geometry.surface,SurfaceEnum) FetchDataToInput(md,inputs,elements,md.geometry.base,BaseEnum) - FetchDataToInput(md,inputs,elements,md.initialization.vx./md.constants.yts,VxEnum) - FetchDataToInput(md,inputs,elements,md.initialization.vy./md.constants.yts,VyEnum) + FetchDataToInput(md,inputs,elements,md.initialization.vx,VxEnum, 1.0/md.constants.yts) + FetchDataToInput(md,inputs,elements,md.initialization.vy,VyEnum, 1.0/md.constants.yts) FetchDataToInput(md,inputs,elements,md.mask.ice_levelset, MaskIceLevelsetEnum) FetchDataToInput(md,inputs,elements,md.mask.ocean_levelset, MaskOceanLevelsetEnum) diff --git a/src/core/elements.jl b/src/core/elements.jl index 34e7240..c7d8f3d 100644 --- a/src/core/elements.jl +++ b/src/core/elements.jl @@ -395,23 +395,23 @@ function IceVolumeAboveFloatation(element::Tria) # {{{ return area*(surface-base+min(rho_water/rho_ice*bed,0.)) end # }}} -function InputCreate(element::Tria,inputs::Inputs,data::Vector{Float64},enum::IssmEnum) #{{{ +function InputCreate(element::Tria,inputs::Inputs,data::Vector{Float64},enum::IssmEnum, scaling::Float64) #{{{ if size(data,1)==inputs.numberofelements - SetTriaInput(inputs,enum,P0Enum,element.sid,data[element.sid]) + SetTriaInput(inputs, enum, P0Enum, element.sid, scaling.*data[element.sid]) elseif size(data,1)==inputs.numberofvertices - SetTriaInput(inputs,enum,P1Enum,element.vertexids,data[element.vertexids]) + SetTriaInput(inputs, enum, P1Enum, element.vertexids, scaling.*data[element.vertexids]) else error("size ",size(data,1)," not supported for ", enum); end return nothing end #}}} -function InputCreate(element::Tria,inputs::Inputs,data::Matrix{Float64},enum::IssmEnum) #{{{ +function InputCreate(element::Tria, inputs::Inputs, data::Matrix{Float64}, enum::IssmEnum, scaling::Float64) #{{{ if size(data,1)==inputs.numberofelements+1 error("not supported yet") elseif size(data,1)==inputs.numberofvertices+1 #Extract time first - times = data[end,:] + times = data[end,:].*(365*24*3600.) #FIXME where should this conversion happen? #Create Transient Input SetTransientInput(inputs, enum, times) @@ -419,7 +419,7 @@ function InputCreate(element::Tria,inputs::Inputs,data::Matrix{Float64},enum::Is #Set values for all time slices for i in 1:length(times) - AddTimeInput(inputs, transientinput, i, P1Enum, element.vertexids, data[element.vertexids]) + AddTimeInput(inputs, transientinput, i, P1Enum, element.vertexids, scaling.*data[element.vertexids]) end else error("size ",size(data,1)," not supported for ", enum); diff --git a/src/core/inputs.jl b/src/core/inputs.jl index acd7a4c..ee3b612 100644 --- a/src/core/inputs.jl +++ b/src/core/inputs.jl @@ -10,6 +10,9 @@ mutable struct TransientInput #{{{ times::Vector{Float64} N::Int64 inputs::Vector{Input} + parameters::Parameters + currentinput::Input + currentstep::Float64 end# }}} #Inputs dataset definition @@ -20,6 +23,18 @@ mutable struct Inputs #{{{ end# }}} #Inputs functions +function Configure(inputs::Inputs, parameters::Parameters) #{{{ + + #Loop over all inputs and find transient inputs + for (key, value) in inputs.lookup + if isa(value, TransientInput) + value.parameters = parameters + end + end + + return nothing + +end#}}} function DuplicateInput(inputs::Inputs, old::IssmEnum, new::IssmEnum)#{{{ #Fetch input that needs to be copied @@ -39,20 +54,20 @@ function GetInput(inputs::Inputs, enum::IssmEnum) #{{{ end #return input - if typeof(inputs.lookup[enum])!=Input error("Input ",enum," is not an Input") end - return inputs.lookup[enum] + input = inputs.lookup[enum] + if isa(input, TransientInput) -end#}}} -function GetTransientInput(inputs::Inputs, enum::IssmEnum) #{{{ + #Find time in parameters + time = FindParam(Float64, input.parameters, TimeEnum) - #Does this input exist - if !haskey(inputs.lookup,enum) - error("Input ",enum," not found") + #Interpolate input if necessary + SetCurrentTimeInput(input, time) + input = input.currentinput end - #return input - if typeof(inputs.lookup[enum])!=TransientInput error("Input ",enum," is not a TransientInput") end - return inputs.lookup[enum] + #Make sure this is a regular input + if ~isa(input, Input) error("Input ",enum," is not an Input") end + return input end#}}} function GetInputAverageValue(input::Input) #{{{ @@ -168,23 +183,36 @@ function SetTriaInput(inputs::Inputs,enum::IssmEnum,interp::IssmEnum,indices::Ve return nothing end#}}} -function AddTimeInput(inputs::Inputs, trinput::TransientInput, index::Int64, interp::IssmEnum,indices::Vector{Int64},values::Vector{Float64}) #{{{ +#TransientInput functions +function GetTransientInput(inputs::Inputs, enum::IssmEnum) #{{{ + + #Does this input exist + if !haskey(inputs.lookup,enum) + error("Input ",enum," not found") + end + + #return input + if typeof(inputs.lookup[enum])!=TransientInput error("Input ",enum," is not a TransientInput") end + return inputs.lookup[enum] + +end#}}} +function AddTimeInput(inputs::Inputs, tr_input::TransientInput, index::Int64, interp::IssmEnum,indices::Vector{Int64},values::Vector{Float64}) #{{{ #Check index - @assert index>0 && index<=trinput.N + @assert index>0 && index<=tr_input.N #Is input already assigned? - if ~isassigned(trinput.inputs, index) + if ~isassigned(tr_input.inputs, index) @assert inputs.numberofelements > 0 if interp==P1Enum - trinput.inputs[index] = Input(trinput.enum, interp,zeros(inputs.numberofvertices), Vector{Float64}(undef,3)) + tr_input.inputs[index] = Input(tr_input.enum, interp,zeros(inputs.numberofvertices), Vector{Float64}(undef,3)) else error("not supported yet") end end #set value - trinput.inputs[index].values[indices] = values + tr_input.inputs[index].values[indices] = values return nothing end#}}} @@ -195,7 +223,7 @@ function SetTransientInput(inputs::Inputs,enum::IssmEnum,times::Vector{Float64}) #it does not exist yet, we need to create it... N = length(times) @assert N>0 - inputs.lookup[enum] = TransientInput(enum, times, length(times), Vector{Input}(undef,N)) + inputs.lookup[enum] = TransientInput(enum, times, length(times), Vector{Input}(undef,N), Parameters(), Input(enum, P0Enum ,zeros(0), Vector{Float64}(undef,0)), -1.) end #Some checks that everything is consistent @@ -205,3 +233,38 @@ function SetTransientInput(inputs::Inputs,enum::IssmEnum,times::Vector{Float64}) return nothing end#}}} +function SetCurrentTimeInput(tr_input::TransientInput, time::Float64) #{{{ + + #Binary search where time is in our series + p = searchsorted(tr_input.times, time) + + println("time is ",time, " and p is ", p) + + if p.stop==0 + #Before our time series + @assert time<=tr_input.times[1] + + #If already processed return + if(tr_input.currentstep==0.) return nothing end + + #Copy first input + tr_input.currentinput = tr_input.inputs[0] + tr_input.currentstep = 0. + + elseif p.start==length(tr_input.times)+1 + #After the end of our time series + @assert time>=tr_input.times[end] + + #If already processed return + if(tr_input.currentstep==Float64(tr_input.N)) return nothing end + + #Copy last input + tr_input.currentinput = tr_input.inputs[end] + tr_input.currentstep = Float64(tr_input.N) + else + #General case... + error("Not implemented yet, see C++ code") + end + + return nothing +end#}}} diff --git a/src/core/modules.jl b/src/core/modules.jl index ab57dd1..f9bec2c 100644 --- a/src/core/modules.jl +++ b/src/core/modules.jl @@ -1,13 +1,13 @@ #Model Processor and Core I/O -function FetchDataToInput(md::model, inputs::Inputs, elements::Vector{Tria}, data::Vector{Float64}, enum::IssmEnum) #{{{ +function FetchDataToInput(md::model, inputs::Inputs, elements::Vector{Tria}, data::Vector{Float64}, enum::IssmEnum, scaling::Float64=1.) #{{{ for i in 1:length(elements) - InputCreate(elements[i],inputs,data,enum) + InputCreate(elements[i], inputs, data, enum, scaling) end return nothing end#}}} -function FetchDataToInput(md::model, inputs::Inputs, elements::Vector{Tria}, data::Matrix{Float64}, enum::IssmEnum) #{{{ +function FetchDataToInput(md::model, inputs::Inputs, elements::Vector{Tria}, data::Matrix{Float64}, enum::IssmEnum, scaling::Float64=1.) #{{{ for i in 1:length(elements) - InputCreate(elements[i],inputs,data,enum) + InputCreate(elements[i], inputs, data, enum, scaling) end return nothing end#}}} @@ -24,7 +24,7 @@ function ModelProcessor(md::model, solutionstring::Symbol) #{{{ CreateElements(elements, md) CreateVertices(vertices, md) CreateParameters(parameters, md, solutionstring) - CreateInputs(inputs,elements, md) + CreateInputs(inputs, elements, md) if solutionstring===:TransientSolution UpdateParameters(TransientAnalysis(), parameters, md) end @@ -185,10 +185,14 @@ end# }}} #Other modules function ConfigureObjectx(elements::Vector{Tria}, nodes::Vector{Node}, vertices::Vector{Vertex}, parameters::Parameters, inputs::Inputs, analysis::Int64) #{{{ + #Configure elements for i in 1:length(elements) Configure(elements[i], nodes, vertices, parameters, inputs, analysis) end + #Configure inputs + Configure(inputs, parameters) + return nothing end# }}} function CreateNodalConstraintsx(nodes::Vector{Node})# {{{