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 26731) @@ -38,4 +38,5 @@ FetchDataToInput(md,inputs,elements,md.basalforcings.groundedice_melting_rate./md.constants.yts,BasalforcingsGroundediceMeltingRateEnum) FetchDataToInput(md,inputs,elements,md.basalforcings.floatingice_melting_rate./md.constants.yts,BasalforcingsFloatingiceMeltingRateEnum) + FetchDataToInput(md,inputs,elements,md.smb.mass_balance./md.constants.yts,SmbMassBalanceEnum) FetchDataToInput(md,inputs,elements,md.mask.ice_levelset, MaskIceLevelsetEnum) FetchDataToInput(md,inputs,elements,md.mask.ocean_levelset, MaskOceanLevelsetEnum) @@ -54,13 +55,18 @@ function Core(analysis::MasstransportAnalysis,femmodel::FemModel)# {{{ - error("not supported yet") + println(" computing mass transport") + SetCurrentConfiguration!(femmodel, analysis) + + InputDuplicatex(femmodel, ThicknessEnum, ThicknessOldEnum) + InputDuplicatex(femmodel, BaseEnum, BaseOldEnum) + InputDuplicatex(femmodel, SurfaceEnum, SurfaceOldEnum) + + solutionsequence_linear(femmodel,analysis) #Save output - RequestedOutputsx(femmodel, [ThicknessEnum, SurfaceEnum, BedEnum]) + RequestedOutputsx(femmodel, [ThicknessEnum, SurfaceEnum, BaseEnum]) end #}}} function CreateKMatrix(analysis::MasstransportAnalysis,element::Tria)# {{{ - - error("not supported") #Internmediaries @@ -70,12 +76,15 @@ 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) - H_input = GetInput(element, ThicknessEnum) - - #Prepare material object - material = Matice(element) - + vx_input = GetInput(element, VxEnum) + vy_input = GetInput(element, VyEnum) + dt = FindParam(element, TimesteppingTimeStepEnum) + stabilization = FindParam(element, MasstransportStabilizationEnum) + + h = CharacteristicLength(element) + #Start integrating gauss = GaussTria(2) @@ -83,41 +92,41 @@ Jdet = JacobianDeterminant(xyz_list, gauss) - NodalFunctionsDerivatives(element,dbasis,xyz_list,gauss) - - H = GetInputValue(H_input, gauss, ig) - mu = ViscositySSA(material, xyz_list, gauss, ig) - + 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[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]) + Ke.values[i ,j] += gauss.weights[ig]*Jdet*basis[i]*basis[j] 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 + + #Advection term + vx = GetInputValue(vx_input, gauss, ig) + vy = GetInputValue(vy_input, gauss, ig) + dvx = GetInputDerivativeValue(vx_input, xyz_list, gauss, ig) + dvy = GetInputDerivativeValue(vy_input, xyz_list, gauss, ig) + for i in 1:numnodes + for j in 1:numnodes + #\phi_i \phi_j \nabla\cdot v + Ke.values[i ,j] += dt*gauss.weights[ig]*Jdet*basis[i]*basis[j]*(dvx[1] + dvy[2]) + #\phi_i v\cdot\nabla\phi_j + Ke.values[i ,j] += dt*gauss.weights[ig]*Jdet*basis[i]*(vx*dbasis[j,1] + vy*dbasis[j,2]) + end + end + + #Stabilization + if(stabilization==0) + #do nothing + elseif (stabilization==1) + vx = GetInputAverageValue(vx_input) + vy = GetInputAverageValue(vy_input) + D = dt*gauss.weights[ig]*Jdet*[h/2*abs(vx) 0; 0 h/2*abs(vy)] + for i in 1:numnodes; for j in 1:numnodes + Ke.values[i ,j] += (dbasis[i,1]*(D[1,1]*dbasis[j,1] + D[1,2]*dbasis[j,2]) + + dbasis[i,2]*(D[2,1]*dbasis[j,1] + D[2,2]*dbasis[j,2])) + end end + else + error("Stabilization ",stabilization, " not supported yet") end end @@ -126,6 +135,4 @@ end #}}} function CreatePVector(analysis::MasstransportAnalysis,element::Tria)# {{{ - - error("not supported") #Internmediaries @@ -138,11 +145,17 @@ #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) + H_input = GetInput(element, ThicknessEnum) + gmb_input = GetInput(element, BasalforcingsGroundediceMeltingRateEnum) + fmb_input = GetInput(element, BasalforcingsFloatingiceMeltingRateEnum) + smb_input = GetInput(element, SmbMassBalanceEnum) + olevelset_input = GetInput(element, MaskOceanLevelsetEnum) + dt = FindParam(element, TimesteppingTimeStepEnum) + stabilization = FindParam(element, MasstransportStabilizationEnum) + + #How much is actually grounded? + phi=GetGroundedPortion(element, xyz_list) #Start integrating - gauss = GaussTria(2) + gauss = GaussTria(3) for ig in 1:gauss.numgauss @@ -150,40 +163,16 @@ NodalFunctions(element, basis, gauss, ig, P1Enum) - H = GetInputValue(H_input, gauss, ig) - ds = GetInputDerivativeValue(s_input, xyz_list, gauss, ig) + smb = GetInputValue(smb_input, gauss, ig) + H = GetInputValue(H_input, gauss, ig) + + #Only apply melt on fully floating cells + if(phi<0.00000001) + mb = GetInputValue(fmb_input, gauss, ig) + else + mb = GetInputValue(gmb_input, gauss, ig) + end 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 + pe.values[i] += gauss.weights[ig]*Jdet*(H + dt*(smb - mb))*basis[i] end end @@ -237,6 +226,6 @@ #Enforce minimum thickness - if(thickness[i]maxstep) maxstep = femmodel.results[i].step end + end + + #Initialize vector now that we know the size + output = Vector{Dict}(undef, maxstep) + for i in 1:maxstep; output[i] = Dict() end + + #Insert results in vector + for i in 1:length(femmodel.results) + result = femmodel.results[i] + step = femmodel.results[i].step + (output[step])[EnumToString(result.enum)] = result.value + end + else + output = Dict() + for i in length(femmodel.results) + result = femmodel.results[i] + output[EnumToString(result.enum)] = result.value + end + end + + md.results[solution] = output end# }}} @@ -182,4 +204,7 @@ return ug +end#}}} +function InputDuplicatex(femmodel::FemModel, oldenum::IssmEnum, newenum::IssmEnum) #{{{ + DuplicateInput(femmodel.inputs, oldenum, newenum) end#}}} function Reducevectorgtofx(ug::IssmVector,nodes::Vector{Node}) #{{{ @@ -281,5 +306,5 @@ #Add to femmodel.results dataset - push!(femmodel.results,result) + push!(femmodel.results, result) else Index: /issm/trunk-jpl/src/jl/core/solutionsequences.jl =================================================================== --- /issm/trunk-jpl/src/jl/core/solutionsequences.jl (revision 26730) +++ /issm/trunk-jpl/src/jl/core/solutionsequences.jl (revision 26731) @@ -1,2 +1,22 @@ +function solutionsequence_linear(femmodel::FemModel,analysis::Analysis) # {{{ + + #Get new matrices + Kff, Kfs, pf = SystemMatricesx(femmodel,analysis) + + #Enforce constraints + ys = CreateNodalConstraintsx(femmodel.nodes) + Reduceloadx!(pf, Kfs, ys) + + #Solve! + uf = Solverx(Kff, pf) + + #Merge uf with ys and update inputs + gsize = NumberOfDofs(femmodel.nodes,GsetEnum) + ug = IssmVector(gsize) + Mergesolutionfromftogx(ug, uf, ys, femmodel.nodes) + + InputUpdateFromSolutionx(analysis, ug, femmodel) + +end# }}} function solutionsequence_nonlinear(femmodel::FemModel,analysis::Analysis,maxiter::Int64,restol::Float64,reltol::Float64,abstol::Float64) # {{{ Index: /issm/trunk-jpl/src/jl/core/solve.jl =================================================================== --- /issm/trunk-jpl/src/jl/core/solve.jl (revision 26730) +++ /issm/trunk-jpl/src/jl/core/solve.jl (revision 26731) @@ -41,5 +41,5 @@ #move results to md - OutputResultsx(femmodel, md) + OutputResultsx(femmodel, md, solutionstring) end# }}} Index: /issm/trunk-jpl/src/jl/core/toolkits.jl =================================================================== --- /issm/trunk-jpl/src/jl/core/toolkits.jl (revision 26730) +++ /issm/trunk-jpl/src/jl/core/toolkits.jl (revision 26731) @@ -125,8 +125,15 @@ end#}}} -function Solverx(A::IssmMatrix,b::IssmVector,xold::IssmVector) #{{{ +function Solverx(A::IssmMatrix, b::IssmVector, xold::IssmVector) #{{{ #Initialize output #x = IssmVector(GetSize(xold)) + + return Solverx(A, b) + +end#}}} +function Solverx(A::IssmMatrix, b::IssmVector) #{{{ + + #Initialize output x = IssmVector(0) @@ -136,4 +143,3 @@ return x - end#}}}