Changeset 22267
- Timestamp:
- 11/18/17 15:19:15 (7 years ago)
- Location:
- issm/trunk-jpl
- Files:
-
- 43 added
- 18 edited
Legend:
- Unmodified
- Added
- Removed
-
issm/trunk-jpl/src/m/classes/SMBgemb.m
r21232 r22267 358 358 359 359 WriteData(fid,prefix,'data',dt,'name','md.smb.dt','format','Double','scale',yts); 360 360 361 361 % Check if smb_dt goes evenly into transient core time step 362 362 if (mod(md.timestepping.time_step,dt) >= 1e-10) -
issm/trunk-jpl/src/m/classes/SMBgradientsela.py
r21469 r22267 9 9 10 10 Usage: 11 SMBgradientsela=SMBgradientsela() ;11 SMBgradientsela=SMBgradientsela() 12 12 """ 13 13 … … 19 19 self.b_min = float('NaN') 20 20 self.requested_outputs = [] 21 self.setdefaultparameters() 21 22 #}}} 22 23 def __repr__(self): # {{{ 23 string=" surface forcings parameters:" 24 string = " surface forcings parameters:" 25 string+= '\n SMB gradients ela parameters:' 24 26 25 string="%s\n%s"%(string,fielddisplay(self,'issmbgradientsela','is smb gradients ela method activated (0 or 1, default is 0)'))26 27 string="%s\n%s"%(string,fielddisplay(self,'ela',' equilibrium line altitude from which deviation is used to calculate smb using the smb gradients ela method [m a.s.l.]')) 27 28 string="%s\n%s"%(string,fielddisplay(self,'b_pos',' vertical smb gradient (dB/dz) above ela')) … … 47 48 return self 48 49 #}}} 50 def setdefaultparameters(self): # {{{ 51 self.b_max=9999. 52 self.b_min=-9999. 53 return self 54 #}}} 49 55 def checkconsistency(self,md,solution,analyses): # {{{ 50 56 … … 56 62 md = checkfield(md,'fieldname','smb.b_min','timeseries',1,'NaN',1,'Inf',1) 57 63 58 md = checkfield(md,'fieldname',' masstransport.requested_outputs','stringrow',1)64 md = checkfield(md,'fieldname','smb.requested_outputs','stringrow',1) 59 65 return md 60 66 # }}} -
issm/trunk-jpl/src/m/classes/amr.py
r22241 r22267 28 28 self.deviatoricerror_threshold = 0. 29 29 self.deviatoricerror_maximum = 0. 30 self.restart=0. 30 31 #set defaults 31 32 self.setdefaultparameters() … … 48 49 string="%s\n%s"%(string,fielddisplay(self,"deviatoricerror_threshold","maximum threshold deviatoricstress error permitted")) 49 50 string="%s\n%s"%(string,fielddisplay(self,"deviatoricerror_maximum","maximum deviatoricstress error permitted")) 51 string="%s\n%s"%(string,fielddisplay(self,'restart',['indicates if ReMesh() will call before first time step'])) 50 52 return string 51 53 #}}} … … 67 69 self.deviatoricerror_threshold = 0 68 70 self.deviatoricerror_maximum = 0 71 self.restart = 0. 69 72 return self 70 73 #}}} … … 84 87 md = checkfield(md,'fieldname','amr.deviatoricerror_threshold','numel',[1],'>=',0,'<=',1,'NaN',1); 85 88 md = checkfield(md,'fieldname','amr.deviatoricerror_maximum','numel',[1],'>=',0,'NaN',1,'Inf',1); 89 md = checkfield(md,'fieldname','amr.restart','numel',[1],'>=',0,'<=',1,'NaN',1) 86 90 return md 87 91 # }}} … … 104 108 WriteData(fid,prefix,'object',self,'fieldname','deviatoricerror_threshold','format','Double'); 105 109 WriteData(fid,prefix,'object',self,'fieldname','deviatoricerror_maximum','format','Double'); 110 WriteData(fid,prefix,'object',self,'class','amr','fieldname','restart','format','Integer') 106 111 # }}} -
issm/trunk-jpl/src/m/classes/basalforcings.py
r21303 r22267 45 45 self.floatingice_melting_rate=np.zeros((md.mesh.numberofvertices)) 46 46 print " no basalforcings.floatingice_melting_rate specified: values set as zero" 47 #if np.all(np.isnan(self.geothermalflux)): 48 #self.geothermalflux=np.zeros((md.mesh.numberofvertices)) 49 #print " no basalforcings.geothermalflux specified: values set as zero" 47 50 48 51 return self -
issm/trunk-jpl/src/m/classes/constants.py
r22122 r22267 12 12 13 13 def __init__(self): # {{{ 14 self.g = 0 15 self.omega = 0 16 self.yts = 0 17 self.referencetemperature = 0 14 self.g = 0. 15 self.omega = 0. 16 self.yts = 0. 17 self.referencetemperature = 0. 18 18 19 19 #set defaults … … 49 49 def checkconsistency(self,md,solution,analyses): # {{{ 50 50 51 md = checkfield(md,'fieldname','constants.g','> ',0,'size',[1])52 md = checkfield(md,'fieldname','constants.omega','>=',0,'size',[1 ])53 md = checkfield(md,'fieldname','constants.yts','>',0,'size',[1 ])54 md = checkfield(md,'fieldname','constants.referencetemperature','size',[1 ])51 md = checkfield(md,'fieldname','constants.g','>=',0,'size',[1,1]) 52 md = checkfield(md,'fieldname','constants.omega','>=',0,'size',[1,1]) 53 md = checkfield(md,'fieldname','constants.yts','>',0,'size',[1,1]) 54 md = checkfield(md,'fieldname','constants.referencetemperature','size',[1,1]) 55 55 56 56 return md -
issm/trunk-jpl/src/m/classes/flowequation.py
r21303 r22267 92 92 md = checkfield(md,'fieldname','flowequation.fe_SSA','values',['P1','P1bubble','P1bubblecondensed','P2','P2bubble']) 93 93 md = checkfield(md,'fieldname','flowequation.fe_HO' ,'values',['P1','P1bubble','P1bubblecondensed','P1xP2','P2xP1','P2','P2bubble','P1xP3','P2xP4']) 94 md = checkfield(md,'fieldname','flowequation.fe_FS' ,'values',['P1P1','P1P1GLS','MINIcondensed','MINI','TaylorHood',' XTaylorHood','OneLayerP4z','CrouzeixRaviart'])94 md = checkfield(md,'fieldname','flowequation.fe_FS' ,'values',['P1P1','P1P1GLS','MINIcondensed','MINI','TaylorHood','LATaylorHood','XTaylorHood','OneLayerP4z','CrouzeixRaviart','LACrouzeixRaviart']) 95 95 md = checkfield(md,'fieldname','flowequation.borderSSA','size',[md.mesh.numberofvertices],'values',[0,1]) 96 96 md = checkfield(md,'fieldname','flowequation.borderHO','size',[md.mesh.numberofvertices],'values',[0,1]) … … 104 104 md = checkfield(md,'fieldname','flowequation.vertex_equation','size',[md.mesh.numberofvertices],'values',[1,2]) 105 105 md = checkfield(md,'fieldname','flowequation.element_equation','size',[md.mesh.numberofelements],'values',[1,2]) 106 elif m.strcmp(md.mesh.domaintype(),'2Dvertical'): 107 md = checkfield(md,'fieldname','flowequation.vertex_equation','size',[md.mesh.numberofvertices],'values',[2,4,5]) 108 md = checkfield(md,'fieldname','flowequation.element_equation','size',[md.mesh.numberofelements],'values',[2,4,5]) 106 109 elif m.strcmp(md.mesh.domaintype(),'3D'): 107 110 md = checkfield(md,'fieldname','flowequation.vertex_equation','size',[md.mesh.numberofvertices],'values',np.arange(0,8+1)) -
issm/trunk-jpl/src/m/classes/mismipbasalforcings.py
r21354 r22267 40 40 def initialize(self,md): # {{{ 41 41 if np.all(np.isnan(self.groundedice_melting_rate)): 42 self.groundedice_melting_rate=np.zeros( md.mesh.numberofvertices)42 self.groundedice_melting_rate=np.zeros((md.mesh.numberofvertices)) 43 43 print ' no basalforcings.groundedice_melting_rate specified: values set as zero' 44 if np.all(np.isnan(self.geothermalflux)): 45 self.geothermalflux=np.zeros((md.mesh.numberofvertices)) 46 print " no basalforcings.geothermalflux specified: values set as zero" 44 47 return self 45 48 #}}} … … 84 87 85 88 floatingice_melting_rate = np.zeros((md.mesh.numberofvertices)) 86 floatingice_melting_rate = md.basalforcings.meltrate_factor*np.tanh((md.geometry.base-md.geometry.bed)/md.basalforcings.threshold_thickness)* np.amax(md.basalforcings.upperdepth_melt-md.geometry.base,0)89 floatingice_melting_rate = md.basalforcings.meltrate_factor*np.tanh((md.geometry.base-md.geometry.bed)/md.basalforcings.threshold_thickness)*(md.basalforcings.upperdepth_melt-md.geometry.base) 87 90 88 91 WriteData(fid,prefix,'name','md.basalforcings.model','data',3,'format','Integer') -
issm/trunk-jpl/src/m/classes/taoinversion.py
r21303 r22267 6 6 from IssmConfig import IssmConfig 7 7 from marshallcostfunctions import marshallcostfunctions 8 from supportedcontrols import * 9 from supportedcostfunctions import * 8 10 9 10 class taoinversion: 11 class taoinversion(object): 11 12 def __init__(self): 12 iscontrol = 0 13 incomplete_adjoint = 0 14 control_parameters = float('NaN') 15 maxsteps = 0 16 maxiter = 0 17 fatol = 0 18 frtol = 0 19 gatol = 0 20 grtol = 0 21 gttol = 0 22 algorithm = '' 23 cost_functions = float('NaN') 24 cost_functions_coefficients = float('NaN') 25 min_parameters = float('NaN') 26 max_parameters = float('NaN') 27 vx_obs = float('NaN') 28 vy_obs = float('NaN') 29 vz_obs = float('NaN') 30 vel_obs = float('NaN') 31 thickness_obs = float('NaN') 32 surface_obs = float('NaN') 13 self.iscontrol = 0 14 self.incomplete_adjoint = 0 15 self.control_parameters = float('NaN') 16 self.maxsteps = 0 17 self.maxiter = 0 18 self.fatol = 0 19 self.frtol = 0 20 self.gatol = 0 21 self.grtol = 0 22 self.gttol = 0 23 self.algorithm = '' 24 self.cost_functions = float('NaN') 25 self.cost_functions_coefficients = float('NaN') 26 self.min_parameters = float('NaN') 27 self.max_parameters = float('NaN') 28 self.vx_obs = float('NaN') 29 self.vy_obs = float('NaN') 30 self.vz_obs = float('NaN') 31 self.vel_obs = float('NaN') 32 self.thickness_obs = float('NaN') 33 self.surface_obs = float('NaN') 34 self.setdefaultparameters() 33 35 34 36 def __repr__(self): … … 98 100 #several responses can be used: 99 101 self.cost_functions=101; 100 101 102 return self 102 103 … … 119 120 120 121 def checkconsistency(self,md,solution,analyses): 121 if not self. control:122 if not self.iscontrol: 122 123 return md 123 124 if not IssmConfig('_HAVE_TAO_')[0]: … … 125 126 126 127 127 num_controls= np. numel(md.inversion.control_parameters)128 num_costfunc= np.size(md.inversion.cost_functions ,2)128 num_controls= np.size(md.inversion.control_parameters) 129 num_costfunc= np.size(md.inversion.cost_functions) 129 130 130 131 md = checkfield(md,'fieldname','inversion.iscontrol','values',[0, 1]) 131 md = checkfield(md,'fieldname','inversion.incomplete_adjoint','values',[0, 132 md = checkfield(md,'fieldname','inversion.incomplete_adjoint','values',[0,1]) 132 133 md = checkfield(md,'fieldname','inversion.control_parameters','cell',1,'values',supportedcontrols()) 133 134 md = checkfield(md,'fieldname','inversion.maxsteps','numel',1,'>=',0) … … 143 144 PETSCMINOR = IssmConfig('_PETSC_MINOR_')[0] 144 145 if(PETSCMAJOR>3 or (PETSCMAJOR==3 and PETSCMINOR>=5)): 145 md = checkfield(md,'fieldname','inversion.algorithm','values', {'blmvm','cg','lmvm'})146 md = checkfield(md,'fieldname','inversion.algorithm','values',['blmvm','cg','lmvm']) 146 147 else: 147 md = checkfield(md,'fieldname','inversion.algorithm','values', {'tao_blmvm','tao_cg','tao_lmvm'})148 md = checkfield(md,'fieldname','inversion.algorithm','values',['tao_blmvm','tao_cg','tao_lmvm']) 148 149 149 150 150 md = checkfield(md,'fieldname','inversion.cost_functions','size', [1,num_costfunc],'values',supportedcostfunctions())151 md = checkfield(md,'fieldname','inversion.cost_functions','size', [num_costfunc],'values',supportedcostfunctions()) 151 152 md = checkfield(md,'fieldname','inversion.cost_functions_coefficients','size',[md.mesh.numberofvertices, num_costfunc],'>=',0) 152 153 md = checkfield(md,'fieldname','inversion.min_parameters','size',[md.mesh.numberofvertices, num_controls]) … … 162 163 md = checkfield(md,'fieldname','inversion.vy_obs','size',[md.mesh.numberofvertices],'NaN',1,'Inf',1) 163 164 164 def marshall(self, md,fid):165 def marshall(self,prefix,md,fid): 165 166 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 167 yts=md.constants.yts; 168 WriteData(fid,prefix,'object',self,'class','inversion','fieldname','iscontrol','format','Boolean') 169 WriteData(fid,prefix,'name','md.inversion.type','data',1,'format','Integer') 170 if not self.iscontrol: 171 return 172 WriteData(fid,prefix,'object',self,'class','inversion','fieldname','incomplete_adjoint','format','Boolean') 173 WriteData(fid,prefix,'object',self,'class','inversion','fieldname','maxsteps','format','Integer') 174 WriteData(fid,prefix,'object',self,'class','inversion','fieldname','maxiter','format','Integer') 175 WriteData(fid,prefix,'object',self,'class','inversion','fieldname','fatol','format','Double') 176 WriteData(fid,prefix,'object',self,'class','inversion','fieldname','frtol','format','Double') 177 WriteData(fid,prefix,'object',self,'class','inversion','fieldname','gatol','format','Double') 178 WriteData(fid,prefix,'object',self,'class','inversion','fieldname','grtol','format','Double') 179 WriteData(fid,prefix,'object',self,'class','inversion','fieldname','gttol','format','Double') 180 WriteData(fid,prefix,'object',self,'class','inversion','fieldname','algorithm','format','String') 181 WriteData(fid,prefix,'object',self,'class','inversion','fieldname','cost_functions_coefficients','format','DoubleMat','mattype',1) 182 WriteData(fid,prefix,'object',self,'class','inversion','fieldname','min_parameters','format','DoubleMat','mattype',3) 183 WriteData(fid,prefix,'object',self,'class','inversion','fieldname','max_parameters','format','DoubleMat','mattype',3) 184 WriteData(fid,prefix,'object',self,'class','inversion','fieldname','vx_obs','format','DoubleMat','mattype',1,'scale',1./yts) 185 WriteData(fid,prefix,'object',self,'class','inversion','fieldname','vy_obs','format','DoubleMat','mattype',1,'scale',1./yts) 186 WriteData(fid,prefix,'object',self,'class','inversion','fieldname','vz_obs','format','DoubleMat','mattype',1,'scale',1./yts) 187 WriteData(fid,prefix,'object',self,'class','inversion','fieldname','thickness_obs','format','DoubleMat','mattype',1) 188 WriteData(fid,prefix,'object',self,'class','inversion','fieldname','surface_obs','format','DoubleMat','mattype',1) 188 189 189 190 num_control_parameters = np.numel(self.control_parameters)191 192 190 #process control parameters 191 num_control_parameters = np.size(self.control_parameters) 192 WriteData(fid,prefix,'object',self,'fieldname','control_parameters','format','StringArray') 193 WriteData(fid,prefix,'data',num_control_parameters,'name','md.inversion.num_control_parameters','format','Integer') 193 194 194 195 num_cost_functions = np.size(self.cost_functions,2)196 197 198 195 #process cost functions 196 num_cost_functions = np.size(self.cost_functions) 197 data= marshallcostfunctions(self.cost_functions) 198 WriteData(fid,prefix,'data',data,'name','md.inversion.cost_functions','format','StringArray') 199 WriteData(fid,prefix,'data',num_cost_functions,'name','md.inversion.num_cost_functions','format','Integer') -
issm/trunk-jpl/src/m/classes/thermal.py
r21542 r22267 101 101 102 102 if 'EnthalpyAnalysis' in analyses and md.thermal.isenthalpy and md.mesh.dimension()==3: 103 pos=np.where(~np.isnan(md.thermal.spctemperature[0:md.mesh.numberofvertices])) 103 TEMP = md.thermal.spctemperature[:-1].flatten(-1) 104 pos=np.where(~np.isnan(TEMP)) 104 105 try: 105 106 spccol=np.size(md.thermal.spctemperature,1) 106 107 except IndexError: 107 108 spccol=1 108 replicate=np.tile(md.geometry.surface-md.mesh.z,(spccol)) 109 control=md.materials.meltingpoint-md.materials.beta*md.materials.rho_ice*md.constants.g*replicate 110 md = checkfield(md,'fieldname','thermal.spctemperature','field',md.thermal.spctemperature[pos],'<=',control[pos],'message',"spctemperature should be below the adjusted melting point") 109 110 replicate=np.tile(md.geometry.surface-md.mesh.z,(spccol)).flatten(-1) 111 112 control=md.materials.meltingpoint-md.materials.beta*md.materials.rho_ice*md.constants.g*replicate+10**-5 113 114 md = checkfield(md,'fieldname','thermal.spctemperature','field',md.thermal.spctemperature.flatten(-1)[pos],'<=',control[pos],'message',"spctemperature should be below the adjusted melting point") 111 115 md = checkfield(md,'fieldname','thermal.isenthalpy','numel',[1],'values',[0,1]) 112 116 md = checkfield(md,'fieldname','thermal.isdynamicbasalspc','numel',[1],'values',[0,1]); -
issm/trunk-jpl/src/m/classes/transient.py
r21802 r22267 12 12 13 13 def __init__(self): # {{{ 14 self.issmb = False14 self.issmb = False 15 15 self.ismasstransport = False 16 16 self.isstressbalance = False … … 23 23 self.ishydrology = False 24 24 self.isslr = False 25 self.iscoupler = False 26 self.amr_frequency = 0 25 27 self.isoceancoupling = False 26 self.iscoupler = False27 amr_frequency = 028 28 self.requested_outputs = [] 29 29 … … 76 76 self.iscoupler = False 77 77 self.amr_frequency = 0 78 79 #default output 80 self.requested_outputs=[] 81 return self 82 #}}} 83 def deactivateall(self):#{{{ 84 self.issmb = False 85 self.ismasstransport = False 86 self.isstressbalance = False 87 self.isthermal = False 88 self.isgroundingline = False 89 self.isgia = False 90 self.isesa = False 91 self.isdamageevolution = False 92 self.ismovingfront = False 93 self.ishydrology = False 94 self.isslr = False 95 self.isoceancoupling = False 96 self.iscoupler = False 97 self.amr_frequency = 0 78 98 79 99 #default output -
issm/trunk-jpl/src/m/consistency/checkfield.py
r22150 r22267 84 84 if options.exist('numel'): 85 85 fieldnumel=options.getfieldvalue('numel') 86 if np.size(field) not in fieldnumel:86 if (type(fieldnumel) == int and np.size(field) != fieldnumel) or (type(fieldnumel) == list and np.size(field) not in fieldnumel): 87 87 if len(fieldnumel)==1: 88 88 md = md.checkmessage(options.getfieldvalue('message',\ … … 101 101 "NaN values found in field '%s'" % fieldname)) 102 102 103 103 104 #check Inf 104 105 if options.getfieldvalue('Inf',0): … … 106 107 md = md.checkmessage(options.getfieldvalue('message',\ 107 108 "Inf values found in field '%s'" % fieldname)) 109 108 110 109 111 #check cell … … 129 131 #check greater 130 132 if options.exist('>='): 131 lowerbound=options.getfieldvalue('>=') 132 if np.any(field<lowerbound): 133 lowerbound = options.getfieldvalue('>=') 134 field2 = np.reshape(field,(np.prod(np.shape(field)),1),order='F').copy() 135 if options.getfieldvalue('timeseries',0): 136 field2 = np.reshape(field[:-1],np.prod(np.shape(field[:-1])),1) 137 138 if options.getfieldvalue('singletimeseries',0): 139 field2 = np.reshape(field[0],np.prod(np.shape(field[0])),1) 140 141 if np.any(field2<lowerbound): 133 142 md = md.checkmessage(options.getfieldvalue('message',\ 134 143 "field '%s' should have values above %d" % (fieldname,lowerbound))) 144 135 145 if options.exist('>'): 136 146 lowerbound=options.getfieldvalue('>') 137 if np.any(field<=lowerbound): 147 field2 = np.reshape(field,(np.prod(np.shape(field)),1),order='F').copy() 148 if options.getfieldvalue('timeseries',0): 149 field2 = np.reshape(field[:-1],np.prod(np.shape(field[:-1])),1) 150 151 if options.getfieldvalue('singletimeseries',0): 152 field2 = np.reshape(field[0],np.prod(np.shape(field[0])),1) 153 154 if np.any(field2<=lowerbound): 138 155 md = md.checkmessage(options.getfieldvalue('message',\ 139 156 "field '%s' should have values above %d" % (fieldname,lowerbound))) … … 142 159 if options.exist('<='): 143 160 upperbound=options.getfieldvalue('<=') 144 if np.any(field>upperbound): 161 field2 = np.reshape(field,(np.prod(np.shape(field)),1),order='F').copy() 162 if options.getfieldvalue('timeseries',0): 163 field2 = np.reshape(field[:-1],np.prod(np.shape(field[:-1])),1) 164 165 if options.getfieldvalue('singletimeseries',0): 166 field2 = np.reshape(field[0],np.prod(np.shape(field[0])),1) 167 168 if np.any(field2>upperbound): 145 169 md = md.checkmessage(options.getfieldvalue('message',\ 146 170 "field '%s' should have values below %d" % (fieldname,upperbound))) 147 171 if options.exist('<'): 148 172 upperbound=options.getfieldvalue('<') 149 if np.any(field>=upperbound): 173 field2 = np.reshape(field,(np.prod(np.shape(field)),1),order='F').copy() 174 if options.getfieldvalue('timeseries',0): 175 field2 = np.reshape(field[:-1],np.prod(np.shape(field[:-1])),1) 176 177 if options.getfieldvalue('singletimeseries',0): 178 field2 = np.reshape(field[0],np.prod(np.shape(field[0])),1) 179 180 if np.any(field2>=upperbound): 150 181 md = md.checkmessage(options.getfieldvalue('message',\ 151 182 "field '%s' should have values below %d" % (fieldname,upperbound))) … … 164 195 #Check forcings (size and times) 165 196 if options.getfieldvalue('timeseries',0): 166 if np.size(field,0)==md.mesh.numberofvertices:197 if np.size(field,0)==md.mesh.numberofvertices or np.size(field,0)==md.mesh.numberofelements: 167 198 if np.ndim(field)>1 and not np.size(field,1)==1: 168 199 md = md.checkmessage(options.getfieldvalue('message',\ 169 200 "field '%s' should have only one column as there are md.mesh.numberofvertices lines" % fieldname)) 170 elif np.size(field,0)==md.mesh.numberofvertices+1 or np.size(field,0)== 2:171 if n ot all(field[-1,:]==np.sort(field[-1,:])):201 elif np.size(field,0)==md.mesh.numberofvertices+1 or np.size(field,0)==md.mesh.numberofelements+1: 202 if np.ndim(field) > 1 and not all(field[-1,:]==np.sort(field[-1,:])): 172 203 md = md.checkmessage(options.getfieldvalue('message',\ 173 204 "field '%s' columns should be sorted chronologically" % fieldname)) 174 if any(field[-1,0:-1]==field[-1,1:]):205 if np.ndim(field) > 1 and any(field[-1,0:-1]==field[-1,1:]): 175 206 md = md.checkmessage(options.getfieldvalue('message',\ 176 207 "field '%s' columns must not contain duplicate timesteps" % fieldname)) … … 198 229 return md 199 230 231 -
issm/trunk-jpl/src/m/dev/devpath.py
r20679 r22267 31 31 #Manual imports for commonly used functions 32 32 from plotmodel import plotmodel 33 from runme import runme 33 34 34 35 #c = get_ipython().config -
issm/trunk-jpl/src/m/mesh/bamg.py
r22216 r22267 1 1 import os.path 2 2 import numpy as np 3 from mesh2d import mesh2d3 from mesh2dvertical import * 4 4 from collections import OrderedDict 5 5 from pairoptions import pairoptions … … 80 80 #Check that file exists 81 81 domainfile=options.getfieldvalue('domain') 82 if not os.path.exists(domainfile): 83 raise IOError("bamg error message: file '%s' not found" % domainfile) 84 domain=expread(domainfile) 82 if type(domainfile) == str: 83 if not os.path.exists(domainfile): 84 raise IOError("bamg error message: file '%s' not found" % domainfile) 85 domain=expread(domainfile) 86 else: 87 domain=domainfile 85 88 86 89 #Build geometry … … 89 92 90 93 #Check that the domain is closed 91 if (domaini ['x'][0]!=domaini['x'][-1] or domaini['y'][0]!=domaini['y'][-1]):94 if (domaini.x[0]!=domaini.x[-1] or domaini.y[0]!=domaini.y[-1]): 92 95 raise RuntimeError("bamg error message: all contours provided in ''domain'' should be closed") 93 96 94 97 #Checks that all holes are INSIDE the principle domain outline 95 98 if i: 96 flags=ContourToNodes(domaini ['x'],domaini['y'],domainfile,0)[0]99 flags=ContourToNodes(domaini.x,domaini.y,domainfile,0)[0] 97 100 if np.any(np.logical_not(flags)): 98 101 raise RuntimeError("bamg error message: All holes should be strictly inside the principal domain") 99 102 100 103 #Add all points to bamg_geometry 101 nods=domaini ['nods']-1 #the domain are closed 0=end102 bamg_geometry.Vertices=np.vstack((bamg_geometry.Vertices,np.vstack((domaini ['x'][0:nods],domaini['y'][0:nods],np.ones((nods)))).T))104 nods=domaini.nods-1 #the domain are closed 0=end 105 bamg_geometry.Vertices=np.vstack((bamg_geometry.Vertices,np.vstack((domaini.x[0:nods],domaini.y[0:nods],np.ones((nods)))).T)) 103 106 bamg_geometry.Edges =np.vstack((bamg_geometry.Edges,np.vstack((np.arange(count+1,count+nods+1),np.hstack((np.arange(count+2,count+nods+1),count+1)),1.*np.ones((nods)))).T)) 104 107 if i: … … 115 118 #update counter 116 119 count+=nods 120 121 if options.getfieldvalue('vertical',0): 122 if np.size(options.getfieldvalue('Markers',[]))!=np.size(bamg_geometry.Edges,0): 123 raise RuntimeError('for 2d vertical mesh, ''Markers'' option is required, and should be of size ' + str(np.size(bamg_geometry.Edges,0))) 124 if np.size(options.getfieldvalue('Markers',[]))==np.size(bamg_geometry.Edges,0): 125 bamg_geometry.Edges[:,2]=options.getfieldvalue('Markers') 117 126 118 127 #take care of rifts … … 323 332 324 333 # plug results onto model 334 if options.getfieldvalue('vertical',0): 335 md.mesh=mesh2dvertical() 336 md.mesh.x=bamgmesh_out['Vertices'][:,0].copy() 337 md.mesh.y=bamgmesh_out['Vertices'][:,1].copy() 338 md.mesh.elements=bamgmesh_out['Triangles'][:,0:3].astype(int) 339 md.mesh.edges=bamgmesh_out['IssmEdges'].astype(int) 340 md.mesh.segments=bamgmesh_out['IssmSegments'][:,0:3].astype(int) 341 md.mesh.segmentmarkers=bamgmesh_out['IssmSegments'][:,3].astype(int) 342 343 #Fill in rest of fields: 344 md.mesh.numberofelements=np.size(md.mesh.elements,axis=0) 345 md.mesh.numberofvertices=np.size(md.mesh.x) 346 md.mesh.numberofedges=np.size(md.mesh.edges,axis=0) 347 md.mesh.vertexonboundary=np.zeros(md.mesh.numberofvertices,bool) 348 md.mesh.vertexonboundary[md.mesh.segments[:,0:2]-1]=True 349 350 #Now, build the connectivity tables for this mesh. Doubled in matlab for some reason 351 md.mesh.vertexonboundary=np.zeros(md.mesh.numberofvertices,) 352 md.mesh.vertexonboundary[md.mesh.segments[:,0:2]-1]=1 353 354 elif options.getfieldvalue('3dsurface',0): 355 md.mesh=mesh3dsurface() 356 md.mesh.x=bamgmesh_out['Vertices'][:,0].copy() 357 md.mesh.y=bamgmesh_out['Vertices'][:,1].copy() 358 md.mesh.z=md.mesh.x 359 md.mesh.z[:]=0 360 md.mesh.elements=bamgmesh_out['Triangles'][:,0:3].astype(int) 361 md.mesh.edges=bamgmesh_out['IssmEdges'].astype(int) 362 md.mesh.segments=bamgmesh_out['IssmSegments'][:,0:3].astype(int) 363 md.mesh.segmentmarkers=bamgmesh_out['IssmSegments'][:,3].astype(int) 364 365 #Fill in rest of fields: 366 md.mesh.numberofelements=np.size(md.mesh.elements,axis=0) 367 md.mesh.numberofvertices=np.size(md.mesh.x) 368 md.mesh.numberofedges=np.size(md.mesh.edges,axis=0) 369 md.mesh.vertexonboundary=np.zeros(md.mesh.numberofvertices,bool) 370 md.mesh.vertexonboundary[md.mesh.segments[:,0:2]-1]=True 371 372 else: 373 md.mesh=mesh2d() 374 md.mesh.x=bamgmesh_out['Vertices'][:,0].copy() 375 md.mesh.y=bamgmesh_out['Vertices'][:,1].copy() 376 md.mesh.elements=bamgmesh_out['Triangles'][:,0:3].astype(int) 377 md.mesh.edges=bamgmesh_out['IssmEdges'].astype(int) 378 md.mesh.segments=bamgmesh_out['IssmSegments'][:,0:3].astype(int) 379 md.mesh.segmentmarkers=bamgmesh_out['IssmSegments'][:,3].astype(int) 380 381 #Fill in rest of fields: 382 md.mesh.numberofelements=np.size(md.mesh.elements,axis=0) 383 md.mesh.numberofvertices=np.size(md.mesh.x) 384 md.mesh.numberofedges=np.size(md.mesh.edges,axis=0) 385 md.mesh.vertexonboundary=np.zeros(md.mesh.numberofvertices,bool) 386 md.mesh.vertexonboundary[md.mesh.segments[:,0:2]-1]=True 387 388 #Bamg private fields 325 389 md.private.bamg=OrderedDict() 326 390 md.private.bamg['mesh']=bamgmesh(bamgmesh_out) 327 391 md.private.bamg['geometry']=bamggeom(bamggeom_out) 328 md.mesh = mesh2d()329 md.mesh.x=bamgmesh_out['Vertices'][:,0].copy()330 md.mesh.y=bamgmesh_out['Vertices'][:,1].copy()331 md.mesh.elements=bamgmesh_out['Triangles'][:,0:3].astype(int)332 md.mesh.edges=bamgmesh_out['IssmEdges'].astype(int)333 md.mesh.segments=bamgmesh_out['IssmSegments'][:,0:3].astype(int)334 md.mesh.segmentmarkers=bamgmesh_out['IssmSegments'][:,3].astype(int)335 336 #Fill in rest of fields:337 md.mesh.numberofelements=np.size(md.mesh.elements,axis=0)338 md.mesh.numberofvertices=np.size(md.mesh.x)339 md.mesh.numberofedges=np.size(md.mesh.edges,axis=0)340 md.mesh.vertexonboundary=np.zeros(md.mesh.numberofvertices,bool)341 md.mesh.vertexonboundary[md.mesh.segments[:,0:2]-1]=True342 392 md.mesh.elementconnectivity=md.private.bamg['mesh'].ElementConnectivity 343 393 md.mesh.elementconnectivity[np.nonzero(np.isnan(md.mesh.elementconnectivity))]=0 -
issm/trunk-jpl/src/m/solve/WriteData.py
r22204 r22267 51 51 yts = options.getfieldvalue('yts') 52 52 if np.ndim(data) > 1: 53 data[-1,:] = data[-1,:]*yts54 else: 55 data[-1] = data[-1]*yts53 data[-1,:] = yts*data[-1,:] 54 else: 55 data[-1] = yts*data[-1] 56 56 57 57 #Step 1: write the enum to identify this record uniquely -
issm/trunk-jpl/src/m/solve/parseresultsfromdisk.py
r21674 r22267 175 175 if fieldname=='BalancethicknessThickeningRate': 176 176 field = field*yts 177 elif fieldname=='Time':178 field = field/yts179 177 elif fieldname=='HydrologyWaterVx': 180 178 field = field*yts … … 191 189 elif fieldname=='BasalforcingsGroundediceMeltingRate': 192 190 field = field*yts 191 elif fieldname=='BasalforcingsFloatingiceMeltingRate': 192 field = field*yts 193 193 elif fieldname=='TotalFloatingBmb': 194 field = field/10.**12 .*yts #(GigaTon/year)194 field = field/10.**12*yts #(GigaTon/year) 195 195 elif fieldname=='TotalGroundedBmb': 196 field = field/10.**12 .*yts #(GigaTon/year)196 field = field/10.**12*yts #(GigaTon/year) 197 197 elif fieldname=='TotalSmb': 198 field = field/10.**12 .*yts #(GigaTon/year)198 field = field/10.**12*yts #(GigaTon/year) 199 199 elif fieldname=='SmbMassBalance': 200 200 field = field*yts 201 elif fieldname=='SmbPrecipitation': 202 field = field*yts 203 elif fieldname=='SmbRunoff': 204 field = field*yts 205 elif fieldname=='SmbEC': 206 field = field*yts 207 elif fieldname=='SmbAccumulation': 208 field = field*yts 209 elif fieldname=='SmbMelt': 210 field = field*yts 211 elif fieldname=='SmbDz_add': 212 field = field*yts 213 elif fieldname=='SmbM_add': 214 field = field*yts 201 215 elif fieldname=='CalvingCalvingrate': 202 216 field = field*yts 217 218 if time !=-9999: 219 time = time/yts 203 220 204 221 saveres=OrderedDict() -
issm/trunk-jpl/test/Par/ISMIPE.par
r21168 r22267 14 14 md.geometry.surface(i)=(1.-coeff)*data(point1,3)+coeff*data(point2,3); 15 15 end 16 16 17 md.geometry.thickness=md.geometry.surface-md.geometry.base; 17 18 md.geometry.thickness(find(~md.geometry.thickness))=0.01; -
issm/trunk-jpl/test/Par/ISMIPE.py
r21409 r22267 13 13 point1=numpy.floor(y/100.) 14 14 point2=numpy.minimum(point1+1,50) 15 coeff=(y-(point1 -1.)*100.)/100.15 coeff=(y-(point1)*100.)/100. 16 16 md.geometry.base[i]=(1.-coeff)*data[point1,1]+coeff*data[point2,1] 17 17 md.geometry.surface[i]=(1.-coeff)*data[point1,2]+coeff*data[point2,2] 18 18 19 md.geometry.thickness=md.geometry.surface-md.geometry.base 19 20 md.geometry.thickness[numpy.nonzero(numpy.logical_not(md.geometry.thickness))]=0.01 -
issm/trunk-jpl/test/Par/SquareThermal.py
r21409 r22267 26 26 print " creating temperatures" 27 27 md.initialization.temperature=(273.-20.)*numpy.ones((md.mesh.numberofvertices)) 28 md.initialization.pressure=numpy.zeros((md.mesh.numberofvertices,)) 29 md.initialization.waterfraction=numpy.zeros((md.mesh.numberofvertices,)) 30 md.initialization.watercolumn=numpy.zeros((md.mesh.numberofvertices,)) 28 31 29 32 print " creating flow law parameter" … … 33 36 print " creating surface mass balance" 34 37 md.smb.mass_balance=numpy.ones((md.mesh.numberofvertices))/md.constants.yts #1m/a 35 md.basalforcings.melting_rate=0.*numpy.ones((md.mesh.numberofvertices))/md.constants.yts #1m/a 38 #md.basalforcings.melting_rate=0.*numpy.ones((md.mesh.numberofvertices))/md.constants.yts #1m/a 39 md.basalforcings.groundedice_melting_rate=0.*numpy.ones((md.mesh.numberofvertices))/md.constants.yts #1m/a 40 md.basalforcings.floatingice_melting_rate=0.*numpy.ones((md.mesh.numberofvertices))/md.constants.yts #1m/a 36 41 37 42 #Deal with boundary conditions: … … 41 46 42 47 print " boundary conditions for thermal model" 43 md.thermal.spctemperature =md.initialization.temperature48 md.thermal.spctemperature[:]=md.initialization.temperature 44 49 md.basalforcings.geothermalflux=numpy.zeros((md.mesh.numberofvertices)) 45 50 md.basalforcings.geothermalflux[numpy.nonzero(md.mask.groundedice_levelset>0.)[0]]=1.*10**-3 #1 mW/m^2
Note:
See TracChangeset
for help on using the changeset viewer.