source: issm/trunk-jpl/test/NightlyRun/test1106.py@ 21060

#Test Name: ISMIPCFS
import numpy
from model import *
from triangle import *
from setmask import *
from parameterize import *
from setflowequation import *
from solve import *
from MatlabFuncs import *
from PythonFuncs import *

"""
This test is a test from the ISMP-HOM Intercomparison project.
Pattyn and Payne 2006
"""

L_list=[5000.,10000.,20000.,40000.,80000.,160000.]
results=[]

for L in L_list:
        md=triangle(model(),"../Exp/Square_%d.exp" % L,L/10.)    #size 3*L 
        md=setmask(md,'','')    #ice sheet test
        md=parameterize(md,'../Par/ISMIPC.py')
        md.friction.coefficient=numpy.sqrt(md.constants.yts*(1000.+1000.*numpy.sin(md.mesh.x.reshape(-1,1)*2.*numpy.pi/L)*numpy.sin(md.mesh.y.reshape(-1,1)*2.*numpy.pi/L)))
        md.extrude(10,1.)

        #Add spc on the borders
        pos=numpy.nonzero(logical_or_n(md.mesh.x==0.,md.mesh.x==numpy.max(md.mesh.x),md.mesh.y==0.,md.mesh.y==numpy.max(md.mesh.y)))
        md.stressbalance.spcvx[pos]=0.
        md.stressbalance.spcvy[pos]=0.
        if   (L==5000.):
                md.stressbalance.spcvx[pos]=15.66
                md.stressbalance.spcvy[pos]=-0.1967
        elif (L==10000.):
                md.stressbalance.spcvx[pos]=16.04
                md.stressbalance.spcvy[pos]=-0.1977
        elif (L==20000.):
                md.stressbalance.spcvx[pos]=16.53
                md.stressbalance.spcvy[pos]=-1.27
        elif (L==40000.):
                md.stressbalance.spcvx[pos]=17.23
                md.stressbalance.spcvy[pos]=-3.17
        elif (L==80000.):
                md.stressbalance.spcvx[pos]=16.68
                md.stressbalance.spcvy[pos]=-2.69
        elif (L==160000.):
                md.stressbalance.spcvx[pos]=16.03
                md.stressbalance.spcvy[pos]=-1.27

        md=setflowequation(md,'FS','all')

        #Compute the stressbalance
        md.cluster=generic('name',oshostname(),'np',8)
        md=solve(md,'Stressbalance')

        #Plot the results and save them
        vx=md.results.StressbalanceSolution.Vx
        vy=md.results.StressbalanceSolution.Vy
        vz=md.results.StressbalanceSolution.Vz
        results.append(md.results.StressbalanceSolution)

#       plotmodel(md,'data',vx,'data',vy,'data',vz,'layer#all',md.mesh.numberoflayers)

#Fields and tolerances to track changes
field_names     =[\
        'Vx5km','Vy5km','Vz5km',\
        'Vx10km','Vy10km','Vz10km',\
        'Vx20km','Vy20km','Vz20km',\
        'Vx40km','Vy40km','Vz40km',\
        'Vx80km','Vy80km','Vz80km',\
        'Vx160km','Vy160km','Vz160km'
]
field_tolerances=[\
        1e-12,1e-12,1e-11,\
        1e-12,1e-12,1e-12,\
        1e-12,1e-12,1e-12,\
        1e-12,1e-12,1e-12,\
        1e-12,1e-12,1e-12,\
        1e-12,1e-11,1e-12,\
]
field_values=[]
for result in results:
        field_values=field_values+[\
                result.Vx,\
                result.Vy,\
                result.Vz,\
                ]