Index: ../trunk-jpl/src/m/materials/paterson.py =================================================================== --- ../trunk-jpl/src/m/materials/paterson.py (revision 13461) +++ ../trunk-jpl/src/m/materials/paterson.py (revision 13462) @@ -1,24 +1,23 @@ -from numpy import * +import numpy def paterson(temperature): + """ + PATERSON - figure out the rigidity of ice for a given temperature - # Local Variables: pos11, pos5, pos10, temperature, pos, T, pos8, pos9, pos6, pos7, pos4, rigidity, pos2, pos3, pos1 - # Function calls: length, zeros, argwhere, paterson, error - #PATERSON - figure out the rigidity of ice for a given temperature - # - # rigidigty (in s^(1/3)Pa) is the flow law paramter in the flow law sigma=B*e(1/3) (Paterson, p97). - # temperature is in Kelvin degrees - # - # Usage: - # rigidity=paterson(temperature) - - pos=argwhere(temperature<0.) - if len(pos): - print 'input temperature should be in Kelvin (positive)' - return [] - + rigidity (in s^(1/3)Pa) is the flow law paramter in the flow law sigma=B*e(1/3) (Paterson, p97). + temperature is in Kelvin degrees + + Usage: + rigidity=paterson(temperature) + """ + + if numpy.any(temperature<0.): + raise RuntimeError("input temperature should be in Kelvin (positive)") + T = temperature-273.15 + #The routine below is equivalent to: + # n=3; T=temperature-273; # %From paterson, # Temp=[0;-2;-5;-10;-15;-20;-25;-30;-35;-40;-45;-50]; @@ -30,21 +29,33 @@ # fittedmodel=fit(Temp,B,'cubicspline'); # rigidity=fittedmodel(temperature); - rigidity=zeros(len(T)) - 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) - 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) - 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) - 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) - 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) - 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) - 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) - 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) - 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) - 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) - 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) + rigidity=numpy.zeros((numpy.size(T,axis=0),1)) + pos1=numpy.nonzero(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) + pos2=numpy.nonzero(numpy.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) + pos3=numpy.nonzero(numpy.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) + pos4=numpy.nonzero(numpy.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) + pos5=numpy.nonzero(numpy.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) + pos6=numpy.nonzero(numpy.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) + pos7=numpy.nonzero(numpy.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) + pos8=numpy.nonzero(numpy.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) + pos9=numpy.nonzero(numpy.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) + pos10=numpy.nonzero(numpy.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) + pos11=numpy.nonzero(-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) #Now make sure that rigidity is positive - pos=argwhere(rigidity<0); rigidity[pos]=1**6 + pos=numpy.nonzero(rigidity<0) + rigidity[pos]=1**6 return rigidity Index: ../trunk-jpl/src/m/materials/paterson.m =================================================================== --- ../trunk-jpl/src/m/materials/paterson.m (revision 13461) +++ ../trunk-jpl/src/m/materials/paterson.m (revision 13462) @@ -1,7 +1,7 @@ function rigidity=paterson(temperature) %PATERSON - figure out the rigidity of ice for a given temperature % -% rigidigty (in s^(1/3)Pa) is the flow law paramter in the flow law sigma=B*e(1/3) (Paterson, p97). +% rigidity (in s^(1/3)Pa) is the flow law paramter in the flow law sigma=B*e(1/3) (Paterson, p97). % temperature is in Kelvin degrees % % Usage: Index: ../trunk-jpl/src/m/solve/parseresultsfromdisk.py =================================================================== --- ../trunk-jpl/src/m/solve/parseresultsfromdisk.py (revision 13461) +++ ../trunk-jpl/src/m/solve/parseresultsfromdisk.py (revision 13462) @@ -51,7 +51,7 @@ if result['step'] in results and \ result['fieldname'] in results[result['step']] and \ not strcmp(result['fieldname'],'SolutionType'): - results[result['step']][result['fieldname']]=numpy.concatenate((results[result['step']][result['fieldname']],result['field']),axis=0) + results[result['step']][result['fieldname']]=numpy.vstack((results[result['step']][result['fieldname']],result['field'])) else: results[result['step']][result['fieldname']]=result['field'] Index: ../trunk-jpl/src/m/classes/thermal.py =================================================================== --- ../trunk-jpl/src/m/classes/thermal.py (revision 13461) +++ ../trunk-jpl/src/m/classes/thermal.py (revision 13462) @@ -71,11 +71,8 @@ md = checkfield(md,'thermal.spctemperature','forcing',1) if EnthalpyAnalysisEnum() in analyses and (md.thermal.isenthalpy or solution==EnthalpySolutionEnum()) and md.mesh.dimension==3: pos=numpy.nonzero(numpy.logical_not(numpy.isnan(md.thermal.spctemperature[0:md.mesh.numberofvertices]))) -# need to know if md.thermal.spctemperature can have multiple columns? -# replicate=numpy.tile(md.geometry.surface-md.mesh.z,(1,numpy.size(md.thermal.spctemperature,axis=1))) -# 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") - replicate=numpy.tile(md.geometry.surface-md.mesh.z,(1)) - 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") + replicate=numpy.tile(md.geometry.surface-md.mesh.z,(1,numpy.size(md.thermal.spctemperature,axis=1))) + 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") md = checkfield(md,'thermal.isenthalpy','numel',[1],'values',[0,1]) return md Index: ../trunk-jpl/src/m/classes/model/model.py =================================================================== --- ../trunk-jpl/src/m/classes/model/model.py (revision 13461) +++ ../trunk-jpl/src/m/classes/model/model.py (revision 13462) @@ -241,7 +241,7 @@ x3d=numpy.concatenate((x3d,md.mesh.x)) y3d=numpy.concatenate((y3d,md.mesh.y)) #nodes are distributed between bed and surface accordingly to the given exponent - z3d=numpy.concatenate((z3d,bed3d+thickness3d*extrusionlist[i])) + z3d=numpy.concatenate((z3d,(bed3d+thickness3d*extrusionlist[i]).reshape(-1))) number_nodes3d=numpy.size(x3d) #number of 3d nodes for the non extruded part of the mesh #Extrude elements @@ -401,7 +401,7 @@ #Put lithostatic pressure if there is an existing pressure if not numpy.any(numpy.isnan(md.initialization.pressure)): - md.initialization.pressure=md.constants.g*md.materials.rho_ice*(md.geometry.surface-md.mesh.z) + md.initialization.pressure=md.constants.g*md.materials.rho_ice*(md.geometry.surface-md.mesh.z.reshape(-1,1)) #special for thermal modeling: md.basalforcings.melting_rate=project3d(md,'vector',md.basalforcings.melting_rate,'type','node','layer',1) Index: ../trunk-jpl/src/m/boundaryconditions/SetIceShelfBC.py =================================================================== --- ../trunk-jpl/src/m/boundaryconditions/SetIceShelfBC.py (revision 13461) +++ ../trunk-jpl/src/m/boundaryconditions/SetIceShelfBC.py (revision 13462) @@ -1,12 +1,12 @@ import os import numpy -from ContourToMesh import ContourToMesh +from ContourToMesh import * 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 + Neumann BC are used on the ice front (an ARGUS contour around the ice front must be given in input) Dirichlet BC are used elsewhere for diagnostic @@ -24,23 +24,23 @@ 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).reshape(-1,1),(md.mesh.x).reshape(-1,1),(md.mesh.y).reshape(-1,1),icefrontfile,'node',2) + nodeinsideicefront=ContourToMesh(md.mesh.elements,md.mesh.x.reshape(-1,1),md.mesh.y.reshape(-1,1),icefrontfile,'node',2) nodeonicefront= numpy.bitwise_and( map(int,md.mesh.vertexonboundary), map(int,nodeinsideicefront[0].ravel()) ) else: - nodeonicefront=numpy.zeros(md.mesh.numberofvertices) + nodeonicefront=numpy.zeros((md.mesh.numberofvertices,1)) # 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=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 numpy.size(md.inversion.vx_obs)==md.mesh.numberofvertices and numpy.size(md.inversion.vy_obs)==md.mesh.numberofvertices: + 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] @@ -55,43 +55,43 @@ 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.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_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.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) + pressureload=numpy.vstack((pressureload,numpy.hstack((pressureload_layer1[:,0:3]+(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.concatenate((pressureload,1*md.mask.elementonfloatingice[pressureload[:,-1].astype('int')-1].reshape((-1,1))),axis=1) + pressureload=numpy.hstack((pressureload,1*md.mask.elementonfloatingice[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) + 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) + 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) + 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) + 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) - md.balancethickness.spcthickness=float('NaN')*numpy.ones(md.mesh.numberofvertices) + 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) + 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: - md.basalforcings.geothermalflux=numpy.zeros(md.mesh.numberofvertices) + md.basalforcings.geothermalflux=numpy.zeros((md.mesh.numberofvertices,1)) else: print ' no thermal boundary conditions created: no observed temperature found' Index: ../trunk-jpl/src/m/parameterization/setflowequation.py =================================================================== --- ../trunk-jpl/src/m/parameterization/setflowequation.py (revision 13461) +++ ../trunk-jpl/src/m/parameterization/setflowequation.py (revision 13462) @@ -93,10 +93,10 @@ fullspcnodes=numpy.logical_or(numpy.logical_not(numpy.isnan(md.diagnostic.spcvx)).astype(int)+ \ numpy.logical_not(numpy.isnan(md.diagnostic.spcvy)).astype(int)+ \ numpy.logical_not(numpy.isnan(md.diagnostic.spcvz)).astype(int)==3, \ - numpy.logical_and(nodeonpattyn,nodeonstokes)).astype(int) #find all the nodes on the boundary of the domain without icefront + numpy.logical_and(nodeonpattyn,nodeonstokes).reshape(-1,1)).astype(int) #find all the nodes on the boundary of the domain without icefront # fullspcelems=double(sum(fullspcnodes(md.mesh.elements),2)==6); %find all the nodes on the boundary of the domain without icefront 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 - stokesflag[numpy.nonzero(fullspcelems)]=0 + stokesflag[numpy.nonzero(fullspcelems.reshape(-1))]=0 nodeonstokes[md.mesh.elements[numpy.nonzero(stokesflag),:].astype(int)-1]=1 #Then complete with NoneApproximation or the other model used if there is no stokes @@ -128,7 +128,7 @@ if numpy.all(numpy.logical_not(numpy.isnan(bordernodes2d))): penalties=numpy.zeros((0,2)) for i in xrange(1,numlayers): - penalties=numpy.concatenate((penalties,numpy.concatenate((bordernodes2d,bordernodes2d+md.mesh.numberofvertices2d*(i)),axis=1)),axis=0) + penalties=numpy.vstack((penalties,numpy.hstack((bordernodes2d,bordernodes2d+md.mesh.numberofvertices2d*(i))))) md.diagnostic.vertex_pairing=penalties elif strcmpi(coupling_method,'tiling'):