CHG: implemented most of transient inputs, need to work on the final …

…scaling for general case

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)

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#}}}

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)
 

src/core/elements.jl

@@ -395,31 +395,31 @@ 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)
 		transientinput = GetTransientInput(inputs, enum)
 
 		#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);

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#}}}

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})# {{{