Ignore:
Timestamp:
11/16/12 08:10:16 (12 years ago)
Author:
Mathieu Morlighem
Message:

merged trunk-jpl and trunk for revision 13974

Location:
issm/trunk
Files:
3 edited

Legend:

Unmodified
Added
Removed
  • issm/trunk

  • issm/trunk/src

  • issm/trunk/src/m/boundaryconditions/SetIceShelfBC.py

    r13395 r13975  
    11import os
    22import numpy
     3from ContourToMesh import *
    34
    45def SetIceShelfBC(md,icefrontfile=''):
     
    67        SETICESHELFBC - Create the boundary conditions for diagnostic and thermal models for a  Ice Shelf with Ice Front
    78
    8            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
    910           must be given in input)
    1011           Dirichlet BC are used elsewhere for diagnostic
     
    2425                if not os.path.exists(icefrontfile):
    2526                        raise IOError("SetIceShelfBC error message: ice front file '%s' not found." % icefrontfile)
    26                 nodeinsideicefront=ContourToMesh(md.mesh.elements,md.mesh.x,md.mesh.y,icefrontfile,'node',2)
    27                 nodeonicefront=double(md.mesh.vertexonboundary and nodeinsideicefront)
     27                [nodeinsideicefront,dum]=ContourToMesh(md.mesh.elements,md.mesh.x.reshape(-1,1),md.mesh.y.reshape(-1,1),icefrontfile,'node',2)
     28                nodeonicefront=numpy.logical_and(md.mesh.vertexonboundary,nodeinsideicefront.reshape(-1)).astype(float)
    2829        else:
    29                 nodeonicefront=numpy.zeros(md.mesh.numberofvertices)
     30                nodeonicefront=numpy.zeros((md.mesh.numberofvertices))
    3031
    3132#       pos=find(md.mesh.vertexonboundary & ~nodeonicefront);
    32         pos=[i for i,(vob,noif) in enumerate(zip(md.mesh.vertexonboundary,nodeonicefront)) if vob and not noif]
    33         md.diagnostic.spcvx=float('NaN')*numpy.ones(md.mesh.numberofvertices)
    34         md.diagnostic.spcvy=float('NaN')*numpy.ones(md.mesh.numberofvertices)
    35         md.diagnostic.spcvz=float('NaN')*numpy.ones(md.mesh.numberofvertices)
     33        pos=numpy.nonzero(numpy.logical_and(md.mesh.vertexonboundary,numpy.logical_not(nodeonicefront)))[0]
     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))
    3637        md.diagnostic.spcvx[pos]=0
    3738        md.diagnostic.spcvy[pos]=0
    3839        md.diagnostic.spcvz[pos]=0
    39         md.diagnostic.referential=float('NaN')*numpy.ones((md.mesh.numberofvertices,6))
     40        md.diagnostic.referential=float('nan')*numpy.ones((md.mesh.numberofvertices,6))
    4041
    4142        #Dirichlet Values
    42         if numpy.size(md.inversion.vx_obs)==md.mesh.numberofvertices and numpy.size(md.inversion.vy_obs)==md.mesh.numberofvertices:
    43                 print '      boundary conditions for diagnostic model: spc set as observed velocities'
     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:
     44                #reshape to rank-2 if necessary to match spc arrays
     45                if numpy.ndim(md.inversion.vx_obs)==1:
     46                        md.inversion.vx_obs=md.inversion.vx_obs.reshape(-1,1)
     47                if numpy.ndim(md.inversion.vy_obs)==1:
     48                        md.inversion.vy_obs=md.inversion.vy_obs.reshape(-1,1)
     49                print "      boundary conditions for diagnostic model: spc set as observed velocities"
    4450                md.diagnostic.spcvx[pos]=md.inversion.vx_obs[pos]
    4551                md.diagnostic.spcvy[pos]=md.inversion.vy_obs[pos]
    4652        else:
    47                 print '      boundary conditions for diagnostic model: spc set as zero'
     53                print "      boundary conditions for diagnostic model: spc set as zero"
    4854
    4955        #segment on Ice Front
    5056        #segment on Neumann (Ice Front)
    5157#       pos=find(nodeonicefront(md.mesh.segments(:,1)) | nodeonicefront(md.mesh.segments(:,2)));
    52         pos=[i for i,(noif1,noif2) in enumerate(zip(nodeonicefront[md.mesh.segments[:,0].astype('int')-1],nodeonicefront[md.mesh.segments[:,1].astype('int')-1])) if noif1 or noif2]
     58        pos=numpy.nonzero(numpy.logical_or(nodeonicefront[md.mesh.segments[:,0].astype(int)-1],nodeonicefront[md.mesh.segments[:,1].astype(int)-1]))[0]
    5359        if   md.mesh.dimension==2:
    5460                pressureload=md.mesh.segments[pos,:]
    5561        elif md.mesh.dimension==3:
    5662#               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)];
    57                 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)
     63                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]))
    5864                pressureload=numpy.zeros((0,5))
    5965                for i in xrange(1,md.mesh.numberoflayers):
    6066#                       pressureload=[pressureload ;pressureload_layer1(:,1:4)+(i-1)*md.mesh.numberofvertices2d pressureload_layer1(:,5)+(i-1)*md.mesh.numberofelements2d ];
    61                         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)
     67                        pressureload=numpy.vstack((pressureload,numpy.hstack((pressureload_layer1[:,0:4]+(i-1)*md.mesh.numberofvertices2d,pressureload_layer1[:,4]+(i-1)*md.mesh.numberofelements2d))))
    6268
    6369        #Add water or air enum depending on the element
    6470#       pressureload=[pressureload 1*md.mask.elementonfloatingice(pressureload(:,end))];
    65         pressureload=numpy.concatenate((pressureload,1*md.mask.elementonfloatingice[pressureload[:,-1].astype('int')-1].reshape((-1,1))),axis=1)
     71        pressureload=numpy.hstack((pressureload,1.*md.mask.elementonfloatingice[pressureload[:,-1].astype('int')-1].reshape(-1,1)))
    6672
    6773        #plug onto model
     
    7076        #Create zeros basalforcings and surfaceforcings
    7177        if numpy.all(numpy.isnan(md.surfaceforcings.precipitation)) and (md.surfaceforcings.ispdd==1):
    72                 md.surfaceforcings.precipitation=numpy.zeros(md.mesh.numberofvertices)
    73                 print '      no surfaceforcings.precipitation specified: values set as zero'
     78                md.surfaceforcings.precipitation=numpy.zeros((md.mesh.numberofvertices,1))
     79                print "      no surfaceforcings.precipitation specified: values set as zero"
    7480        if numpy.all(numpy.isnan(md.surfaceforcings.mass_balance)) and (md.surfaceforcings.ispdd==0):
    75                 md.surfaceforcings.mass_balance=numpy.zeros(md.mesh.numberofvertices)
    76                 print '      no surfaceforcings.mass_balance specified: values set as zero'
     81                md.surfaceforcings.mass_balance=numpy.zeros((md.mesh.numberofvertices,1))
     82                print "      no surfaceforcings.mass_balance specified: values set as zero"
    7783        if numpy.all(numpy.isnan(md.basalforcings.melting_rate)):
    78                 md.basalforcings.melting_rate=numpy.zeros(md.mesh.numberofvertices)
    79                 print '      no basalforcings.melting_rate specified: values set as zero'
     84                md.basalforcings.melting_rate=numpy.zeros((md.mesh.numberofvertices,1))
     85                print "      no basalforcings.melting_rate specified: values set as zero"
    8086        if numpy.all(numpy.isnan(md.balancethickness.thickening_rate)):
    81                 md.balancethickness.thickening_rate=numpy.zeros(md.mesh.numberofvertices)
    82                 print '      no balancethickness.thickening_rate specified: values set as zero'
     87                md.balancethickness.thickening_rate=numpy.zeros((md.mesh.numberofvertices,1))
     88                print "      no balancethickness.thickening_rate specified: values set as zero"
    8389
    84         md.prognostic.spcthickness=float('NaN')*numpy.ones(md.mesh.numberofvertices)
    85         md.balancethickness.spcthickness=float('NaN')*numpy.ones(md.mesh.numberofvertices)
     90        md.prognostic.spcthickness=float('nan')*numpy.ones((md.mesh.numberofvertices,1))
     91        md.balancethickness.spcthickness=float('nan')*numpy.ones((md.mesh.numberofvertices,1))
    8692
    87         if numpy.size(md.initialization.temperature)==md.mesh.numberofvertices:
    88                 md.thermal.spctemperature=float('NaN')*numpy.ones(md.mesh.numberofvertices)
     93        if isinstance(md.initialization.temperature,numpy.ndarray) and numpy.size(md.initialization.temperature,axis=0)==md.mesh.numberofvertices:
     94                md.thermal.spctemperature=float('nan')*numpy.ones((md.mesh.numberofvertices,1))
    8995#               pos=find(md.mesh.vertexonsurface); md.thermal.spctemperature(pos)=md.initialization.temperature(pos); %impose observed temperature on surface
    90                 pos=[i for i,vos in enumerate(md.mesh.vertexonsurface) if vos]
    91                 md.thermal.spctemperature[pos]=md.initialization.temperature[pos]    # impose observed temperature on surface
    92                 if not numpy.size(md.basalforcings.geothermalflux)==md.mesh.numberofvertices:
    93                         md.basalforcings.geothermalflux=numpy.zeros(md.mesh.numberofvertices)
     96                pos=numpy.nonzero(md.mesh.vertexonsurface)[0]
     97                md.thermal.spctemperature[pos]=md.initialization.temperature[pos]    #impose observed temperature on surface
     98                if not isinstance(md.basalforcings.geothermalflux,numpy.ndarray) or not numpy.size(md.basalforcings.geothermalflux,axis=0)==md.mesh.numberofvertices:
     99                        md.basalforcings.geothermalflux=numpy.zeros((md.mesh.numberofvertices,1))
    94100        else:
    95                 print '      no thermal boundary conditions created: no observed temperature found'
     101                print "      no thermal boundary conditions created: no observed temperature found"
    96102
    97103        return md
Note: See TracChangeset for help on using the changeset viewer.