From 81a3f6808c17b8f4907bd0bbddf4f746cf87f6ad Mon Sep 17 00:00:00 2001 From: Cheng Gong Date: Mon, 25 Mar 2024 15:48:34 -0400 Subject: [PATCH 01/15] test optm --- src/core/control.jl | 9 +++++---- src/core/costfunctions.jl | 2 +- 2 files changed, 6 insertions(+), 5 deletions(-) diff --git a/src/core/control.jl b/src/core/control.jl index 844bea2..db78889 100644 --- a/src/core/control.jl +++ b/src/core/control.jl @@ -1,9 +1,10 @@ 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) #{{{ @@ -28,7 +29,7 @@ function computeGradient(md::model, femmodel::FemModel) #{{{ 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_∂α)) + @time autodiff(Enzyme.Reverse, costfunction, Duplicated(α, ∂J_∂α), Duplicated(femmodel,dfemmodel)) #Put gradient in results InputUpdateFromVectorx(femmodel, ∂J_∂α, GradientEnum, VertexSIdEnum) diff --git a/src/core/costfunctions.jl b/src/core/costfunctions.jl index 38591ba..a872389 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 From ae8299a3fdc6c4fe6a0f35e61d1f2b15dadd3e31 Mon Sep 17 00:00:00 2001 From: Cheng Gong Date: Mon, 25 Mar 2024 15:49:32 -0400 Subject: [PATCH 02/15] add test case --- test/testoptimization.jl | 30 ++++++++++++++++++++++++++++++ 1 file changed, 30 insertions(+) create mode 100755 test/testoptimization.jl diff --git a/test/testoptimization.jl b/test/testoptimization.jl new file mode 100755 index 0000000..3e26a57 --- /dev/null +++ b/test/testoptimization.jl @@ -0,0 +1,30 @@ +module enzymeDiff + +using dJUICE +using MAT +using Test +using Enzyme + +Enzyme.Compiler.RunAttributor[] = false + +#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.independent_string = "FrictionCoefficient" + +md = solve(md, :sb) + +# compute gradient by finite differences at each node +addJ = md.results["StressbalanceSolution"]["Gradient"] + +end From eb3ec423ebd77818b26427350743dec0bd3b024d Mon Sep 17 00:00:00 2001 From: Cheng Gong Date: Mon, 25 Mar 2024 16:27:19 -0400 Subject: [PATCH 03/15] change the order in calling costfunction, fix tests --- test/testad.jl | 7 ++++--- test/testad2.jl | 7 ++++--- 2 files changed, 8 insertions(+), 6 deletions(-) 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 From b7a7db020afa3802a057ee53c290fbb1b61a6e4a Mon Sep 17 00:00:00 2001 From: Cheng Gong Date: Mon, 25 Mar 2024 17:28:21 -0400 Subject: [PATCH 04/15] fix optimization problem --- src/core/control.jl | 7 ++++--- test/testoptimization.jl | 4 +++- 2 files changed, 7 insertions(+), 4 deletions(-) diff --git a/src/core/control.jl b/src/core/control.jl index db78889..29cc220 100644 --- a/src/core/control.jl +++ b/src/core/control.jl @@ -9,15 +9,16 @@ 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 n = length(α) 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 diff --git a/test/testoptimization.jl b/test/testoptimization.jl index 3e26a57..a72608d 100755 --- a/test/testoptimization.jl +++ b/test/testoptimization.jl @@ -20,11 +20,13 @@ 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 = solve(md, :sb) # compute gradient by finite differences at each node -addJ = md.results["StressbalanceSolution"]["Gradient"] +#addJ = md.results["StressbalanceSolution"]["Gradient"] end From 9fceae683dde9997b330c5345d57a0d76c8f2110 Mon Sep 17 00:00:00 2001 From: Cheng Gong Date: Tue, 26 Mar 2024 11:00:59 -0400 Subject: [PATCH 05/15] update computeGradient --- src/core/control.jl | 15 ++++++++++----- src/core/costfunctions.jl | 4 ++-- 2 files changed, 12 insertions(+), 7 deletions(-) diff --git a/src/core/control.jl b/src/core/control.jl index 29cc220..9d31402 100644 --- a/src/core/control.jl +++ b/src/core/control.jl @@ -11,7 +11,10 @@ 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()) @@ -25,14 +28,16 @@ function computeGradient(md::model, femmodel::FemModel) #{{{ α = 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]) +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 - #println("CALLING AUTODIFF, prepare to die...") @time autodiff(Enzyme.Reverse, costfunction, Duplicated(α, ∂J_∂α), Duplicated(femmodel,dfemmodel)) - - #Put gradient in results - InputUpdateFromVectorx(femmodel, ∂J_∂α, GradientEnum, VertexSIdEnum) - RequestedOutputsx(femmodel, [GradientEnum]) end#}}} diff --git a/src/core/costfunctions.jl b/src/core/costfunctions.jl index a872389..0dc6268 100644 --- a/src/core/costfunctions.jl +++ b/src/core/costfunctions.jl @@ -9,9 +9,9 @@ function costfunction(α::Vector{Float64}, femmodel::FemModel) #{{{ 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) From bf47e05cd4a1e7972de728b69361e0bed5dc36a9 Mon Sep 17 00:00:00 2001 From: Cheng Gong Date: Tue, 26 Mar 2024 16:01:34 -0400 Subject: [PATCH 06/15] test case of LBFGS --- test/testoptimization.jl | 16 +++++++++++++++- 1 file changed, 15 insertions(+), 1 deletion(-) diff --git a/test/testoptimization.jl b/test/testoptimization.jl index a72608d..0a8eb9a 100755 --- a/test/testoptimization.jl +++ b/test/testoptimization.jl @@ -6,6 +6,7 @@ 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 @@ -24,7 +25,20 @@ 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 = solve(md, :sb) +α = 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"] From f0cd102977d846a3410b2d43098e307fd09a6342 Mon Sep 17 00:00:00 2001 From: Cheng Gong Date: Tue, 2 Apr 2024 15:31:50 -0600 Subject: [PATCH 07/15] ADD: calving class --- src/usr/classes.jl | 18 ++++++++++++++++-- 1 file changed, 16 insertions(+), 2 deletions(-) diff --git a/src/usr/classes.jl b/src/usr/classes.jl index 982a3bf..35642c1 100644 --- a/src/usr/classes.jl +++ b/src/usr/classes.jl @@ -264,6 +264,18 @@ function Base.show(io::IO, this::Inversion)# {{{ IssmStructDisp(io, this) end# }}} # }}} +#Calving{{{ +mutable struct Calving + law::Int64 + calvingrate::Vector{Float64} +end +function Calving() #{{{ + return Calving( 0, Vector{Float64}(undef,0)) +end# }}} +function Base.show(io::IO, this::Calving)# {{{ + IssmStructDisp(io, this) +end# }}} +# }}} #Model structure mutable struct model{Mesh<:AbstractMesh, Friction<:AbstractFriction} @@ -282,18 +294,19 @@ mutable struct model{Mesh<:AbstractMesh, Friction<:AbstractFriction} masstransport::Masstransport transient::Transient inversion::Inversion + calving::Calving end function model() #{{{ return model( Mesh2dTriangle(), Geometry(), Mask(), Materials(), Initialization(),Stressbalance(), Constants(), Dict(), BuddFriction(), Basalforcings(), SMBforcings(), Timestepping(), - Masstransport(), Transient(), Inversion()) + Masstransport(), Transient(), Inversion(), Calving()) end#}}} function model(md::model; mesh::AbstractMesh=md.mesh, friction::AbstractFriction=md.friction) #{{{ 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) end#}}} function model(matmd::Dict,verbose::Bool=true) #{{{ @@ -371,6 +384,7 @@ 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) "calving" @printf "%19s: %-26s -- %s\n" "results" typeof(md.results) "model results" end# }}} From ac74282d34fc8d03a7e4970e1602e9e91ced5aaf Mon Sep 17 00:00:00 2001 From: Cheng Gong Date: Tue, 2 Apr 2024 15:36:45 -0600 Subject: [PATCH 08/15] add levelset --- src/usr/classes.jl | 17 +++++++++++++++-- 1 file changed, 15 insertions(+), 2 deletions(-) diff --git a/src/usr/classes.jl b/src/usr/classes.jl index 35642c1..22dc564 100644 --- a/src/usr/classes.jl +++ b/src/usr/classes.jl @@ -276,6 +276,17 @@ function Base.show(io::IO, this::Calving)# {{{ IssmStructDisp(io, this) end# }}} # }}} +#Levelset{{{ +mutable struct Levelset + spclevelset::Vector{Float64} +end +function Levelset() #{{{ + return Levelset(Vector{Float64}(undef,0)) +end# }}} +function Base.show(io::IO, this::Levelset)# {{{ + IssmStructDisp(io, this) +end# }}} +# }}} #Model structure mutable struct model{Mesh<:AbstractMesh, Friction<:AbstractFriction} @@ -295,18 +306,19 @@ mutable struct model{Mesh<:AbstractMesh, Friction<:AbstractFriction} transient::Transient inversion::Inversion calving::Calving + levelset::Levelset end function model() #{{{ return model( Mesh2dTriangle(), Geometry(), Mask(), Materials(), Initialization(),Stressbalance(), Constants(), Dict(), BuddFriction(), Basalforcings(), SMBforcings(), Timestepping(), - Masstransport(), Transient(), Inversion(), Calving()) + Masstransport(), Transient(), Inversion(), Calving(), Levelset()) end#}}} function model(md::model; mesh::AbstractMesh=md.mesh, friction::AbstractFriction=md.friction) #{{{ 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.calving) + md.masstransport, md.transient, md.inversion, md.calving, md.levelset) end#}}} function model(matmd::Dict,verbose::Bool=true) #{{{ @@ -385,6 +397,7 @@ function Base.show(io::IO, md::model)# {{{ @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) "calving" + @printf "%19s: %-26s -- %s\n" "levelset" typeof(md.levelset) "Level-set parameters" @printf "%19s: %-26s -- %s\n" "results" typeof(md.results) "model results" end# }}} From 496a6457d15614eeaaab20a88c73d57175fea38c Mon Sep 17 00:00:00 2001 From: Cheng Gong Date: Tue, 2 Apr 2024 16:03:40 -0600 Subject: [PATCH 09/15] CHG: use Abstract type for calving --- src/usr/classes.jl | 18 +++++++++--------- 1 file changed, 9 insertions(+), 9 deletions(-) diff --git a/src/usr/classes.jl b/src/usr/classes.jl index 22dc564..340e880 100644 --- a/src/usr/classes.jl +++ b/src/usr/classes.jl @@ -264,15 +264,15 @@ function Base.show(io::IO, this::Inversion)# {{{ IssmStructDisp(io, this) end# }}} # }}} -#Calving{{{ -mutable struct Calving - law::Int64 +#Calving {{{ +abstract type AbstractCalving end +mutable struct DefaultCalving <: AbstractCalving calvingrate::Vector{Float64} end -function Calving() #{{{ - return Calving( 0, Vector{Float64}(undef,0)) +function DefaultCalving() #{{{ + return DefaultCalving(Vector{Float64}(undef,0)) end# }}} -function Base.show(io::IO, this::Calving)# {{{ +function Base.show(io::IO, this::DefaultCalving)# {{{ IssmStructDisp(io, this) end# }}} # }}} @@ -289,7 +289,7 @@ 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 @@ -312,9 +312,9 @@ function model() #{{{ return model( Mesh2dTriangle(), Geometry(), Mask(), Materials(), Initialization(),Stressbalance(), Constants(), Dict(), BuddFriction(), Basalforcings(), SMBforcings(), Timestepping(), - Masstransport(), Transient(), Inversion(), Calving(), Levelset()) + Masstransport(), Transient(), Inversion(), DefaultCalving(), Levelset()) 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, From 675f0be6281c0a2174563c860408e1f775f8d4af Mon Sep 17 00:00:00 2001 From: Cheng Gong Date: Tue, 2 Apr 2024 16:20:05 -0600 Subject: [PATCH 10/15] ADD: md.frontalforcings --- src/usr/classes.jl | 24 ++++++++++++++++++++---- 1 file changed, 20 insertions(+), 4 deletions(-) diff --git a/src/usr/classes.jl b/src/usr/classes.jl index 340e880..15bf3ce 100644 --- a/src/usr/classes.jl +++ b/src/usr/classes.jl @@ -287,6 +287,18 @@ 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, Calving<:AbstractCalving} @@ -307,18 +319,21 @@ mutable struct model{Mesh<:AbstractMesh, Friction<:AbstractFriction, Calving<:Ab 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(), DefaultCalving(), Levelset()) + Masstransport(), Transient(), Inversion(), DefaultCalving(), + Levelset(), Frontalforcings()) end#}}} 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.calving, md.levelset) + md.masstransport, md.transient, md.inversion, md.calving, + md.levelset, md.frontalforcings) end#}}} function model(matmd::Dict,verbose::Bool=true) #{{{ @@ -396,8 +411,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) "calving" - @printf "%19s: %-26s -- %s\n" "levelset" typeof(md.levelset) "Level-set parameters" + @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# }}} From 6a8af851fb6c020767a8f152a63469b691bc8a52 Mon Sep 17 00:00:00 2001 From: Cheng Gong Date: Tue, 2 Apr 2024 21:41:35 -0700 Subject: [PATCH 11/15] add levelsetAnalysis --- src/core/analyses/levelsetanalysis.jl | 232 ++++++++++++++++++++++++++ src/core/elements.jl | 39 +++++ src/usr/classes.jl | 6 +- 3 files changed, 276 insertions(+), 1 deletion(-) create mode 100644 src/core/analyses/levelsetanalysis.jl diff --git a/src/core/analyses/levelsetanalysis.jl b/src/core/analyses/levelsetanalysis.jl new file mode 100644 index 0000000..8ba8c04 --- /dev/null +++ b/src/core/analyses/levelsetanalysis.jl @@ -0,0 +1,232 @@ +#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(spcvx[i]) + push!(constraints,Constraint(count,i,1,spclevelset)) + count+=1 + end + end + + return nothing +end#}}} +function CreateNodes(analysis::LevelsetAnalysis,nodes::Vector{Node},md::model) #{{{ + + numdof = 2 + 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)# {{{ + + #Activate formulation + SetCurrentConfiguration!(femmodel, analysis) + + #Call solution sequence to compute new speeds + println(" call computational core:"); + solutionsequence_linear(femmodel,analysis) + + 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#}}} diff --git a/src/core/elements.jl b/src/core/elements.jl index 060682f..09b2fb2 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, value, 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 diff --git a/src/usr/classes.jl b/src/usr/classes.jl index 15bf3ce..dd48626 100644 --- a/src/usr/classes.jl +++ b/src/usr/classes.jl @@ -279,9 +279,13 @@ 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)) + return Levelset(Vector{Float64}(undef,0), 1, 10, 1, 1.0e12) end# }}} function Base.show(io::IO, this::Levelset)# {{{ IssmStructDisp(io, this) From 4dbf2bf45e6390847f43762a9e9702a6a2fa8184 Mon Sep 17 00:00:00 2001 From: Cheng Gong Date: Fri, 5 Apr 2024 10:25:24 -0700 Subject: [PATCH 12/15] add Levelset, but singlar K matrix --- src/core/analyses/levelsetanalysis.jl | 14 ++++- src/core/analyses/transientanalysis.jl | 4 ++ src/core/elements.jl | 59 ++++++++++++++++++++ src/core/modules.jl | 5 +- src/core/solve.jl | 1 + test/test807.jl | 74 ++++++++++++++++++++++++++ 6 files changed, 154 insertions(+), 3 deletions(-) create mode 100755 test/test807.jl diff --git a/src/core/analyses/levelsetanalysis.jl b/src/core/analyses/levelsetanalysis.jl index 8ba8c04..9cf530a 100644 --- a/src/core/analyses/levelsetanalysis.jl +++ b/src/core/analyses/levelsetanalysis.jl @@ -10,8 +10,8 @@ function CreateConstraints(analysis::LevelsetAnalysis,constraints::Vector{Constr count = 1 for i in 1:md.mesh.numberofvertices - if ~isnan(spcvx[i]) - push!(constraints,Constraint(count,i,1,spclevelset)) + if ~isnan(spclevelset[i]) + push!(constraints,Constraint(count,i,1,spclevelset[i])) count+=1 end end @@ -78,6 +78,8 @@ end#}}} #Finite Element Analysis function Core(analysis::LevelsetAnalysis,femmodel::FemModel)# {{{ + # moving front + MovingFrontalVel(femmodel) #Activate formulation SetCurrentConfiguration!(femmodel, analysis) @@ -230,3 +232,11 @@ 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/elements.jl b/src/core/elements.jl index 09b2fb2..abd1619 100644 --- a/src/core/elements.jl +++ b/src/core/elements.jl @@ -412,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..80b06f5 100644 --- a/src/core/modules.jl +++ b/src/core/modules.jl @@ -27,7 +27,7 @@ function ModelProcessor(md::model, solutionstring::Symbol) #{{{ if solutionstring===:StressbalanceSolution analyses = Analysis[StressbalanceAnalysis()] elseif solutionstring===:TransientSolution - analyses = Analysis[StressbalanceAnalysis(), MasstransportAnalysis()] + analyses = Analysis[StressbalanceAnalysis(), MasstransportAnalysis(), LevelsetAnalysis()] else error(solutionstring, " not supported by ModelProcessor") end @@ -117,6 +117,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/test/test807.jl b/test/test807.jl new file mode 100755 index 0000000..d812662 --- /dev/null +++ b/test/test807.jl @@ -0,0 +1,74 @@ +using dJUICE + +md = model() +md = triangle(md,issmdir()*"/test/Exp/Square.exp",50000.) +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 = ((x .- alpha*Lx).>0) .- ((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) From 0f7778d4b7d610a406d889d724b36ef1590a0f20 Mon Sep 17 00:00:00 2001 From: Cheng Gong Date: Sun, 7 Apr 2024 21:23:28 -0400 Subject: [PATCH 13/15] add levelset, need to compare with ISSM matrix --- src/core/analyses/levelsetanalysis.jl | 5 ++++- src/core/elements.jl | 2 +- test/test807.jl | 4 ++-- 3 files changed, 7 insertions(+), 4 deletions(-) diff --git a/src/core/analyses/levelsetanalysis.jl b/src/core/analyses/levelsetanalysis.jl index 9cf530a..b0c3a45 100644 --- a/src/core/analyses/levelsetanalysis.jl +++ b/src/core/analyses/levelsetanalysis.jl @@ -20,7 +20,7 @@ function CreateConstraints(analysis::LevelsetAnalysis,constraints::Vector{Constr end#}}} function CreateNodes(analysis::LevelsetAnalysis,nodes::Vector{Node},md::model) #{{{ - numdof = 2 + 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 @@ -87,6 +87,9 @@ function Core(analysis::LevelsetAnalysis,femmodel::FemModel)# {{{ println(" call computational core:"); solutionsequence_linear(femmodel,analysis) + # TODO: add reinitialization + # save + RequestedOutputsx(femmodel, [MaskIceLevelsetEnum]) return nothing end #}}} function CreateKMatrix(analysis::LevelsetAnalysis,element::Tria)# {{{ diff --git a/src/core/elements.jl b/src/core/elements.jl index abd1619..1514889 100644 --- a/src/core/elements.jl +++ b/src/core/elements.jl @@ -247,7 +247,7 @@ function InputUpdateFromSolutionOneDof(element::Tria, ug::Vector{Float64},enum:: for i in 1:numdof values[i]=ug[doflist[i]] end #Add back to Mask - AddInput(element, enum, value, P1Enum) + AddInput(element, enum, values, P1Enum) return nothing diff --git a/test/test807.jl b/test/test807.jl index d812662..8a09b18 100755 --- a/test/test807.jl +++ b/test/test807.jl @@ -1,7 +1,7 @@ using dJUICE md = model() -md = triangle(md,issmdir()*"/test/Exp/Square.exp",50000.) +md = triangle(md,issmdir()*"/test/Exp/Square.exp",200000.) md = setmask(md,"all","") #Geometry @@ -52,7 +52,7 @@ md.basalforcings.groundedice_melting_rate=5*ones(md.mesh.numberofvertices) # mask Lx = xmax - xmin alpha = 2.0/3.0 -md.mask.ice_levelset = ((x .- alpha*Lx).>0) .- ((x .- alpha*Lx).<0) +md.mask.ice_levelset = ((md.mesh.x .- alpha*Lx).>0) .- ((md.mesh.x .- alpha*Lx).<0) md.levelset.kill_icebergs = 0 # time stepping From a3f6ed3f1a65f2c81b8247635e5eb88bc1b20265 Mon Sep 17 00:00:00 2001 From: Cheng Gong Date: Mon, 8 Apr 2024 09:26:00 -0400 Subject: [PATCH 14/15] BUG: set default `md.transient.ismovingfront=false` --- src/core/modules.jl | 6 +++++- src/usr/classes.jl | 2 +- test/runtests.jl | 2 +- 3 files changed, 7 insertions(+), 3 deletions(-) diff --git a/src/core/modules.jl b/src/core/modules.jl index 80b06f5..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(), LevelsetAnalysis()] + if md.transient.ismovingfront + analyses = Analysis[StressbalanceAnalysis(), MasstransportAnalysis(), LevelsetAnalysis()] + else + analyses = Analysis[StressbalanceAnalysis(), MasstransportAnalysis()] + end else error(solutionstring, " not supported by ModelProcessor") end diff --git a/src/usr/classes.jl b/src/usr/classes.jl index dd48626..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) 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") From 193fc68485382cdd3fd38cf3b30300a1b84ba652 Mon Sep 17 00:00:00 2001 From: Cheng Gong Date: Mon, 8 Apr 2024 11:53:23 -0400 Subject: [PATCH 15/15] add data file for test807 --- test/Data/SquareShelf.arch | Bin 0 -> 25767 bytes 1 file changed, 0 insertions(+), 0 deletions(-) create mode 100644 test/Data/SquareShelf.arch diff --git a/test/Data/SquareShelf.arch b/test/Data/SquareShelf.arch new file mode 100644 index 0000000000000000000000000000000000000000..a397911a54937c9ad99a493e0c6afec348fda17d GIT binary patch literal 25767 zcmai*30O_r8~;y+WG*TdaheB3gk*ddk&={>3OUVlohA}$3)f7U=Z-`&W=u#sWk|+r z7Va_66e9Eg*lWGpKhOX7{Lk~`v%lZ>UE|(+oxRsRS0a(LkVqtIs=r2)Q&nP@>Zq~G z@Mf+Ge?D;N%ZPKQwP*2%v-AvGWov-JnolIy&G=6SJ%V0Z>^37+`ZW}}4 z!^N3hCA~ngV%tGGpJBi?y;py_E+5ufeQRf4XAjFle2Z4Q)k69U-6j`@8NrS#i>`I= z zkBf#qJNrVcjgjKJk0&1(__KVqsRQq$cC$x6ZPmQe4WLx;q+g` zdn@_miXR)6zn{iWzmf6w{gqTcN#5Chce6@z%6j^y_yV_12b% z^>MX32VOh>j?Z{`?aQ6?rJ|o2JmX#R&hFs+}i zjSDY8vV27e*P#TWdN*Er;mr>i$jxY8c(NNGF<3W!PDnAd8(98j%fsvZ@^vLQmp0Re z;Rmh`pMEHVpVcL?kLS@ie)gn$W*-8U!U+AwD~tP@@a`Ladgyo4;FE{f+^$K=;^W73 zlU#mj#%Gr4Y*fp=#JfdqF8FK~$S-%9{jANgt^9bCr2N!pW&DDt#qE?M_4!#V-QP^s z%utnIT(-7wK5xAH>Dl}8E`z`8n7I{a*76DG+o&z@?7?s9u_3Kv&-LOwK%e3F$RTIY zn%MjMnhr%UWY)uhrDr?yz3L`fm0f5Kwk5Ydui5_@Jo=|7gY6A@pHhR!($QmiuXU^R z*L^<4dkyQbaC!YzUiQ*q(6b3Ud2ZkhF8k19KBCMJ9&Ct&Aibd0OHZ8PV?#D+YTwwz z2V`EI@T^w`AL6BSkY;q@J!_4&b(~_vCj|{E9lBk@CwBR^ZhqAOK6Ua_J+IdleDq#Z z^WykKKJt)(;mGb@{FqM4pC1#-~Q;f2j4* z=hFs6YhOr7ELU+w6^5Gn|xAf$cIlGw({~vo2x9L)ZN1)gD1@#XdVb#WN{FBLDePbGU`Tiq z6ng$&CCJlt`&6xY4k`29HZ3da3&{f(j@HW9fae9$CiZE&;NBgWl2(8HV9qN@U1i&L z%tjqPV{H2Ucrubt-{^h(R4&iY469S$GDnT)#_Stv)%O@bGb8=KxFb{eDK|cZjJ8>X z_3^3I9W=N`ANVxby!!gNOMJ$Y%%T2@E&L3Pn2old_wi{%zeilH*eB*E&MlHc$t6)o zG-;+$M9f7UcRtU2edjbj$?0nA)9$~(kmJ^ zKX3+?gTSOub2H-~!Zc4!ZqbDHFg1 zPKC*@%=kIEUHI5rrHyxvkn?^ER!`NK70yTPDoA-~yp&J+GPu{__eFea?4eOr*`Z=R z9xvk8U;n8G!~SWq?!t=hFtI+p#Q%B;gr)bKW!ftWBFonBiG z6Rae=b|=;GGpqB*!q_~1%7e{)?|lg7_i=L$t=@PaMw{GiIpt;?zq&lOmF~vtSU(>! z=BaM0@tgUCBUy`ly}R;jW9;p!Z3Foo$v+YYn+iFqm3t^4jTY>Yu_j-N7d6=mkW!~{=rI`^CJd*1`G%7A=J6);)go<_sXd+gZw zwTJ%X`~u8JLfr8-hx-nn2+8w{1G}GU3MmmE$MnAJ0jVp%%Uyj3BrS8j|LMV5h}d1N zEL<<%pJMr_R%rnbZ`wo5f+Mz1LZ(5~;bDq(%XL8B^mXE(u?dhCymq9$mNi7wT~%-8 z+X3W10-mp&X#=shx<)tswFUy)UiVMg)LpF4|GItDDUUb^XqDOaLmO4OR;JrEdZhAE z2hPuEacKl3z75C@)i8v3=OOy*nnyvxm|ME;7Jj^UhJSsBUBCGLCg&&Z+0jL;uQz$P zeeE$XKF;LRiV9g*eo0-o6Lo&w_(Q$qy1+Mu10(R7C4RNlFgABXyBQMh!3$DV62j(MJSiittBB3Jh#kX&I%|*?3?)k&O@&W8#v8r9y{YJoRp8L1a&P#A4 zt9j|#%T2|+vYJUX%F{Vep3ywO^F~W3zB2XWtwZaf%yW?Am?b$-(fHsBhtP>oF@DnC zV)-Q~|96?8=9I;7@?7zV4^<^lZZ+`c>xl|DnLaZ1#_TXSniCb6`ne16^Oi@{ly-xY zgTl6d&bbR^Ir3xMn(9EQQ_Ex1TXcZZ{d(#d^CN+8dt=@>=h;x!C$Tznm#STFU(?`@ z^>A`&&o-K0e?nQCCk4Y&XThP2;}30WpTO}JLmhh@IRi&$o2GU0(SpNoE;Y`x^M&Gf zrH+nC)8K^u(-vEs4H3t8;#6PRvg4znu(bR83ClyEWVJ_a`nYd!^up`xE$95;VBV|E zk#7RU{O$g{gBB4n?-ccTPMW0Q*lgnq?*NhHXDS z)<v;lK_*WurwuziR&dipT`m)p^Cp?P>%HwHC~;yQu-YE3e)&{S^*d;(wf|+;j$V z)?ckSSFs8T3e(LyjC%sR+dO}_u+<~jVSGNU^{+=_eVa3E#|<#hggm*?Tsym3(O;3C zcP3`XpKr^JkF;xADorU8HQ)YoikuayefS#v%bJ34@3?#`A?b`OP$e;)d* z&20r`Z#3&Gwj037I)(SaIuj_rarW^tZX%SAY3J|tG6Tv>zm&KH&WDOVeSPmWS^^ad zdTK24n+6s7>+-h#vm8!7Uh%TkmS<4W)-=N*Z5y1d81d<5(GRFt%D+v!Wdjw`@v;c- z<52D!Q9g3UMJV6eHuU<4csP0Icu?i7YABzxdx_oY^-#{ek{W2bz{w3svefP_P#$1y zP_(fM$_L~eKi;wvoYZi#ebMd#oaFmh9*;_gvOJCcRnYo zy){}izHlAN`;;EK-N{y*hqCy|y0iD#!%6%6+sjSQLs?1Lp1vVhp)C9Q&0RbCLRqZ(I$fR}&YvTsrV@EzyA{kmLrf9!tPM}3Iueo9#Sb@A!rP)dv+U+z)@r9UedO*~f# zCv}<}t0!xqO!@CYS?+u|YVj+7GCC1|`y?H+yHzgpxshmwk7|rOP>OVR602wdrT~!rDDWihc_o!m@GeB>Sbo zuv}-PLqe`<-KQVO%#;TX?)*X2sr!3fsWq(;b>`f`rG8%1V7AGGPtomE*IToT>u-f^ zfwI&W=bHxq2YbC^`#XC5fZVI;{6fhxv0Q1NF#mFYb=dd%{qBrP)%U=9BeySiO<;f6 z%5m4s{yZgYd3#}%UiZ&%D9`pBd3+4E?JRglnz_Sv`^3rFnU%1osH%Lfo*%4h=;5?BFDqQB%bB_><{NzXh!W-u&MF;TtaiKZ@(*7uE}k_r7s`q>~0E zxAqUsT`?F+*5`*`x)u+|!Y#W`Z2B3FJL}Dw+~yq|y0Rs7w!>D~@BEgJ`~GJa=HUE` z{dMbW;b3I+>xHYUp|F!IdUst1vE5}p4|Pw}4Hor=Ng4mFyjcXxMjD0nh|GsArc10c8($U|;x# zrZFeiz=;!2tfs&GCDy+=rh{qaM+NMh;+f^2tqUic=eyk*e-(I(tcqTdD}J=jKDiZE|1#NMvHSoOOf>kX z-?EF4<7BewZ2oXDFE=^wsm{!)kagxttHS1*uxiZEi`{Rwg*@Yt&wFke2GcWEk4i*ev(99$A_fwAWfw^gEW7Sq|f|<38oh=iNK}ML>_xO!k7)OFUeC5RN zE5uzB$m*Uj@bj&qkT=Y+$v;o0K*8ki2CWVLffdS-c4zwixj$0niYLb%GGIwSYR96H z*I`!qniKc6j9|;1p?7a&RKr?s`07F4q1X-*)}B(n%sv$e1zlP-`r6|&th~}RWf{zb zy_Q9DEevl!VR)nIJJeRgTBlzTMefd64hdV9jZE9sb`9hQ=S7|T=O$$D>DH`#L<}r^ zy*_VCeK%O|(&6W;PyJwHaqdd{xaF{K)U%fhnw)}z+h$$XKmJ6u&ICBL>HKHi$lJJ~ zXGn1s2gC@yF{Y)SVp*ymPo-Rqv} z`qP-1J>f+*6jvSJ>T~6a*uUJj)90O1y}fPUu-3-xTR7~Wd&aS|@*?c-f4g9K_ua5; z&9;85lv6WC@p*MR6_ zIn*ZYE(>|q2ZRelERVjC^*kYt+Vq>AppjCI8_3{xcT$RI%CIn|}YX146ila79w+z2NFOrxYdi&()vp2-}C8-;D*p--;o*45yFy~LbucA6GZ(c|& zJJ0tVGG#cizM7=@-M1^TYCdY4Eg=@|{Z^b%A49Ac|2lN@bdZ>Db2Mt!JbO8@@#q{L)yRuj&6r$pCgy`;uaxcxCuYec^U=2iVtjKXmlM>2m`jOrQbjhg ze6U{GdqH1f9T?JLLATq)#>TkS<(Ln|YJ=yj6AAN)Wf$Y6#yO3M<(b?sjWJ_hD8{@MvP0t-te%WPs>*4Kvv&4ND}}W-XRHeLtR9y{dR< zW4er3-;ZlkKPy=^k3LPyy-SIW(;+v_mQrH%xoKt7^{t4-8h>KDwK1`3mp7&JPgC*y zTdjC=@mb%O#Hv&?tWT7V*p7{DrG3s>w?F!6f_5hzC9&0b5P4@+${+veh}ceDwh>zk zGT5lm3)1z_JInj^LkM0^czv1-k9BYBoK2+CPpb{y*Ak-@9i0XZb|pre43y6!dlRYe zyX@JXW~6ieRWJXlVWi`Q3hk?Fr;*O>UYyW0d_!~}$%O05sP&r0r=8ak19lzjvFoi% z?a%8o2OJ=sO6wlwl{_UnS2{KK+#EvmM^#KRoY#lw?su61UsH(QrM)igJu->5zs9J7 zrq76h?`3(a(MO^+Xnebl+gcIrjss&VjpB&bwcVee+)yLZ0UOla%eD~XeFHz|p710c z&wugRl6R1_r}xwEKK=fX_TO)eEV9bG2LGj*JSZkV%o#C&7jGdq{Wc9&Ig9v zCECL)4Yv=uO-y4>FPPKvI5C~*m>*XgN{l8<%U+f6mKZPVQ4*%rgGfKQM7`^pLQKuC zX?>gel9=XiTi{jePRx6)b6jh8ju^lD`pdV`XJV`|<6O1$Ixz|wGjGxzS0cSNV9Kn2 z97*S7OWiwT9Ef2?_v;H1^oT{TmU8!b$BEg*T%pnHf z23Mc#kWEbX%PZcsaVN&^HxxQun-JsmS4SoIWRfn!Ht2-(ze7yt_SH-;Hz6Hkdpz}0 zA4qiO@8I2Dm=e7L$D{RwTaga)Up|g)ltemsE{&h!+mmQU)w>0sI88c?-8W_XCUw%j zAU*O{v?||QrzE@2ZPKRszggvT@<{s|7nYb6xD(B7-cLBqwWLYG)^BbFJZU-pT*2W2 zPBCu(_7i}V47bj)8<qxM8|Jd5KeQIG~sqz;e ze)@15G4NX$Y@Vt{bl=XJSzVS&noqWMb=Ec`(v^oFEV#Ch;QFj}J9LEI(G1e!cE5q+ z$9j;a=8tw=(&#|iTF$z2HE_c3wtjzzWW=D}XD0iT_B*TRrDisU`cD()8ck_VjOxD* z?QFl0NGE>b*6H{WYebGuczM|50Fty`(yZ?r!rq#5(M;`$-&6C#WQ?|;1S z4M+XxcXh}eVt8oLyPfTG#PvXbsk+y_!c?N;)B>)muJaa7hwGNc3*T)*J)M5duO4tn zbsx@H)O*A$Ytm`pBE8thg~afq|C-dPElEfHm$$=bI}w9l>KWz^E#YgsTcI1hSChs@ zzsAG2+3<5V_q!gbTQ_uc?+2$z(*T#J$p_zw=hNoP-9_iqInwpqT>m?M*NLX7KfF&e zC+d%6zsr%vblrejY2beSo7$xLcJ%=FgfsACbJddtUqVQ$&ckD)F1{okhcxHz_*W8> zAlGZV#}tY6Hw))q-3#eM)Ruio8+-mLY0kABy7ia?X>@+gSM8PUez)=a1GQeu-f=G=Frv+Pwaa>;vuWhmzcNX5x=6?4B zKk!`9ZN=aGA*qq9jOxjA!J+M%Ub4hb~_SI)c_qhu=^ORkv>u`9o;m2P41r_kUbkNFnKk|s0`r+KP*0Jz2V6Zgu zu5i!)t`nY}_z)L;bQ-+0kBR*EWE9l9AA!o~1;n(|!PRYb_mgHvf_!f%WTb_L-t^%4 znM5!3_NLW3IYjF3zj%7@EYc+_+1$jUm)O6@+P^+-*S!VLAN)I4bM$(sJw0s6xocDZ zo_AtxJL|k%d8dxW$K3 zTm(@69mjg{q4RHITwGjnaND(i_d#seR@b$G#(b+`M+M)e!v zeyc%FDT_Z=5{svYip~}KN&iA0kb6KT0etsvTpBRS)Ux z;j3dt_SC*DcZrL-?x5$6_r9kQ{mDDuUf45Ee0VeP-acpmc}hB*&1}8iWgjs(wmh=C zR%>FoF7v>MkGF}I>q9#S(;%YTsYTWf%X>sSaB1P3Z#hJ>U$%62ekRe%{gAg({r&HA z_WMJc1sm*I8(c_QZ~S$nkA|bTsI}<8*)?%}D1Iq)n*4s#Oy^TXzhhqc%NN$-BBaG$ zdfGqOJfZ0tN76n*GyBc2v*LAw=ilj#bBJY3|A$1OdN_ROcE8>wnV9w2n>laVVq#-< zqqNGkL~I}LJM%u7hq`?YCZ?Bd@4DSeC#I$|?fuo=h{@}Ov@}CAVsduTG-t(VVia>> zP40iwiBZaqNsp?_h)Gy8Nu$InV%k-I%b)-yF_sw2Dm=27Sm_$tUUT|REViAg4vDKF z7CIVRzUp))9f!^+YBNHgbdc_O(=Mlo7_@k?toK|C(kVLD%-`cO(f;_k<-{rtqSt1g zi+vMyqT|!`w8G^8(f!iIUCsRxX*Y&%wcvRo=~5GWxSR8T#7KVV=+^*Mzt$^r%L43) ze#^XenaWTSN>b5V@C}7 z<9Z@qcdOl{OY7>HfpeDQZ(IZQ0q5j+`z1dspiR6iPVt{7|)M+(Ya(QGT zY3x13E+OTKnE#^3;BK4in~3^tI6u-VmlMDI@k{^Pu#NK*O*DyddD@~y)z-xHgvs~X z-bSQ-=+*EWEu`?JspKtp@*CDeOm=SYeo_}gbX&iuSiZUm(FtrGzTSKh(QbF3;$qqw zalO#Yp8cQhmKQ|3CdR9qeGt)Gnb-1c<~5=wb(j^^4~Tl>4X3Pj7eH;_IR&4e$HVJ{ z$AP||^5OI2mX&hOLTt~p*gLmnP!%zo?l^6g=SyNTtaGP%H|`SC&Eso^^vowlj}*%< zuTLOmsj(}SXHAHO$;m|-Z=R_7`K-iY<_ThKm>5<4GLdLat>2NEA^vWHXz5uQUVQAX zdZ^7Z9n(01v?yre^7GeZY!{xja=1JpELhajNoOg?&#Lyy6vOLKCN_yP&%8786U*6ZeLeAZ=3rtYOR7A6sxz@Unzuh_ z#yB0(9LQWI2gR=%t{@_cb?mF9wCb!3j64JI>wSij%r1AabUw&%5 zCe32%ynmjaN?PnqGOoM7LR`N!e)R1fZgdo4!_&voo6^Q=IUDMGENR^M|m9#2CVMt?XY*+81h`O|2~@t$q4^wN$ewZCX~T zCpS{9pQabAY)gR{H@~v;m#aQ8S+#JHZDJo{Hh*H-)0?V?Q?>~xK6hwFj89WNw3~W& zAG*1p2Rs|P_*akU9YmTlq{=tEni%`Xxdc_jkrsDj-L`KSMYL`W{?at3is-IAe14&+ zqxigFar5_MmV{U>Fqrrx`8ctDFz4jx<@?0z44=$3?B^hs_rKy;4~`e>M?I94Cyc{6 z;Q-i=6kr7NozQ2*a-*5{_=`?VJ2T>WkTD$~5k6u&9!!h%GL8A{KXQN@^O4xT3-j@M z!FJJ)R|pB>c!qFP<8>jEQJrE(X)Qnk9zb_Sco^QAi<80%MT8 zZ~!=NcaCao2gfy1uq6=B;^O=;js0_Hd0vcTc;Ntsa+EkP*w2y758`PY`-$}hDXEWf ztPk%~9EXpPrWzMV^Tz-s`kp-1qXqn#PcNYUSmt~2RF7lYlcO4a>>uW1+?(Zx%V>TW zPWKpWuX4@)Zz3FqZQQ0@lOD%0|iI`%Lj7I=NXv}7Ni8kxQxa}{p~BLAE;!MHPpCXV0}SS z7UwCGIm%!KCC-<&ASHm~M_lc;}mQp`NfJA)$ zzk`%1{9k2$4}80+f?D%4nwl&PP8~ zLGwacUW}CL5Tuf7jAMJz%nwn}I1=-)e5jPhF+Y*zqaVuAcqC7W^@K^O#(1cbYV;!& zOf%v<1uL1)F)IGj%7%WplQcr*cVx?3k*fnIVjB1RF^UD0Vzv9@B5TNCv87Fd7#|w}E z<|QkskM+cHRL95~`cbGk=JSlesGv;dC~@79vw2Kr@pu`H#{=cxeiHL#JdI;{7ePv3 zG~@99q2+A;y^I7?1o|cbOQeuDbc>}LInGlwM8^`h#lo3is0TPgB zQJtiqi~-67O-geCs`pq%ylxU%Ik}SN$rY3dQbrF(yT1rbYf{Qy71Do(Idbqr~et zk!f6yaJ|I)73VdfhL%gpV)b2S^q`CtpgMm_s^u(>>t7N_{X`Z|Wbq_T=JS-1e-YP} zMCN1rkxH76^@j;k0?AV9C&?(WJ=94U2O2jPpt_!9REu$*YV2pCf@v8gUY|Hm*dMVx zu)NERaz-U1PZ`aaR6}*LCZhmVzm!zt{hq{Aoh(2VXZ1!(sUM}F3}sCApgLMcnGBS1 zQpy+^Ww3%W5Gdn#%4AkQwinOh$%QP=^5O)kI{z$=I>nCW$)%LZmnjDekbnkbij?}v zHIym9s7X1Lqx9e@Qv^r=II7(kwV5BnIEEJv0O!S?`9m1-y>AFd^I~O`F$zkYuXsWJ zS--H|fz0pA=*tNQ;`?hq=J#Oi`4{^z?S};60D3b`r%wU{n8tk6-I$N@W5UuOIxr70K#xVU#iA^^5(nX7#$W zINm3&%*Q;WOeO?CE~gBXGRi2uILa|V=`KJ5$dRZy8pm-CL;o*g9O=X2UW_&!^#!iA%?SpgF$2@Ew?`s?v`d+O3SXK__#fSO9 zLUj`G0^tAxnUC6=`C%N5d$V|ul4`7HtRN+T@o*W9BRzrYF^pnA1X}_)pFS*(^@`+# zL>R~NMY8;HjP@*!>j&;%5$5Cli2E|W%=c&IG49I=0dd{x#q#hu&Y!1o@Vht^|81gNAt%2rLjHi&jcau&%PGM4jjKEPHC_uFz229H+ zgQb*szXUQL^ROQ{t{@>@wJtE?{U}l?1XTaSIVI*JvEE>e1I@$p6eOb>=QoI>S|+8$ ze&BOZ5Knz63b+Z^m%O2*H*>IAerNH~{P)wif}^ z$9l0o@j77TCUGo}XH+&Mjx#|d3fJoVXouDCxJp`?C{ zlrmC5nOM`1!5r05j4=wPffDBn`;GTWvK`CgD6#)Io+t(NqXkI7k`Xmt?>JsuUyU)3 z6Al3LVwFtGC{f2lLqA$cb+VK)i4n&o&KF1H@NpuCwu|!+ z$FVq)XPQwiNY(p(7WLzqM&kQvoRY@zJsF>8#F z55~fVj1i;+kl0=l$MUlNqJn7ws{O0~QD2gj)`7+qN2 z0Hz5ewlkP%oQFY7_hIbI=m)|9jN>Tja3$itSwF#+0Qy+of%!=6M~alzlUzgT%kr>& zY`5oMeh;RHaKZuL_<;G?4`6ZZ51t2{PiL0zPYZ$*%Xeb_2&SbV{~Q2eh(GGe}0FE{lM!6>mSbYhB4y&Vt?^{8q2$} zc05`BNJbA{H~_p(upihj>?f8FWch(oT5g;G)jkj-mdE->Gaviw0fGTIA2^@x9QASD z{w|N@JXqckAzcEPkK@4c1xkek!22**Nj2WDV}QjuN*|`Nzg|4`y*NtD$9~{=v7g?+ z@)^f)RQoVud>o4hDrp>jY;UY!tKP$S8pnQ*<*3Gd)Zq&1hw+r?i~SbtKff!(ar!Xd z3urk!k3Kxr*l(;yB<6E8&zGn46Yyt!LXG=5c)#L(Dr#09*B@MuaG%JJqxJhU`my}+ zg00#wVwB5-1Muf5@%o=2NY%a(M|~`h^Xp%a95hOp`nlU1k0!|1Ii#qQFGMC^M>aI z8KdyOot(v zkL_VQAqvUwkGMa9*G&lb+y1k!7A;7XmI?>r@84p1s>OOZs-py`_RnNghXQ4Yk}{5` zjFwV{Dk#IGlvp0~Llx8yl~N{gOv@<41gM_tl~iMUkqV|6F)va^eQXcs59c+E<$Loi z4;Yya<{D}|FX2k+V?0ujfA+1ycp4Aq7?qR}GNz@Bz{oRllo5>h{UB#LoTu?{C8hAM z>ik8}1Hf{~a0Si7cCjARK^%*-JaL^7Y}N0}WHcVhQQ~+LfoeRjQ7j(G{7@y0M=&NT zsE(3SisuWcPPj~o<>Xln^J1h_h~b5-#9+><7G6DHlu2fryBinOpE29VI2E4f%%g-8XqM<0d3PfBA{DSPim6V}O<2oJ9d|bbA9g72+ zkLz5Jlxkd;;y9|Yf7q`$CH0ZXg$*^X3q1r|wSHWtar7OTp90huwUp`<0jhm7MyzKf J$KoLH{{SMY*AV~! literal 0 HcmV?d00001