source: issm/oecreview/Archive/25834-26739/ISSM-26712-26713.diff

../trunk-jpl/src/jl/md/parameterization.jl

@@ -16 +16 @@
 """
 function setmask(md::model,floatingicename::String,groundedicename::String)
 
-        elementonfloatingice = FlagElements( md, floatingicename);
-        elementongroundedice = FlagElements( md, groundedicename);
+        elementonfloatingice = FlagElements( md, floatingicename)
+        elementongroundedice = FlagElements( md, groundedicename)
 
         elementonfloatingice = convert( Array{Float64}, (elementonfloatingice.>0) .&  (elementongroundedice.==0.))
         elementongroundedice = convert( Array{Float64}, elementonfloatingice.==0.)
@@ -44 +44 @@
         elseif region == "all"
                 flags = ones(md.mesh.numberofelements)
         else
-                flags = InContour(md.mesh.x,md.mesh.y,region);
+                xcenter = md.mesh.x[md.mesh.elements]*[1;1;1]/3
+                ycenter = md.mesh.y[md.mesh.elements]*[1;1;1]/3
+                flags = ContourToNodes(xcenter, ycenter, region, 2.)
         end
 
         return flags

../trunk-jpl/src/jl/test101.jl

@@ -25 +25 @@
 y     = ISSM.archread(ISSM.issmdir()*"/test/Data/SquareShelfConstrained.arch","y")
 vx    = ISSM.archread(ISSM.issmdir()*"/test/Data/SquareShelfConstrained.arch","vx")
 vy    = ISSM.archread(ISSM.issmdir()*"/test/Data/SquareShelfConstrained.arch","vy")
-index = ISSM.archread(ISSM.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)
+index = Int.(ISSM.archread(ISSM.issmdir()*"/test/Data/SquareShelfConstrained.arch","index"))
+md.initialization.vx=ISSM.InterpFromMeshToMesh2d(index,x,y,vx,md.mesh.x,md.mesh.y,0.0)
+md.initialization.vy=ISSM.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.stressbalance.restol=0.05
 md.stressbalance.reltol=0.05

../trunk-jpl/src/jl/test501.jl

@@ -5 +5 @@
 md = ISSM.model()
 md = ISSM.triangle(md,"../../test/Exp/Pig.exp",10000.)
 md = ISSM.setmask(md,"../../test/Exp/PigShelves.exp","../../test/Exp/PigIslands.exp")
-error("not implemented yet")
 
-#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 .+20
-md.geometry.surface   = md.geometry.base+md.geometry.thickness
-md.geometry.bed       = md.geometry.base
+#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])
 
-#Initial velocity
-#x     = ISSM.archread(ISSM.issmdir()*"/test/Data/SquareShelfConstrained.arch","x")
-#y     = ISSM.archread(ISSM.issmdir()*"/test/Data/SquareShelfConstrained.arch","y")
-#vx    = ISSM.archread(ISSM.issmdir()*"/test/Data/SquareShelfConstrained.arch","vx")
-#vy    = ISSM.archread(ISSM.issmdir()*"/test/Data/SquareShelfConstrained.arch","vy")
-#index = ISSM.archread(ISSM.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.friction.coefficient=50*ones(md.mesh.numberofvertices)
 
 md.stressbalance.restol=0.05
-md.stressbalance.reltol=0.05
+md.stressbalance.reltol=1.0
 md.stressbalance.abstol=NaN
 
 #Boundary conditions