source: issm/oecreview/Archive/16554-17801/ISSM-17723-17724.diff@ 17802

Last change on this file since 17802 was 17802, checked in by Mathieu Morlighem, 11 years ago

Added archives
File size: 10.4 KB
RevLine 
[17802]1Index: ../trunk-jpl/src/m/classes/model.py
2===================================================================
3--- ../trunk-jpl/src/m/classes/model.py (revision 17723)
4+++ ../trunk-jpl/src/m/classes/model.py (revision 17724)
5@@ -41,15 +41,16 @@
6 from radaroverlay import radaroverlay
7 from miscellaneous import miscellaneous
8 from private import private
9-from EnumDefinitions import *
10 from mumpsoptions import mumpsoptions
11 from iluasmoptions import iluasmoptions
12 from project3d import project3d
13+from project2d import project2d
14 from FlagElements import FlagElements
15 from NodeConnectivity import NodeConnectivity
16 from ElementConnectivity import ElementConnectivity
17 from contourenvelope import contourenvelope
18 import MatlabFuncs as m
19+from DepthAverage import DepthAverage
20 #}}}
21
22 class model(object):
23@@ -174,7 +175,7 @@
24 return string
25 # }}}
26 def checkmessage(self,string): # {{{
27- print ("model not consistent: %s" % string)
28+ print "model not consistent: ", string
29 self.private.isconsistent=False
30 return self
31 # }}}
32@@ -327,7 +328,7 @@
33
34 #Edges
35 if m.strcmp(md.mesh.domaintype(),'2Dhorizontal'):
36- if numpy.ndim(md2.mesh.edges)>1 and numpy.size(md2.mesh.edges,axis=1)>1: #do not use ~isnan because there are some NaNs...
37+ if numpy.ndim(md2.mesh.edges)>1 and numpy.size(md2.mesh.edges,axis=1)>1: #do not use ~isnan because there are some npy.nans...
38 #renumber first two columns
39 pos=numpy.nonzero(md2.mesh.edges[:,3]!=-1)[0]
40 md2.mesh.edges[: ,0]=Pnode[md2.mesh.edges[:,0]-1]
41@@ -393,8 +394,8 @@
42 md2.stressbalance.spcvx[nodestoflag2]=md2.inversion.vx_obs[nodestoflag2]
43 md2.stressbalance.spcvy[nodestoflag2]=md2.inversion.vy_obs[nodestoflag2]
44 else:
45- md2.stressbalance.spcvx[nodestoflag2]=float('NaN')
46- md2.stressbalance.spcvy[nodestoflag2]=float('NaN')
47+ md2.stressbalance.spcvx[nodestoflag2]=npy.nan
48+ md2.stressbalance.spcvy[nodestoflag2]=npy.nan
49 print "\n!! extract warning: spc values should be checked !!\n\n"
50 #put 0 for vz
51 md2.stressbalance.spcvz[nodestoflag2]=0
52@@ -639,9 +640,9 @@
53 md.stressbalance.spcvx=project3d(md,'vector',md.stressbalance.spcvx,'type','node')
54 md.stressbalance.spcvy=project3d(md,'vector',md.stressbalance.spcvy,'type','node')
55 md.stressbalance.spcvz=project3d(md,'vector',md.stressbalance.spcvz,'type','node')
56- md.thermal.spctemperature=project3d(md,'vector',md.thermal.spctemperature,'type','node','layer',md.mesh.numberoflayers,'padding',float('NaN'))
57+ md.thermal.spctemperature=project3d(md,'vector',md.thermal.spctemperature,'type','node','layer',md.mesh.numberoflayers,'padding',numpy.nan)
58 if isinstance(md.initialization.temperature,numpy.ndarray) and numpy.size(md.initialization.temperature,axis=0)==md.mesh.numberofvertices:
59- md.thermal.spctemperature=float('nan')*numpy.ones((md.mesh.numberofvertices,1))
60+ md.thermal.spctemperature=numpy.nan*numpy.ones((md.mesh.numberofvertices,1))
61 if hasattr(md.mesh,'vertexonsurface'):
62 pos=numpy.nonzero(md.mesh.vertexonsurface)[0]
63 md.thermal.spctemperature[pos]=md.initialization.temperature[pos] #impose observed temperature on surface
64@@ -702,3 +703,125 @@
65
66 return md
67 # }}}
68+ def collapse(md): #{{{
69+ '''
70+ collapses a 3d mesh into a 2d mesh
71+
72+ This routine collapses a 3d model into a 2d model and collapses all
73+ the fileds of the 3d model by taking their depth-averaged values
74+
75+ Usage:
76+ md=collapse(md)
77+ '''
78+
79+ #Check that the model is really a 3d model
80+ if md.mesh.domaintype().lower() != '3d':
81+ raise StandardError("only a 3D model can be collapsed")
82+
83+ #drag is limited to nodes that are on the bedrock.
84+ md.friction.coefficient=project2d(md,md.friction.coefficient,1)
85+
86+ #p and q (same deal, except for element that are on the bedrock: )
87+ md.friction.p=project2d(md,md.friction.p,1)
88+ md.friction.q=project2d(md,md.friction.q,1)
89+
90+ #observations
91+ if not numpy.isnan(md.inversion.vx_obs).all(): md.inversion.vx_obs=project2d(md,md.inversion.vx_obs,md.mesh.numberoflayers)
92+ if not numpy.isnan(md.inversion.vy_obs).all(): md.inversion.vy_obs=project2d(md,md.inversion.vy_obs,md.mesh.numberoflayers)
93+ if not numpy.isnan(md.inversion.vel_obs).all(): md.inversion.vel_obs=project2d(md,md.inversion.vel_obs,md.mesh.numberoflayers)
94+ if not numpy.isnan(md.inversion.cost_functions_coefficients).all(): md.inversion.cost_functions_coefficients=project2d(md,md.inversion.cost_functions_coefficients,md.mesh.numberoflayers)
95+ if md.inversion.min_parameters.size>1: md.inversion.min_parameters=project2d(md,md.inversion.min_parameters,md.mesh.numberoflayers)
96+ if md.inversion.max_parameters.size>1: md.inversion.max_parameters=project2d(md,md.inversion.max_parameters,md.mesh.numberoflayers)
97+ if not numpy.isnan(md.surfaceforcings.mass_balance).all():
98+ md.surfaceforcings.mass_balance=project2d(md,md.surfaceforcings.mass_balance,md.mesh.numberoflayers)
99+
100+ if not numpy.isnan(md.balancethickness.thickening_rate).all(): md.balancethickness.thickening_rate=project2d(md,md.balancethickness.thickening_rate,md.mesh.numberoflayers)
101+
102+ #results
103+ if not numpy.isnan(md.initialization.vx).all(): md.initialization.vx=DepthAverage(md,md.initialization.vx)
104+ if not numpy.isnan(md.initialization.vy).all(): md.initialization.vy=DepthAverage(md,md.initialization.vy)
105+ if not numpy.isnan(md.initialization.vz).all(): md.initialization.vz=DepthAverage(md,md.initialization.vz)
106+ if not numpy.isnan(md.initialization.vel).all(): md.initialization.vel=DepthAverage(md,md.initialization.vel)
107+ if not numpy.isnan(md.initialization.temperature).all(): md.initialization.temperature=DepthAverage(md,md.initialization.temperature)
108+
109+ #gia
110+ if not numpy.isnan(md.gia.mantle_viscosity).all(): md.gia.mantle_viscosity=project2d(md,md.gia.mantle_viscosity,1)
111+ if not numpy.isnan(md.gia.lithosphere_thickness).all(): md.gia.lithosphere_thickness=project2d(md,md.gia.lithosphere_thickness,1)
112+
113+ #elementstype
114+ if not numpy.isnan(md.flowequation.element_equation).all():
115+ md.flowequation.element_equation=project2d(md,md.flowequation.element_equation,1)
116+ md.flowequation.vertex_equation=project2d(md,md.flowequation.vertex_equation,1)
117+ md.flowequation.borderSSA=project2d(md,md.flowequation.borderSSA,1)
118+ md.flowequation.borderHO=project2d(md,md.flowequation.borderHO,1)
119+ md.flowequation.borderFS=project2d(md,md.flowequation.borderFS,1)
120+
121+ #boundary conditions
122+ md.stressbalance.spcvx=project2d(md,md.stressbalance.spcvx,md.mesh.numberoflayers)
123+ md.stressbalance.spcvy=project2d(md,md.stressbalance.spcvy,md.mesh.numberoflayers)
124+ md.stressbalance.spcvz=project2d(md,md.stressbalance.spcvz,md.mesh.numberoflayers)
125+ md.stressbalance.referential=project2d(md,md.stressbalance.referential,md.mesh.numberoflayers)
126+ md.stressbalance.loadingforce=project2d(md,md.stressbalance.loadingforce,md.mesh.numberoflayers)
127+ md.masstransport.spcthickness=project2d(md,md.masstransport.spcthickness,md.mesh.numberoflayers)
128+ md.damage.spcdamage=project2d(md,md.damage.spcdamage,md.mesh.numberoflayers-1)
129+ md.thermal.spctemperature=project2d(md,md.thermal.spctemperature,md.mesh.numberoflayers-1)
130+
131+ #materials
132+ md.materials.rheology_B=DepthAverage(md,md.materials.rheology_B)
133+ md.materials.rheology_n=project2d(md,md.materials.rheology_n,1)
134+
135+ #damage:
136+ md.damage.D=DepthAverage(md,md.damage.D)
137+
138+ #special for thermal modeling:
139+ md.basalforcings.melting_rate=project2d(md,md.basalforcings.melting_rate,1)
140+ md.basalforcings.geothermalflux=project2d(md,md.basalforcings.geothermalflux,1) #bedrock only gets geothermal flux
141+
142+ #update of connectivity matrix
143+ md.mesh.average_vertex_connectivity=25
144+
145+ #Collapse the mesh
146+ nodes2d=md.mesh.numberofvertices2d
147+ elements2d=md.mesh.numberofelements2d
148+
149+ #parameters
150+ md.geometry.surface=project2d(md,md.geometry.surface,1)
151+ md.geometry.thickness=project2d(md,md.geometry.thickness,1)
152+ md.geometry.base=project2d(md,md.geometry.base,1)
153+ md.geometry.bed=project2d(md,md.geometry.bed,1)
154+ md.mesh.vertexonboundary=project2d(md,md.mesh.vertexonboundary,1)
155+ md.mesh.elementconnectivity=project2d(md,md.mesh.elementconnectivity,1)
156+ md.mask.groundedice_levelset=project2d(md,md.mask.groundedice_levelset,1)
157+ md.mask.ice_levelset=project2d(md,md.mask.ice_levelset,1)
158+
159+ #lat long
160+ if md.mesh.lat.size==md.mesh.numberofvertices: md.mesh.lat=project2d(md,md.mesh.lat,1)
161+ if md.mesh.long.size==md.mesh.numberofvertices: md.mesh.long=project2d(md,md.mesh.long,1)
162+
163+ #Initialize with the 2d mesh
164+ mesh=mesh2d()
165+ mesh.x=md.mesh.x2d
166+ mesh.y=md.mesh.y2d
167+ mesh.z=numpy.zeros_like(md.mesh.x2d)
168+ mesh.numberofvertices=md.mesh.numberofvertices2d
169+ mesh.numberofelements=md.mesh.numberofelements2d
170+ mesh.elements=md.mesh.elements2d
171+ md.mesh=mesh
172+
173+ #Keep a trace of lower and upper nodes
174+ md.mesh.lowervertex=numpy.nan
175+ md.mesh.uppervertex=numpy.nan
176+ md.mesh.lowerelements=numpy.nan
177+ md.mesh.upperelements=numpy.nan
178+
179+ #Remove old mesh
180+ md.mesh.x2d=numpy.nan
181+ md.mesh.y2d=numpy.nan
182+ md.mesh.elements2d=numpy.nan
183+ md.mesh.numberofelements2d=md.mesh.numberofelements
184+ md.mesh.numberofvertices2d=md.mesh.numberofvertices
185+ md.mesh.numberoflayers=0
186+
187+ return md
188+
189+#}}}
190Index: ../trunk-jpl/src/m/classes/model.m
191===================================================================
192--- ../trunk-jpl/src/m/classes/model.m (revision 17723)
193+++ ../trunk-jpl/src/m/classes/model.m (revision 17724)
194@@ -159,7 +159,7 @@
195 error('collapse error message: only 3d mesh can be collapsed')
196 end
197
198- %Start with changing alle the fields from the 3d mesh
199+ %Start with changing all the fields from the 3d mesh
200
201 %drag is limited to nodes that are on the bedrock.
202 md.friction.coefficient=project2d(md,md.friction.coefficient,1);
203@@ -243,13 +243,14 @@
204 if numel(md.mesh.long)==md.mesh.numberofvertices, md.mesh.long=project2d(md,md.mesh.long,1); end
205
206 %Initialize with the 2d mesh
207- md.mesh=mesh2d();
208- md.mesh.x=md.mesh.x2d;
209- md.mesh.y=md.mesh.y2d;
210- md.mesh.z=zeros(size(md.mesh.x2d));
211- md.mesh.numberofvertices=md.mesh.numberofvertices2d;
212- md.mesh.numberofelements=md.mesh.numberofelements2d;
213- md.mesh.elements=md.mesh.elements2d;
214+ mesh=mesh2d();
215+ mesh.x=md.mesh.x2d;
216+ mesh.y=md.mesh.y2d;
217+ mesh.z=zeros(size(md.mesh.x2d));
218+ mesh.numberofvertices=md.mesh.numberofvertices2d;
219+ mesh.numberofelements=md.mesh.numberofelements2d;
220+ mesh.elements=md.mesh.elements2d;
221+ md.mesh=mesh;
222
223 %Keep a trace of lower and upper nodes
224 md.mesh.lowervertex=NaN;
Note: See TracBrowser for help on using the repository browser.