Index: /u/astrid-r1b/morlighe/issmuci/trunk-jpl/../trunk-jpl/src/m/utils/BC/SetIceShelfBC.py =================================================================== --- /u/astrid-r1b/morlighe/issmuci/trunk-jpl/../trunk-jpl/src/m/utils/BC/SetIceShelfBC.py (revision 12853) +++ /u/astrid-r1b/morlighe/issmuci/trunk-jpl/../trunk-jpl/src/m/utils/BC/SetIceShelfBC.py (revision 12854) @@ -1,91 +1,98 @@ -from numpy import * +import os +import numpy + def SetIceShelfBC(md,icefrontfile=''): - #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 - # must be given in input) - # Dirichlet BC are used elsewhere for diagnostic - # - # Usage: - # md=SetIceShelfBC(md,varargin) - # - # Example: - # md=SetIceShelfBC(md); - # md=SetIceShelfBC(md,'Front.exp'); - # - # See also: SETICESHEETBC, SETMARINEICESHEETBC + """ + 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 + must be given in input) + Dirichlet BC are used elsewhere for diagnostic + + Usage: + md=SetIceShelfBC(md,varargin) + + Example: + md=SetIceShelfBC(md); + md=SetIceShelfBC(md,'Front.exp'); + + See also: SETICESHEETBC, SETMARINEICESHEETBC + """ #node on Dirichlet (boundary and ~icefront) - if not icefrontfile: - nodeonicefront=zeros(md.mesh.numberofvertices) + if icefrontfile: + if not os.path.exists(icefrontfile): + raise IOError("SetIceShelfBC error message: ice front file '%s' not found." % icefrontfile) + nodeinsideicefront=ContourToMesh(md.mesh.elements,md.mesh.x,md.mesh.y,icefrontfile,'node',2) + nodeonicefront=double(md.mesh.vertexonboundary and nodeinsideicefront) else: - if not os.path.isfile(icefrontfile): - print 'SetIceShelfBC error message: ice front file '+icefrontfile+ ' not found' - return [] - nodeinsideicefront=ContourToMesh(md.mesh.elements,md.mesh.x,md.mesh.y,icefrontfile,'node',2) - nodeonicefront=double(md.mesh.vertexonboundary.astype(bool) and nodeinsideicefront.astype(bool)) - - pos=argwhere(logical_and(md.mesh.vertexonboundary.astype(bool),~nodeonicefront.astype(bool))) - md.diagnostic.spcvx=NaN*ones(md.mesh.numberofvertices); - md.diagnostic.spcvy=NaN*ones(md.mesh.numberofvertices); - md.diagnostic.spcvz=NaN*ones(md.mesh.numberofvertices); - md.diagnostic.spcvx[pos]=0; - md.diagnostic.spcvy[pos]=0; - md.diagnostic.spcvz[pos]=0; - md.diagnostic.referential=float('NaN')*ones((md.mesh.numberofvertices,6),float); + nodeonicefront=numpy.zeros(md.mesh.numberofvertices) +# 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[pos]=0 + md.diagnostic.spcvy[pos]=0 + md.diagnostic.spcvz[pos]=0 + md.diagnostic.referential=float('NaN')*numpy.ones((md.mesh.numberofvertices,6)) + #Dirichlet Values - if ~isnan(md.inversion.vx_obs) and ~isnan(md.inversion.vy_obs): - if (len(md.inversion.vx_obs)==md.mesh.numberofvertices) and (len(md.inversion.vy_obs)==md.mesh.numberofvertices): - print ' boundary conditions for diagnostic model: spc set as observed velocities' - md.diagnostic.spcvx[pos]=md.inversion.vx_obs[pos] - md.diagnostic.spcvy[pos]=md.inversion.vy_obs[pos] + if numpy.size(md.inversion.vx_obs)==md.mesh.numberofvertices and numpy.size(md.inversion.vy_obs)==md.mesh.numberofvertices: + print ' boundary conditions for diagnostic model: spc set as observed velocities' + md.diagnostic.spcvx[pos]=md.inversion.vx_obs[pos] + md.diagnostic.spcvy[pos]=md.inversion.vy_obs[pos] else: print ' boundary conditions for diagnostic model: spc set as zero' #segment on Ice Front #segment on Neumann (Ice Front) - segs1=md.mesh.segments[:,0].astype(int)-1 - segs2=md.mesh.segments[:,1].astype(int)-1 - - pos=argwhere(logical_and(nodeonicefront[segs1].astype(bool),nodeonicefront[segs2].astype(bool))) - if (md.mesh.dimension==2): - pressureload=md.mesh.segments[pos,:]; +# pos=find(nodeonicefront(md.mesh.segments(:,1)) | nodeonicefront(md.mesh.segments(:,2))); + 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] + if md.mesh.dimension==2: + pressureload=md.mesh.segments[pos,:] elif md.mesh.dimension==3: - pressureload_layer=[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=[]; - for i in range(1,md.mesh.numberoflayers-1): - pressureload=[pressureload ,pressureload_layer1[:,1:4]+(i-1)*md.mesh.numberofvertices2d, pressureload_layer1[:,5]+(i-1)*md.mesh.numberofelements2d ]; +# 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=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) #Add water or air enum depending on the element - pressureload=[pressureload, 1*md.mask.elementonfloatingice(pressureload[:,end])]; +# 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) #plug onto model - md.diagnostic.icefront=pressureload; + md.diagnostic.icefront=pressureload #Create zeros basalforcings and surfaceforcings - if isnan(md.surfaceforcings.precipitation): - md.surfaceforcings.precipitation=zeros(md.mesh.numberofvertices); + if numpy.isnan(md.surfaceforcings.precipitation).all(): + md.surfaceforcings.precipitation=numpy.zeros(md.mesh.numberofvertices) print ' no surfaceforcings.precipitation specified: values set as zero' - if isnan(md.surfaceforcings.mass_balance): - md.surfaceforcings.mass_balance=zeros(md.mesh.numberofvertices); + if numpy.isnan(md.surfaceforcings.mass_balance).all(): + md.surfaceforcings.mass_balance=numpy.zeros(md.mesh.numberofvertices) print ' no surfaceforcings.mass_balance specified: values set as zero' - if isnan(md.basalforcings.melting_rate): - md.basalforcings.melting_rate=zeros(md.mesh.numberofvertices) + if numpy.isnan(md.basalforcings.melting_rate).all(): + md.basalforcings.melting_rate=numpy.zeros(md.mesh.numberofvertices) print ' no basalforcings.melting_rate specified: values set as zero' - if isnan(md.balancethickness.thickening_rate): - md.balancethickness.thickening_rate=zeros(md.mesh.numberofvertices); + if numpy.isnan(md.balancethickness.thickening_rate).all(): + md.balancethickness.thickening_rate=numpy.zeros(md.mesh.numberofvertices) print ' no balancethickness.thickening_rate specified: values set as zero' - md.prognostic.spcthickness=NaN*ones(md.mesh.numberofvertices); - md.balancethickness.spcthickness=NaN*ones(md.mesh.numberofvertices); + md.prognostic.spcthickness=float('NaN')*numpy.ones(md.mesh.numberofvertices) + md.balancethickness.spcthickness=float('NaN')*numpy.ones(md.mesh.numberofvertices) - if (len(md.initialization.temperature)==md.mesh.numberofvertices): - md.thermal.spctemperature=NaN*ones(md.mesh.numberofvertices); - pos=argwhere(md.mesh.vertexonsurface); md.thermal.spctemperature[pos]=md.initialization.temperature[pos]; #impose observed temperature on surface - if (len(md.basalforcings.geothermalflux) !=md.mesh.numberofvertices): - md.basalforcings.geothermalflux=zeros(md.mesh.numberofvertices); + if numpy.size(md.initialization.temperature)==md.mesh.numberofvertices: + md.thermal.spctemperature=float('NaN')*numpy.ones(md.mesh.numberofvertices) +# 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) else: print ' no thermal boundary conditions created: no observed temperature found' return md +