1 | #Test Name: ISMIPCFS
|
---|
2 | import numpy
|
---|
3 | from model import *
|
---|
4 | from triangle import *
|
---|
5 | from setmask import *
|
---|
6 | from parameterize import *
|
---|
7 | from setflowequation import *
|
---|
8 | from EnumDefinitions import *
|
---|
9 | from solve import *
|
---|
10 | from MatlabFuncs import *
|
---|
11 | from PythonFuncs import *
|
---|
12 |
|
---|
13 | """
|
---|
14 | This test is a test from the ISMP-HOM Intercomparison project.
|
---|
15 | Pattyn and Payne 2006
|
---|
16 | """
|
---|
17 |
|
---|
18 | L_list=[5000.,10000.,20000.,40000.,80000.,160000.]
|
---|
19 | results=[]
|
---|
20 |
|
---|
21 | for L in L_list:
|
---|
22 | md=triangle(model(),"../Exp/Square_%d.exp" % L,L/10.) #size 3*L
|
---|
23 | md=setmask(md,'','') #ice sheet test
|
---|
24 | md=parameterize(md,'../Par/ISMIPC.py')
|
---|
25 | 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)))
|
---|
26 | md.extrude(10,1.)
|
---|
27 |
|
---|
28 | #Add spc on the borders
|
---|
29 | 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)))
|
---|
30 | md.stressbalance.spcvx[pos]=0.
|
---|
31 | md.stressbalance.spcvy[pos]=0.
|
---|
32 | if (L==5000.):
|
---|
33 | md.stressbalance.spcvx[pos]=15.66
|
---|
34 | md.stressbalance.spcvy[pos]=-0.1967
|
---|
35 | elif (L==10000.):
|
---|
36 | md.stressbalance.spcvx[pos]=16.04
|
---|
37 | md.stressbalance.spcvy[pos]=-0.1977
|
---|
38 | elif (L==20000.):
|
---|
39 | md.stressbalance.spcvx[pos]=16.53
|
---|
40 | md.stressbalance.spcvy[pos]=-1.27
|
---|
41 | elif (L==40000.):
|
---|
42 | md.stressbalance.spcvx[pos]=17.23
|
---|
43 | md.stressbalance.spcvy[pos]=-3.17
|
---|
44 | elif (L==80000.):
|
---|
45 | md.stressbalance.spcvx[pos]=16.68
|
---|
46 | md.stressbalance.spcvy[pos]=-2.69
|
---|
47 | elif (L==160000.):
|
---|
48 | md.stressbalance.spcvx[pos]=16.03
|
---|
49 | md.stressbalance.spcvy[pos]=-1.27
|
---|
50 |
|
---|
51 | md=setflowequation(md,'FS','all')
|
---|
52 |
|
---|
53 | #Compute the stressbalance
|
---|
54 | md.cluster=generic('name',oshostname(),'np',8)
|
---|
55 | md=solve(md,'StressbalanceSolution')
|
---|
56 |
|
---|
57 | #Plot the results and save them
|
---|
58 | vx=md.results.StressbalanceSolution.Vx
|
---|
59 | vy=md.results.StressbalanceSolution.Vy
|
---|
60 | vz=md.results.StressbalanceSolution.Vz
|
---|
61 | results.append(md.results.StressbalanceSolution)
|
---|
62 |
|
---|
63 | # plotmodel(md,'data',vx,'data',vy,'data',vz,'layer#all',md.mesh.numberoflayers)
|
---|
64 |
|
---|
65 | #Fields and tolerances to track changes
|
---|
66 | field_names =[\
|
---|
67 | 'Vx5km','Vy5km','Vz5km',\
|
---|
68 | 'Vx10km','Vy10km','Vz10km',\
|
---|
69 | 'Vx20km','Vy20km','Vz20km',\
|
---|
70 | 'Vx40km','Vy40km','Vz40km',\
|
---|
71 | 'Vx80km','Vy80km','Vz80km',\
|
---|
72 | 'Vx160km','Vy160km','Vz160km'
|
---|
73 | ]
|
---|
74 | field_tolerances=[\
|
---|
75 | 1e-12,1e-12,1e-11,\
|
---|
76 | 1e-12,1e-12,1e-12,\
|
---|
77 | 1e-12,1e-12,1e-12,\
|
---|
78 | 1e-12,1e-12,1e-12,\
|
---|
79 | 1e-12,1e-12,1e-12,\
|
---|
80 | 1e-12,1e-11,1e-12,\
|
---|
81 | ]
|
---|
82 | field_values=[]
|
---|
83 | for result in results:
|
---|
84 | field_values=field_values+[\
|
---|
85 | result.Vx,\
|
---|
86 | result.Vy,\
|
---|
87 | result.Vz,\
|
---|
88 | ]
|
---|