Index: ../trunk-jpl/src/m/geometry/FlagElements.py =================================================================== --- ../trunk-jpl/src/m/geometry/FlagElements.py (revision 13469) +++ ../trunk-jpl/src/m/geometry/FlagElements.py (revision 13470) @@ -40,13 +40,13 @@ if not os.path.exists(region): if len(region)>3 and not strcmp(region[-4:],'.exp'): raise IOError("Error: File 'region' not found!" % region) - raise RuntimeError("FlagElements -- basinzoom not yet converted.") + raise RuntimeError("FlagElements.py calling basinzoom.py is not complete.") xlim,ylim=basinzoom('basin',region) flag_nodes=numpy.logical_and(numpy.logical_and(md.mesh.xxlim[0]),numpy.logical_and(md.mesh.yylim[0])).astype(float) flag=numpy.prod(flag_nodes[md.mesh.elements],axis=1) else: #ok, flag elements - [flag,fnone]=ContourToMesh(md.mesh.elements[:,0:3],md.mesh.x.reshape(-1,1),md.mesh.y.reshape(-1,1),region,'element',1) + [flag,dum]=ContourToMesh(md.mesh.elements[:,0:3],md.mesh.x.reshape(-1,1),md.mesh.y.reshape(-1,1),region,'element',1) if invert: flag=numpy.logical_not(flag) Index: ../trunk-jpl/src/m/mesh/bamg.py =================================================================== --- ../trunk-jpl/src/m/mesh/bamg.py (revision 13469) +++ ../trunk-jpl/src/m/mesh/bamg.py (revision 13470) @@ -120,7 +120,7 @@ raise RuntimeError("one rift has all its points outside of the domain outline") elif numpy.any(numpy.logical_not(flags)): - raise RuntimeError("bamg.m for rifts is not complete.") + raise RuntimeError("bamg.py for rifts is not complete.") #We LOTS of work to do print("Rift tip outside of or on the domain has been detected and is being processed...") """ @@ -346,7 +346,7 @@ def processgeometry(geom,tol,outline): # {{{ - raise RuntimeError("bamg.m/processgeometry is not complete.") + raise RuntimeError("bamg.py/processgeometry is not complete.") #Deal with edges print("Checking Edge crossing...") i=0 Index: ../trunk-jpl/src/m/classes/model/model.py =================================================================== --- ../trunk-jpl/src/m/classes/model/model.py (revision 13469) +++ ../trunk-jpl/src/m/classes/model/model.py (revision 13470) @@ -166,6 +166,7 @@ def checkmessage(self,string): # {{{ print ("model not consistent: %s" % string) self.private.isconsistent=False + return self # }}} def extrude(md,*args): # {{{ Index: ../trunk-jpl/src/m/boundaryconditions/SetMarineIceSheetBC.m =================================================================== --- ../trunk-jpl/src/m/boundaryconditions/SetMarineIceSheetBC.m (revision 13469) +++ ../trunk-jpl/src/m/boundaryconditions/SetMarineIceSheetBC.m (revision 13470) @@ -52,7 +52,7 @@ end md.hydrology.spcwatercolumn=zeros(md.mesh.numberofvertices,2); -pos=find(md.mesh.vertexonboundary); +pos=find(md.mesh.vertexonboundary); md.hydrology.spcwatercolumn(pos,1)=1; %segment on Neumann (Ice Front) @@ -100,7 +100,7 @@ pos=find(md.mesh.vertexonsurface); md.thermal.spctemperature(pos)=md.initialization.temperature(pos); %impose observed temperature on surface if (length(md.basalforcings.geothermalflux)~=md.mesh.numberofvertices), md.basalforcings.geothermalflux=zeros(md.mesh.numberofvertices,1); - md.basalforcings.geothermalflux(find(md.mask.vertexongroundedice))=50*10^-3; %50mW/m2 + md.basalforcings.geothermalflux(find(md.mask.vertexongroundedice))=50.*10.^-3; %50mW/m2 end else disp(' no thermal boundary conditions created: no observed temperature found'); Index: ../trunk-jpl/src/m/boundaryconditions/SetIceShelfBC.py =================================================================== --- ../trunk-jpl/src/m/boundaryconditions/SetIceShelfBC.py (revision 13469) +++ ../trunk-jpl/src/m/boundaryconditions/SetIceShelfBC.py (revision 13470) @@ -25,32 +25,32 @@ if not os.path.exists(icefrontfile): raise IOError("SetIceShelfBC error message: ice front file '%s' not found." % icefrontfile) [nodeinsideicefront,dum]=ContourToMesh(md.mesh.elements,md.mesh.x.reshape(-1,1),md.mesh.y.reshape(-1,1),icefrontfile,'node',2) - nodeonicefront=numpy.logical_and(md.mesh.vertexonboundary,nodeinsideicefront).astype(float) + nodeonicefront=numpy.logical_and(md.mesh.vertexonboundary,nodeinsideicefront.reshape(-1)).astype(float) else: - nodeonicefront=numpy.zeros((md.mesh.numberofvertices,1)) + 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,1)) - md.diagnostic.spcvy=float('NaN')*numpy.ones((md.mesh.numberofvertices,1)) - md.diagnostic.spcvz=float('NaN')*numpy.ones((md.mesh.numberofvertices,1)) + pos=numpy.nonzero(numpy.logical_and(md.mesh.vertexonboundary,numpy.logical_not(nodeonicefront)))[0] + 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 md.diagnostic.spcvz[pos]=0 - md.diagnostic.referential=float('NaN')*numpy.ones((md.mesh.numberofvertices,6)) + md.diagnostic.referential=float('nan')*numpy.ones((md.mesh.numberofvertices,6)) #Dirichlet Values 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' + 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' + print " boundary conditions for diagnostic model: spc set as zero" #segment on Ice Front #segment on Neumann (Ice Front) # 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] + pos=numpy.nonzero(numpy.logical_or(nodeonicefront[md.mesh.segments[:,0].astype(int)-1],nodeonicefront[md.mesh.segments[:,1].astype(int)-1]))[0] if md.mesh.dimension==2: pressureload=md.mesh.segments[pos,:] elif md.mesh.dimension==3: @@ -59,11 +59,11 @@ 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.vstack((pressureload,numpy.hstack((pressureload_layer1[:,0:3]+(i-1)*md.mesh.numberofvertices2d,pressureload_layer1[:,4]+(i-1)*md.mesh.numberofelements2d)))) + pressureload=numpy.vstack((pressureload,numpy.hstack((pressureload_layer1[:,0:4]+(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.hstack((pressureload,1*md.mask.elementonfloatingice[pressureload[:,-1].astype('int')-1].reshape((-1,1)))) + pressureload=numpy.hstack((pressureload,1.*md.mask.elementonfloatingice[pressureload[:,-1].astype('int')-1].reshape(-1,1))) #plug onto model md.diagnostic.icefront=pressureload @@ -71,29 +71,29 @@ #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,1)) - print ' no surfaceforcings.precipitation specified: values set as zero' + 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,1)) - print ' no surfaceforcings.mass_balance specified: values set as zero' + 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,1)) - print ' no basalforcings.melting_rate specified: values set as zero' + 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,1)) - print ' no balancethickness.thickening_rate specified: values set as zero' + print " no balancethickness.thickening_rate specified: values set as zero" - md.prognostic.spcthickness=float('NaN')*numpy.ones((md.mesh.numberofvertices,1)) - md.balancethickness.spcthickness=float('NaN')*numpy.ones((md.mesh.numberofvertices,1)) + 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,1)) + if isinstance(md.initialization.temperature,numpy.ndarray) and numpy.size(md.initialization.temperature,axis=0)==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: + pos=numpy.nonzero(md.mesh.vertexonsurface)[0] + md.thermal.spctemperature[pos]=md.initialization.temperature[pos] #impose observed temperature on surface + if not isinstance(md.basalforcings.geothermalflux,numpy.ndarray) or not numpy.size(md.basalforcings.geothermalflux,axis=0)==md.mesh.numberofvertices: md.basalforcings.geothermalflux=numpy.zeros((md.mesh.numberofvertices,1)) else: - print ' no thermal boundary conditions created: no observed temperature found' + print " no thermal boundary conditions created: no observed temperature found" return md Index: ../trunk-jpl/src/m/boundaryconditions/SetMarineIceSheetBC.py =================================================================== --- ../trunk-jpl/src/m/boundaryconditions/SetMarineIceSheetBC.py (revision 0) +++ ../trunk-jpl/src/m/boundaryconditions/SetMarineIceSheetBC.py (revision 13470) @@ -0,0 +1,111 @@ +import os +import numpy +from ContourToMesh import * + +def SetMarineIceSheetBC(md,icefrontfile=''): + """ + SETICEMARINESHEETBC - Create the boundary conditions for diagnostic and thermal models for a Marine Ice Sheet with Ice Front + + Neumann BC are used on the ice front (an ARGUS contour around the ice front + can be given in input, or it will be deduced as onfloatingice & onboundary) + Dirichlet BC are used elsewhere for diagnostic + + Usage: + md=SetMarineIceSheetBC(md,icefrontfile) + md=SetMarineIceSheetBC(md) + + Example: + md=SetMarineIceSheetBC(md,'Front.exp') + md=SetMarineIceSheetBC(md) + + See also: SETICESHELFBC, SETMARINEICESHEETBC + """ + + #node on Dirichlet (boundary and ~icefront) + if icefrontfile: + #User provided Front.exp, use it + if not os.path.exists(icefrontfile): + raise IOError("SetMarineIceSheetBC error message: ice front file '%s' not found." % icefrontfile) + [nodeinsideicefront,dum]=ContourToMesh(md.mesh.elements,md.mesh.x.reshape(-1,1),md.mesh.y.reshape(-1,1),icefrontfile,'node',2) + vertexonicefront=numpy.logical_and(md.mesh.vertexonboundary,nodeinsideicefront.reshape(-1)).astype(float) + else: + #Guess where the ice front is + vertexonfloatingice=numpy.zeros((md.mesh.numberofvertices)) + vertexonfloatingice[md.mesh.elements[numpy.nonzero(md.mask.elementonfloatingice),:].astype(int)-1]=1 + vertexonicefront=numpy.logical_and(md.mesh.vertexonboundary,vertexonfloatingice).astype(float) + +# pos=find(md.mesh.vertexonboundary & ~vertexonicefront); + pos=numpy.nonzero(numpy.logical_and(md.mesh.vertexonboundary,numpy.logical_not(vertexonicefront)))[0] + if not numpy.size(pos): + print("SetMarineIceSheetBC warning: ice front all around the glacier, no dirichlet found. Dirichlet must be added manually.") + + 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 + md.diagnostic.spcvz[pos]=0 + md.diagnostic.referential=float('nan')*numpy.ones((md.mesh.numberofvertices,6)) + + #Dirichlet Values + 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] + md.diagnostic.spcvy[pos]=md.inversion.vy_obs[pos] + else: + print(" boundary conditions for diagnostic model: spc set as zero") + + md.hydrology.spcwatercolumn=numpy.zeros((md.mesh.numberofvertices,2)) + pos=numpy.nonzero(md.mesh.vertexonboundary)[0] + md.hydrology.spcwatercolumn[pos,0]=1 + + #segment on Neumann (Ice Front) +# pos=find(vertexonicefront(md.mesh.segments(:,1)) | vertexonicefront(md.mesh.segments(:,2))); + pos=numpy.nonzero(numpy.logical_or(vertexonicefront[md.mesh.segments[:,0].astype(int)-1],vertexonicefront[md.mesh.segments[:,1].astype(int)-1]))[0] + if md.mesh.dimension==2: + pressureload=md.mesh.segments[pos,:] + 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.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.vstack((pressureload,numpy.hstack((pressureload_layer1[:,0:4]+(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))+ 0*md.mask.elementongroundedice(pressureload(:,end))]; + pressureload=numpy.hstack((pressureload,1.*md.mask.elementonfloatingice[pressureload[:,-1].astype('int')-1].reshape(-1,1)+0.*md.mask.elementongroundedice[pressureload[:,-1].astype('int')-1].reshape(-1,1))) + + #plug onto model + md.diagnostic.icefront=pressureload + + #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,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,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,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,1)) + print(" no balancethickness.thickening_rate specified: values set as zero") + + md.prognostic.spcthickness=float('nan')*numpy.ones((md.mesh.numberofvertices,1)) + md.balancethickness.spcthickness=float('nan')*numpy.ones((md.mesh.numberofvertices,1)) + + if isinstance(md.initialization.temperature,numpy.ndarray) and numpy.size(md.initialization.temperature,axis=0)==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=numpy.nonzero(md.mesh.vertexonsurface)[0] + md.thermal.spctemperature[pos]=md.initialization.temperature[pos] #impose observed temperature on surface + if not isinstance(md.basalforcings.geothermalflux,numpy.ndarray) or not numpy.size(md.basalforcings.geothermalflux,axis=0)==md.mesh.numberofvertices: + md.basalforcings.geothermalflux=numpy.zeros((md.mesh.numberofvertices,1)) + md.basalforcings.geothermalflux[numpy.nonzero(md.mask.vertexongroundedice)]=50.*10.**-3; #50mW/m2 + else: + print(" no thermal boundary conditions created: no observed temperature found"); + + return md +