Changeset 13462
- Timestamp:
- 09/27/12 09:44:58 (12 years ago)
- Location:
- issm/trunk-jpl/src/m
- Files:
-
- 7 edited
Legend:
- Unmodified
- Added
- Removed
-
issm/trunk-jpl/src/m/boundaryconditions/SetIceShelfBC.py
r13455 r13462 1 1 import os 2 2 import numpy 3 from ContourToMesh import ContourToMesh3 from ContourToMesh import * 4 4 5 5 def SetIceShelfBC(md,icefrontfile=''): … … 7 7 SETICESHELFBC - Create the boundary conditions for diagnostic and thermal models for a Ice Shelf with Ice Front 8 8 9 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 10 10 must be given in input) 11 11 Dirichlet BC are used elsewhere for diagnostic … … 25 25 if not os.path.exists(icefrontfile): 26 26 raise IOError("SetIceShelfBC error message: ice front file '%s' not found." % icefrontfile) 27 nodeinsideicefront=ContourToMesh( (md.mesh.elements).reshape(-1,1),(md.mesh.x).reshape(-1,1),(md.mesh.y).reshape(-1,1),icefrontfile,'node',2)27 nodeinsideicefront=ContourToMesh(md.mesh.elements,md.mesh.x.reshape(-1,1),md.mesh.y.reshape(-1,1),icefrontfile,'node',2) 28 28 nodeonicefront= numpy.bitwise_and( map(int,md.mesh.vertexonboundary), map(int,nodeinsideicefront[0].ravel()) ) 29 29 else: 30 nodeonicefront=numpy.zeros( md.mesh.numberofvertices)30 nodeonicefront=numpy.zeros((md.mesh.numberofvertices,1)) 31 31 32 32 # pos=find(md.mesh.vertexonboundary & ~nodeonicefront); 33 33 pos=[i for i,(vob,noif) in enumerate(zip(md.mesh.vertexonboundary,nodeonicefront)) if vob and not noif] 34 md.diagnostic.spcvx=float('NaN')*numpy.ones( md.mesh.numberofvertices)35 md.diagnostic.spcvy=float('NaN')*numpy.ones( md.mesh.numberofvertices)36 md.diagnostic.spcvz=float('NaN')*numpy.ones( md.mesh.numberofvertices)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)) 37 37 md.diagnostic.spcvx[pos]=0 38 38 md.diagnostic.spcvy[pos]=0 … … 41 41 42 42 #Dirichlet Values 43 if numpy.size(md.inversion.vx_obs)==md.mesh.numberofvertices and numpy.size(md.inversion.vy_obs)==md.mesh.numberofvertices: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 44 print ' boundary conditions for diagnostic model: spc set as observed velocities' 45 45 md.diagnostic.spcvx[pos]=md.inversion.vx_obs[pos] … … 56 56 elif md.mesh.dimension==3: 57 57 # 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)]; 58 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)58 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])) 59 59 pressureload=numpy.zeros((0,5)) 60 60 for i in xrange(1,md.mesh.numberoflayers): 61 61 # pressureload=[pressureload ;pressureload_layer1(:,1:4)+(i-1)*md.mesh.numberofvertices2d pressureload_layer1(:,5)+(i-1)*md.mesh.numberofelements2d ]; 62 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)62 pressureload=numpy.vstack((pressureload,numpy.hstack((pressureload_layer1[:,0:3]+(i-1)*md.mesh.numberofvertices2d,pressureload_layer1[:,4]+(i-1)*md.mesh.numberofelements2d)))) 63 63 64 64 #Add water or air enum depending on the element 65 65 # pressureload=[pressureload 1*md.mask.elementonfloatingice(pressureload(:,end))]; 66 pressureload=numpy. concatenate((pressureload,1*md.mask.elementonfloatingice[pressureload[:,-1].astype('int')-1].reshape((-1,1))),axis=1)66 pressureload=numpy.hstack((pressureload,1*md.mask.elementonfloatingice[pressureload[:,-1].astype('int')-1].reshape((-1,1)))) 67 67 68 68 #plug onto model … … 71 71 #Create zeros basalforcings and surfaceforcings 72 72 if numpy.all(numpy.isnan(md.surfaceforcings.precipitation)) and (md.surfaceforcings.ispdd==1): 73 md.surfaceforcings.precipitation=numpy.zeros( md.mesh.numberofvertices)73 md.surfaceforcings.precipitation=numpy.zeros((md.mesh.numberofvertices,1)) 74 74 print ' no surfaceforcings.precipitation specified: values set as zero' 75 75 if numpy.all(numpy.isnan(md.surfaceforcings.mass_balance)) and (md.surfaceforcings.ispdd==0): 76 md.surfaceforcings.mass_balance=numpy.zeros( md.mesh.numberofvertices)76 md.surfaceforcings.mass_balance=numpy.zeros((md.mesh.numberofvertices,1)) 77 77 print ' no surfaceforcings.mass_balance specified: values set as zero' 78 78 if numpy.all(numpy.isnan(md.basalforcings.melting_rate)): 79 md.basalforcings.melting_rate=numpy.zeros( md.mesh.numberofvertices)79 md.basalforcings.melting_rate=numpy.zeros((md.mesh.numberofvertices,1)) 80 80 print ' no basalforcings.melting_rate specified: values set as zero' 81 81 if numpy.all(numpy.isnan(md.balancethickness.thickening_rate)): 82 md.balancethickness.thickening_rate=numpy.zeros( md.mesh.numberofvertices)82 md.balancethickness.thickening_rate=numpy.zeros((md.mesh.numberofvertices,1)) 83 83 print ' no balancethickness.thickening_rate specified: values set as zero' 84 84 85 md.prognostic.spcthickness=float('NaN')*numpy.ones( md.mesh.numberofvertices)86 md.balancethickness.spcthickness=float('NaN')*numpy.ones( md.mesh.numberofvertices)85 md.prognostic.spcthickness=float('NaN')*numpy.ones((md.mesh.numberofvertices,1)) 86 md.balancethickness.spcthickness=float('NaN')*numpy.ones((md.mesh.numberofvertices,1)) 87 87 88 88 if numpy.size(md.initialization.temperature)==md.mesh.numberofvertices: 89 md.thermal.spctemperature=float('NaN')*numpy.ones( md.mesh.numberofvertices)89 md.thermal.spctemperature=float('NaN')*numpy.ones((md.mesh.numberofvertices,1)) 90 90 # pos=find(md.mesh.vertexonsurface); md.thermal.spctemperature(pos)=md.initialization.temperature(pos); %impose observed temperature on surface 91 91 pos=[i for i,vos in enumerate(md.mesh.vertexonsurface) if vos] 92 92 md.thermal.spctemperature[pos]=md.initialization.temperature[pos] # impose observed temperature on surface 93 93 if not numpy.size(md.basalforcings.geothermalflux)==md.mesh.numberofvertices: 94 md.basalforcings.geothermalflux=numpy.zeros( md.mesh.numberofvertices)94 md.basalforcings.geothermalflux=numpy.zeros((md.mesh.numberofvertices,1)) 95 95 else: 96 96 print ' no thermal boundary conditions created: no observed temperature found' -
issm/trunk-jpl/src/m/classes/model/model.py
r13360 r13462 242 242 y3d=numpy.concatenate((y3d,md.mesh.y)) 243 243 #nodes are distributed between bed and surface accordingly to the given exponent 244 z3d=numpy.concatenate((z3d, bed3d+thickness3d*extrusionlist[i]))244 z3d=numpy.concatenate((z3d,(bed3d+thickness3d*extrusionlist[i]).reshape(-1))) 245 245 number_nodes3d=numpy.size(x3d) #number of 3d nodes for the non extruded part of the mesh 246 246 … … 402 402 #Put lithostatic pressure if there is an existing pressure 403 403 if not numpy.any(numpy.isnan(md.initialization.pressure)): 404 md.initialization.pressure=md.constants.g*md.materials.rho_ice*(md.geometry.surface-md.mesh.z )404 md.initialization.pressure=md.constants.g*md.materials.rho_ice*(md.geometry.surface-md.mesh.z.reshape(-1,1)) 405 405 406 406 #special for thermal modeling: -
issm/trunk-jpl/src/m/classes/thermal.py
r13439 r13462 72 72 if EnthalpyAnalysisEnum() in analyses and (md.thermal.isenthalpy or solution==EnthalpySolutionEnum()) and md.mesh.dimension==3: 73 73 pos=numpy.nonzero(numpy.logical_not(numpy.isnan(md.thermal.spctemperature[0:md.mesh.numberofvertices]))) 74 # need to know if md.thermal.spctemperature can have multiple columns? 75 # replicate=numpy.tile(md.geometry.surface-md.mesh.z,(1,numpy.size(md.thermal.spctemperature,axis=1))) 76 # 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") 77 replicate=numpy.tile(md.geometry.surface-md.mesh.z,(1)) 78 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") 74 replicate=numpy.tile(md.geometry.surface-md.mesh.z,(1,numpy.size(md.thermal.spctemperature,axis=1))) 75 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") 79 76 md = checkfield(md,'thermal.isenthalpy','numel',[1],'values',[0,1]) 80 77 -
issm/trunk-jpl/src/m/materials/paterson.m
r13011 r13462 2 2 %PATERSON - figure out the rigidity of ice for a given temperature 3 3 % 4 % rigidi gty (in s^(1/3)Pa) is the flow law paramter in the flow law sigma=B*e(1/3) (Paterson, p97).4 % rigidity (in s^(1/3)Pa) is the flow law paramter in the flow law sigma=B*e(1/3) (Paterson, p97). 5 5 % temperature is in Kelvin degrees 6 6 % -
issm/trunk-jpl/src/m/materials/paterson.py
r13011 r13462 1 from numpy import * 1 import numpy 2 2 3 3 def paterson(temperature): 4 """ 5 PATERSON - figure out the rigidity of ice for a given temperature 4 6 5 # Local Variables: pos11, pos5, pos10, temperature, pos, T, pos8, pos9, pos6, pos7, pos4, rigidity, pos2, pos3, pos1 6 # Function calls: length, zeros, argwhere, paterson, error 7 #PATERSON - figure out the rigidity of ice for a given temperature 8 # 9 # rigidigty (in s^(1/3)Pa) is the flow law paramter in the flow law sigma=B*e(1/3) (Paterson, p97). 10 # temperature is in Kelvin degrees 11 # 12 # Usage: 13 # rigidity=paterson(temperature) 14 15 pos=argwhere(temperature<0.) 16 if len(pos): 17 print 'input temperature should be in Kelvin (positive)' 18 return [] 19 7 rigidity (in s^(1/3)Pa) is the flow law paramter in the flow law sigma=B*e(1/3) (Paterson, p97). 8 temperature is in Kelvin degrees 9 10 Usage: 11 rigidity=paterson(temperature) 12 """ 13 14 if numpy.any(temperature<0.): 15 raise RuntimeError("input temperature should be in Kelvin (positive)") 16 20 17 T = temperature-273.15 18 21 19 #The routine below is equivalent to: 20 22 21 # n=3; T=temperature-273; 23 22 # %From paterson, … … 31 30 # rigidity=fittedmodel(temperature); 32 31 33 rigidity=zeros(len(T)) 34 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) 35 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) 36 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) 37 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) 38 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) 39 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) 40 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) 41 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) 42 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) 43 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) 44 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) 32 rigidity=numpy.zeros((numpy.size(T,axis=0),1)) 33 pos1=numpy.nonzero(T<=-45) 34 rigidity[pos1]=10**8*(-0.000292866376675*(T[pos1]+50)**3+ 0.011672640664130*(T[pos1]+50)**2 -0.325004442485481*(T[pos1]+50)+ 6.524779401948101) 35 pos2=numpy.nonzero(numpy.logical_and(-45<=T,T<-40)) 36 rigidity[pos2]=10**8*(-0.000292866376675*(T[pos2]+45)**3+ 0.007279645014004*(T[pos2]+45)**2 -0.230243014094813*(T[pos2]+45)+ 5.154964909039554) 37 pos3=numpy.nonzero(numpy.logical_and(-40<=T,T<-35)) 38 rigidity[pos3]=10**8*(0.000072737147457*(T[pos3]+40)**3+ 0.002886649363879*(T[pos3]+40)**2 -0.179411542205399*(T[pos3]+40)+ 4.149132666831214) 39 pos4=numpy.nonzero(numpy.logical_and(-35<=T,T<-30)) 40 rigidity[pos4]=10**8*(-0.000086144770023*(T[pos4]+35)**3+ 0.003977706575736*(T[pos4]+35)**2 -0.145089762507325*(T[pos4]+35)+ 3.333333333333331) 41 pos5=numpy.nonzero(numpy.logical_and(-30<=T,T<-25)) 42 rigidity[pos5]=10**8*(-0.000043984685769*(T[pos5]+30)**3+ 0.002685535025386*(T[pos5]+30)**2 -0.111773554501713*(T[pos5]+30)+ 2.696559088937191) 43 pos6=numpy.nonzero(numpy.logical_and(-25<=T,T<-20)) 44 rigidity[pos6]=10**8*(-0.000029799523463*(T[pos6]+25)**3+ 0.002025764738854*(T[pos6]+25)**2 -0.088217055680511*(T[pos6]+25)+ 2.199331606342181) 45 pos7=numpy.nonzero(numpy.logical_and(-20<=T,T<-15)) 46 rigidity[pos7]=10**8*(0.000136920904777*(T[pos7]+20)**3+ 0.001578771886910*(T[pos7]+20)**2 -0.070194372551690*(T[pos7]+20)+ 1.805165505978111) 47 pos8=numpy.nonzero(numpy.logical_and(-15<=T,T<-10)) 48 rigidity[pos8]=10**8*(-0.000899763781026*(T[pos8]+15)**3+ 0.003632585458564*(T[pos8]+15)**2 -0.044137585824322*(T[pos8]+15)+ 1.510778053489523) 49 pos9=numpy.nonzero(numpy.logical_and(-10<=T,T<-5)) 50 rigidity[pos9]=10**8*(0.001676964325070*(T[pos9]+10)**3- 0.009863871256831*(T[pos9]+10)**2 -0.075294014815659*(T[pos9]+10)+ 1.268434288203714) 51 pos10=numpy.nonzero(numpy.logical_and(-5<=T,T<-2)) 52 rigidity[pos10]=10**8*(-0.003748937622487*(T[pos10]+5)**3+0.015290593619213*(T[pos10]+5)**2 -0.048160403003748*(T[pos10]+5)+ 0.854987973338348) 53 pos11=numpy.nonzero(-2<=T) 54 rigidity[pos11]=10**8*(-0.003748937622488*(T[pos11]+2)**3-0.018449844983174*(T[pos11]+2)**2 -0.057638157095631*(T[pos11]+2)+ 0.746900791092860) 45 55 46 56 #Now make sure that rigidity is positive 47 pos=argwhere(rigidity<0); rigidity[pos]=1**6 57 pos=numpy.nonzero(rigidity<0) 58 rigidity[pos]=1**6 48 59 49 60 return rigidity -
issm/trunk-jpl/src/m/parameterization/setflowequation.py
r13244 r13462 94 94 numpy.logical_not(numpy.isnan(md.diagnostic.spcvy)).astype(int)+ \ 95 95 numpy.logical_not(numpy.isnan(md.diagnostic.spcvz)).astype(int)==3, \ 96 numpy.logical_and(nodeonpattyn,nodeonstokes) ).astype(int) #find all the nodes on the boundary of the domain without icefront96 numpy.logical_and(nodeonpattyn,nodeonstokes).reshape(-1,1)).astype(int) #find all the nodes on the boundary of the domain without icefront 97 97 # fullspcelems=double(sum(fullspcnodes(md.mesh.elements),2)==6); %find all the nodes on the boundary of the domain without icefront 98 98 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 99 stokesflag[numpy.nonzero(fullspcelems )]=099 stokesflag[numpy.nonzero(fullspcelems.reshape(-1))]=0 100 100 nodeonstokes[md.mesh.elements[numpy.nonzero(stokesflag),:].astype(int)-1]=1 101 101 … … 129 129 penalties=numpy.zeros((0,2)) 130 130 for i in xrange(1,numlayers): 131 penalties=numpy. concatenate((penalties,numpy.concatenate((bordernodes2d,bordernodes2d+md.mesh.numberofvertices2d*(i)),axis=1)),axis=0)131 penalties=numpy.vstack((penalties,numpy.hstack((bordernodes2d,bordernodes2d+md.mesh.numberofvertices2d*(i))))) 132 132 md.diagnostic.vertex_pairing=penalties 133 133 -
issm/trunk-jpl/src/m/solve/parseresultsfromdisk.py
r13116 r13462 52 52 result['fieldname'] in results[result['step']] and \ 53 53 not strcmp(result['fieldname'],'SolutionType'): 54 results[result['step']][result['fieldname']]=numpy. concatenate((results[result['step']][result['fieldname']],result['field']),axis=0)54 results[result['step']][result['fieldname']]=numpy.vstack((results[result['step']][result['fieldname']],result['field'])) 55 55 else: 56 56 results[result['step']][result['fieldname']]=result['field']
Note:
See TracChangeset
for help on using the changeset viewer.