Changeset 13975 for issm/trunk/src/m/boundaryconditions/SetIceShelfBC.py
- Timestamp:
- 11/16/12 08:10:16 (12 years ago)
- Location:
- issm/trunk
- Files:
-
- 3 edited
Legend:
- Unmodified
- Added
- Removed
-
issm/trunk
- Property svn:mergeinfo changed
/issm/trunk-jpl merged: 13397-13398,13401,13407-13582,13584-13974
- Property svn:mergeinfo changed
-
issm/trunk/src
- Property svn:mergeinfo changed
-
issm/trunk/src/m/boundaryconditions/SetIceShelfBC.py
r13395 r13975 1 1 import os 2 2 import numpy 3 from ContourToMesh import * 3 4 4 5 def SetIceShelfBC(md,icefrontfile=''): … … 6 7 SETICESHELFBC - Create the boundary conditions for diagnostic and thermal models for a Ice Shelf with Ice Front 7 8 8 Neumann BC are used on the ice front (an A NRGUS contour around the ice front9 Neumann BC are used on the ice front (an ARGUS contour around the ice front 9 10 must be given in input) 10 11 Dirichlet BC are used elsewhere for diagnostic … … 24 25 if not os.path.exists(icefrontfile): 25 26 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) 28 29 else: 29 nodeonicefront=numpy.zeros( md.mesh.numberofvertices)30 nodeonicefront=numpy.zeros((md.mesh.numberofvertices)) 30 31 31 32 # 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)) 36 37 md.diagnostic.spcvx[pos]=0 37 38 md.diagnostic.spcvy[pos]=0 38 39 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)) 40 41 41 42 #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" 44 50 md.diagnostic.spcvx[pos]=md.inversion.vx_obs[pos] 45 51 md.diagnostic.spcvy[pos]=md.inversion.vy_obs[pos] 46 52 else: 47 print ' boundary conditions for diagnostic model: spc set as zero'53 print " boundary conditions for diagnostic model: spc set as zero" 48 54 49 55 #segment on Ice Front 50 56 #segment on Neumann (Ice Front) 51 57 # 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] 53 59 if md.mesh.dimension==2: 54 60 pressureload=md.mesh.segments[pos,:] 55 61 elif md.mesh.dimension==3: 56 62 # 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])) 58 64 pressureload=numpy.zeros((0,5)) 59 65 for i in xrange(1,md.mesh.numberoflayers): 60 66 # 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)))) 62 68 63 69 #Add water or air enum depending on the element 64 70 # 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))) 66 72 67 73 #plug onto model … … 70 76 #Create zeros basalforcings and surfaceforcings 71 77 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" 74 80 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" 77 83 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" 80 86 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" 83 89 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)) 86 92 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)) 89 95 # 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] # 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)) 94 100 else: 95 print ' no thermal boundary conditions created: no observed temperature found'101 print " no thermal boundary conditions created: no observed temperature found" 96 102 97 103 return md
Note:
See TracChangeset
for help on using the changeset viewer.