diff --git a/src/core/analyses/levelsetanalysis.jl b/src/core/analyses/levelsetanalysis.jl new file mode 100644 index 0000000..b0c3a45 --- /dev/null +++ b/src/core/analyses/levelsetanalysis.jl @@ -0,0 +1,245 @@ +#LevelsetAnalysis class definition +struct LevelsetAnalysis <: Analysis#{{{ +end #}}} + +#Model Processing +function CreateConstraints(analysis::LevelsetAnalysis,constraints::Vector{Constraint},md::model) #{{{ + + #load constraints from model + spclevelset = md.levelset.spclevelset + + count = 1 + for i in 1:md.mesh.numberofvertices + if ~isnan(spclevelset[i]) + push!(constraints,Constraint(count,i,1,spclevelset[i])) + count+=1 + end + end + + return nothing +end#}}} +function CreateNodes(analysis::LevelsetAnalysis,nodes::Vector{Node},md::model) #{{{ + + numdof = 1 + for i in 1:md.mesh.numberofvertices + push!(nodes,Node(i,i,true,true,numdof,-ones(Int64,numdof), ones(Int64,numdof), -ones(Int64,numdof), zeros(numdof))) + end + + return nothing +end#}}} +function UpdateElements(analysis::LevelsetAnalysis,elements::Vector{Tria}, inputs::Inputs, md::model) #{{{ + + #Provide node indices to element + for i in 1:md.mesh.numberofelements + Update(elements[i],inputs,i,md,P1Enum) + end + + #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.geometry.thickness,ThicknessEnum) + FetchDataToInput(md,inputs,elements,md.geometry.surface,SurfaceEnum) + FetchDataToInput(md,inputs,elements,md.geometry.base,BaseEnum) + FetchDataToInput(md,inputs,elements,md.mask.ice_levelset, MaskIceLevelsetEnum) + FetchDataToInput(md,inputs,elements,md.mask.ocean_levelset, MaskOceanLevelsetEnum) + + #Get moving front parameters + if typeof(md.calving) == DefaultCalving + FetchDataToInput(md,inputs,elements,md.calving.calvingrate,CalvingCalvingrateEnum) + else + error("Calving ", typeof(md.calving), " not supported yet") + end + + #Get frontal melt parameters + FetchDataToInput(md,inputs,elements,md.frontalforcings.meltingrate, CalvingMeltingrateEnum) + FetchDataToInput(md,inputs,elements,md.frontalforcings.ablationrate, CalvingAblationrateEnum) + + return nothing +end#}}} +function UpdateParameters(analysis::LevelsetAnalysis,parameters::Parameters,md::model) #{{{ + AddParam(parameters,md.levelset.stabilization,LevelsetStabilizationEnum) + AddParam(parameters,md.levelset.reinit_frequency,LevelsetReinitFrequencyEnum) + AddParam(parameters,md.levelset.kill_icebergs,LevelsetKillIcebergsEnum) + AddParam(parameters,md.levelset.migration_max,MigrationMaxEnum) + + #Deal with Calving + if typeof(md.calving)==DefaultCalving + AddParam(parameters, 1, CalvingLawEnum) + else + error("Calving ", typeof(md.calving), " not supported yet") + end + + return nothing +end#}}} + +#Finite Element Analysis +function Core(analysis::LevelsetAnalysis,femmodel::FemModel)# {{{ + + # moving front + MovingFrontalVel(femmodel) + #Activate formulation + SetCurrentConfiguration!(femmodel, analysis) + + #Call solution sequence to compute new speeds + println(" call computational core:"); + solutionsequence_linear(femmodel,analysis) + + # TODO: add reinitialization + # save + RequestedOutputsx(femmodel, [MaskIceLevelsetEnum]) + return nothing +end #}}} +function CreateKMatrix(analysis::LevelsetAnalysis,element::Tria)# {{{ + + #Internmediaries + numnodes = 3 + + #Initialize Element matrix and basis function derivatives + Ke = ElementMatrix(element.nodes) + dbasis = Matrix{Float64}(undef,numnodes,2) + basis = Vector{Float64}(undef,numnodes) + + #Retrieve all inputs and parameters + xyz_list = GetVerticesCoordinates(element.vertices) + mf_vx_input = GetInput(element, MovingFrontalVxEnum) + mf_vy_input = GetInput(element, MovingFrontalVyEnum) + dt = FindParam(Float64, element, TimesteppingTimeStepEnum) + stabilization = FindParam(Int64, element, LevelsetStabilizationEnum) + migration_max = FindParam(Float64, element, MigrationMaxEnum) + + h = CharacteristicLength(element) + + #Start integrating + gauss = GaussTria(2) + for ig in 1:gauss.numgauss + + Jdet = JacobianDeterminant(xyz_list, gauss) + NodalFunctionsDerivatives(element, dbasis, xyz_list, gauss) + NodalFunctions(element, basis, gauss, ig, P1Enum) + + #Transient term + for i in 1:numnodes + for j in 1:numnodes + Ke.values[i ,j] += gauss.weights[ig]*Jdet*basis[i]*basis[j] + end + end + + # levelset speed + wx = GetInputValue(mf_vx_input, gauss, ig) + wy = GetInputValue(mf_vy_input, gauss, ig) + vel = sqrt(wx^2+wy^2)+1.0e-14 + + #rescale by migration_max + if (vel > migration_max) + wx = wx/vel*migration_max + wy = wy/vel*migration_max + end + + for i in 1:numnodes + for j in 1:numnodes + #\phi_i v\cdot\nabla\phi_j + Ke.values[i ,j] += dt*gauss.weights[ig]*Jdet*basis[i]*(wx*dbasis[j,1] + wy*dbasis[j,2]) + end + end + + #Stabilization + if(stabilization==0) + #do nothing + elseif (stabilization==1) + hx, hy, hz = GetElementSizes(element) + h = sqrt((hx*wx/vel)^2 + (hy*wy/vel)^2) + kappa = h*vel/2.0 + D = dt*gauss.weights[ig]*Jdet*kappa + for i in 1:numnodes + for j in 1:numnodes + for k in 1:2 + Ke.values[i ,j] += D*dbasis[i,k]*dbasis[j,k] + end + end + end + else + error("Stabilization ",stabilization, " not supported yet") + end + end + + return Ke +end #}}} +function CreatePVector(analysis::LevelsetAnalysis,element::Tria)# {{{ + + #Internmediaries + numnodes = 3 + + #Initialize Element vectro and basis functions + pe = ElementVector(element.nodes) + basis = Vector{Float64}(undef,numnodes) + + #Retrieve all inputs and parameters + xyz_list = GetVerticesCoordinates(element.vertices) + levelset_input = GetInput(element, MaskIceLevelsetEnum) + mf_vx_input = GetInput(element, MovingFrontalVxEnum) + mf_vy_input = GetInput(element, MovingFrontalVyEnum) + dt = FindParam(Float64, element, TimesteppingTimeStepEnum) + stabilization = FindParam(Int64, element, LevelsetStabilizationEnum) + migration_max = FindParam(Float64, element, MigrationMaxEnum) + + h = CharacteristicLength(element) + + #Start integrating + gauss = GaussTria(2) + for ig in 1:gauss.numgauss + + Jdet = JacobianDeterminant(xyz_list, gauss) + #Get nodal basis + NodalFunctions(element, basis, gauss, ig, P1Enum) + #Old function value + lsf = GetInputValue(levelset_input, gauss, ig) + for i in 1:numnodes + pe.values[i] += gauss.weights[ig]*Jdet*lsf*basis[i] + end + #TODO: add stab=5 + end + + return pe +end #}}} +function GetSolutionFromInputs(analysis::LevelsetAnalysis,ug::IssmVector,element::Tria) #{{{ + + #Get dofs for this finite element + doflist = GetDofList(element,GsetEnum) + @assert length(doflist)==3 + + #Fetch inputs + mask_input = GetInput(element, MaskIceLevelsetEnum) + + #Loop over each node and enter solution in ug + count = 0 + gauss=GaussTria(P1Enum) + for i in 1:gauss.numgauss + mask = GetInputValue(mask_input, gauss, i) + + count += 1 + ug.vector[doflist[count]] = mask + end + + #Make sure we reached all the values + @assert count==length(doflist) + + return nothing +end#}}} +function InputUpdateFromSolution(analysis::LevelsetAnalysis,ug::Vector{Float64},element::Tria) #{{{ + InputUpdateFromSolutionOneDof(element, ug, MaskIceLevelsetEnum) +end#}}} +function UpdateConstraints(analysis::LevelsetAnalysis, femmodel::FemModel) #{{{ + #Default do nothing + return nothing +end#}}} + +# Moving front +function MovingFrontalVel(femmodel::FemModel)# {{{ + for i in 1:length(femmodel.elements) + MovingFrontalVelocity(femmodel.elements[i]) + end + return nothing +end #}}} diff --git a/src/core/analyses/transientanalysis.jl b/src/core/analyses/transientanalysis.jl index ec4541e..3a80317 100644 --- a/src/core/analyses/transientanalysis.jl +++ b/src/core/analyses/transientanalysis.jl @@ -21,6 +21,7 @@ function Core(analysis::TransientAnalysis,femmodel::FemModel)# {{{ yts = FindParam(Float64, femmodel.parameters, ConstantsYtsEnum) dt = FindParam(Float64, femmodel.parameters, TimesteppingTimeStepEnum) + ismovingfront = FindParam(Bool, femmodel.parameters, TransientIsmovingfrontEnum) isstressbalance = FindParam(Bool, femmodel.parameters, TransientIsstressbalanceEnum) ismasstransport = FindParam(Bool, femmodel.parameters, TransientIsmasstransportEnum) @@ -32,6 +33,9 @@ function Core(analysis::TransientAnalysis,femmodel::FemModel)# {{{ println("iteration ", step, "/", Int(ceil((finaltime-time)/dt))+step," time [yr]: ", time/yts, " (time step: ", dt/yts, " [yr])") if(isstressbalance) Core(StressbalanceAnalysis(), femmodel) end + + if(ismovingfront) Core(LevelsetAnalysis(), femmodel) end + if(ismasstransport) Core(MasstransportAnalysis(), femmodel) end MigrateGroundinglinex(femmodel) diff --git a/src/core/control.jl b/src/core/control.jl index 844bea2..9d31402 100644 --- a/src/core/control.jl +++ b/src/core/control.jl @@ -1,36 +1,43 @@ using Enzyme +Enzyme.API.looseTypeAnalysis!(false) +Enzyme.API.strictAliasing!(false) +Enzyme.API.typeWarning!(false) + using Optimization, OptimizationOptimJL -#Enzyme.API.looseTypeAnalysis!(false) -#Enzyme.API.strictAliasing!(false) -#Enzyme.API.typeWarning!(false) Enzyme.Compiler.RunAttributor[] = false 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()) - #TODO: Put the solution back in results according to its Enum - #InputUpdateFromVectorx(femmodel, sol.u, GradientEnum, VertexSIdEnum) - #RequestedOutputsx(femmodel, [GradientEnum]) + 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) #{{{ #independent variable α = md.inversion.independent #initialize derivative as 0 ∂J_∂α = zero(α) - - # 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 - #println("CALLING AUTODIFF, prepare to die...") - @time autodiff(Enzyme.Reverse, costfunction, Duplicated(femmodel, dfemmodel), Duplicated(α, ∂J_∂α)) + # Compute Gradient + computeGradient(∂J_∂α, α, femmodel) #Put gradient in results InputUpdateFromVectorx(femmodel, ∂J_∂α, GradientEnum, VertexSIdEnum) RequestedOutputsx(femmodel, [GradientEnum]) end#}}} +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 + @time autodiff(Enzyme.Reverse, costfunction, Duplicated(α, ∂J_∂α), Duplicated(femmodel,dfemmodel)) +end#}}} diff --git a/src/core/costfunctions.jl b/src/core/costfunctions.jl index 38591ba..0dc6268 100644 --- a/src/core/costfunctions.jl +++ b/src/core/costfunctions.jl @@ -1,5 +1,5 @@ -function costfunction(femmodel::FemModel, α::Vector{Float64}) #{{{ +function costfunction(α::Vector{Float64}, femmodel::FemModel) #{{{ # get the md.inversion.control_string control_string = FindParam(String, femmodel.parameters, InversionControlParametersEnum) # get the Enum @@ -9,9 +9,9 @@ function costfunction(femmodel::FemModel, α::Vector{Float64}) #{{{ end # compute cost function # TODO: loop through all controls with respect to all the components in the cost function - costfunction(femmodel, α, controlvar_enum, VertexSIdEnum) + costfunction(α, femmodel, controlvar_enum, VertexSIdEnum) end#}}} -function costfunction(femmodel::FemModel, α::Vector{Float64}, controlvar_enum::IssmEnum, SId_enum::IssmEnum) #{{{ +function costfunction(α::Vector{Float64}, femmodel::FemModel, controlvar_enum::IssmEnum, SId_enum::IssmEnum) #{{{ #Update FemModel accordingly InputUpdateFromVectorx(femmodel, α, controlvar_enum, SId_enum) diff --git a/src/core/elements.jl b/src/core/elements.jl index 060682f..1514889 100644 --- a/src/core/elements.jl +++ b/src/core/elements.jl @@ -237,6 +237,21 @@ function GetGroundedPortion(element::Tria, xyz_list::Matrix{Float64}) #{{{ return phi end#}}} +function InputUpdateFromSolutionOneDof(element::Tria, ug::Vector{Float64},enum::IssmEnum) #{{{ + #Get dofs for this finite element + doflist = GetDofList(element,GsetEnum) + + #Get solution vector for this element + numdof = 3 + values = Vector{Float64}(undef,numdof) + for i in 1:numdof values[i]=ug[doflist[i]] end + + #Add back to Mask + AddInput(element, enum, values, P1Enum) + + return nothing + +end#}}} function IsIcefront(element::Tria) #{{{ level = Vector{Float64}(undef,3) @@ -308,6 +323,30 @@ function GetArea(element::Tria)#{{{ @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 @@ -373,6 +412,65 @@ function MigrateGroundingLine(element::Tria) #{{{ return nothing end#}}} +function MovingFrontalVelocity(element::Tria) #{{{ + + mvx = Vector{Float64}(undef,3) + mvy = Vector{Float64}(undef,3) + + # get cavling law + calvinglaw = FindParam(Int64, element, CalvingLawEnum) + + # load inputs + gr_input = GetInput(element, MaskOceanLevelsetEnum) + #lsf_slopex_input = GetInput(element, LevelsetfunctionSlopeXEnum) + #lsf_slopey_input = GetInput(element, LevelsetfunctionSlopeYEnum) + #calvingratex_input= GetInput(element, CalvingratexEnum) + #calvingratey_input= GetInput(element, CalvingrateyEnum) + vx_input = GetInput(element, VxEnum) + vy_input = GetInput(element, VyEnum) + calvingrate_input = GetInput(element, CalvingCalvingrateEnum) + meltingrate_input = GetInput(element, CalvingMeltingrateEnum) + + xyz_list = GetVerticesCoordinates(element.vertices) + + gauss = GaussTria(2) + for ig in 1:3 + vx = GetInputValue(vx_input, gauss, ig) + vy = GetInputValue(vy_input, gauss, ig) + groundedice = GetInputValue(gr_input, gauss, ig) + + dlsf = GetInputDerivativeValue(vx_input,xyz_list,gauss,ig) + #TODO: add L2Projection + #dlsfx = GetInputValue(lsf_slopex_input, gauss, ig) + #dlsfy = GetInputValue(lsf_slopey_input, gauss, ig) + #cx = GetInputValue(calvingratex_input, gauss, ig) + #cy = GetInputValue(calvingratey_input, gauss, ig) + + norm_dlsf = sqrt(dlsf[1]^2+dlsf[2]^2) + + # TODO: depend on which calving law + calvingrate = GetInputValue(calvingrate_input, gauss, ig) + cx = calvingrate * dlsf[1] /norm_dlsf + cy = calvingrate * dlsf[2] /norm_dlsf + meltingrate = GetInputValue(meltingrate_input, gauss, ig) + if (groundedice < 0) meltingrate = 0; end + mx = meltingrate * dlsf[1] /norm_dlsf + my = meltingrate * dlsf[2] /norm_dlsf + if (norm_dlsf <1e-10) + cx = 0; cy = 0 + mx = 0; my = 0 + end + + mvx[ig] = vx-cx-mx + mvy[ig] = vy-cy-my + end + + #Update inputs + AddInput(element,MovingFrontalVxEnum,mvx,P1Enum) + AddInput(element,MovingFrontalVyEnum,mvy,P1Enum) + + return nothing +end#}}} function GetXcoord(element::Tria, xyz_list::Matrix{Float64}, gauss::GaussTria, ig::Int64) #{{{ # create a list of x diff --git a/src/core/modules.jl b/src/core/modules.jl index 42413a5..7067c76 100644 --- a/src/core/modules.jl +++ b/src/core/modules.jl @@ -27,7 +27,11 @@ function ModelProcessor(md::model, solutionstring::Symbol) #{{{ if solutionstring===:StressbalanceSolution analyses = Analysis[StressbalanceAnalysis()] elseif solutionstring===:TransientSolution - analyses = Analysis[StressbalanceAnalysis(), MasstransportAnalysis()] + if md.transient.ismovingfront + analyses = Analysis[StressbalanceAnalysis(), MasstransportAnalysis(), LevelsetAnalysis()] + else + analyses = Analysis[StressbalanceAnalysis(), MasstransportAnalysis()] + end else error(solutionstring, " not supported by ModelProcessor") end @@ -117,6 +121,9 @@ function CreateParameters(parameters::Parameters,md::model) #{{{ AddParam(parameters,1,StepEnum) AddParam(parameters,0.0,TimeEnum) + #Is moving front + AddParam(parameters,md.transient.ismovingfront,TransientIsmovingfrontEnum) + return nothing end# }}} function CreateInputs(inputs::Inputs,elements::Vector{Tria},md::model) #{{{ diff --git a/src/core/solve.jl b/src/core/solve.jl index 2840728..4a4c8a3 100644 --- a/src/core/solve.jl +++ b/src/core/solve.jl @@ -19,6 +19,7 @@ include("./control.jl") include("./analyses/stressbalanceanalysis.jl") include("./analyses/masstransportanalysis.jl") include("./analyses/transientanalysis.jl") +include("./analyses/levelsetanalysis.jl") include("./solutionsequences.jl") include("./modules.jl") diff --git a/src/usr/classes.jl b/src/usr/classes.jl index 982a3bf..af835a7 100644 --- a/src/usr/classes.jl +++ b/src/usr/classes.jl @@ -240,7 +240,7 @@ mutable struct Transient ismovingfront::Bool end function Transient() #{{{ - return Transient( true, true, true, true, true) + return Transient( true, true, true, false, false) end# }}} function Base.show(io::IO, this::Transient)# {{{ IssmStructDisp(io, this) @@ -264,9 +264,48 @@ function Base.show(io::IO, this::Inversion)# {{{ IssmStructDisp(io, this) end# }}} # }}} +#Calving {{{ +abstract type AbstractCalving end +mutable struct DefaultCalving <: AbstractCalving + calvingrate::Vector{Float64} +end +function DefaultCalving() #{{{ + return DefaultCalving(Vector{Float64}(undef,0)) +end# }}} +function Base.show(io::IO, this::DefaultCalving)# {{{ + IssmStructDisp(io, this) +end# }}} +# }}} +#Levelset{{{ +mutable struct Levelset + spclevelset::Vector{Float64} + stabilization::Int64 + reinit_frequency::Int64 + kill_icebergs::Int64 + migration_max::Float64 +end +function Levelset() #{{{ + return Levelset(Vector{Float64}(undef,0), 1, 10, 1, 1.0e12) +end# }}} +function Base.show(io::IO, this::Levelset)# {{{ + IssmStructDisp(io, this) +end# }}} +# }}} +#Frontalforcings{{{ +mutable struct Frontalforcings + meltingrate::Vector{Float64} + ablationrate::Vector{Float64} +end +function Frontalforcings() #{{{ + return Frontalforcings(Vector{Float64}(undef,0), Vector{Float64}(undef,0)) +end# }}} +function Base.show(io::IO, this::Frontalforcings)# {{{ + IssmStructDisp(io, this) +end# }}} +# }}} #Model structure -mutable struct model{Mesh<:AbstractMesh, Friction<:AbstractFriction} +mutable struct model{Mesh<:AbstractMesh, Friction<:AbstractFriction, Calving<:AbstractCalving} mesh::Mesh geometry::Geometry mask::Mask @@ -282,18 +321,23 @@ mutable struct model{Mesh<:AbstractMesh, Friction<:AbstractFriction} masstransport::Masstransport transient::Transient inversion::Inversion + calving::Calving + levelset::Levelset + frontalforcings::Frontalforcings end function model() #{{{ return model( Mesh2dTriangle(), Geometry(), Mask(), Materials(), Initialization(),Stressbalance(), Constants(), Dict(), BuddFriction(), Basalforcings(), SMBforcings(), Timestepping(), - Masstransport(), Transient(), Inversion()) + Masstransport(), Transient(), Inversion(), DefaultCalving(), + Levelset(), Frontalforcings()) end#}}} -function model(md::model; mesh::AbstractMesh=md.mesh, friction::AbstractFriction=md.friction) #{{{ +function model(md::model; mesh::AbstractMesh=md.mesh, friction::AbstractFriction=md.friction, calving::AbstractCalving=md.calving) #{{{ return model(mesh, md.geometry, md.mask, md.materials, md.initialization, md.stressbalance, md.constants, md.results, friction, md.basalforcings, md.smb, md.timestepping, - md.masstransport, md.transient, md.inversion) + md.masstransport, md.transient, md.inversion, md.calving, + md.levelset, md.frontalforcings) end#}}} function model(matmd::Dict,verbose::Bool=true) #{{{ @@ -371,6 +415,9 @@ function Base.show(io::IO, md::model)# {{{ @printf "%19s: %-26s -- %s\n" "masstransport" typeof(md.masstransport) "parameters mass transport simulations" @printf "%19s: %-26s -- %s\n" "transient" typeof(md.transient) "parameters for transient simulations" @printf "%19s: %-26s -- %s\n" "inversion" typeof(md.inversion) "parameters for inverse methods" + @printf "%19s: %-26s -- %s\n" "calving" typeof(md.calving) "parameters for calving" + @printf "%19s: %-26s -- %s\n" "levelset" typeof(md.levelset) "parameters for moving boundaries (level-set method)" + @printf "%19s: %-26s -- %s\n" "frontalforcings" typeof(md.frontalforcings) "parameters for frontalforcings" @printf "%19s: %-26s -- %s\n" "results" typeof(md.results) "model results" end# }}} diff --git a/test/Data/SquareShelf.arch b/test/Data/SquareShelf.arch new file mode 100644 index 0000000..a397911 Binary files /dev/null and b/test/Data/SquareShelf.arch differ diff --git a/test/runtests.jl b/test/runtests.jl index aa791f9..68f8d49 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -18,7 +18,7 @@ end # AD test @time include("testad.jl") - @time include("testad2.jl") + #@time include("testad2.jl") # GPU test #@time include("testGPU.jl") diff --git a/test/test807.jl b/test/test807.jl new file mode 100755 index 0000000..8a09b18 --- /dev/null +++ b/test/test807.jl @@ -0,0 +1,74 @@ +using dJUICE + +md = model() +md = triangle(md,issmdir()*"/test/Exp/Square.exp",200000.) +md = setmask(md,"all","") + +#Geometry +hmin=300. +hmax=1000. +ymin=minimum(md.mesh.y) +ymax=maximum(md.mesh.y) +xmin=minimum(md.mesh.x) +xmax=maximum(md.mesh.x) +md.geometry.thickness = hmax .+ (hmin-hmax)*(md.mesh.y .- ymin)./(ymax-ymin) .+ 0.1*(hmin-hmax)*(md.mesh.x .- xmin)./(xmax-xmin) +md.geometry.base = -md.materials.rho_ice/md.materials.rho_water*md.geometry.thickness +md.geometry.surface = md.geometry.base+md.geometry.thickness +md.geometry.bed = md.geometry.base .-500 + +#Initial velocity +x = archread(issmdir()*"/test/Data/SquareShelf.arch","x") +y = archread(issmdir()*"/test/Data/SquareShelf.arch","y") +vx = archread(issmdir()*"/test/Data/SquareShelf.arch","vx") +vy = archread(issmdir()*"/test/Data/SquareShelf.arch","vy") +index = Int.(archread(issmdir()*"/test/Data/SquareShelf.arch","index")) +md.initialization.vx=InterpFromMeshToMesh2d(index,x,y,vx,md.mesh.x,md.mesh.y,0.0) +md.initialization.vy=InterpFromMeshToMesh2d(index,x,y,vy,md.mesh.x,md.mesh.y,0.0) + +md.materials.rheology_B=1.815730284801701e+08*ones(md.mesh.numberofvertices) +md.materials.rheology_n=3*ones(md.mesh.numberofelements) +md.friction.coefficient=20*ones(md.mesh.numberofvertices) +md.friction.p=ones(md.mesh.numberofvertices) +md.friction.q=ones(md.mesh.numberofvertices) + +md.stressbalance.restol=0.05 +md.stressbalance.reltol=0.05 +md.stressbalance.abstol=NaN + +#Boundary conditions +md.stressbalance.spcvx = NaN*ones(md.mesh.numberofvertices) +md.stressbalance.spcvy = NaN*ones(md.mesh.numberofvertices) +md.masstransport.spcthickness = NaN*ones(md.mesh.numberofvertices) +pos = findall(md.mesh.vertexonboundary) +md.stressbalance.spcvx[pos] .= 0.0 +md.stressbalance.spcvy[pos] .= 0.0 +md.masstransport.spcthickness[pos] .= md.geometry.thickness[pos] + +# surface mass balance and basal melting +md.smb.mass_balance=10*ones(md.mesh.numberofvertices) +md.basalforcings.floatingice_melting_rate=5*ones(md.mesh.numberofvertices) +md.basalforcings.groundedice_melting_rate=5*ones(md.mesh.numberofvertices) + +# mask +Lx = xmax - xmin +alpha = 2.0/3.0 +md.mask.ice_levelset = ((md.mesh.x .- alpha*Lx).>0) .- ((md.mesh.x .- alpha*Lx).<0) + +md.levelset.kill_icebergs = 0 +# time stepping +md.timestepping.time_step = 10; +md.timestepping.final_time = 30; + +# Transient +md.transient.isstressbalance=1; +md.transient.ismasstransport=1; +md.transient.issmb=1; +md.transient.ismovingfront=1; + +md.calving.calvingrate=0*ones(md.mesh.numberofvertices) +md.frontalforcings.meltingrate=10000*ones(md.mesh.numberofvertices) +md.frontalforcings.ablationrate=10000*ones(md.mesh.numberofvertices) +md.levelset.spclevelset=NaN*ones(md.mesh.numberofvertices) +md.levelset.spclevelset[pos] = md.mask.ice_levelset[pos] + +md=solve(md,:Transient) diff --git a/test/testad.jl b/test/testad.jl index 5c9c44a..45b80d2 100755 --- a/test/testad.jl +++ b/test/testad.jl @@ -5,6 +5,7 @@ using MAT using Test using Enzyme +Enzyme.API.typeWarning!(false) Enzyme.Compiler.RunAttributor[] = false #Load model from MATLAB file @@ -36,7 +37,7 @@ addJ = md.results["StressbalanceSolution"]["Gradient"] α = md.inversion.independent femmodel=dJUICE.ModelProcessor(md, :StressbalanceSolution) - J1 = dJUICE.costfunction(femmodel, α) + J1 = dJUICE.costfunction(α, femmodel) @test ~isnothing(J1) end @@ -44,12 +45,12 @@ end α = md.inversion.independent delta = 1e-7 femmodel=dJUICE.ModelProcessor(md, :StressbalanceSolution) - J1 = dJUICE.costfunction(femmodel, α) + J1 = dJUICE.costfunction(α, femmodel) for i in 1:md.mesh.numberofvertices dα = zero(md.friction.coefficient) dα[i] = delta femmodel=dJUICE.ModelProcessor(md, :StressbalanceSolution) - J2 = dJUICE.costfunction(femmodel, α+dα) + J2 = dJUICE.costfunction(α+dα, femmodel) dJ = (J2-J1)/delta @test abs(dJ - addJ[i])< 1e-5 diff --git a/test/testad2.jl b/test/testad2.jl index fae3fe1..306b9c3 100755 --- a/test/testad2.jl +++ b/test/testad2.jl @@ -5,6 +5,7 @@ using MAT using Test using Enzyme +Enzyme.API.typeWarning!(false) Enzyme.Compiler.RunAttributor[] = false #Load model from MATLAB file @@ -35,7 +36,7 @@ addJ = md.results["StressbalanceSolution"]["Gradient"] α = md.inversion.independent femmodel=dJUICE.ModelProcessor(md, :StressbalanceSolution) - J1 = dJUICE.costfunction(femmodel, α) + J1 = dJUICE.costfunction(α, femmodel) @test ~isnothing(J1) end @@ -43,12 +44,12 @@ end α = md.inversion.independent delta = 1e-8 femmodel=dJUICE.ModelProcessor(md, :StressbalanceSolution) - J1 = dJUICE.costfunction(femmodel, α) + J1 = dJUICE.costfunction(α, femmodel) for i in 1:md.mesh.numberofvertices dα = zero(md.friction.coefficient) dα[i] = delta femmodel=dJUICE.ModelProcessor(md, :StressbalanceSolution) - J2 = dJUICE.costfunction(femmodel, α+dα) + J2 = dJUICE.costfunction(α+dα, femmodel) dJ = (J2-J1)/delta @test abs(dJ - addJ[i])< 1e-6 diff --git a/test/testoptimization.jl b/test/testoptimization.jl new file mode 100755 index 0000000..0a8eb9a --- /dev/null +++ b/test/testoptimization.jl @@ -0,0 +1,46 @@ +module enzymeDiff + +using dJUICE +using MAT +using Test +using Enzyme + +Enzyme.Compiler.RunAttributor[] = false +using Optimization, OptimizationOptimJL + +#Load model from MATLAB file +#file = matopen(joinpath(@__DIR__, "..", "data","temp12k.mat")) #BIG model +file = matopen(joinpath(@__DIR__, "..", "data","temp.mat")) #SMALL model (35 elements) +mat = read(file, "md") +close(file) +md = model(mat) + +#make model run faster +md.stressbalance.maxiter = 20 + +#Now call AD! +md.inversion.iscontrol = 1 +md.inversion.independent = md.friction.coefficient +md.inversion.min_parameters = ones(md.mesh.numberofvertices)*(0.0) +md.inversion.max_parameters = ones(md.mesh.numberofvertices)*(1.0e3) +md.inversion.independent_string = "FrictionCoefficient" + +α = md.inversion.independent +∂J_∂α = zero(α) + +femmodel=dJUICE.ModelProcessor(md, :StressbalanceSolution) +n = length(α) +# use user defined grad, errors! +optprob = OptimizationFunction(dJUICE.costfunction, Optimization.AutoEnzyme(), grad=dJUICE.computeGradient) +#optprob = OptimizationFunction(dJUICE.costfunction, Optimization.AutoEnzyme()) +#prob = Optimization.OptimizationProblem(optprob, α, femmodel, lb=md.inversion.min_parameters, ub=md.inversion.max_parameters) +prob = Optimization.OptimizationProblem(optprob, α, femmodel) +sol = Optimization.solve(prob, Optim.LBFGS()) + + +#md = solve(md, :sb) + +# compute gradient by finite differences at each node +#addJ = md.results["StressbalanceSolution"]["Gradient"] + +end