Changeset 13462


Ignore:
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:
7 edited

Legend:

Unmodified
Added
Removed
  • issm/trunk-jpl/src/m/boundaryconditions/SetIceShelfBC.py

    r13455 r13462  
    11import os
    22import numpy
    3 from ContourToMesh import ContourToMesh
     3from ContourToMesh import *
    44
    55def SetIceShelfBC(md,icefrontfile=''):
     
    77        SETICESHELFBC - Create the boundary conditions for diagnostic and thermal models for a  Ice Shelf with Ice Front
    88
    9            Neumann BC are used on the ice front (an ANRGUS contour around the ice front
     9           Neumann BC are used on the ice front (an ARGUS contour around the ice front
    1010           must be given in input)
    1111           Dirichlet BC are used elsewhere for diagnostic
     
    2525                if not os.path.exists(icefrontfile):
    2626                        raise IOError("SetIceShelfBC error message: ice front file '%s' not found." % icefrontfile)
    27                 nodeinsideicefront=ContourToMesh((md.mesh.elements).reshape(-1,1),(md.mesh.x).reshape(-1,1),(md.mesh.y).reshape(-1,1),icefrontfile,'node',2)
     27                nodeinsideicefront=ContourToMesh(md.mesh.elements,md.mesh.x.reshape(-1,1),md.mesh.y.reshape(-1,1),icefrontfile,'node',2)
    2828                nodeonicefront= numpy.bitwise_and( map(int,md.mesh.vertexonboundary), map(int,nodeinsideicefront[0].ravel()) )
    2929        else:
    30                 nodeonicefront=numpy.zeros(md.mesh.numberofvertices)
     30                nodeonicefront=numpy.zeros((md.mesh.numberofvertices,1))
    3131
    3232#       pos=find(md.mesh.vertexonboundary & ~nodeonicefront);
    3333        pos=[i for i,(vob,noif) in enumerate(zip(md.mesh.vertexonboundary,nodeonicefront)) if vob and not noif]
    34         md.diagnostic.spcvx=float('NaN')*numpy.ones(md.mesh.numberofvertices)
    35         md.diagnostic.spcvy=float('NaN')*numpy.ones(md.mesh.numberofvertices)
    36         md.diagnostic.spcvz=float('NaN')*numpy.ones(md.mesh.numberofvertices)
     34        md.diagnostic.spcvx=float('NaN')*numpy.ones((md.mesh.numberofvertices,1))
     35        md.diagnostic.spcvy=float('NaN')*numpy.ones((md.mesh.numberofvertices,1))
     36        md.diagnostic.spcvz=float('NaN')*numpy.ones((md.mesh.numberofvertices,1))
    3737        md.diagnostic.spcvx[pos]=0
    3838        md.diagnostic.spcvy[pos]=0
     
    4141
    4242        #Dirichlet Values
    43         if numpy.size(md.inversion.vx_obs)==md.mesh.numberofvertices and numpy.size(md.inversion.vy_obs)==md.mesh.numberofvertices:
     43        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:
    4444                print '      boundary conditions for diagnostic model: spc set as observed velocities'
    4545                md.diagnostic.spcvx[pos]=md.inversion.vx_obs[pos]
     
    5656        elif md.mesh.dimension==3:
    5757#               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)];
    58                 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)
     58                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]))
    5959                pressureload=numpy.zeros((0,5))
    6060                for i in xrange(1,md.mesh.numberoflayers):
    6161#                       pressureload=[pressureload ;pressureload_layer1(:,1:4)+(i-1)*md.mesh.numberofvertices2d pressureload_layer1(:,5)+(i-1)*md.mesh.numberofelements2d ];
    62                         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)
     62                        pressureload=numpy.vstack((pressureload,numpy.hstack((pressureload_layer1[:,0:3]+(i-1)*md.mesh.numberofvertices2d,pressureload_layer1[:,4]+(i-1)*md.mesh.numberofelements2d))))
    6363
    6464        #Add water or air enum depending on the element
    6565#       pressureload=[pressureload 1*md.mask.elementonfloatingice(pressureload(:,end))];
    66         pressureload=numpy.concatenate((pressureload,1*md.mask.elementonfloatingice[pressureload[:,-1].astype('int')-1].reshape((-1,1))),axis=1)
     66        pressureload=numpy.hstack((pressureload,1*md.mask.elementonfloatingice[pressureload[:,-1].astype('int')-1].reshape((-1,1))))
    6767
    6868        #plug onto model
     
    7171        #Create zeros basalforcings and surfaceforcings
    7272        if numpy.all(numpy.isnan(md.surfaceforcings.precipitation)) and (md.surfaceforcings.ispdd==1):
    73                 md.surfaceforcings.precipitation=numpy.zeros(md.mesh.numberofvertices)
     73                md.surfaceforcings.precipitation=numpy.zeros((md.mesh.numberofvertices,1))
    7474                print '      no surfaceforcings.precipitation specified: values set as zero'
    7575        if numpy.all(numpy.isnan(md.surfaceforcings.mass_balance)) and (md.surfaceforcings.ispdd==0):
    76                 md.surfaceforcings.mass_balance=numpy.zeros(md.mesh.numberofvertices)
     76                md.surfaceforcings.mass_balance=numpy.zeros((md.mesh.numberofvertices,1))
    7777                print '      no surfaceforcings.mass_balance specified: values set as zero'
    7878        if numpy.all(numpy.isnan(md.basalforcings.melting_rate)):
    79                 md.basalforcings.melting_rate=numpy.zeros(md.mesh.numberofvertices)
     79                md.basalforcings.melting_rate=numpy.zeros((md.mesh.numberofvertices,1))
    8080                print '      no basalforcings.melting_rate specified: values set as zero'
    8181        if numpy.all(numpy.isnan(md.balancethickness.thickening_rate)):
    82                 md.balancethickness.thickening_rate=numpy.zeros(md.mesh.numberofvertices)
     82                md.balancethickness.thickening_rate=numpy.zeros((md.mesh.numberofvertices,1))
    8383                print '      no balancethickness.thickening_rate specified: values set as zero'
    8484
    85         md.prognostic.spcthickness=float('NaN')*numpy.ones(md.mesh.numberofvertices)
    86         md.balancethickness.spcthickness=float('NaN')*numpy.ones(md.mesh.numberofvertices)
     85        md.prognostic.spcthickness=float('NaN')*numpy.ones((md.mesh.numberofvertices,1))
     86        md.balancethickness.spcthickness=float('NaN')*numpy.ones((md.mesh.numberofvertices,1))
    8787
    8888        if numpy.size(md.initialization.temperature)==md.mesh.numberofvertices:
    89                 md.thermal.spctemperature=float('NaN')*numpy.ones(md.mesh.numberofvertices)
     89                md.thermal.spctemperature=float('NaN')*numpy.ones((md.mesh.numberofvertices,1))
    9090#               pos=find(md.mesh.vertexonsurface); md.thermal.spctemperature(pos)=md.initialization.temperature(pos); %impose observed temperature on surface
    9191                pos=[i for i,vos in enumerate(md.mesh.vertexonsurface) if vos]
    9292                md.thermal.spctemperature[pos]=md.initialization.temperature[pos]    # impose observed temperature on surface
    9393                if not numpy.size(md.basalforcings.geothermalflux)==md.mesh.numberofvertices:
    94                         md.basalforcings.geothermalflux=numpy.zeros(md.mesh.numberofvertices)
     94                        md.basalforcings.geothermalflux=numpy.zeros((md.mesh.numberofvertices,1))
    9595        else:
    9696                print '      no thermal boundary conditions created: no observed temperature found'
  • issm/trunk-jpl/src/m/classes/model/model.py

    r13360 r13462  
    242242                        y3d=numpy.concatenate((y3d,md.mesh.y))
    243243                        #nodes are distributed between bed and surface accordingly to the given exponent
    244                         z3d=numpy.concatenate((z3d,bed3d+thickness3d*extrusionlist[i]))
     244                        z3d=numpy.concatenate((z3d,(bed3d+thickness3d*extrusionlist[i]).reshape(-1)))
    245245                number_nodes3d=numpy.size(x3d)    #number of 3d nodes for the non extruded part of the mesh
    246246
     
    402402                #Put lithostatic pressure if there is an existing pressure
    403403                if not numpy.any(numpy.isnan(md.initialization.pressure)):
    404                         md.initialization.pressure=md.constants.g*md.materials.rho_ice*(md.geometry.surface-md.mesh.z)
     404                        md.initialization.pressure=md.constants.g*md.materials.rho_ice*(md.geometry.surface-md.mesh.z.reshape(-1,1))
    405405
    406406                #special for thermal modeling:
  • issm/trunk-jpl/src/m/classes/thermal.py

    r13439 r13462  
    7272                if EnthalpyAnalysisEnum() in analyses and (md.thermal.isenthalpy or solution==EnthalpySolutionEnum()) and md.mesh.dimension==3:
    7373                        pos=numpy.nonzero(numpy.logical_not(numpy.isnan(md.thermal.spctemperature[0:md.mesh.numberofvertices])))
    74 #                       need to know if md.thermal.spctemperature can have multiple columns?
    75 #                       replicate=numpy.tile(md.geometry.surface-md.mesh.z,(1,numpy.size(md.thermal.spctemperature,axis=1)))
    76 #                       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")
    77                         replicate=numpy.tile(md.geometry.surface-md.mesh.z,(1))
    78                         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")
     74                        replicate=numpy.tile(md.geometry.surface-md.mesh.z,(1,numpy.size(md.thermal.spctemperature,axis=1)))
     75                        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")
    7976                        md = checkfield(md,'thermal.isenthalpy','numel',[1],'values',[0,1])
    8077
  • issm/trunk-jpl/src/m/materials/paterson.m

    r13011 r13462  
    22%PATERSON - figure out the rigidity of ice for a given temperature
    33%
    4 %   rigidigty (in s^(1/3)Pa) is the flow law paramter in the flow law sigma=B*e(1/3) (Paterson, p97).
     4%   rigidity (in s^(1/3)Pa) is the flow law paramter in the flow law sigma=B*e(1/3) (Paterson, p97).
    55%   temperature is in Kelvin degrees
    66%
  • issm/trunk-jpl/src/m/materials/paterson.py

    r13011 r13462  
    1 from numpy import *
     1import numpy
    22
    33def paterson(temperature):
     4        """
     5        PATERSON - figure out the rigidity of ice for a given temperature
    46
    5     # Local Variables: pos11, pos5, pos10, temperature, pos, T, pos8, pos9, pos6, pos7, pos4, rigidity, pos2, pos3, pos1
    6     # Function calls: length, zeros, argwhere, paterson, error
    7     #PATERSON - figure out the rigidity of ice for a given temperature
    8     #
    9     #   rigidigty (in s^(1/3)Pa) is the flow law paramter in the flow law sigma=B*e(1/3) (Paterson, p97).
    10     #   temperature is in Kelvin degrees
    11     #
    12     #   Usage:
    13     #      rigidity=paterson(temperature)
    14    
    15         pos=argwhere(temperature<0.)
    16         if len(pos):
    17                 print 'input temperature should be in Kelvin (positive)'
    18                 return []
    19    
     7           rigidity (in s^(1/3)Pa) is the flow law paramter in the flow law sigma=B*e(1/3) (Paterson, p97).
     8           temperature is in Kelvin degrees
     9
     10           Usage:
     11              rigidity=paterson(temperature)
     12        """
     13       
     14        if numpy.any(temperature<0.):
     15                raise RuntimeError("input temperature should be in Kelvin (positive)")
     16       
    2017        T = temperature-273.15
     18
    2119        #The routine below is equivalent to:
     20
    2221        # n=3; T=temperature-273;
    2322        # %From paterson,
     
    3130        # rigidity=fittedmodel(temperature);
    3231
    33         rigidity=zeros(len(T))
    34         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)
    35         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)
    36         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)
    37         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)
    38         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)
    39         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)
    40         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)
    41         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)
    42         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)
    43         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)
    44         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)
     32        rigidity=numpy.zeros((numpy.size(T,axis=0),1))
     33        pos1=numpy.nonzero(T<=-45)
     34        rigidity[pos1]=10**8*(-0.000292866376675*(T[pos1]+50)**3+ 0.011672640664130*(T[pos1]+50)**2  -0.325004442485481*(T[pos1]+50)+  6.524779401948101)
     35        pos2=numpy.nonzero(numpy.logical_and(-45<=T,T<-40))
     36        rigidity[pos2]=10**8*(-0.000292866376675*(T[pos2]+45)**3+ 0.007279645014004*(T[pos2]+45)**2  -0.230243014094813*(T[pos2]+45)+  5.154964909039554)
     37        pos3=numpy.nonzero(numpy.logical_and(-40<=T,T<-35))
     38        rigidity[pos3]=10**8*(0.000072737147457*(T[pos3]+40)**3+  0.002886649363879*(T[pos3]+40)**2  -0.179411542205399*(T[pos3]+40)+  4.149132666831214)
     39        pos4=numpy.nonzero(numpy.logical_and(-35<=T,T<-30))
     40        rigidity[pos4]=10**8*(-0.000086144770023*(T[pos4]+35)**3+ 0.003977706575736*(T[pos4]+35)**2  -0.145089762507325*(T[pos4]+35)+  3.333333333333331)
     41        pos5=numpy.nonzero(numpy.logical_and(-30<=T,T<-25))
     42        rigidity[pos5]=10**8*(-0.000043984685769*(T[pos5]+30)**3+ 0.002685535025386*(T[pos5]+30)**2  -0.111773554501713*(T[pos5]+30)+  2.696559088937191)
     43        pos6=numpy.nonzero(numpy.logical_and(-25<=T,T<-20))
     44        rigidity[pos6]=10**8*(-0.000029799523463*(T[pos6]+25)**3+ 0.002025764738854*(T[pos6]+25)**2  -0.088217055680511*(T[pos6]+25)+  2.199331606342181)
     45        pos7=numpy.nonzero(numpy.logical_and(-20<=T,T<-15))
     46        rigidity[pos7]=10**8*(0.000136920904777*(T[pos7]+20)**3+  0.001578771886910*(T[pos7]+20)**2  -0.070194372551690*(T[pos7]+20)+  1.805165505978111)
     47        pos8=numpy.nonzero(numpy.logical_and(-15<=T,T<-10))
     48        rigidity[pos8]=10**8*(-0.000899763781026*(T[pos8]+15)**3+ 0.003632585458564*(T[pos8]+15)**2  -0.044137585824322*(T[pos8]+15)+  1.510778053489523)
     49        pos9=numpy.nonzero(numpy.logical_and(-10<=T,T<-5))
     50        rigidity[pos9]=10**8*(0.001676964325070*(T[pos9]+10)**3-  0.009863871256831*(T[pos9]+10)**2  -0.075294014815659*(T[pos9]+10)+  1.268434288203714)
     51        pos10=numpy.nonzero(numpy.logical_and(-5<=T,T<-2))
     52        rigidity[pos10]=10**8*(-0.003748937622487*(T[pos10]+5)**3+0.015290593619213*(T[pos10]+5)**2  -0.048160403003748*(T[pos10]+5)+  0.854987973338348)
     53        pos11=numpy.nonzero(-2<=T)
     54        rigidity[pos11]=10**8*(-0.003748937622488*(T[pos11]+2)**3-0.018449844983174*(T[pos11]+2)**2  -0.057638157095631*(T[pos11]+2)+  0.746900791092860)
    4555
    4656        #Now make sure that rigidity is positive
    47         pos=argwhere(rigidity<0);        rigidity[pos]=1**6
     57        pos=numpy.nonzero(rigidity<0)
     58        rigidity[pos]=1**6
    4859
    4960        return rigidity
  • issm/trunk-jpl/src/m/parameterization/setflowequation.py

    r13244 r13462  
    9494                                              numpy.logical_not(numpy.isnan(md.diagnostic.spcvy)).astype(int)+ \
    9595                                              numpy.logical_not(numpy.isnan(md.diagnostic.spcvz)).astype(int)==3, \
    96                                               numpy.logical_and(nodeonpattyn,nodeonstokes)).astype(int)    #find all the nodes on the boundary of the domain without icefront
     96                                              numpy.logical_and(nodeonpattyn,nodeonstokes).reshape(-1,1)).astype(int)    #find all the nodes on the boundary of the domain without icefront
    9797#               fullspcelems=double(sum(fullspcnodes(md.mesh.elements),2)==6);         %find all the nodes on the boundary of the domain without icefront
    9898                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
    99                 stokesflag[numpy.nonzero(fullspcelems)]=0
     99                stokesflag[numpy.nonzero(fullspcelems.reshape(-1))]=0
    100100                nodeonstokes[md.mesh.elements[numpy.nonzero(stokesflag),:].astype(int)-1]=1
    101101
     
    129129                        penalties=numpy.zeros((0,2))
    130130                        for     i in xrange(1,numlayers):
    131                                 penalties=numpy.concatenate((penalties,numpy.concatenate((bordernodes2d,bordernodes2d+md.mesh.numberofvertices2d*(i)),axis=1)),axis=0)
     131                                penalties=numpy.vstack((penalties,numpy.hstack((bordernodes2d,bordernodes2d+md.mesh.numberofvertices2d*(i)))))
    132132                        md.diagnostic.vertex_pairing=penalties
    133133
  • issm/trunk-jpl/src/m/solve/parseresultsfromdisk.py

    r13116 r13462  
    5252                   result['fieldname'] in results[result['step']] and \
    5353                   not strcmp(result['fieldname'],'SolutionType'):
    54                         results[result['step']][result['fieldname']]=numpy.concatenate((results[result['step']][result['fieldname']],result['field']),axis=0)
     54                        results[result['step']][result['fieldname']]=numpy.vstack((results[result['step']][result['fieldname']],result['field']))
    5555                else:
    5656                        results[result['step']][result['fieldname']]=result['field']
Note: See TracChangeset for help on using the changeset viewer.