source: issm/trunk-jpl/src/jl/test501.jl@ 26714

#!/Applications/Julia-1.6.app/Contents/Resources/julia/bin/julia
include("issm.jl")
using .ISSM

md = ISSM.model()
md = ISSM.triangle(md,"../../test/Exp/Pig.exp",10000.)
md = ISSM.setmask(md,"../../test/Exp/PigShelves.exp","../../test/Exp/PigIslands.exp")

#Initial velocity and geometry
x         = ISSM.archread(ISSM.issmdir()*"/test/Data/Pig.arch","x")
y         = ISSM.archread(ISSM.issmdir()*"/test/Data/Pig.arch","y")
vx_obs    = ISSM.archread(ISSM.issmdir()*"/test/Data/Pig.arch","vx_obs")
vy_obs    = ISSM.archread(ISSM.issmdir()*"/test/Data/Pig.arch","vy_obs")
index     = Int.(ISSM.archread(ISSM.issmdir()*"/test/Data/Pig.arch","index"))
surface   = ISSM.archread(ISSM.issmdir()*"/test/Data/Pig.arch","surface")
thickness = ISSM.archread(ISSM.issmdir()*"/test/Data/Pig.arch","thickness")
bed       = ISSM.archread(ISSM.issmdir()*"/test/Data/Pig.arch","bed")
md.initialization.vx=ISSM.InterpFromMeshToMesh2d(index, x, y, vx_obs, md.mesh.x, md.mesh.y, 0.0)
md.initialization.vy=ISSM.InterpFromMeshToMesh2d(index, x, y, vy_obs, md.mesh.x, md.mesh.y, 0.0)
md.geometry.surface = ISSM.InterpFromMeshToMesh2d(index, x, y, surface, md.mesh.x, md.mesh.y, 0.0)
md.geometry.thickness = ISSM.InterpFromMeshToMesh2d(index, x, y, thickness, md.mesh.x, md.mesh.y, 0.0)
md.geometry.base=md.geometry.surface .- md.geometry.thickness
md.geometry.bed =md.geometry.base
pos = findall(md.mask.ocean_levelset.<0)
md.geometry.bed[pos] = ISSM.InterpFromMeshToMesh2d(index, x, y, bed, md.mesh.x[pos], md.mesh.y[pos])

md.materials.rheology_B=1.815730284801701e+08*ones(md.mesh.numberofvertices)
md.materials.rheology_n=3*ones(md.mesh.numberofelements)
md.friction.coefficient=50*ones(md.mesh.numberofvertices)

md.stressbalance.restol=0.05
md.stressbalance.reltol=1.0
md.stressbalance.abstol=NaN

#Boundary conditions
pos = findall(vec(sum(Int64.(md.mask.ocean_levelset[md.mesh.elements].<0), dims=2)) .> 0.0)
vertexonfloatingice=zeros(md.mesh.numberofvertices)
vertexonfloatingice[md.mesh.elements[pos,:]] .= 1
nodefront=(md.mesh.vertexonboundary .& (vertexonfloatingice.>0))
md.mask.ice_levelset[findall(nodefront)] .= 0

md.stressbalance.spcvx = NaN*ones(md.mesh.numberofvertices)
md.stressbalance.spcvy = NaN*ones(md.mesh.numberofvertices)
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=ISSM.solve(md,"Stressbalance")