Changeset 13462

Timestamp:

09/27/12 09:44:58 (12 years ago)

Author:

jschierm

Message:

CHG: Changed all data 1-d vectors to 2-d arrays to allow multiple columns.

Location:

issm/trunk-jpl/src/m

Files:

• boundaryconditions/SetIceShelfBC.py (modified) (6 diffs)
• classes/model/model.py (modified) (2 diffs)
• classes/thermal.py (modified) (1 diff)
• materials/paterson.m (modified) (1 diff)
• materials/paterson.py (modified) (2 diffs)
• parameterization/setflowequation.py (modified) (2 diffs)
• solve/parseresultsfromdisk.py (modified) (1 diff)

issm/trunk-jpl/src/m/boundaryconditions/SetIceShelfBC.py

@@ -1 +1 @@
 import os
 import numpy
-from ContourToMesh import ContourToMesh
+from ContourToMesh import *
 
 def SetIceShelfBC(md,icefrontfile=''):
@@ -7 +7 @@
         SETICESHELFBC - Create the boundary conditions for diagnostic and thermal models for a  Ice Shelf with Ice Front
 
-           Neumann BC are used on the ice front (an ANRGUS contour around the ice front
+           Neumann BC are used on the ice front (an ARGUS contour around the ice front
            must be given in input)
            Dirichlet BC are used elsewhere for diagnostic
@@ -25 +25 @@
                 if not os.path.exists(icefrontfile):
                         raise IOError("SetIceShelfBC error message: ice front file '%s' not found." % icefrontfile)
-                nodeinsideicefront=ContourToMesh((md.mesh.elements).reshape(-1,1),(md.mesh.x).reshape(-1,1),(md.mesh.y).reshape(-1,1),icefrontfile,'node',2)
+                nodeinsideicefront=ContourToMesh(md.mesh.elements,md.mesh.x.reshape(-1,1),md.mesh.y.reshape(-1,1),icefrontfile,'node',2)
                 nodeonicefront= numpy.bitwise_and( map(int,md.mesh.vertexonboundary), map(int,nodeinsideicefront[0].ravel()) )
         else:
-                nodeonicefront=numpy.zeros(md.mesh.numberofvertices)
+                nodeonicefront=numpy.zeros((md.mesh.numberofvertices,1))
 
 #       pos=find(md.mesh.vertexonboundary & ~nodeonicefront);
         pos=[i for i,(vob,noif) in enumerate(zip(md.mesh.vertexonboundary,nodeonicefront)) if vob and not noif]
-        md.diagnostic.spcvx=float('NaN')*numpy.ones(md.mesh.numberofvertices)
-        md.diagnostic.spcvy=float('NaN')*numpy.ones(md.mesh.numberofvertices)
-        md.diagnostic.spcvz=float('NaN')*numpy.ones(md.mesh.numberofvertices)
+        md.diagnostic.spcvx=float('NaN')*numpy.ones((md.mesh.numberofvertices,1))
+        md.diagnostic.spcvy=float('NaN')*numpy.ones((md.mesh.numberofvertices,1))
+        md.diagnostic.spcvz=float('NaN')*numpy.ones((md.mesh.numberofvertices,1))
         md.diagnostic.spcvx[pos]=0
         md.diagnostic.spcvy[pos]=0
@@ -41 +41 @@
 
         #Dirichlet Values
-        if numpy.size(md.inversion.vx_obs)==md.mesh.numberofvertices and numpy.size(md.inversion.vy_obs)==md.mesh.numberofvertices:
+        if isinstance(md.inversion.vx_obs,numpy.ndarray) and numpy.size(md.inversion.vx_obs,axis=0)==md.mesh.numberofvertices and isinstance(md.inversion.vy_obs,numpy.ndarray) and numpy.size(md.inversion.vy_obs,axis=0)==md.mesh.numberofvertices:
                 print '      boundary conditions for diagnostic model: spc set as observed velocities'
                 md.diagnostic.spcvx[pos]=md.inversion.vx_obs[pos]
@@ -56 +56 @@
         elif md.mesh.dimension==3:
 #               pressureload_layer1=[md.mesh.segments(pos,1:2)  md.mesh.segments(pos,2)+md.mesh.numberofvertices2d  md.mesh.segments(pos,1)+md.mesh.numberofvertices2d  md.mesh.segments(pos,3)];
-                pressureload_layer1=numpy.concatenate((md.mesh.segments[pos,0:2],md.mesh.segments[pos,1]+md.mesh.numberofvertices2d,md.mesh.segments[pos,0]+md.mesh.numberofvertices2d,md.mesh.segments[pos,2]),axis=1)
+                pressureload_layer1=numpy.hstack((md.mesh.segments[pos,0:2],md.mesh.segments[pos,1]+md.mesh.numberofvertices2d,md.mesh.segments[pos,0]+md.mesh.numberofvertices2d,md.mesh.segments[pos,2]))
                 pressureload=numpy.zeros((0,5))
                 for i in xrange(1,md.mesh.numberoflayers):
 #                       pressureload=[pressureload ;pressureload_layer1(:,1:4)+(i-1)*md.mesh.numberofvertices2d pressureload_layer1(:,5)+(i-1)*md.mesh.numberofelements2d ];
-                        pressureload=numpy.concatenate((pressureload,numpy.concatenate((pressureload_layer1[:,0:3]+(i-1)*md.mesh.numberofvertices2d,pressureload_layer1[:,4]+(i-1)*md.mesh.numberofelements2d),axis=1)),axis=0)
+                        pressureload=numpy.vstack((pressureload,numpy.hstack((pressureload_layer1[:,0:3]+(i-1)*md.mesh.numberofvertices2d,pressureload_layer1[:,4]+(i-1)*md.mesh.numberofelements2d))))
 
         #Add water or air enum depending on the element
 #       pressureload=[pressureload 1*md.mask.elementonfloatingice(pressureload(:,end))];
-        pressureload=numpy.concatenate((pressureload,1*md.mask.elementonfloatingice[pressureload[:,-1].astype('int')-1].reshape((-1,1))),axis=1)
+        pressureload=numpy.hstack((pressureload,1*md.mask.elementonfloatingice[pressureload[:,-1].astype('int')-1].reshape((-1,1))))
 
         #plug onto model
@@ -71 +71 @@
         #Create zeros basalforcings and surfaceforcings
         if numpy.all(numpy.isnan(md.surfaceforcings.precipitation)) and (md.surfaceforcings.ispdd==1):
-                md.surfaceforcings.precipitation=numpy.zeros(md.mesh.numberofvertices)
+                md.surfaceforcings.precipitation=numpy.zeros((md.mesh.numberofvertices,1))
                 print '      no surfaceforcings.precipitation specified: values set as zero'
         if numpy.all(numpy.isnan(md.surfaceforcings.mass_balance)) and (md.surfaceforcings.ispdd==0):
-                md.surfaceforcings.mass_balance=numpy.zeros(md.mesh.numberofvertices)
+                md.surfaceforcings.mass_balance=numpy.zeros((md.mesh.numberofvertices,1))
                 print '      no surfaceforcings.mass_balance specified: values set as zero'
         if numpy.all(numpy.isnan(md.basalforcings.melting_rate)):
-                md.basalforcings.melting_rate=numpy.zeros(md.mesh.numberofvertices)
+                md.basalforcings.melting_rate=numpy.zeros((md.mesh.numberofvertices,1))
                 print '      no basalforcings.melting_rate specified: values set as zero'
         if numpy.all(numpy.isnan(md.balancethickness.thickening_rate)):
-                md.balancethickness.thickening_rate=numpy.zeros(md.mesh.numberofvertices)
+                md.balancethickness.thickening_rate=numpy.zeros((md.mesh.numberofvertices,1))
                 print '      no balancethickness.thickening_rate specified: values set as zero'
 
-        md.prognostic.spcthickness=float('NaN')*numpy.ones(md.mesh.numberofvertices)
-        md.balancethickness.spcthickness=float('NaN')*numpy.ones(md.mesh.numberofvertices)
+        md.prognostic.spcthickness=float('NaN')*numpy.ones((md.mesh.numberofvertices,1))
+        md.balancethickness.spcthickness=float('NaN')*numpy.ones((md.mesh.numberofvertices,1))
 
         if numpy.size(md.initialization.temperature)==md.mesh.numberofvertices:
-                md.thermal.spctemperature=float('NaN')*numpy.ones(md.mesh.numberofvertices)
+                md.thermal.spctemperature=float('NaN')*numpy.ones((md.mesh.numberofvertices,1))
 #               pos=find(md.mesh.vertexonsurface); md.thermal.spctemperature(pos)=md.initialization.temperature(pos); %impose observed temperature on surface
                 pos=[i for i,vos in enumerate(md.mesh.vertexonsurface) if vos]
                 md.thermal.spctemperature[pos]=md.initialization.temperature[pos]    # impose observed temperature on surface
                 if not numpy.size(md.basalforcings.geothermalflux)==md.mesh.numberofvertices:
-                        md.basalforcings.geothermalflux=numpy.zeros(md.mesh.numberofvertices)
+                        md.basalforcings.geothermalflux=numpy.zeros((md.mesh.numberofvertices,1))
         else:
                 print '      no thermal boundary conditions created: no observed temperature found'

issm/trunk-jpl/src/m/classes/model/model.py

@@ -242 +242 @@
                         y3d=numpy.concatenate((y3d,md.mesh.y))
                         #nodes are distributed between bed and surface accordingly to the given exponent
-                        z3d=numpy.concatenate((z3d,bed3d+thickness3d*extrusionlist[i]))
+                        z3d=numpy.concatenate((z3d,(bed3d+thickness3d*extrusionlist[i]).reshape(-1)))
                 number_nodes3d=numpy.size(x3d)    #number of 3d nodes for the non extruded part of the mesh
 
@@ -402 +402 @@
                 #Put lithostatic pressure if there is an existing pressure
                 if not numpy.any(numpy.isnan(md.initialization.pressure)):
-                        md.initialization.pressure=md.constants.g*md.materials.rho_ice*(md.geometry.surface-md.mesh.z)
+                        md.initialization.pressure=md.constants.g*md.materials.rho_ice*(md.geometry.surface-md.mesh.z.reshape(-1,1))
 
                 #special for thermal modeling:

issm/trunk-jpl/src/m/classes/thermal.py

@@ -72 +72 @@
                 if EnthalpyAnalysisEnum() in analyses and (md.thermal.isenthalpy or solution==EnthalpySolutionEnum()) and md.mesh.dimension==3:
                         pos=numpy.nonzero(numpy.logical_not(numpy.isnan(md.thermal.spctemperature[0:md.mesh.numberofvertices])))
-#                       need to know if md.thermal.spctemperature can have multiple columns?
-#                       replicate=numpy.tile(md.geometry.surface-md.mesh.z,(1,numpy.size(md.thermal.spctemperature,axis=1)))
-#                       md = checkfield(md,'thermal.spctemperature[numpy.nonzero(numpy.logical_not(numpy.isnan(md.thermal.spctemperature[0:md.mesh.numberofvertices,:])))]','<',md.materials.meltingpoint-md.materials.beta*md.materials.rho_ice*md.constants.g*replicate[pos],'message',"spctemperature should be below the adjusted melting point")
-                        replicate=numpy.tile(md.geometry.surface-md.mesh.z,(1))
-                        md = checkfield(md,'thermal.spctemperature[numpy.nonzero(numpy.logical_not(numpy.isnan(md.thermal.spctemperature[0:md.mesh.numberofvertices])))]','<',md.materials.meltingpoint-md.materials.beta*md.materials.rho_ice*md.constants.g*replicate[pos],'message',"spctemperature should be below the adjusted melting point")
+                        replicate=numpy.tile(md.geometry.surface-md.mesh.z,(1,numpy.size(md.thermal.spctemperature,axis=1)))
+                        md = checkfield(md,'thermal.spctemperature[numpy.nonzero(numpy.logical_not(numpy.isnan(md.thermal.spctemperature[0:md.mesh.numberofvertices,:])))]','<',md.materials.meltingpoint-md.materials.beta*md.materials.rho_ice*md.constants.g*replicate[pos],'message',"spctemperature should be below the adjusted melting point")
                         md = checkfield(md,'thermal.isenthalpy','numel',[1],'values',[0,1])
 

issm/trunk-jpl/src/m/materials/paterson.m

@@ -2 +2 @@
 %PATERSON - figure out the rigidity of ice for a given temperature
 %
-%   rigidigty (in s^(1/3)Pa) is the flow law paramter in the flow law sigma=B*e(1/3) (Paterson, p97). 
+%   rigidity (in s^(1/3)Pa) is the flow law paramter in the flow law sigma=B*e(1/3) (Paterson, p97). 
 %   temperature is in Kelvin degrees
 %

issm/trunk-jpl/src/m/materials/paterson.py

@@ -1 @@
-from numpy import *
+import numpy
 
 def paterson(temperature):
+        """
+        PATERSON - figure out the rigidity of ice for a given temperature
 
-    # Local Variables: pos11, pos5, pos10, temperature, pos, T, pos8, pos9, pos6, pos7, pos4, rigidity, pos2, pos3, pos1
-    # Function calls: length, zeros, argwhere, paterson, error
-    #PATERSON - figure out the rigidity of ice for a given temperature
-    #
-    #   rigidigty (in s^(1/3)Pa) is the flow law paramter in the flow law sigma=B*e(1/3) (Paterson, p97). 
-    #   temperature is in Kelvin degrees
-    #
-    #   Usage:
-    #      rigidity=paterson(temperature)
-    
-        pos=argwhere(temperature<0.)
-        if len(pos):
-                print 'input temperature should be in Kelvin (positive)'
-                return []
-    
+           rigidity (in s^(1/3)Pa) is the flow law paramter in the flow law sigma=B*e(1/3) (Paterson, p97). 
+           temperature is in Kelvin degrees
+
+           Usage:
+              rigidity=paterson(temperature)
+        """
+        
+        if numpy.any(temperature<0.):
+                raise RuntimeError("input temperature should be in Kelvin (positive)")
+        
         T = temperature-273.15
+
         #The routine below is equivalent to:
+
         # n=3; T=temperature-273;
         # %From paterson,
@@ -31 +30 @@
         # rigidity=fittedmodel(temperature);
 
-        rigidity=zeros(len(T))
-        pos1=argwhere(T<=-45);           rigidity[pos1]=10**8*(-0.000292866376675*(T[pos1]+50)**3+ 0.011672640664130*(T[pos1]+50)**2  -0.325004442485481*(T[pos1]+50)+  6.524779401948101)
-        pos2=argwhere(logical_and(-45<=T,T<-40));   rigidity[pos2]=10**8*(-0.000292866376675*(T[pos2]+45)**3+ 0.007279645014004*(T[pos2]+45)**2  -0.230243014094813*(T[pos2]+45)+  5.154964909039554)
-        pos3=argwhere(logical_and(-40<=T,T<-35));   rigidity[pos3]=10**8*(0.000072737147457*(T[pos3]+40)**3+  0.002886649363879*(T[pos3]+40)**2  -0.179411542205399*(T[pos3]+40)+  4.149132666831214)
-        pos4=argwhere(logical_and(-35<=T,T<-30));   rigidity[pos4]=10**8*(-0.000086144770023*(T[pos4]+35)**3+ 0.003977706575736*(T[pos4]+35)**2  -0.145089762507325*(T[pos4]+35)+  3.333333333333331)
-        pos5=argwhere(logical_and(-30<=T,T<-25));   rigidity[pos5]=10**8*(-0.000043984685769*(T[pos5]+30)**3+ 0.002685535025386*(T[pos5]+30)**2  -0.111773554501713*(T[pos5]+30)+  2.696559088937191)
-        pos6=argwhere(logical_and(-25<=T,T<-20));   rigidity[pos6]=10**8*(-0.000029799523463*(T[pos6]+25)**3+ 0.002025764738854*(T[pos6]+25)**2  -0.088217055680511*(T[pos6]+25)+  2.199331606342181)
-        pos7=argwhere(logical_and(-20<=T,T<-15));   rigidity[pos7]=10**8*(0.000136920904777*(T[pos7]+20)**3+  0.001578771886910*(T[pos7]+20)**2  -0.070194372551690*(T[pos7]+20)+  1.805165505978111)
-        pos8=argwhere(logical_and(-15<=T,T<-10));   rigidity[pos8]=10**8*(-0.000899763781026*(T[pos8]+15)**3+ 0.003632585458564*(T[pos8]+15)**2  -0.044137585824322*(T[pos8]+15)+  1.510778053489523)
-        pos9=argwhere(logical_and(-10<=T,T<-5));    rigidity[pos9]=10**8*(0.001676964325070*(T[pos9]+10)**3-  0.009863871256831*(T[pos9]+10)**2  -0.075294014815659*(T[pos9]+10)+  1.268434288203714)
-        pos10=argwhere(logical_and(-5<=T,T<-2));    rigidity[pos10]=10**8*(-0.003748937622487*(T[pos10]+5)**3+0.015290593619213*(T[pos10]+5)**2  -0.048160403003748*(T[pos10]+5)+  0.854987973338348)
-        pos11=argwhere(-2<=T);           rigidity[pos11]=10**8*(-0.003748937622488*(T[pos11]+2)**3-0.018449844983174*(T[pos11]+2)**2  -0.057638157095631*(T[pos11]+2)+  0.746900791092860)
+        rigidity=numpy.zeros((numpy.size(T,axis=0),1))
+        pos1=numpy.nonzero(T<=-45)
+        rigidity[pos1]=10**8*(-0.000292866376675*(T[pos1]+50)**3+ 0.011672640664130*(T[pos1]+50)**2  -0.325004442485481*(T[pos1]+50)+  6.524779401948101)
+        pos2=numpy.nonzero(numpy.logical_and(-45<=T,T<-40))
+        rigidity[pos2]=10**8*(-0.000292866376675*(T[pos2]+45)**3+ 0.007279645014004*(T[pos2]+45)**2  -0.230243014094813*(T[pos2]+45)+  5.154964909039554)
+        pos3=numpy.nonzero(numpy.logical_and(-40<=T,T<-35))
+        rigidity[pos3]=10**8*(0.000072737147457*(T[pos3]+40)**3+  0.002886649363879*(T[pos3]+40)**2  -0.179411542205399*(T[pos3]+40)+  4.149132666831214)
+        pos4=numpy.nonzero(numpy.logical_and(-35<=T,T<-30))
+        rigidity[pos4]=10**8*(-0.000086144770023*(T[pos4]+35)**3+ 0.003977706575736*(T[pos4]+35)**2  -0.145089762507325*(T[pos4]+35)+  3.333333333333331)
+        pos5=numpy.nonzero(numpy.logical_and(-30<=T,T<-25))
+        rigidity[pos5]=10**8*(-0.000043984685769*(T[pos5]+30)**3+ 0.002685535025386*(T[pos5]+30)**2  -0.111773554501713*(T[pos5]+30)+  2.696559088937191)
+        pos6=numpy.nonzero(numpy.logical_and(-25<=T,T<-20))
+        rigidity[pos6]=10**8*(-0.000029799523463*(T[pos6]+25)**3+ 0.002025764738854*(T[pos6]+25)**2  -0.088217055680511*(T[pos6]+25)+  2.199331606342181)
+        pos7=numpy.nonzero(numpy.logical_and(-20<=T,T<-15))
+        rigidity[pos7]=10**8*(0.000136920904777*(T[pos7]+20)**3+  0.001578771886910*(T[pos7]+20)**2  -0.070194372551690*(T[pos7]+20)+  1.805165505978111)
+        pos8=numpy.nonzero(numpy.logical_and(-15<=T,T<-10))
+        rigidity[pos8]=10**8*(-0.000899763781026*(T[pos8]+15)**3+ 0.003632585458564*(T[pos8]+15)**2  -0.044137585824322*(T[pos8]+15)+  1.510778053489523)
+        pos9=numpy.nonzero(numpy.logical_and(-10<=T,T<-5))
+        rigidity[pos9]=10**8*(0.001676964325070*(T[pos9]+10)**3-  0.009863871256831*(T[pos9]+10)**2  -0.075294014815659*(T[pos9]+10)+  1.268434288203714)
+        pos10=numpy.nonzero(numpy.logical_and(-5<=T,T<-2))
+        rigidity[pos10]=10**8*(-0.003748937622487*(T[pos10]+5)**3+0.015290593619213*(T[pos10]+5)**2  -0.048160403003748*(T[pos10]+5)+  0.854987973338348)
+        pos11=numpy.nonzero(-2<=T)
+        rigidity[pos11]=10**8*(-0.003748937622488*(T[pos11]+2)**3-0.018449844983174*(T[pos11]+2)**2  -0.057638157095631*(T[pos11]+2)+  0.746900791092860)
 
         #Now make sure that rigidity is positive
-        pos=argwhere(rigidity<0);        rigidity[pos]=1**6 
+        pos=numpy.nonzero(rigidity<0)
+        rigidity[pos]=1**6 
 
         return rigidity

issm/trunk-jpl/src/m/parameterization/setflowequation.py

@@ -94 +94 @@
                                               numpy.logical_not(numpy.isnan(md.diagnostic.spcvy)).astype(int)+ \
                                               numpy.logical_not(numpy.isnan(md.diagnostic.spcvz)).astype(int)==3, \
-                                              numpy.logical_and(nodeonpattyn,nodeonstokes)).astype(int)    #find all the nodes on the boundary of the domain without icefront
+                                              numpy.logical_and(nodeonpattyn,nodeonstokes).reshape(-1,1)).astype(int)    #find all the nodes on the boundary of the domain without icefront
 #               fullspcelems=double(sum(fullspcnodes(md.mesh.elements),2)==6);         %find all the nodes on the boundary of the domain without icefront
                 fullspcelems=(numpy.sum(fullspcnodes[md.mesh.elements.astype(int)-1],axis=1)==6).astype(int)    #find all the nodes on the boundary of the domain without icefront
-                stokesflag[numpy.nonzero(fullspcelems)]=0
+                stokesflag[numpy.nonzero(fullspcelems.reshape(-1))]=0
                 nodeonstokes[md.mesh.elements[numpy.nonzero(stokesflag),:].astype(int)-1]=1
 
@@ -129 +129 @@
                         penalties=numpy.zeros((0,2))
                         for     i in xrange(1,numlayers):
-                                penalties=numpy.concatenate((penalties,numpy.concatenate((bordernodes2d,bordernodes2d+md.mesh.numberofvertices2d*(i)),axis=1)),axis=0)
+                                penalties=numpy.vstack((penalties,numpy.hstack((bordernodes2d,bordernodes2d+md.mesh.numberofvertices2d*(i)))))
                         md.diagnostic.vertex_pairing=penalties
 

issm/trunk-jpl/src/m/solve/parseresultsfromdisk.py

@@ -52 +52 @@
                    result['fieldname'] in results[result['step']] and \
                    not strcmp(result['fieldname'],'SolutionType'):
-                        results[result['step']][result['fieldname']]=numpy.concatenate((results[result['step']][result['fieldname']],result['field']),axis=0)
+                        results[result['step']][result['fieldname']]=numpy.vstack((results[result['step']][result['fieldname']],result['field']))
                 else:
                         results[result['step']][result['fieldname']]=result['field']