From 849a8bf97314ae13f81a6f731f6810cefde1e3ee Mon Sep 17 00:00:00 2001 From: Cheng Gong Date: Thu, 6 Jun 2024 09:08:46 -0400 Subject: [PATCH] add IceVolume and IceVolumeAboveFloatation to the cost function list. Add md.inversion.dependent_string to choose different cost functions --- src/core/control.jl | 83 +++-- src/core/costfunctions.jl | 29 -- src/core/elements.jl | 625 +++++++++++++++++++++++--------------- src/core/femmodel.jl | 4 + src/core/gauss.jl | 165 ++++++++-- src/core/inputs.jl | 143 ++++----- src/core/modules.jl | 38 ++- src/core/parameters.jl | 13 + src/core/solve.jl | 6 +- src/usr/classes.jl | 3 +- 10 files changed, 698 insertions(+), 411 deletions(-) diff --git a/src/core/control.jl b/src/core/control.jl index 06a51e0..19e8549 100644 --- a/src/core/control.jl +++ b/src/core/control.jl @@ -1,43 +1,72 @@ using Enzyme -Enzyme.API.looseTypeAnalysis!(false) -Enzyme.API.strictAliasing!(false) Enzyme.API.typeWarning!(false) Enzyme.Compiler.RunAttributor[] = false using Optimization, OptimizationOptimJL - -function Control_Core(md::model, femmodel::FemModel) #{{{ - # solve for optimization - # TODO: just a first try, need to add all the features - α = md.inversion.independent - ∂J_∂α = zero(α) - n = length(α) - # use user defined grad, errors! - #optprob = OptimizationFunction(costfunction, Optimization.AutoEnzyme(), grad=computeGradient(∂J_∂α, α, femmodel)) - optprob = OptimizationFunction(costfunction, Optimization.AutoEnzyme()) - prob = Optimization.OptimizationProblem(optprob, α, femmodel, lb=md.inversion.min_parameters, ub=md.inversion.max_parameters) - sol = Optimization.solve(prob, Optim.LBFGS()) - - independent_enum = StringToEnum(md.inversion.independent_string) - InputUpdateFromVectorx(femmodel, sol.u, independent_enum, VertexSIdEnum) - RequestedOutputsx(femmodel, [independent_enum]) -end#}}} -function computeGradient(md::model, femmodel::FemModel) #{{{ +function Control_Core(md::model, femmodel::FemModel, solution::Symbol) #{{{ #independent variable α = md.inversion.independent #initialize derivative as 0 ∂J_∂α = zero(α) - # Compute Gradient - computeGradient(∂J_∂α, α, femmodel) - - #Put gradient in results - InputUpdateFromVectorx(femmodel, ∂J_∂α, GradientEnum, VertexSIdEnum) - RequestedOutputsx(femmodel, [GradientEnum]) + if solution ===:grad + # only compute the gradient + ComputeGradient(∂J_∂α, α, femmodel) + #Put gradient in results + InputUpdateFromVectorx(femmodel, ∂J_∂α, GradientEnum, VertexSIdEnum) + RequestedOutputsx(femmodel, [GradientEnum]) + else + # optimization + # use user defined grad, errors! + #optprob = OptimizationFunction(costfunction, Optimization.AutoEnzyme(), grad=computeGradient(∂J_∂α, α, femmodel)) + optprob = OptimizationFunction(costfunction, Optimization.AutoEnzyme()) + prob = Optimization.OptimizationProblem(optprob, α, femmodel, lb=md.inversion.min_parameters, ub=md.inversion.max_parameters) + sol = Optimization.solve(prob, Optim.LBFGS()) + independent_enum = StringToEnum(md.inversion.independent_string) + InputUpdateFromVectorx(femmodel, sol.u, independent_enum, VertexSIdEnum) + RequestedOutputsx(femmodel, [independent_enum]) + end end#}}} -function computeGradient(∂J_∂α::Vector{Float64}, α::Vector{Float64}, femmodel::FemModel) #{{{ +function ComputeGradient(∂J_∂α::Vector{Float64}, α::Vector{Float64}, femmodel::FemModel) #{{{ # zero ALL depth of the model, make sure we get correct gradient dfemmodel = Enzyme.Compiler.make_zero(Base.Core.Typeof(femmodel), IdDict(), femmodel) # compute the gradient autodiff(Enzyme.Reverse, costfunction, Duplicated(α, ∂J_∂α), Duplicated(femmodel,dfemmodel)) end#}}} +function costfunction(α::Vector{Float64}, femmodel::FemModel) #{{{ + # get the md.inversion.control_string + control_string = FindParam(String, femmodel.parameters, InversionControlParametersEnum) + # get the Enum + controlvar_enum = StringToEnum(control_string) + if isnothing(controlvar_enum) + error(control_string, " is not defined in DJUICE, therefore the derivative with respect to ", control_string, " is meaningless") + end + + # get the cost function list from md.inversion.dependent_string + cost_list = FindParam(Vector{String}, femmodel.parameters, InversionCostFunctionsEnum) + cost_enum_list = map(StringToEnum, cost_list) + + # compute cost function + # TODO: loop through all controls with respect to all the components in the cost function + CostFunctionx(femmodel, α, controlvar_enum, VertexSIdEnum, cost_enum_list) +end#}}} +function CostFunctionx(femmodel::FemModel, α::Vector{Float64}, controlvar_enum::IssmEnum, SId_enum::IssmEnum, cost_enum_list::Vector{IssmEnum}) #{{{ + #Update FemModel accordingly + InputUpdateFromVectorx(femmodel, α, controlvar_enum, SId_enum) + + #solve PDE + analysis = StressbalanceAnalysis() + Core(analysis, femmodel) + + #update all values of the cost functions + RequestedOutputsx(femmodel, cost_enum_list) + + #Compute cost function + J = 0.0 + for i in 1:length(cost_enum_list) + J += femmodel.results[end-i+1].value + end + + #return cost function + return J +end#}}} diff --git a/src/core/costfunctions.jl b/src/core/costfunctions.jl index 8e1fbd0..825d447 100644 --- a/src/core/costfunctions.jl +++ b/src/core/costfunctions.jl @@ -1,33 +1,4 @@ -function costfunction(α::Vector{Float64}, femmodel::FemModel) #{{{ - # get the md.inversion.control_string - control_string = FindParam(String, femmodel.parameters, InversionControlParametersEnum) - # get the Enum - controlvar_enum = StringToEnum(control_string) - if isnothing(controlvar_enum) - error(control_string, " is not defined in DJUICE, therefore the derivative with respect to ", control_string, " is meaningless") - end - # compute cost function - # TODO: loop through all controls with respect to all the components in the cost function - costfunction(α, femmodel, controlvar_enum, VertexSIdEnum) -end#}}} -function costfunction(α::Vector{Float64}, femmodel::FemModel, controlvar_enum::IssmEnum, SId_enum::IssmEnum) #{{{ - #Update FemModel accordingly - InputUpdateFromVectorx(femmodel, α, controlvar_enum, SId_enum) - - #solve PDE - analysis = StressbalanceAnalysis() - Core(analysis, femmodel) - - #Compute cost function - J = SurfaceAbsVelMisfitx(femmodel) - - #J += ControlVariableAbsGradientx(femmodel, α, controlvar_enum) - - #return cost function - return J -end#}}} - # The misfit functions function SurfaceAbsVelMisfitx(femmodel::FemModel) #{{{ diff --git a/src/core/elements.jl b/src/core/elements.jl index 1514889..ce08350 100644 --- a/src/core/elements.jl +++ b/src/core/elements.jl @@ -28,17 +28,6 @@ function Tria(sid::Int64, pid::Int64, vertexids::Vector{Int64}) #{{{ end #}}} #Element functions -function InputCreate(element::Tria,inputs::Inputs,data::Vector{Float64},enum::IssmEnum) #{{{ - if size(data,1)==inputs.numberofelements - SetTriaInput(inputs,enum,P0Enum,element.sid,data[element.sid]) - elseif size(data,1)==inputs.numberofvertices - SetTriaInput(inputs,enum,P1Enum,element.vertexids,data[element.vertexids]) - else - error("size ",size(data,1)," not supported for ", enum); - end - - return nothing -end #}}} function AddInput(element::Tria,inputenum::IssmEnum,data::Vector{Float64},interpolation::IssmEnum) #{{{ if interpolation==P1Enum @assert length(data)==3 @@ -49,40 +38,6 @@ function AddInput(element::Tria,inputenum::IssmEnum,data::Vector{Float64},interp return nothing end #}}} -function InputUpdateFromVector(element::Tria, vector::Vector{Float64}, enum::IssmEnum, layout::IssmEnum) #{{{ - - lidlist = element.vertexids - data = Vector{Float64}(undef, 3) - - if(layout==VertexSIdEnum) - for i in 1:3 - data[i] = vector[element.vertices[i].sid] - @assert isfinite(data[i]) - end - SetTriaInput(element.inputs, enum, P1Enum, lidlist, data) - else - error("layout ", layout, " not supported yet"); - end - - return nothing -end #}}} -function Update(element::Tria, inputs::Inputs, index::Int64, md::model, finiteelement::IssmEnum) #{{{ - - if finiteelement==P1Enum - numnodes = 3 - nodeids = Vector{Int64}(undef,numnodes) - nodeids[1] = md.mesh.elements[index,1] - nodeids[2] = md.mesh.elements[index,2] - nodeids[3] = md.mesh.elements[index,3] - - push!(element.nodes_ids_list, nodeids) - push!(element.nodes_list, Vector{Node}(undef, numnodes)) - else - error("not supported yet") - end - - return nothing -end #}}} function Configure(element::Tria,nodes::Vector{Node},vertices::Vector{Vertex},parameters::Parameters,inputs::Inputs,index::Int64) # {{{ #Configure vertices @@ -104,6 +59,26 @@ function Configure(element::Tria,nodes::Vector{Node},vertices::Vector{Vertex},pa return nothing end # }}} +function CharacteristicLength(element::Tria) #{{{ + + return sqrt(2*GetArea(element)) +end#}}} +function FindParam(::Type{T}, element::Tria, enum::IssmEnum) where T # {{{ + + return FindParam(T, element.parameters, enum) + +end # }}} +function GetArea(element::Tria)#{{{ + + #Get xyz list + xyz_list = GetVerticesCoordinates(element.vertices) + x1 = xyz_list[1,1]; y1 = xyz_list[1,2] + x2 = xyz_list[2,1]; y2 = xyz_list[2,2] + x3 = xyz_list[3,1]; y3 = xyz_list[3,2] + + @assert x2*y3 - y2*x3 + x1*y2 - y1*x2 + x3*y1 - y3*x1>0 + return (x2*y3 - y2*x3 + x1*y2 - y1*x2 + x3*y1 - y3*x1)/2 +end#}}} function GetDofList(element::Tria,setenum::IssmEnum) # {{{ #Define number of nodes @@ -126,6 +101,174 @@ function GetDofList(element::Tria,setenum::IssmEnum) # {{{ return doflist end # }}} +function GetElementSizes(element::Tria)#{{{ + + #Get xyz list + xyz_list = GetVerticesCoordinates(element.vertices) + x1 = xyz_list[1,1]; y1 = xyz_list[1,2] + x2 = xyz_list[2,1]; y2 = xyz_list[2,2] + x3 = xyz_list[3,1]; y3 = xyz_list[3,2] + + xmin=xyz_list[1,1]; xmax=xyz_list[1,1]; + ymin=xyz_list[1,2]; ymax=xyz_list[1,2]; + + for i in [2,3] + if(xyz_list[i,1]xmax) xmax=xyz_list[i,1] end + if(xyz_list[i,2]ymax) ymax=xyz_list[i,2] end + end + + hx = xmax-xmin + hy = ymax-ymin + hz = 0. + + return hx, hy, hz +end#}}} +function GetFractionGeometry(element::Tria, gl::Vector{Float64})#{{{ + trapezeisnegative = true + # Be sure that values are not zero + if(gl[1]==0.); gl[1]=gl[1]+eps(Float64) end + if(gl[2]==0.); gl[2]=gl[2]+eps(Float64) end + if(gl[3]==0.); gl[3]=gl[3]+eps(Float64) end + + # Check that not all nodes are positive or negative: + if(gl[1]>0. && gl[2]>0. && gl[3]>0.) # All positive + point=1 + f1=1. + f2=1. + elseif (gl[1]<0 && gl[2]<0 && gl[3]<0) #All negative + point=1 + f1=0. + f2=0. + else + if(gl[1]*gl[2]*gl[3]<0) + trapezeisnegative=false # no matter what configuration, there has to be two positive vertices, which means the trapeze is positive. + end + + if (gl[1]*gl[2]>0) # Nodes 1 and 2 are similar, so points must be found on segment 1-3 and 2-3 + point=3 + f1=gl[3]/(gl[3]-gl[1]) + f2=gl[3]/(gl[3]-gl[2]) + elseif (gl[2]*gl[3]>0) # Nodes 2 and 3 are similar, so points must be found on segment 1-2 and 1-3 + point=1 + f1=gl[1]/(gl[1]-gl[2]) + f2=gl[1]/(gl[1]-gl[3]) + elseif (gl[1]*gl[3]>0) # Nodes 1 and 3 are similar, so points must be found on segment 2-1 and 2-3 + point=2 + f1=gl[2]/(gl[2]-gl[3]) + f2=gl[2]/(gl[2]-gl[1]) + else + error("case not possible"); + end + end + if (trapezeisnegative) + phi=1-f1*f2 + else + phi=f1*f2 + end + + # Compute the weights + gauss = GaussTria(point, f1, f2, !trapezeisnegative, 2) + numnodes = 3 + + total_weight = 0.0 + weights = Vector{Float64}(undef, numnodes) + loadweights_g = Vector{Float64}(undef,numnodes) + + for ig in 1:gauss.numgauss + NodalFunctions(element, loadweights_g, gauss, ig, P1Enum) + for i in 1:numnodes + weights[i] += loadweights_g[i]*gauss.weights[ig] + end + total_weight += gauss.weights[ig] + end + # Normalizing to phi such that weights provide coefficients for integration over subelement (for averaging:phi*weights) + if (total_weight>0.) + for i in 1:numnodes + weights[i] /= (total_weight/phi); + end + else + for i in 1:numnodes + weights[i]=0. + end + end + + return weights, phi +end#}}} +function GetGroundedPortion(element::Tria, xyz_list::Matrix{Float64}) #{{{ + + level = Vector{Float64}(undef,3) + GetInputListOnVertices!(element, level, MaskOceanLevelsetEnum) + + #Be sure that values are not zero + epsilon = 1.e-15 + for i in 1:3 + if(level[i]==0.) level[i]=level[i]+epsilon end + end + + if level[1]>0 && level[2]>0 && level[3]>0 + #Completely grounded + phi = 1.0 + elseif level[1]<0 && level[2]<0 && level[3]<0 + #Completely floating + phi = 0.0 + else + #Partially floating, + if(level[1]*level[2]>0) #Nodes 0 and 1 are similar, so points must be found on segment 0-2 and 1-2 + s1=level[3]/(level[3]-level[2]); + s2=level[3]/(level[3]-level[1]); + elseif(level[2]*level[3]>0) #Nodes 1 and 2 are similar, so points must be found on segment 0-1 and 0-2 + s1=level[1]/(level[1]-level[2]); + s2=level[1]/(level[1]-level[3]); + elseif(level[1]*level[3]>0) #Nodes 0 and 2 are similar, so points must be found on segment 1-0 and 1-2 + s1=level[2]/(level[2]-level[1]); + s2=level[2]/(level[2]-level[3]); + else + error("not supposed to be here...") + end + + if(level[1]*level[2]*level[3]>0) + phi = s1*s2 + else + phi = (1-s1*s2) + end + end + + return phi +end#}}} +function GetIcefrontCoordinates!(element::Tria, xyz_front::Matrix{Float64}, xyz_list::Matrix{Float64}, levelsetenum::IssmEnum) #{{{ + + #Intermediaries + level = Vector{Float64}(undef,3) + indicesfront = Vector{Int64}(undef,3) + + #Recover value of levelset for all vertices + GetInputListOnVertices!(element, level, levelsetenum) + + #Get nodes where there is no ice + num_frontnodes = 0 + for i in 1:3 + if(level[i]>=0.) + num_frontnodes += 1 + indicesfront[num_frontnodes] = i + end + end + @assert num_frontnodes==2 + + #Arrange order of frontnodes such that they are oriented counterclockwise + NUMVERTICES = 3 + if((NUMVERTICES+indicesfront[1]-indicesfront[2])%NUMVERTICES != NUMVERTICES-1) + index=indicesfront[1] + indicesfront[1]=indicesfront[2] + indicesfront[2]=index + end + + #Return nodes + xyz_front[1,:]=xyz_list[indicesfront[1],:] + xyz_front[2,:]=xyz_list[indicesfront[2],:] + return nothing +end#}}} function GetInput(element::Tria,enum::IssmEnum) # {{{ input = GetInput(element.inputs,enum) @@ -146,6 +289,19 @@ function GetInputListOnNodes!(element::Tria, vector::Vector{Float64}, enum::Issm return nothing end # }}} +function GetInputListOnVertices!(element::Tria, vector::Vector{Float64}, enum::IssmEnum) # {{{ + + #Get Input first + input = GetInput(element, enum) + + #Get value at each vertex (i.e. P1 Nodes) + gauss=GaussTria(P1Enum) + for i in 1:gauss.numgauss + vector[i] = GetInputValue(input, gauss, i) + end + + return nothing +end # }}} function GetInputValue(element::Tria, vector::Vector{Float64}, gauss::GaussTria, ig::Int64, enum::IssmEnum) # {{{ # Allocate basis functions @@ -164,79 +320,105 @@ function GetInputValue(element::Tria, vector::Vector{Float64}, gauss::GaussTria, return value end # }}} -function GetInputListOnVertices!(element::Tria, vector::Vector{Float64}, enum::IssmEnum) # {{{ +function GetXcoord(element::Tria, xyz_list::Matrix{Float64}, gauss::GaussTria, ig::Int64) #{{{ - #Get Input first - input = GetInput(element, enum) + # create a list of x + x_list = Vector{Float64}(undef,3) - #Get value at each vertex (i.e. P1 Nodes) - gauss=GaussTria(P1Enum) - for i in 1:gauss.numgauss - vector[i] = GetInputValue(input, gauss, i) + for i in 1:3 + x_list[i] = xyz_list[i,1] end - return nothing -end # }}} -function FindParam(::Type{T}, element::Tria, enum::IssmEnum) where T # {{{ + # Get value at gauss point + return GetInputValue(element, x_list, gauss, ig, P1Enum); +end#}}} +function GetYcoord(element::Tria, xyz_list::Matrix{Float64}, gauss::GaussTria, ig::Int64) #{{{ - return FindParam(T, element.parameters, enum) + # create a list of y + y_list = Vector{Float64}(undef,3) -end # }}} -function InputServe!(element::Tria,input::Input) # {{{ + for i in 1:3 + y_list[i] = xyz_list[i,2] + end - if input.interp==P0Enum - input.element_values[1] = input.values[element.sid] - elseif input.interp==P1Enum - for i in 1:3 - input.element_values[i] = input.values[element.vertices[i].sid] - end + # Get value at gauss point + return GetInputValue(element, y_list, gauss, ig, P1Enum); +end#}}} +function IceVolume(element::Tria) # {{{ + + if (!IsIceInElement(element)); return 0.0 end + + lsf = Vector{Float64}(undef,3) + GetInputListOnVertices!(element, lsf, MaskIceLevelsetEnum) + + # partially ice-covered element + if (lsf[1]*lsf[2]<=0 || lsf[1]*lsf[3]<=0 || lsf[2]*lsf[3]<=0) + surfaces = Vector{Float64}(undef,3) + Hice = Vector{Float64}(undef,3) + bases = Vector{Float64}(undef,3) + GetInputListOnVertices!(element, surfaces, SurfaceEnum) + GetInputListOnVertices!(element, bases, BaseEnum) + weights, phi = GetFractionGeometry(element, lsf) + Hice = surfaces - bases + Haverage = sum(weights.*Hice)/phi + area_basetot = GetArea(element) + area_base = phi*area_basetot else - error("interpolation ",input.interp," not supported yet") + area_base = GetArea(element) + surface_input = GetInput(element, SurfaceEnum) + base_input = GetInput(element, BaseEnum) + + surface = GetInputAverageValue(surface_input) + base = GetInputAverageValue(base_input) + Haverage = surface - base end - return nothing + return area_base*Haverage end # }}} -function GetGroundedPortion(element::Tria, xyz_list::Matrix{Float64}) #{{{ +function IceVolumeAboveFloatation(element::Tria) # {{{ - level = Vector{Float64}(undef,3) - GetInputListOnVertices!(element, level, MaskOceanLevelsetEnum) + if (!IsIceInElement(element) || IsAllFloating(element)); return 0.0 end - #Be sure that values are not zero - epsilon = 1.e-15 - for i in 1:3 - if(level[i]==0.) level[i]=level[i]+epsilon end + rho_ice = FindParam(Float64, element, MaterialsRhoIceEnum) + rho_water = FindParam(Float64, element, MaterialsRhoSeawaterEnum) + + area = GetArea(element) + + surface_input = GetInput(element, SurfaceEnum) + base_input = GetInput(element, BaseEnum) + bed_input = GetInput(element, BedEnum) + + surface = GetInputAverageValue(surface_input) + base = GetInputAverageValue(base_input) + bed = GetInputAverageValue(bed_input) + + return area*(surface-base+min(rho_water/rho_ice*bed,0.)) +end # }}} +function InputCreate(element::Tria,inputs::Inputs,data::Vector{Float64},enum::IssmEnum) #{{{ + if size(data,1)==inputs.numberofelements + SetTriaInput(inputs,enum,P0Enum,element.sid,data[element.sid]) + elseif size(data,1)==inputs.numberofvertices + SetTriaInput(inputs,enum,P1Enum,element.vertexids,data[element.vertexids]) + else + error("size ",size(data,1)," not supported for ", enum); end - if level[1]>0 && level[2]>0 && level[3]>0 - #Completely grounded - phi = 1.0 - elseif level[1]<0 && level[2]<0 && level[3]<0 - #Completely floating - phi = 0.0 - else - #Partially floating, - if(level[1]*level[2]>0) #Nodes 0 and 1 are similar, so points must be found on segment 0-2 and 1-2 - s1=level[3]/(level[3]-level[2]); - s2=level[3]/(level[3]-level[1]); - elseif(level[2]*level[3]>0) #Nodes 1 and 2 are similar, so points must be found on segment 0-1 and 0-2 - s1=level[1]/(level[1]-level[2]); - s2=level[1]/(level[1]-level[3]); - elseif(level[1]*level[3]>0) #Nodes 0 and 2 are similar, so points must be found on segment 1-0 and 1-2 - s1=level[2]/(level[2]-level[1]); - s2=level[2]/(level[2]-level[3]); - else - error("not supposed to be here...") - end + return nothing +end #}}} +function InputServe!(element::Tria,input::Input) # {{{ - if(level[1]*level[2]*level[3]>0) - phi = s1*s2 - else - phi = (1-s1*s2) + if input.interp==P0Enum + input.element_values[1] = input.values[element.sid] + elseif input.interp==P1Enum + for i in 1:3 + input.element_values[i] = input.values[element.vertices[i].sid] end + else + error("interpolation ",input.interp," not supported yet") end - return phi -end#}}} + return nothing +end # }}} function InputUpdateFromSolutionOneDof(element::Tria, ug::Vector{Float64},enum::IssmEnum) #{{{ #Get dofs for this finite element doflist = GetDofList(element,GsetEnum) @@ -252,6 +434,33 @@ function InputUpdateFromSolutionOneDof(element::Tria, ug::Vector{Float64},enum:: return nothing end#}}} +function InputUpdateFromVector(element::Tria, vector::Vector{Float64}, enum::IssmEnum, layout::IssmEnum) #{{{ + + lidlist = element.vertexids + data = Vector{Float64}(undef, 3) + + if(layout==VertexSIdEnum) + for i in 1:3 + data[i] = vector[element.vertices[i].sid] + @assert isfinite(data[i]) + end + SetTriaInput(element.inputs, enum, P1Enum, lidlist, data) + else + error("layout ", layout, " not supported yet"); + end + + return nothing +end #}}} +function IsAllFloating(element::Tria) #{{{ + input = GetInput(element, MaskOceanLevelsetEnum) + + # by default use subelment migration + if GetInputMax() <= 0. + return true; + else + return false + end +end#}}} function IsIcefront(element::Tria) #{{{ level = Vector{Float64}(undef,3) @@ -280,90 +489,6 @@ function IsIceInElement(element::Tria) #{{{ end end#}}} -function GetIcefrontCoordinates!(element::Tria, xyz_front::Matrix{Float64}, xyz_list::Matrix{Float64}, levelsetenum::IssmEnum) #{{{ - - #Intermediaries - level = Vector{Float64}(undef,3) - indicesfront = Vector{Int64}(undef,3) - - #Recover value of levelset for all vertices - GetInputListOnVertices!(element, level, levelsetenum) - - #Get nodes where there is no ice - num_frontnodes = 0 - for i in 1:3 - if(level[i]>=0.) - num_frontnodes += 1 - indicesfront[num_frontnodes] = i - end - end - @assert num_frontnodes==2 - - #Arrange order of frontnodes such that they are oriented counterclockwise - NUMVERTICES = 3 - if((NUMVERTICES+indicesfront[1]-indicesfront[2])%NUMVERTICES != NUMVERTICES-1) - index=indicesfront[1] - indicesfront[1]=indicesfront[2] - indicesfront[2]=index - end - - #Return nodes - xyz_front[1,:]=xyz_list[indicesfront[1],:] - xyz_front[2,:]=xyz_list[indicesfront[2],:] - return nothing -end#}}} -function GetArea(element::Tria)#{{{ - - #Get xyz list - xyz_list = GetVerticesCoordinates(element.vertices) - x1 = xyz_list[1,1]; y1 = xyz_list[1,2] - x2 = xyz_list[2,1]; y2 = xyz_list[2,2] - x3 = xyz_list[3,1]; y3 = xyz_list[3,2] - - @assert x2*y3 - y2*x3 + x1*y2 - y1*x2 + x3*y1 - y3*x1>0 - return (x2*y3 - y2*x3 + x1*y2 - y1*x2 + x3*y1 - y3*x1)/2 -end#}}} -function GetElementSizes(element::Tria)#{{{ - - #Get xyz list - xyz_list = GetVerticesCoordinates(element.vertices) - x1 = xyz_list[1,1]; y1 = xyz_list[1,2] - x2 = xyz_list[2,1]; y2 = xyz_list[2,2] - x3 = xyz_list[3,1]; y3 = xyz_list[3,2] - - xmin=xyz_list[1,1]; xmax=xyz_list[1,1]; - ymin=xyz_list[1,2]; ymax=xyz_list[1,2]; - - for i in [2,3] - if(xyz_list[i,1]xmax) xmax=xyz_list[i,1] end - if(xyz_list[i,2]ymax) ymax=xyz_list[i,2] end - end - - hx = xmax-xmin - hy = ymax-ymin - hz = 0. - - return hx, hy, hz -end#}}} -function NormalSection(element::Tria, xyz_front::Matrix{Float64}) #{{{ - - #Build output pointing vector - nx = xyz_front[2,2] - xyz_front[1,2] - ny = -xyz_front[2,1] + xyz_front[1,1] - - #normalize - norm = sqrt(nx^2 + ny^2) - nx = nx/norm - ny = ny/norm - - return nx, ny -end#}}} -function CharacteristicLength(element::Tria) #{{{ - - return sqrt(2*GetArea(element)) -end#}}} function MigrateGroundingLine(element::Tria) #{{{ h = Vector{Float64}(undef,3) @@ -471,56 +596,69 @@ function MovingFrontalVelocity(element::Tria) #{{{ return nothing end#}}} -function GetXcoord(element::Tria, xyz_list::Matrix{Float64}, gauss::GaussTria, ig::Int64) #{{{ +function NormalSection(element::Tria, xyz_front::Matrix{Float64}) #{{{ - # create a list of x - x_list = Vector{Float64}(undef,3) + #Build output pointing vector + nx = xyz_front[2,2] - xyz_front[1,2] + ny = -xyz_front[2,1] + xyz_front[1,1] - for i in 1:3 - x_list[i] = xyz_list[i,1] - end + #normalize + norm = sqrt(nx^2 + ny^2) + nx = nx/norm + ny = ny/norm - # Get value at gauss point - return GetInputValue(element, x_list, gauss, ig, P1Enum); + return nx, ny end#}}} -function GetYcoord(element::Tria, xyz_list::Matrix{Float64}, gauss::GaussTria, ig::Int64) #{{{ +function Update(element::Tria, inputs::Inputs, index::Int64, md::model, finiteelement::IssmEnum) #{{{ - # create a list of y - y_list = Vector{Float64}(undef,3) + if finiteelement==P1Enum + numnodes = 3 + nodeids = Vector{Int64}(undef,numnodes) + nodeids[1] = md.mesh.elements[index,1] + nodeids[2] = md.mesh.elements[index,2] + nodeids[3] = md.mesh.elements[index,3] - for i in 1:3 - y_list[i] = xyz_list[i,2] + push!(element.nodes_ids_list, nodeids) + push!(element.nodes_list, Vector{Node}(undef, numnodes)) + else + error("not supported yet") end - # Get value at gauss point - return GetInputValue(element, y_list, gauss, ig, P1Enum); -end#}}} + return nothing +end #}}} #Finite Element stuff -function JacobianDeterminant(xyz_list::Matrix{Float64}, gauss::GaussTria) #{{{ +function GaussTria(element::Tria, xyz_list::Matrix{Float64}, xyz_list_front::Matrix{Float64}, order::Int64) #{{{ - #Get Jacobian Matrix - J = Jacobian(xyz_list) + area_coordinates = Matrix{Float64}(undef,2,3) + GetAreaCoordinates!(element, area_coordinates, xyz_list_front, xyz_list) - #Get its determinant - Jdet = Matrix2x2Determinant(J) + return GaussTria(area_coordinates, order) +end# }}} +function GetAreaCoordinates!(element::Tria, area_coordinates::Matrix{Float64}, xyz_zero::Matrix{Float64}, xyz_list::Matrix{Float64})#{{{ - #check and return - if(Jdet<0) error("negative Jacobian Determinant") end - return Jdet + numpoints = size(area_coordinates,1) + area = GetArea(element) -end#}}} -function JacobianDeterminantSurface(xyz_list::Matrix{Float64}, gauss::GaussTria) #{{{ + #Copy original xyz_list + xyz_bis=copy(xyz_list) + for i in 1:numpoints + for j in 1:3 - x1 = xyz_list[1,1]; y1 = xyz_list[1,2] - x2 = xyz_list[2,1]; y2 = xyz_list[2,2] - Jdet = .5*sqrt((x2-x1)^2 + (y2-y1)^2) + #Change appropriate line + xyz_bis[j,:] = xyz_zero[i,:] - #check and return - if(Jdet<0) error("negative Jacobian Determinant") end - return Jdet + #Compute area fraction + area_portion=abs(xyz_bis[2,1]*xyz_bis[3,2] - xyz_bis[2,2]*xyz_bis[3,1] + xyz_bis[1,1]*xyz_bis[2,2] - xyz_bis[1,2]*xyz_bis[2,1] + xyz_bis[3,1]*xyz_bis[1,2] - xyz_bis[3,2]*xyz_bis[1,1])/2 + area_coordinates[i,j] = area_portion/area -end#}}} + #reinitialize xyz_list + xyz_bis[j,:] = xyz_list[j,:] + end + end + + return nothing +end #}}} function Jacobian(xyz_list::Matrix{Float64}) #{{{ J = Matrix{Float64}(undef,2,2) @@ -539,6 +677,30 @@ function Jacobian(xyz_list::Matrix{Float64}) #{{{ return J end#}}} +function JacobianDeterminant(xyz_list::Matrix{Float64}, gauss::GaussTria) #{{{ + + #Get Jacobian Matrix + J = Jacobian(xyz_list) + + #Get its determinant + Jdet = Matrix2x2Determinant(J) + + #check and return + if(Jdet<0) error("negative Jacobian Determinant") end + return Jdet + +end#}}} +function JacobianDeterminantSurface(xyz_list::Matrix{Float64}, gauss::GaussTria) #{{{ + + x1 = xyz_list[1,1]; y1 = xyz_list[1,2] + x2 = xyz_list[2,1]; y2 = xyz_list[2,2] + Jdet = .5*sqrt((x2-x1)^2 + (y2-y1)^2) + + #check and return + if(Jdet<0) error("negative Jacobian Determinant") end + return Jdet + +end#}}} function JacobianInvert(xyz_list::Matrix{Float64}, gauss::GaussTria) #{{{ #Get Jacobian matrix @@ -616,34 +778,3 @@ function NumberofNodesTria(finiteelement) #{{{ return nothing end#}}} -function GaussTria(element::Tria, xyz_list::Matrix{Float64}, xyz_list_front::Matrix{Float64}, order::Int64) #{{{ - - area_coordinates = Matrix{Float64}(undef,2,3) - GetAreaCoordinates!(element, area_coordinates, xyz_list_front, xyz_list) - - return GaussTria(area_coordinates, order) -end# }}} -function GetAreaCoordinates!(element::Tria, area_coordinates::Matrix{Float64}, xyz_zero::Matrix{Float64}, xyz_list::Matrix{Float64})#{{{ - - numpoints = size(area_coordinates,1) - area = GetArea(element) - - #Copy original xyz_list - xyz_bis=copy(xyz_list) - for i in 1:numpoints - for j in 1:3 - - #Change appropriate line - xyz_bis[j,:] = xyz_zero[i,:] - - #Compute area fraction - area_portion=abs(xyz_bis[2,1]*xyz_bis[3,2] - xyz_bis[2,2]*xyz_bis[3,1] + xyz_bis[1,1]*xyz_bis[2,2] - xyz_bis[1,2]*xyz_bis[2,1] + xyz_bis[3,1]*xyz_bis[1,2] - xyz_bis[3,2]*xyz_bis[1,1])/2 - area_coordinates[i,j] = area_portion/area - - #reinitialize xyz_list - xyz_bis[j,:] = xyz_list[j,:] - end - end - - return nothing -end #}}} diff --git a/src/core/femmodel.jl b/src/core/femmodel.jl index cd839a5..ee63760 100644 --- a/src/core/femmodel.jl +++ b/src/core/femmodel.jl @@ -39,3 +39,7 @@ function SetCurrentConfiguration!(femmodel::FemModel, analysis::Analysis) #{{{ return nothing end#}}} +function AddResult!(femmodel::FemModel, result::Result) #{{{ + # TODO: maybe need to check if the result already exist, then overwrite + push!(femmodel.results, result) +end#}}} diff --git a/src/core/gauss.jl b/src/core/gauss.jl index 8648935..c11004e 100644 --- a/src/core/gauss.jl +++ b/src/core/gauss.jl @@ -20,34 +20,7 @@ end# }}} #Gauss constructor function GaussTria(order::Int64) #{{{ - - #=Gauss quadrature points for the triangle. - Higher-order points from D.A. Dunavant, "High Degree Efficient - Symmetrical Gaussian Quadrature Rules for the Triangle", IJNME, - Vol. 21, pp. 1129-1148 (1985), as transcribed for Probe rules3.=# - - if(order==1) - npoints = 1 - weights = [1.732050807568877] - coords1 = [0.333333333333333] - coords2 = [0.333333333333333] - coords3 = [0.333333333333333] - elseif(order==2) - npoints = 3 - weights = [0.577350269189625; 0.577350269189625; 0.577350269189625] - coords1 = [0.666666666666667; 0.166666666666667; 0.166666666666667] - coords2 = [0.166666666666667; 0.666666666666667; 0.166666666666667] - coords3 = [0.166666666666667; 0.166666666666667; 0.666666666666667] - elseif(order==3) - npoints = 4 - weights = [-0.974278579257493; 0.902109795608790; 0.902109795608790; 0.902109795608790] - coords1 = [ 0.333333333333333; 0.600000000000000; 0.200000000000000; 0.200000000000000] - coords2 = [ 0.333333333333333; 0.200000000000000; 0.600000000000000; 0.200000000000000] - coords3 = [ 0.333333333333333; 0.200000000000000; 0.200000000000000; 0.600000000000000] - else - error("order ",order," not supported yet"); - end - + npoints, weights, coords1, coords2, coords3 = GaussLegendreTria(order) return GaussTria(npoints,weights,coords1,coords2,coords3) end# }}} function GaussTria(finiteelement::IssmEnum) #{{{ @@ -113,3 +86,139 @@ function GaussTria(area_coordinates::Matrix{Float64}, order::Int64) #{{{ return GaussTria(npoints, weights, coords1, coords2, coords3) end# }}} +function GaussTria(index::Int64, r1::Float64, r2::Float64, mainlyfloating::Bool, order::Int64) #{{{ + xy_list = Matrix{Float64}(undef,3,2) + + if (mainlyfloating) + # Get gauss points + numgauss, weights, coords1, coords2, coords3 = GaussLegendreTria(order) + + xy_list[1,1]=0.; xy_list[1,2]=0.; + xy_list[2,1]=r1; xy_list[2,2]=0.; + xy_list[3,1]=0.; xy_list[3,2]=r2; + + for ii in 1:numgauss + x = coords1[ii]*xy_list[1,1] + coords2[ii]*xy_list[2,1] + coords3[ii]*xy_list[3,1]; + y = coords1[ii]*xy_list[1,2] + coords2[ii]*xy_list[2,2] + coords3[ii]*xy_list[3,2]; + + if (index==1) + coords1[ii] = 1.0-x-y + coords2[ii] = x + coords3[ii] = y + elseif (index==2) + coords1[ii] = y + coords2[ii] = 1.0-x-y + coords3[ii] = x + elseif (index==3) + coords1[ii] = x + coords2[ii] = y + coords3[ii] = 1.0-x-y + else + error("index ", index, " not supported yet") + end + weights[ii] = weights[ii]*r1*r2 + end + else + # Double number of gauss points + gauss1 = GaussTria(order) + xy_list[1,1]=r1; xy_list[1,2]=0.; + xy_list[2,1]=0.; xy_list[2,2]=1.; + xy_list[3,1]=0.; xy_list[3,2]=r2; + + for ii in 1:gauss1.numgauss + x = gauss1.coords1[ii]*xy_list[1,1] + gauss1.coords2[ii]*xy_list[2,1] + gauss1.coords3[ii]*xy_list[3,1]; + y = gauss1.coords1[ii]*xy_list[1,2] + gauss1.coords2[ii]*xy_list[2,2] + gauss1.coords3[ii]*xy_list[3,2]; + + if (index==1) + gauss1.coords1[ii] = 1.0-x-y + gauss1.coords2[ii] = x + gauss1.coords3[ii] = y + elseif (index==2) + gauss1.coords1[ii] = y + gauss1.coords2[ii] = 1.0-x-y + gauss1.coords3[ii] = x + elseif (index==3) + gauss1.coords1[ii] = x + gauss1.coords2[ii] = y + gauss1.coords3[ii] = 1.0-x-y + else + error("index ", index, " not supported yet") + end + gauss1.weights[ii] = gauss1.weights[ii]*r1*(1.0-r2) + end + + gauss2 = GaussTria(order) + xy_list[1,1]=r1; xy_list[1,2]=0.; + xy_list[2,1]=1.; xy_list[2,2]=0.; + xy_list[3,1]=0.; xy_list[3,2]=1.; + + for ii in 1:gauss2.numgauss + x = gauss2.coords1[ii]*xy_list[1,1] + gauss2.coords2[ii]*xy_list[2,1] + gauss2.coords3[ii]*xy_list[3,1]; + y = gauss2.coords1[ii]*xy_list[1,2] + gauss2.coords2[ii]*xy_list[2,2] + gauss2.coords3[ii]*xy_list[3,2]; + + if (index==1) + gauss2.coords1[ii] = 1.0-x-y + gauss2.coords2[ii] = x + gauss2.coords3[ii] = y + elseif (index==2) + gauss2.coords1[ii] = y + gauss2.coords2[ii] = 1.0-x-y + gauss2.coords3[ii] = x + elseif (index==3) + gauss2.coords1[ii] = x + gauss2.coords2[ii] = y + gauss2.coords3[ii] = 1.0-x-y + else + error("index ", index, " not supported yet") + end + gauss2.weights[ii] = gauss2.weights[ii]*(1.0-r1) + end + + numgauss = gauss1.numgauss + gauss2.numgauss + weights = vcat(gauss1.weights, gauss2.weights) + coords1 = vcat(gauss1.coords1, gauss2.coords1) + coords2 = vcat(gauss1.coords2, gauss2.coords2) + coords3 = vcat(gauss1.coords3, gauss2.coords3) + end + + return GaussTria(numgauss, weights, coords1, coords2, coords3) +end# }}} + +#Numerics +function GaussLegendreTria(order::Int64) #{{{ + + #=Gauss quadrature points for the triangle. + Higher-order points from D.A. Dunavant, "High Degree Efficient + Symmetrical Gaussian Quadrature Rules for the Triangle", IJNME, + Vol. 21, pp. 1129-1148 (1985), as transcribed for Probe rules3.=# + + if(order==1) + npoints = 1 + weights = [1.732050807568877] + coords1 = [0.333333333333333] + coords2 = [0.333333333333333] + coords3 = [0.333333333333333] + elseif(order==2) + npoints = 3 + weights = [0.577350269189625; 0.577350269189625; 0.577350269189625] + coords1 = [0.666666666666667; 0.166666666666667; 0.166666666666667] + coords2 = [0.166666666666667; 0.666666666666667; 0.166666666666667] + coords3 = [0.166666666666667; 0.166666666666667; 0.666666666666667] + elseif(order==3) + npoints = 4 + weights = [-0.974278579257493; 0.902109795608790; 0.902109795608790; 0.902109795608790] + coords1 = [ 0.333333333333333; 0.600000000000000; 0.200000000000000; 0.200000000000000] + coords2 = [ 0.333333333333333; 0.200000000000000; 0.600000000000000; 0.200000000000000] + coords3 = [ 0.333333333333333; 0.200000000000000; 0.200000000000000; 0.600000000000000] + elseif(order==4) + npoints = 5 + weights = [0.386908262797819; 0.386908262797819; 0.386908262797819; 0.190442006391807; 0.190442006391807; 0.190442006391807] + coords1 = [0.108103018168070; 0.445948490915965; 0.445948490915965; 0.816847572980459; 0.091576213509771; 0.091576213509771] + coords2 = [0.445948490915965; 0.108103018168070; 0.445948490915965; 0.091576213509771; 0.816847572980459; 0.091576213509771] + coords3 = [0.445948490915965; 0.445948490915965; 0.108103018168070; 0.091576213509771; 0.091576213509771; 0.816847572980459] + else + error("order ",order," not supported yet"); + end + + return npoints, weights, coords1, coords2, coords3 +end# }}} diff --git a/src/core/inputs.jl b/src/core/inputs.jl index 049e9e6..2e496fe 100644 --- a/src/core/inputs.jl +++ b/src/core/inputs.jl @@ -14,64 +14,27 @@ mutable struct Inputs #{{{ end# }}} #Inputs functions -function GetInput(inputs::Inputs,enum::IssmEnum) #{{{ - - #Does this input exist - if !haskey(inputs.lookup,enum) - error("Input ",enum," not found") - end - - #return input - return inputs.lookup[enum] +function DuplicateInput(inputs::Inputs, old::IssmEnum, new::IssmEnum)#{{{ -end#}}} -function SetTriaInput(inputs::Inputs,enum::IssmEnum,interp::IssmEnum,index::Int64,value::Float64) #{{{ + #Fetch input that needs to be copied + oldinput = inputs.lookup[old] - #Does this input exist - if !haskey(inputs.lookup,enum) - #it does not exist yet, we need to create it... - @assert inputs.numberofelements > 0 - if interp==P0Enum - inputs.lookup[enum] = Input(enum,interp,zeros(inputs.numberofelements),Vector{Float64}(undef,1)) - elseif interp==P1Enum - inputs.lookup[enum] = Input(enum,interp,zeros(inputs.numberofvertices),Vector{Float64}(undef,3)) - else - error("not supported yet") - end + if typeof(oldinput)==Input + inputs.lookup[new] = Input(new, oldinput.interp, copy(oldinput.values), copy(oldinput.element_values)) end - #Get this input and check type - input::Input = inputs.lookup[enum] - if typeof(input)!=Input error("input type not consistent") end - if interp!=input.interp error("input interpolations not consistent") end - - #set value - input.values[index] = value - return nothing end#}}} -function SetTriaInput(inputs::Inputs,enum::IssmEnum,interp::IssmEnum,indices::Vector{Int64},values::Vector{Float64}) #{{{ +function GetInput(inputs::Inputs,enum::IssmEnum) #{{{ #Does this input exist if !haskey(inputs.lookup,enum) - #it does not exist yet, we need to create it... - @assert inputs.numberofvertices>0 - if interp==P1Enum - inputs.lookup[enum] = Input(enum,interp,zeros(inputs.numberofvertices),Vector{Float64}(undef,3)) - else - error("not supported yet") - end + error("Input ",enum," not found") end - #Get this input and check type - input::Input = inputs.lookup[enum] - if typeof(input)!=Input error("input type not consistent") end - if interp!=input.interp error("input interpolations not consistent") end - - #set value - input.values[indices] = values + #return input + return inputs.lookup[enum] - return nothing end#}}} function GetInputAverageValue(input::Input) #{{{ @@ -84,24 +47,6 @@ function GetInputAverageValue(input::Input) #{{{ return value/numnodes -end#}}} -function GetInputMin(input::Input) #{{{ - - return minimum(input.element_values) - -end#}}} -function GetInputValue(input::Input,gauss::GaussTria,i::Int64) #{{{ - - if input.interp==P0Enum - return input.element_values[1] - elseif input.interp==P1Enum - value = input.element_values[1]*gauss.coords1[i] + input.element_values[2]*gauss.coords2[i] + input.element_values[3]*gauss.coords3[i] - else - error("not implemented yet") - end - - return value - end#}}} function GetInputDerivativeValue(input::Input,xyz_list::Matrix{Float64},gauss::GaussTria,i::Int64) #{{{ @@ -132,14 +77,74 @@ function GetInputDerivativeValue(input::Input,xyz_list::Matrix{Float64},gauss::G return dp end#}}} -function DuplicateInput(inputs::Inputs, old::IssmEnum, new::IssmEnum)#{{{ +function GetInputMax(input::Input) #{{{ - #Fetch input that needs to be copied - oldinput = inputs.lookup[old] + return maximum(input.element_values) - if typeof(oldinput)==Input - inputs.lookup[new] = Input(new, oldinput.interp, copy(oldinput.values), copy(oldinput.element_values)) +end#}}} +function GetInputMin(input::Input) #{{{ + + return minimum(input.element_values) + +end#}}} +function GetInputValue(input::Input,gauss::GaussTria,i::Int64) #{{{ + + if input.interp==P0Enum + return input.element_values[1] + elseif input.interp==P1Enum + value = input.element_values[1]*gauss.coords1[i] + input.element_values[2]*gauss.coords2[i] + input.element_values[3]*gauss.coords3[i] + else + error("not implemented yet") + end + + return value + +end#}}} +function SetTriaInput(inputs::Inputs,enum::IssmEnum,interp::IssmEnum,index::Int64,value::Float64) #{{{ + + #Does this input exist + if !haskey(inputs.lookup,enum) + #it does not exist yet, we need to create it... + @assert inputs.numberofelements > 0 + if interp==P0Enum + inputs.lookup[enum] = Input(enum,interp,zeros(inputs.numberofelements),Vector{Float64}(undef,1)) + elseif interp==P1Enum + inputs.lookup[enum] = Input(enum,interp,zeros(inputs.numberofvertices),Vector{Float64}(undef,3)) + else + error("not supported yet") + end end + #Get this input and check type + input::Input = inputs.lookup[enum] + if typeof(input)!=Input error("input type not consistent") end + if interp!=input.interp error("input interpolations not consistent") end + + #set value + input.values[index] = value + + return nothing +end#}}} +function SetTriaInput(inputs::Inputs,enum::IssmEnum,interp::IssmEnum,indices::Vector{Int64},values::Vector{Float64}) #{{{ + + #Does this input exist + if !haskey(inputs.lookup,enum) + #it does not exist yet, we need to create it... + @assert inputs.numberofvertices>0 + if interp==P1Enum + inputs.lookup[enum] = Input(enum,interp,zeros(inputs.numberofvertices),Vector{Float64}(undef,3)) + else + error("not supported yet") + end + end + + #Get this input and check type + input::Input = inputs.lookup[enum] + if typeof(input)!=Input error("input type not consistent") end + if interp!=input.interp error("input interpolations not consistent") end + + #set value + input.values[indices] = values + return nothing end#}}} diff --git a/src/core/modules.jl b/src/core/modules.jl index 7b9a447..85bb579 100644 --- a/src/core/modules.jl +++ b/src/core/modules.jl @@ -60,6 +60,7 @@ function ModelProcessor(md::model, solutionstring::Symbol) #{{{ FetchDataToInput(md, inputs, elements, md.inversion.vx_obs./md.constants.yts,VxObsEnum) FetchDataToInput(md, inputs, elements, md.inversion.vy_obs./md.constants.yts,VyObsEnum) AddParam(parameters, md.inversion.independent_string, InversionControlParametersEnum) + AddParam(parameters, md.inversion.dependent_string, InversionCostFunctionsEnum) end #Build FemModel @@ -136,7 +137,6 @@ function CreateInputs(inputs::Inputs,elements::Vector{Tria},md::model) #{{{ end# }}} function OutputResultsx(femmodel::FemModel, md::model, solutionkey::Symbol)# {{{ - if solutionkey===:TransientSolution #Compute maximum number of steps @@ -346,15 +346,43 @@ function RequestedOutputsx(femmodel::FemModel,outputlist::Vector{IssmEnum})# {{{ input_copy = Vector{Float64}(undef, length(input.values)) copyto!(input_copy, input.values) result = Result(outputlist[i], step, time, input_copy) - - #Add to femmodel.results dataset - push!(femmodel.results, result) else - println("Output ",outputlist[i]," not supported yet") + if outputlist[i] == IceVolumeEnum + double_result = IceVolumex(femmodel) + elseif outputlist[i] == IceVolumeAboveFloatationEnum + double_result = IceVolumeAboveFloatationx(femmodel) + elseif outputlist[i] == SurfaceAbsVelMisfitEnum + double_result = SurfaceAbsVelMisfitx(femmodel) + elseif outputlist[i] == DragCoefficientAbsGradientEnum #TODO: define a new Enum + double_result = ControlVariableAbsGradientx(femmodel) + else + error("Output ",outputlist[i]," not supported yet") + end + + result = Result(outputlist[i], step, time, double_result) end + #Add to femmodel.results dataset + AddResult!(femmodel, result) +# push!(femmodel.results, result) end return nothing end# }}} +function IceVolumeAboveFloatationx(femmodel::FemModel)# {{{ + + total_ice_volume_af = 0.0 + for i=1:length(femmodel.elements) + total_ice_volume_af += IceVolumeAboveFloatation(femmodel.elements[i]) + end + return total_ice_volume_af +end# }}} +function IceVolumex(femmodel::FemModel)# {{{ + + total_ice_volume = 0.0 + for i=1:length(femmodel.elements) + total_ice_volume += IceVolume(femmodel.elements[i]) + end + return total_ice_volume +end# }}} function MigrateGroundinglinex(femmodel::FemModel)# {{{ for i=1:length(femmodel.elements) diff --git a/src/core/parameters.jl b/src/core/parameters.jl index ac90697..2b1be2a 100644 --- a/src/core/parameters.jl +++ b/src/core/parameters.jl @@ -16,6 +16,10 @@ struct StringParam <: Parameter #{{{ enum::IssmEnum value::String end# }}} +mutable struct StringArrayParam <: Parameter #{{{ + enum::IssmEnum + value::Vector{String} +end# }}} mutable struct FluxChainParam <: Parameter #{{{ enum::IssmEnum value::Vector{Flux.Chain{}} @@ -47,6 +51,9 @@ end#}}} function GetParameterValue(param::StringParam) #{{{ return param.value::String end#}}} +function GetParameterValue(param::StringArrayParam) #{{{ + return param.value::Vector{String} +end#}}} function GetParameterValue(param::FluxChainParam) #{{{ return param.value::Vector{Flux.Chain{}} end#}}} @@ -79,6 +86,12 @@ function AddParam(parameters::Parameters,value::String, enum::IssmEnum) #{{{ return nothing end#}}} +function AddParam(parameters::Parameters,value::Vector{String}, enum::IssmEnum) #{{{ + + parameters.lookup[enum] = StringArrayParam(enum,value) + + return nothing +end#}}} function AddParam(parameters::Parameters,value::Vector{Flux.Chain{}}, enum::IssmEnum) #{{{ parameters.lookup[enum] = FluxChainParam(enum,value) diff --git a/src/core/solve.jl b/src/core/solve.jl index 4a4c8a3..d627414 100644 --- a/src/core/solve.jl +++ b/src/core/solve.jl @@ -42,11 +42,7 @@ function solve(md::model, solution::Symbol) #{{{ #Solve (FIXME: to be improved later...) if (md.inversion.iscontrol) # solve inverse problem - if solution===:grad - computeGradient(md, femmodel) - else - Control_Core(md, femmodel) - end + Control_Core(md, femmodel, solution) else # otherwise forward problem if(solutionkey===:StressbalanceSolution) analysis = StressbalanceAnalysis() diff --git a/src/usr/classes.jl b/src/usr/classes.jl index af835a7..7de726c 100644 --- a/src/usr/classes.jl +++ b/src/usr/classes.jl @@ -256,9 +256,10 @@ mutable struct Inversion independent::Vector{Float64} maxiter::Int64 independent_string::String + dependent_string::Vector{String} end function Inversion() #{{{ - return Inversion( false, Vector{Float64}(undef,0), Vector{Float64}(undef,0), Vector{Float64}(undef,0), Vector{Float64}(undef,0), Vector{Float64}(undef,0), 0, "Friction") + return Inversion( false, Vector{Float64}(undef,0), Vector{Float64}(undef,0), Vector{Float64}(undef,0), Vector{Float64}(undef,0), Vector{Float64}(undef,0), 0, "Friction", Vector{String}(undef,0)) end# }}} function Base.show(io::IO, this::Inversion)# {{{ IssmStructDisp(io, this)