Index: /issm/trunk-jpl/src/jl/core/analyses/masstransportanalysis.jl =================================================================== --- /issm/trunk-jpl/src/jl/core/analyses/masstransportanalysis.jl (revision 26730) +++ /issm/trunk-jpl/src/jl/core/analyses/masstransportanalysis.jl (revision 26730) @@ -0,0 +1,272 @@ +#MasstransportAnalysis class definition +struct MasstransportAnalysis <: Analysis#{{{ +end #}}} + +#Model Processing +function CreateConstraints(analysis::MasstransportAnalysis,constraints::Vector{Constraint},md::model) #{{{ + + #load constraints from model + spcthickness = md.masstransport.spcthickness + + count = 1 + for i in 1:md.mesh.numberofvertices + if ~isnan(spcthickness[i]) + push!(constraints,Constraint(count,i,1,spcthickness[i])) + count+=1 + end + end + +end#}}} +function CreateNodes(analysis::MasstransportAnalysis,nodes::Vector{Node},md::model) #{{{ + + numdof = 1 + for i in 1:md.mesh.numberofvertices + push!(nodes,Node(i,i,true,numdof,-ones(Int64,numdof), ones(Int64,numdof), -ones(Int64,numdof), zeros(numdof))) + end +end#}}} +function UpdateElements(analysis::MasstransportAnalysis,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.geometry.thickness,ThicknessEnum) + FetchDataToInput(md,inputs,elements,md.geometry.surface,SurfaceEnum) + FetchDataToInput(md,inputs,elements,md.geometry.base,BaseEnum) + FetchDataToInput(md,inputs,elements,md.basalforcings.groundedice_melting_rate./md.constants.yts,BasalforcingsGroundediceMeltingRateEnum) + FetchDataToInput(md,inputs,elements,md.basalforcings.floatingice_melting_rate./md.constants.yts,BasalforcingsFloatingiceMeltingRateEnum) + FetchDataToInput(md,inputs,elements,md.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) + +end#}}} +function UpdateParameters(analysis::MasstransportAnalysis,parameters::Parameters,md::model) #{{{ + + AddParam(parameters, md.masstransport.min_thickness, MasstransportMinThicknessEnum) + AddParam(parameters, md.masstransport.stabilization, MasstransportStabilizationEnum) + +end#}}} + +#Finite Element Analysis +function Core(analysis::MasstransportAnalysis,femmodel::FemModel)# {{{ + + error("not supported yet") + + #Save output + RequestedOutputsx(femmodel, [ThicknessEnum, SurfaceEnum, BedEnum]) + +end #}}} +function CreateKMatrix(analysis::MasstransportAnalysis,element::Tria)# {{{ + + error("not supported") + + #Internmediaries + numnodes = 3 + + #Initialize Element matrix and basis function derivatives + Ke = ElementMatrix(element.nodes) + dbasis = Matrix{Float64}(undef,numnodes,2) + + #Retrieve all inputs and parameters + xyz_list = GetVerticesCoordinates(element.vertices) + H_input = GetInput(element, ThicknessEnum) + + #Prepare material object + material = Matice(element) + + #Start integrating + gauss = GaussTria(2) + for ig in 1:gauss.numgauss + + Jdet = JacobianDeterminant(xyz_list, gauss) + NodalFunctionsDerivatives(element,dbasis,xyz_list,gauss) + + H = GetInputValue(H_input, gauss, ig) + mu = ViscositySSA(material, xyz_list, gauss, ig) + + for i in 1:numnodes + for j in 1:numnodes + Ke.values[2*i-1,2*j-1] += gauss.weights[ig]*Jdet*mu*H*(4*dbasis[j,1]*dbasis[i,1] + dbasis[j,2]*dbasis[i,2]) + Ke.values[2*i-1,2*j ] += gauss.weights[ig]*Jdet*mu*H*(2*dbasis[j,2]*dbasis[i,1] + dbasis[j,1]*dbasis[i,2]) + Ke.values[2*i ,2*j-1] += gauss.weights[ig]*Jdet*mu*H*(2*dbasis[j,1]*dbasis[i,2] + dbasis[j,2]*dbasis[i,1]) + Ke.values[2*i ,2*j ] += gauss.weights[ig]*Jdet*mu*H*(4*dbasis[j,2]*dbasis[i,2] + dbasis[j,1]*dbasis[i,1]) + end + end + end + + #Add basal friction + phi=GetGroundedPortion(element, xyz_list) + + if(phi>0) + basis = Vector{Float64}(undef,numnodes) + friction = CoreFriction(element) + + #Start integrating + gauss = GaussTria(2) + for ig in 1:gauss.numgauss + + Jdet = JacobianDeterminant(xyz_list, gauss) + NodalFunctions(element, basis, gauss, ig, P1Enum) + + alpha2 = Alpha2(friction, gauss, ig) + + for i in 1:numnodes + for j in 1:numnodes + Ke.values[2*i-1,2*j-1] += gauss.weights[ig]*Jdet*phi*alpha2*basis[i]*basis[j] + Ke.values[2*i ,2*j ] += gauss.weights[ig]*Jdet*phi*alpha2*basis[i]*basis[j] + end + end + end + end + + return Ke +end #}}} +function CreatePVector(analysis::MasstransportAnalysis,element::Tria)# {{{ + + error("not supported") + + #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) + H_input = GetInput(element, ThicknessEnum) + s_input = GetInput(element, SurfaceEnum) + rho_ice = FindParam(element, MaterialsRhoIceEnum) + g = FindParam(element, ConstantsGEnum) + + #Start integrating + gauss = GaussTria(2) + for ig in 1:gauss.numgauss + + Jdet = JacobianDeterminant(xyz_list, gauss) + NodalFunctions(element, basis, gauss, ig, P1Enum) + + H = GetInputValue(H_input, gauss, ig) + ds = GetInputDerivativeValue(s_input, xyz_list, gauss, ig) + + for i in 1:numnodes + pe.values[2*i-1] += -gauss.weights[ig]*Jdet*rho_ice*g*H*ds[1]*basis[i] + pe.values[2*i ] += -gauss.weights[ig]*Jdet*rho_ice*g*H*ds[2]*basis[i] + end + end + + if(IsIcefront(element)) + + #Get additional parameters and inputs + b_input = GetInput(element, BaseEnum) + rho_water = FindParam(element, MaterialsRhoSeawaterEnum) + + #Get normal and ice front coordinates + xyz_list_front = Matrix{Float64}(undef,2,3) + GetIcefrontCoordinates!(element, xyz_list_front, xyz_list, MaskIceLevelsetEnum) + nx, ny = NormalSection(element, xyz_list_front) + + gauss = GaussTria(element, xyz_list, xyz_list_front, 3) + for ig in 1:gauss.numgauss + + Jdet = JacobianDeterminantSurface(xyz_list_front, gauss) + NodalFunctions(element, basis, gauss, ig, P1Enum) + + H = GetInputValue(H_input, gauss, ig) + b = GetInputValue(b_input, gauss, ig) + sl = 0 + + term = 0.5*g*rho_ice*H^2 + 0.5*g*rho_water*(min(0, H+b-sl)^2 - min(0, b-sl)^2) + + for i in 1:numnodes + pe.values[2*i-1] += gauss.weights[ig]*Jdet*term*nx*basis[i] + pe.values[2*i ] += gauss.weights[ig]*Jdet*term*ny*basis[i] + end + end + end + + return pe +end #}}} +function GetSolutionFromInputs(analysis::MasstransportAnalysis,ug::IssmVector,element::Tria) #{{{ + + #Get dofs for this finite element + doflist = GetDofList(element,GsetEnum) + @assert length(doflist)==3 + + #Fetch inputs + thickness_input = GetInput(element, ThicknessEnum) + + #Loop over each node and enter solution in ug + count = 0 + gauss=GaussTria(P1Enum) + for i in 1:gauss.numgauss + thickness = GetInputValue(thickness_input, gauss, i) + + count += 1 + ug.vector[doflist[count]] = thickness + end + + #Make sure we reached all the values + @assert count==length(doflist) + +end#}}} +function InputUpdateFromSolution(analysis::MasstransportAnalysis,ug::Vector{Float64},element::Tria) #{{{ + + #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 + + #Get some parameters + rho_water = FindParam(element, MaterialsRhoSeawaterEnum) + rho_ice = FindParam(element, MaterialsRhoIceEnum) + H_min = FindParam(element, MasstransportMinThicknessEnum) + + #Now split solution vector into x and y components + numnodes = 3 + thickness = Vector{Float64}(undef,numnodes) + for i in 1:numnodes + thickness[i]=values[i] + @assert isfinite(thickness[i]) + + #Enforce minimum thickness + if(thickness[i]0.) + #this is grounded ice: just add thickness to base. + newsurface[i] = bed[i]+newthickness[i]; #surface = bed + newthickness + newbase[i] = bed[i]; #new base at new bed + else + #this is an ice shelf: hydrostatic equilibrium + newsurface[i] = newthickness[i]*(1-rho_ice/rho_water) + sealevel[i] + newbase[i] = newthickness[i]*(-rho_ice/rho_water) + sealevel[i] + end + end + + AddInput(element, SurfaceEnum, newsurface, P1Enum) + AddInput(element, BaseEnum, newbase, P1Enum) +end#}}} Index: /issm/trunk-jpl/src/jl/core/analyses/stressbalanceanalysis.jl =================================================================== --- /issm/trunk-jpl/src/jl/core/analyses/stressbalanceanalysis.jl (revision 26729) +++ /issm/trunk-jpl/src/jl/core/analyses/stressbalanceanalysis.jl (revision 26730) @@ -33,5 +33,5 @@ #Provide node indices to element - for i in 1:md.mesh.numberofvertices + for i in 1:md.mesh.numberofelements Update(elements[i],inputs,i,md,P1Enum) end @@ -72,4 +72,7 @@ function Core(analysis::StressbalanceAnalysis,femmodel::FemModel)# {{{ + #Set current analysis to Stressnalance + SetCurrentConfiguration!(femmodel, analysis) + #Fetch parameters relevant to solution sequence maxiter = FindParam(femmodel.parameters,StressbalanceMaxiterEnum) @@ -78,5 +81,9 @@ abstol = FindParam(femmodel.parameters,StressbalanceAbstolEnum) + #Call solution sequence to compute new speeds solutionsequence_nonlinear(femmodel,analysis,maxiter,restol,reltol,abstol) + + #Save output + RequestedOutputsx(femmodel, [VxEnum,VyEnum,VelEnum]) end #}}} Index: /issm/trunk-jpl/src/jl/core/analyses/transientanalysis.jl =================================================================== --- /issm/trunk-jpl/src/jl/core/analyses/transientanalysis.jl (revision 26730) +++ /issm/trunk-jpl/src/jl/core/analyses/transientanalysis.jl (revision 26730) @@ -0,0 +1,38 @@ +#TransientAnalysis class definition +struct TransientAnalysis <: Analysis#{{{ +end #}}} + +function UpdateParameters(analysis::TransientAnalysis,parameters::Parameters,md::model) #{{{ + + AddParam(parameters, md.constants.yts, ConstantsYtsEnum) + AddParam(parameters, md.timestepping.start_time/md.constants.yts, TimeEnum) + AddParam(parameters, md.timestepping.final_time/md.constants.yts, TimesteppingFinalTimeEnum) + AddParam(parameters, md.timestepping.time_step/md.constants.yts, TimesteppingTimeStepEnum) + AddParam(parameters, md.transient.isstressbalance, TransientIsstressbalanceEnum) + AddParam(parameters, md.transient.ismasstransport, TransientIsmasstransportEnum) + +end#}}} +function Core(analysis::TransientAnalysis,femmodel::FemModel)# {{{ + + step = FindParam(femmodel.parameters, StepEnum) + time = FindParam(femmodel.parameters, TimeEnum) + finaltime = FindParam(femmodel.parameters, TimesteppingFinalTimeEnum) + yts = FindParam(femmodel.parameters, ConstantsYtsEnum) + dt = FindParam(femmodel.parameters, TimesteppingTimeStepEnum) + + isstressbalance = FindParam(femmodel.parameters, TransientIsstressbalanceEnum) + ismasstransport = FindParam(femmodel.parameters, TransientIsmasstransportEnum) + + while(time < finaltime - (yts*eps(Float64))) #make sure we run up to finaltime. + + step+=1 + time+=dt + AddParam(femmodel.parameters, time, TimeEnum); + AddParam(femmodel.parameters, step, StepEnum); + println("iteration ", step, "/", ceil((finaltime-time)/dt)+step," time [yr]: ", time/yts, " (time step: ", dt/yts, ")") + + if(isstressbalance) Core(StressbalanceAnalysis(), femmodel) end + if(ismasstransport) Core(MasstransportAnalysis(), femmodel) end + end + +end #}}} Index: /issm/trunk-jpl/src/jl/core/constraints.jl =================================================================== --- /issm/trunk-jpl/src/jl/core/constraints.jl (revision 26729) +++ /issm/trunk-jpl/src/jl/core/constraints.jl (revision 26730) @@ -14,5 +14,5 @@ #Apply Constraint - ApplyConstraint(node,constraint.dof,constraint.value) + ApplyConstraint(node, constraint.dof, constraint.value) end# }}} Index: /issm/trunk-jpl/src/jl/core/elements.jl =================================================================== --- /issm/trunk-jpl/src/jl/core/elements.jl (revision 26729) +++ /issm/trunk-jpl/src/jl/core/elements.jl (revision 26730) @@ -3,13 +3,28 @@ sid::Int64 pid::Int64 + #vertexids::Int64[3] #vertices::Vertex[3] vertexids::Vector{Int64} vertices::Vector{Vertex} - nodeids::Vector{Int64} + nodes::Vector{Node} + nodes_list::Vector{Vector{Node}} + nodes_ids_list::Vector{Vector{Int64}} + parameters::Parameters inputs::Inputs end# }}} +function Tria(sid::Int64, pid::Int64, vertexids::Vector{Int64}) #{{{ + + #This is the default constructor, at this point we don't have much information + tempparams = Parameters(Dict{IssmEnum,Parameter}()) + tempinputs = Inputs(-1,-1,Dict{IssmEnum,Input}()) + return Tria(sid, pid, + vertexids, Vector{Vertex}(undef,3), + Vector{Node}(undef,0), Vector{Vector{Node}}(undef,0), Vector{Vector{Int64}}(undef,0), + tempparams, tempinputs) + +end #}}} #Element functions @@ -31,17 +46,20 @@ end end #}}} -function Update(element::Tria,inputs::Inputs,index::Int64,md::model,finiteelement::IssmEnum) #{{{ +function Update(element::Tria, inputs::Inputs, index::Int64, md::model, finiteelement::IssmEnum) #{{{ if finiteelement==P1Enum numnodes = 3 - element.nodeids = Vector{Int64}(undef,numnodes) - element.nodeids[1] = md.mesh.elements[index,1] - element.nodeids[2] = md.mesh.elements[index,2] - element.nodeids[3] = md.mesh.elements[index,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 end #}}} -function Configure(element::Tria,nodes::Vector{Node},vertices::Vector{Vertex},parameters::Parameters,inputs::Inputs) # {{{ +function Configure(element::Tria,nodes::Vector{Node},vertices::Vector{Vertex},parameters::Parameters,inputs::Inputs,index::Int64) # {{{ #Configure vertices @@ -51,9 +69,12 @@ #Configure nodes (assuming P1 finite elements) - for i in 1:3 - element.nodes[i] = nodes[element.nodeids[i]] + nodes_list = element.nodes_list[index] + nodes_ids_list = element.nodes_ids_list[index] + for i in 1:3 + nodes_list[i] = nodes[nodes_ids_list[i]] end #Point to real datasets + element.nodes = element.nodes_list[index] element.parameters = parameters element.inputs = inputs Index: /issm/trunk-jpl/src/jl/core/femmodel.jl =================================================================== --- /issm/trunk-jpl/src/jl/core/femmodel.jl (revision 26729) +++ /issm/trunk-jpl/src/jl/core/femmodel.jl (revision 26730) @@ -1,14 +1,40 @@ #femmodel class definition mutable struct FemModel #{{{ + analyses::Vector{Analysis} elements::Vector{Tria} vertices::Vector{Vertex} + nodes::Vector{Node} + nodes_list::Vector{Vector{Node}} + parameters::Parameters inputs::Inputs + constraints::Vector{Constraint} + constraints_list::Vector{Vector{Constraint}} + #loads::Vector{Loads} + results::Vector{Result} end#}}} #femmodel functions +function SetCurrentConfiguration!(femmodel::FemModel, analysis::Analysis) #{{{ + #Find the index of this analysis + index = -1 + for i in 1:length(femmodel.analyses) + if(typeof(femmodel.analyses[i]) == typeof(analysis)) index = i end + end + if(index<1) error("Could not find analysis ",analysis, " in femmodel") end + + #Plug right nodes onto element + for i in 1:length(femmodel.elements) + femmodel.elements[i].nodes = femmodel.elements[i].nodes_list[index] + end + + #Plug in nodes and other datasets + femmodel.nodes = femmodel.nodes_list[index] + femmodel.constraints = femmodel.constraints_list[index] + +end#}}} Index: /issm/trunk-jpl/src/jl/core/modules.jl =================================================================== --- /issm/trunk-jpl/src/jl/core/modules.jl (revision 26729) +++ /issm/trunk-jpl/src/jl/core/modules.jl (revision 26730) @@ -5,37 +5,60 @@ end end#}}} -function ModelProcessor(md::model) #{{{ +function ModelProcessor(md::model, solutionstring::String) #{{{ #Initialize structures elements = Vector{Tria}(undef,0) vertices = Vector{Vertex}(undef,0) - nodes = Vector{Node}(undef,0) - constraints = Vector{Constraint}(undef,0) results = Vector{Result}(undef,0) parameters = Parameters(Dict{IssmEnum,Parameter}()) - inputs = Inputs(md.mesh.numberofelements,md.mesh.numberofvertices,Dict{IssmEnum,Input}()) + inputs = Inputs(md.mesh.numberofelements, md.mesh.numberofvertices, Dict{IssmEnum,Input}()) #Create elements, vertices and materials (independent of the analysis) - CreateElements(elements,md) - CreateVertices(vertices,md) - CreateParameters(parameters,md) - CreateInputs(inputs,elements,md) + CreateElements(elements, md) + CreateVertices(vertices, md) + CreateParameters(parameters, md) + CreateInputs(inputs,elements, md) + if solutionstring=="TransientSolution" + UpdateParameters(TransientAnalysis(), parameters, md) + end #Now create analysis specific data structure - analysis = StressbalanceAnalysis() - UpdateParameters(analysis,parameters,md) - CreateNodes(analysis,nodes,md) - UpdateElements(analysis,elements,inputs,md) - CreateConstraints(analysis,constraints,md) - - #Configure objects - ConfigureObjectx(elements,nodes,vertices,parameters,inputs) - - #Constrain and Number nodes - SpcNodesx(nodes,constraints,parameters) - NodesDofx(nodes,parameters) + if solutionstring=="StressbalanceSolution" + analyses = [StressbalanceAnalysis()] + elseif solutionstring=="TransientSolution" + analyses = [StressbalanceAnalysis(), MasstransportAnalysis()] + else + error(solutionstring, " not supported by ModelProcessor") + end + + #Initialize analysis specific datasets + numanalyses = length(analyses) + nodes = Vector{Vector{Node}}(undef,numanalyses) + constraints = Vector{Vector{Constraint}}(undef,numanalyses) + for i in 1:numanalyses + analysis = analyses[i] + println(" creating datasets for analysis ", typeof(analysis)) + nodes[i] = Vector{Node}(undef,0) + constraints[i] = Vector{Constraint}(undef,0) + + UpdateParameters(analysis, parameters, md) + UpdateElements(analysis, elements, inputs, md) + CreateNodes(analysis, nodes[i], md) + CreateConstraints(analysis, constraints[i], md) + + #Configure objects + ConfigureObjectx(elements, nodes[i], vertices, parameters, inputs, i) + + #Constrain and Number nodes + SpcNodesx(nodes[i], constraints[i], parameters) + NodesDofx(nodes[i], parameters) + end #Build FemModel - return FemModel(elements, vertices, nodes, parameters, inputs, constraints, results) + return FemModel(analyses, elements, vertices, + Vector{Node}(undef,0), nodes, + parameters, inputs, + Vector{Constraint}(undef,0), constraints, + results) end# }}} function CreateElements(elements::Vector{Tria},md::model) #{{{ @@ -44,6 +67,5 @@ @assert length(elements)==0 - tempparams = Parameters(Dict{IssmEnum,Parameter}()) - tempinputs = Inputs(-1,-1,Dict{IssmEnum,Input}()) + count = 0 @@ -55,5 +77,5 @@ #Call constructor and add to dataset elements - push!(elements,Tria(i,count,vertexids,Vector{Vertex}(undef,3),nodeids,Vector{Node}(undef,3),tempparams,tempinputs)) + push!(elements,Tria(i,count, vertexids)) end @@ -85,4 +107,5 @@ AddParam(parameters,1,StepEnum) AddParam(parameters,0.0,TimeEnum) + end# }}} function CreateInputs(inputs::Inputs,elements::Vector{Tria},md::model) #{{{ @@ -104,8 +127,8 @@ #Other modules -function ConfigureObjectx(elements::Vector{Tria},nodes::Vector{Node},vertices::Vector{Vertex},parameters::Parameters,inputs::Inputs) #{{{ +function ConfigureObjectx(elements::Vector{Tria}, nodes::Vector{Node}, vertices::Vector{Vertex}, parameters::Parameters, inputs::Inputs, analysis::Int64) #{{{ for i in 1:length(elements) - Configure(elements[i],nodes,vertices,parameters,inputs) + Configure(elements[i], nodes, vertices, parameters, inputs, analysis) end Index: /issm/trunk-jpl/src/jl/core/parameters.jl =================================================================== --- /issm/trunk-jpl/src/jl/core/parameters.jl (revision 26729) +++ /issm/trunk-jpl/src/jl/core/parameters.jl (revision 26730) @@ -8,4 +8,8 @@ enum::IssmEnum value::Int64 +end# }}} +struct BoolParam <: Parameter #{{{ + enum::IssmEnum + value::Bool end# }}} @@ -22,4 +26,7 @@ return param.value end#}}} +function GetParameterValue(param::BoolParam) #{{{ + return param.value +end#}}} #Parameters functions @@ -29,12 +36,12 @@ end#}}} -function AddParam(parameters::Parameters,value::Int64,enum::IssmEnum) #{{{ +function AddParam(parameters::Parameters,value::Int64, enum::IssmEnum) #{{{ parameters.lookup[enum] = IntParam(enum,value) end#}}} -function AddParam(parameters::Parameters,value::Int64,enum::IssmEnum) #{{{ +function AddParam(parameters::Parameters,value::Bool, enum::IssmEnum) #{{{ - parameters.lookup[enum] = IntParam(enum,value) + parameters.lookup[enum] = BoolParam(enum,value) end#}}} Index: /issm/trunk-jpl/src/jl/core/solve.jl =================================================================== --- /issm/trunk-jpl/src/jl/core/solve.jl (revision 26729) +++ /issm/trunk-jpl/src/jl/core/solve.jl (revision 26730) @@ -12,7 +12,12 @@ include("./matice.jl") include("./friction.jl") +include("./analyses/analysis.jl") include("./femmodel.jl") -include("./analyses/analysis.jl") + +#All analyses include("./analyses/stressbalanceanalysis.jl") +include("./analyses/masstransportanalysis.jl") +include("./analyses/transientanalysis.jl") + include("./solutionsequences.jl") include("./modules.jl") @@ -20,16 +25,18 @@ include("./utils.jl") -function IssmCore(md::model) #{{{ - +function IssmCore(md::model, solutionstring::String) #{{{ #Construct FemModel - femmodel=ModelProcessor(md) + femmodel=ModelProcessor(md, solutionstring) - #Solve - analysis = StressbalanceAnalysis() + #Solve FIXME: to be improved later... + if(solutionstring=="StressbalanceSolution") + analysis = StressbalanceAnalysis() + elseif (solutionstring=="TransientSolution") + analysis = TransientAnalysis() + else + error("not supported") + end Core(analysis, femmodel) - - #Save output - RequestedOutputsx(femmodel, [VxEnum,VyEnum,VelEnum]) #move results to md Index: /issm/trunk-jpl/src/jl/tests/test.jl =================================================================== --- /issm/trunk-jpl/src/jl/tests/test.jl (revision 26729) +++ /issm/trunk-jpl/src/jl/tests/test.jl (revision 26730) @@ -5,5 +5,6 @@ #Load model from MATLAB file -file = matopen("./Pig.mat") +#file = matopen("./Pig.mat") +file = matopen("/Users/mmorligh/Desktop/issmuci/trunk-jpl/test/NightlyRun/temp.mat") mat = read(file, "md") close(file) Index: /issm/trunk-jpl/src/jl/tests/test101.jl =================================================================== --- /issm/trunk-jpl/src/jl/tests/test101.jl (revision 26729) +++ /issm/trunk-jpl/src/jl/tests/test101.jl (revision 26730) @@ -25,6 +25,6 @@ vy = archread(issmdir()*"/test/Data/SquareShelfConstrained.arch","vy") index = Int.(archread(issmdir()*"/test/Data/SquareShelfConstrained.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.initialization.vx=0 .*InterpFromMeshToMesh2d(index,x,y,vx,md.mesh.x,md.mesh.y,0.0) +md.initialization.vy=0 .*InterpFromMeshToMesh2d(index,x,y,vy,md.mesh.x,md.mesh.y,0.0) md.materials.rheology_B=1.815730284801701e+08*ones(md.mesh.numberofvertices) Index: /issm/trunk-jpl/src/jl/tests/test208.jl =================================================================== --- /issm/trunk-jpl/src/jl/tests/test208.jl (revision 26730) +++ /issm/trunk-jpl/src/jl/tests/test208.jl (revision 26730) @@ -0,0 +1,59 @@ +#!/Applications/Julia-1.6.app/Contents/Resources/julia/bin/julia +include("../issm.jl") +using .ISSM + +md = model() +md = triangle(md,issmdir()*"/test/Exp/Square.exp", 150000.) +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 .-10 + +#Initial velocity +x = archread(issmdir()*"/test/Data/SquareShelfConstrained.arch","x") +y = archread(issmdir()*"/test/Data/SquareShelfConstrained.arch","y") +vx = archread(issmdir()*"/test/Data/SquareShelfConstrained.arch","vx") +vy = archread(issmdir()*"/test/Data/SquareShelfConstrained.arch","vy") +index = archread(issmdir()*"/test/Data/SquareShelfConstrained.arch","index") +md.initialization.vx=zeros(md.mesh.numberofvertices)#InterpFromMeshToMesh2d(index,x,y,vx,md.mesh.x,md.mesh.y) +md.initialization.vy=zeros(md.mesh.numberofvertices)#InterpFromMeshToMesh2d(index,x,y,vy,md.mesh.x,md.mesh.y) + +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.stressbalance.restol=0.1 +md.stressbalance.reltol=0.02 +md.stressbalance.abstol=NaN +md.timestepping.start_time = 0.0 +md.timestepping.final_time = 3.0 +md.timestepping.time_step = 1.0 + +#Boundary conditions +nodefront=ContourToNodes(md.mesh.x,md.mesh.y,issmdir()*"/test/Exp/SquareFront.exp",2.0) .& md.mesh.vertexonboundary +md.stressbalance.spcvx = NaN*ones(md.mesh.numberofvertices) +md.stressbalance.spcvy = NaN*ones(md.mesh.numberofvertices) +pos = findall(md.mesh.vertexonboundary .& .~nodefront) +md.mask.ice_levelset[findall(nodefront)] .= 0 + +segmentsfront=md.mask.ice_levelset[md.mesh.segments[:,1:2]]==0 +segments = findall(vec(sum(Int64.(md.mask.ice_levelset[md.mesh.segments[:,1:2]].==0), dims=2)) .!=2) +pos=md.mesh.segments[segments,1:2] +md.stressbalance.spcvx[pos] .= 0.0 +md.stressbalance.spcvy[pos] .= 0.0 + +md.smb.mass_balance = zeros(md.mesh.numberofvertices) +md.basalforcings.groundedice_melting_rate = zeros(md.mesh.numberofvertices) +md.basalforcings.floatingice_melting_rate = ones(md.mesh.numberofvertices) +md.masstransport.spcthickness = NaN*ones(md.mesh.numberofvertices) + +md=solve(md,"Transient") Index: /issm/trunk-jpl/src/jl/usr/classes.jl =================================================================== --- /issm/trunk-jpl/src/jl/usr/classes.jl (revision 26729) +++ /issm/trunk-jpl/src/jl/usr/classes.jl (revision 26730) @@ -142,4 +142,69 @@ end # }}} +#Basalforcings {{{ +mutable struct Basalforcings + groundedice_melting_rate::Vector{Float64} + floatingice_melting_rate::Vector{Float64} +end +function Basalforcings() #{{{ + return Basalforcings( Vector{Float64}(undef,0), Vector{Float64}(undef,0)) +end# }}} +function Base.show(io::IO, this::Basalforcings)# {{{ + IssmStructDisp(io, this) +end# }}} +# }}} +#Surfaceforcings {{{ +mutable struct SMBforcings + mass_balance::Vector{Float64} +end +function SMBforcings() #{{{ + return SMBforcings( Vector{Float64}(undef,0)) +end# }}} +function Base.show(io::IO, this::SMBforcings)# {{{ + IssmStructDisp(io, this) +end# }}} +# }}} +#Timestepping{{{ +abstract type AbstractTimestepping end +mutable struct Timestepping <: AbstractTimestepping + start_time::Float64 + final_time::Float64 + time_step::Float64 +end +function Timestepping() #{{{ + return Timestepping( 0., 0., 0.) +end# }}} +function Base.show(io::IO, this::Timestepping)# {{{ + IssmStructDisp(io, this) +end# }}} +# }}} +#Masstransport {{{ +mutable struct Masstransport + spcthickness::Vector{Float64} + min_thickness::Float64 + stabilization::Int64 +end +function Masstransport() #{{{ + return Masstransport( Vector{Float64}(undef,0), 10.0, 1) +end# }}} +function Base.show(io::IO, this::Masstransport)# {{{ + IssmStructDisp(io, this) +end# }}} +# }}} +#Transient {{{ +mutable struct Transient + issmb::Bool + ismasstransport::Bool + isstressbalance::Bool + isgroundingline::Bool + ismovingfront::Bool +end +function Transient() #{{{ + return Transient( true, true, true, true, true) +end# }}} +function Base.show(io::IO, this::Transient)# {{{ + IssmStructDisp(io, this) +end# }}} +# }}} #Model structure @@ -154,7 +219,15 @@ results::Dict friction::AbstractFriction + basalforcings::Basalforcings + smb::SMBforcings + timestepping::Timestepping + masstransport::Masstransport + transient::Transient end function model() #{{{ - return model( Mesh2dTriangle(), Geometry(), Mask(), Materials(), Initialization(),Stressbalance(), Constants(), Dict(), BuddFriction()) + return model( Mesh2dTriangle(), Geometry(), Mask(), Materials(), + Initialization(),Stressbalance(), Constants(), Dict(), + BuddFriction(), Basalforcings(), SMBforcings(), Timestepping(), + Masstransport(), Transient()) end#}}} function model(matmd::Dict) #{{{ @@ -185,4 +258,7 @@ setfield!(mdfield, Symbol(name2), Int64(value_matlab)) + elseif typeof(value_matlab)==Float64 && typeof(value_julia)==Bool + setfield!(mdfield, Symbol(name2), Bool(value_matlab)) + elseif typeof(value_matlab)==Matrix{Float64} && typeof(value_julia)==Vector{Float64} if(size(value_matlab,2)!=1) error("only one column expected") end @@ -218,9 +294,12 @@ @printf "%19s: %-26s -- %s\n" "materials" typeof(md.materials) "material properties" @printf "%19s: %-26s -- %s\n" "initialization" typeof(md.initialization) "initial state" - @printf "%19s: %-26s -- %s\n" "stressbalance" typeof(md.stressbalance) "stress balance parameters" @printf "%19s: %-26s -- %s\n" "constants" typeof(md.constants) "physical constants" @printf "%19s: %-26s -- %s\n" "friction" typeof(md.friction) "basal friction" - -end# }}} - - + @printf "%19s: %-26s -- %s\n" "basalforcings" typeof(md.basalforcings) "basal forcings" + @printf "%19s: %-26s -- %s\n" "smb" typeof(md.smb) "surface mass balance" + @printf "%19s: %-26s -- %s\n" "timestepping" typeof(md.timestepping) "time stepping for transient simulations" + @printf "%19s: %-26s -- %s\n" "stressbalance" typeof(md.stressbalance) "parameters stress balance simulations" + @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" + +end# }}} Index: /issm/trunk-jpl/src/jl/usr/utils.jl =================================================================== --- /issm/trunk-jpl/src/jl/usr/utils.jl (revision 26729) +++ /issm/trunk-jpl/src/jl/usr/utils.jl (revision 26730) @@ -207,9 +207,11 @@ if solution=="sb" || solution=="Stressbalance" solutionstring = "StressbalanceSolution" + elseif solution=="tr" || solution=="Transient" + solutionstring = "TransientSolution" else error("solutionstring "*solution*" not supported!"); end - IssmCore(md) + IssmCore(md, solutionstring) return md