Index: /issm/trunk-jpl/src/py3/archive/arch.py =================================================================== --- /issm/trunk-jpl/src/py3/archive/arch.py (revision 23708) +++ /issm/trunk-jpl/src/py3/archive/arch.py (revision 23709) @@ -95,12 +95,12 @@ print('Source file: ') - print('\t{0}'.format(filename)) + print(('\t{0}'.format(filename))) print('Variables: ') result=read_field(fid) while result: - print('\t{0}'.format(result['field_name'])) - print('\t\tSize:\t\t{0}'.format(result['size'])) - print('\t\tDatatype:\t{0}'.format(result['data_type'])) + print(('\t{0}'.format(result['field_name']))) + print(('\t\tSize:\t\t{0}'.format(result['size']))) + print(('\t\tDatatype:\t{0}'.format(result['data_type']))) # go to next result result=read_field(fid) Index: /issm/trunk-jpl/src/py3/array/MatlabArray.py =================================================================== --- /issm/trunk-jpl/src/py3/array/MatlabArray.py (revision 23708) +++ /issm/trunk-jpl/src/py3/array/MatlabArray.py (revision 23709) @@ -29,5 +29,5 @@ ''' if type(ain) != type(aval): - print(allequal.__doc__) + print((allequal.__doc__)) raise RuntimeError("ain and aval must be of the same type") Index: /issm/trunk-jpl/src/py3/classes/SMBgemb.py =================================================================== --- /issm/trunk-jpl/src/py3/classes/SMBgemb.py (revision 23708) +++ /issm/trunk-jpl/src/py3/classes/SMBgemb.py (revision 23709) @@ -14,6 +14,6 @@ def __init__(self): # {{{ - #each one of these properties is a transient forcing to the GEMB model, loaded from meteorological data derived - #from an automatic weather stations (AWS). Each property is therefore a matrix, of size (numberofvertices x number + #each one of these properties is a transient forcing to the GEMB model, loaded from meteorological data derived + #from an automatic weather stations (AWS). Each property is therefore a matrix, of size (numberofvertices x number #of time steps. ) @@ -27,7 +27,7 @@ #ismelt #isdensification - #isturbulentflux - - #inputs: + #isturbulentflux + + #inputs: Ta = float('NaN') #2 m air temperature, in Kelvin V = float('NaN') #wind speed (m/s-1) @@ -37,7 +37,6 @@ eAir = float('NaN') #screen level vapor pressure [Pa] pAir = float('NaN') #surface pressure [Pa] - Tmean = float('NaN') #mean annual temperature [K] - Vmean = float('NaN') #mean annual wind velocity [m s-1] + Vmean = float('NaN') #mean annual wind velocity [m s-1] C = float('NaN') #mean annual snow accumulation [kg m-2 yr-1] Tz = float('NaN') #height above ground at which temperature (T) was sampled [m] @@ -47,5 +46,5 @@ aValue = float('NaN') #Albedo forcing at every element. Used only if aIdx == 0. teValue = float('NaN') #Outward longwave radiation thermal emissivity forcing at every element (default in code is 1) - + # Initialization of snow properties Dzini = float('NaN') #cell depth (m) @@ -60,34 +59,30 @@ Sizeini = float('NaN') #Number of layers - #settings: + #settings: aIdx = float('NaN') #method for calculating albedo and subsurface absorption (default is 1) - # 0: direct input from aValue parameter - # 1: effective grain radius [Gardner & Sharp, 2009] - # 2: effective grain radius [Brun et al., 2009] - # 3: density and cloud amount [Greuell & Konzelmann, 1994] - # 4: exponential time decay & wetness [Bougamont & Bamber, 2005] - # 5: ingest MODIS mode, direct input from aValue parameter applied to surface ice only - + # 0: direct input from aValue parameter + # 1: effective grain radius [Gardner & Sharp, 2009] + # 2: effective grain radius [Brun et al., 2009] + # 3: density and cloud amount [Greuell & Konzelmann, 1994] + # 4: exponential time decay & wetness [Bougamont & Bamber, 2005] + # 5: ingest MODIS mode, direct input from aValue parameter applied to surface ice only swIdx = float('NaN') # apply all SW to top grid cell (0) or allow SW to penetrate surface (1) (default 1) - denIdx = float('NaN') #densification model to use (default is 2): - # 1 = emperical model of Herron and Langway (1980) - # 2 = semi-emperical model of Anthern et al. (2010) - # 3 = DO NOT USE: physical model from Appix B of Anthern et al. (2010) - # 4 = DO NOT USE: emperical model of Li and Zwally (2004) - # 5 = DO NOT USE: modified emperical model (4) by Helsen et al. (2008) - # 6 = Antarctica semi-emperical model of Ligtenberg et al. (2011) - # 7 = Greenland semi-emperical model of Kuipers Munneke et al. (2015) - - dsnowIdx = float('NaN') #model for fresh snow accumulation density (default is 1): - # 0 = Original GEMB value, 150 kg/m^3 - # 1 = Antarctica value of fresh snow density, 350 kg/m^3 - # 2 = Greenland value of fresh snow density, 315 kg/m^3, Fausto et al. (2008) - # 3 = Antarctica model of Kaspers et al. (2004) - # 4 = Greenland model of Kuipers Munneke et al. (2015) - + # 1 = emperical model of Herron and Langway (1980) + # 2 = semi-emperical model of Anthern et al. (2010) + # 3 = DO NOT USE: physical model from Appix B of Anthern et al. (2010) + # 4 = DO NOT USE: emperical model of Li and Zwally (2004) + # 5 = DO NOT USE: modified emperical model (4) by Helsen et al. (2008) + # 6 = Antarctica semi-emperical model of Ligtenberg et al. (2011) + # 7 = Greenland semi-emperical model of Kuipers Munneke et al. (2015) + dsnowIdx = float('NaN') #model for fresh snow accumulation density (default is 1): + # 0 = Original GEMB value, 150 kg/m^3 + # 1 = Antarctica value of fresh snow density, 350 kg/m^3 + # 2 = Greenland value of fresh snow density, 315 kg/m^3, Fausto et al. (2008) + # 3 = Antarctica model of Kaspers et al. (2004) + # 4 = Greenland model of Kuipers Munneke et al. (2015) zTop = float('NaN') # depth over which grid length is constant at the top of the snopack (default 10) [m] - dzTop = float('NaN') # initial top vertical grid spacing (default .05) [m] - dzMin = float('NaN') # initial min vertical allowable grid spacing (default dzMin/2) [m] + dzTop = float('NaN') # initial top vertical grid spacing (default .05) [m] + dzMin = float('NaN') # initial min vertical allowable grid spacing (default dzMin/2) [m] zY = float('NaN') # strech grid cells bellow top_z by a [top_dz * y ^ (cells bellow top_z)] zMax = float('NaN') #initial max model depth (default is min(thickness,500)) [m] @@ -95,17 +90,17 @@ outputFreq = float('NaN') #output frequency in days (default is monthly, 30) - #specific albedo parameters: - #Method 1 and 2: + #specific albedo parameters: + #Method 1 and 2: aSnow = float('NaN') # new snow albedo (0.64 - 0.89) aIce = float('NaN') # range 0.27-0.58 for old snow - #Method 3: Radiation Correction Factors -> only used for met station data and Greuell & Konzelmann, 1994 albedo + #Method 3: Radiation Correction Factors -> only used for met station data and Greuell & Konzelmann, 1994 albedo cldFrac = float('NaN') # average cloud amount - #Method 4: additonal tuning parameters albedo as a funtion of age and water content (Bougamont et al., 2005) - t0wet = float('NaN') # time scale for wet snow (15-21.9) - t0dry = float('NaN') # warm snow timescale (30) - K = float('NaN') # time scale temperature coef. (7) + #Method 4: additonal tuning parameters albedo as a funtion of age and water content (Bougamont et al., 2005) + t0wet = float('NaN') # time scale for wet snow (15-21.9) + t0dry = float('NaN') # warm snow timescale (30) + K = float('NaN') # time scale temperature coef. (7) adThresh = float('NaN') # Apply aIdx method to all areas with densities below this value, - # or else apply direct input value from aValue, allowing albedo to be altered. - # Default value is rho water (1023 kg m-3). + # or else apply direct input value from aValue, allowing albedo to be altered. + # Default value is rho water (1023 kg m-3). #densities: @@ -114,10 +109,10 @@ #thermo: ThermoDeltaTScaling = float('NaN') #scaling factor to multiply the thermal diffusion timestep (delta t) - + requested_outputs = [] - #Several fields are missing from the standard GEMB model, which are capture intrinsically by ISSM. - #dateN: that's the last row of the above fields. - #dt: included in dateN. Not an input. + #Several fields are missing from the standard GEMB model, which are capture intrinsically by ISSM. + #dateN: that's the last row of the above fields. + #dt: included in dateN. Not an input. #elev: this is taken from the ISSM surface itself. @@ -128,6 +123,5 @@ #string = "#s\n#s"%(string,fielddisplay(self,'requested_outputs','additional outputs requested')) string = ' surface forcings for SMB GEMB model :' - - string="%s\n%s"%(string,fielddisplay(self,'issmbgradients','is smb gradients method activated (0 or 1, default is 0)')) + string = "%s\n%s"%(string,fielddisplay(self,'issmbgradients','is smb gradients method activated (0 or 1, default is 0)')) string = "%s\n%s"%(string,fielddisplay(self,'isgraingrowth','run grain growth module (default true)')) string = "%s\n%s"%(string,fielddisplay(self,'isalbedo','run albedo module (default true)')) @@ -147,5 +141,5 @@ string = "%s\n%s"%(string,fielddisplay(self,'Tmean','mean annual temperature [K]')) string = "%s\n%s"%(string,fielddisplay(self,'C','mean annual snow accumulation [kg m-2 yr-1]')) - string = "%s\n%s"%(string,fielddisplay(self,'Vmean','mean annual temperature [m s-1] (default 10 m/s)')) + string = "%s\n%s"%(string,fielddisplay(self,'Vmean','mean annual temperature [m s-1] (default 10 m/s)')) string = "%s\n%s"%(string,fielddisplay(self,'Tz','height above ground at which temperature (T) was sampled [m]')) string = "%s\n%s"%(string,fielddisplay(self,'Vz','height above ground at which wind (V) eas sampled [m]')) @@ -161,12 +155,10 @@ string = "%s\n%s"%(string,fielddisplay(self,'adThresh','Apply aIdx method to all areas with densities below this value,','or else apply direct input value from aValue, allowing albedo to be altered.')) string = "%s\n%s"%(string,fielddisplay(self,'aIdx',['method for calculating albedo and subsurface absorption (default is 1)', - '0: direct input from aValue parameter', - '1: effective grain radius [Gardner & Sharp, 2009]', - '2: effective grain radius [Brun et al., 2009]', - '3: density and cloud amount [Greuell & Konzelmann, 1994]', - '4: exponential time decay & wetness [Bougamont & Bamber, 2005]'])) - + '0: direct input from aValue parameter', + '1: effective grain radius [Gardner & Sharp, 2009]', + '2: effective grain radius [Brun et al., 2009]', + '3: density and cloud amount [Greuell & Konzelmann, 1994]', + '4: exponential time decay & wetness [Bougamont & Bamber, 2005]'])) string = "%s\n%s"%(string,fielddisplay(self,'teValue','Outward longwave radiation thermal emissivity forcing at every element (default in code is 1)')) - #snow properties init string = "%s\n%s"%(string,fielddisplay(self,'Dzini','Initial cell depth when restart [m]')) @@ -180,6 +172,6 @@ string = "%s\n%s"%(string,fielddisplay(self,'Tini','Initial snow temperature when restart [K]')) string = "%s\n%s"%(string,fielddisplay(self,'Sizeini','Initial number of layers when restart [K]')) - - #additional albedo parameters: + + #additional albedo parameters: if type(self.aIdx) == list or type(self.aIdx) == type(np.array([1,2])) and (self.aIdx == [1,2] or (1 in self.aIdx and 2 in self.aIdx)): string = "%s\n%s"%(string,fielddisplay(self,'aSnow','new snow albedo (0.64 - 0.89)')) @@ -195,22 +187,20 @@ string = "%s\n%s"%(string,fielddisplay(self,'t0dry','warm snow timescale (30) [d]')) string = "%s\n%s"%(string,fielddisplay(self,'K','time scale temperature coef. (7) [d]')) - + string = "%s\n%s"%(string,fielddisplay(self,'swIdx','apply all SW to top grid cell (0) or allow SW to penetrate surface (1) [default 1]')) string = "%s\n%s"%(string,fielddisplay(self,'denIdx',['densification model to use (default is 2):', - '1 = emperical model of Herron and Langway (1980)', - '2 = semi-emperical model of Anthern et al. (2010)', - '3 = DO NOT USE: physical model from Appix B of Anthern et al. (2010)', - '4 = DO NOT USE: emperical model of Li and Zwally (2004)', - '5 = DO NOT USE: modified emperical model (4) by Helsen et al. (2008)', - '6 = Antarctica semi-emperical model of Ligtenberg et al. (2011)', - '7 = Greenland semi-emperical model of Kuipers Munneke et al. (2015)'])) - - string = "%s\n%s"%(string,fielddisplay(self,'dsnowIdx',['model for fresh snow accumulation density (default is 1):', - '0 = Original GEMB value, 150 kg/m^3', - '1 = Antarctica value of fresh snow density, 350 kg/m^3', - '2 = Greenland value of fresh snow density, 315 kg/m^3, Fausto et al. (2008)', - '3 = Antarctica model of Kaspers et al. (2004), Make sure to set Vmean accurately', - '4 = Greenland model of Kuipers Munneke et al. (2015)'])); - + '1 = emperical model of Herron and Langway (1980)', + '2 = semi-emperical model of Anthern et al. (2010)', + '3 = DO NOT USE: physical model from Appix B of Anthern et al. (2010)', + '4 = DO NOT USE: emperical model of Li and Zwally (2004)', + '5 = DO NOT USE: modified emperical model (4) by Helsen et al. (2008)', + '6 = Antarctica semi-emperical model of Ligtenberg et al. (2011)', + '7 = Greenland semi-emperical model of Kuipers Munneke et al. (2015)'])) + string = "%s\n%s"%(string,fielddisplay(self,'dsnowIdx',['model for fresh snow accumulation density (default is 1):', + '0 = Original GEMB value, 150 kg/m^3', + '1 = Antarctica value of fresh snow density, 350 kg/m^3', + '2 = Greenland value of fresh snow density, 315 kg/m^3, Fausto et al. (2008)', + '3 = Antarctica model of Kaspers et al. (2004), Make sure to set Vmean accurately', + '4 = Greenland model of Kuipers Munneke et al. (2015)'])); string = "%s\n%s"%(string,fielddisplay(self,'requested_outputs','additional outputs requested')) return string @@ -247,9 +237,9 @@ self.isdensification = 1 self.isturbulentflux = 1 - + self.aIdx = 1 self.swIdx = 1 self.denIdx = 2 - self.dsnowIdx = 1 + self.dsnowIdx = 1 self.zTop = 10*np.ones((mesh.numberofelements,)) self.dzTop = .05* np.ones((mesh.numberofelements,)) @@ -258,6 +248,6 @@ self.ThermoDeltaTScaling = 1/11.0 - self.Vmean = 10*np.ones((mesh.numberofelements,)) - + self.Vmean = 10*np.ones((mesh.numberofelements,)) + self.zMax = 250*np.ones((mesh.numberofelements,)) self.zMin = 130*np.ones((mesh.numberofelements,)) @@ -268,5 +258,5 @@ self.aSnow = 0.85 self.aIce = 0.48 - self.cldFrac = 0.1 + self.cldFrac = 0.1 self.t0wet = 15 self.t0dry = 30 @@ -276,7 +266,7 @@ self.teValue = np.ones((mesh.numberofelements,)); self.aValue = self.aSnow*np.ones(mesh.numberofelements,); - + self.Dzini = 0.05*np.ones((mesh.numberofelements,2)) - self.Dini = 910.0*np.ones((mesh.numberofelements,2)) + self.Dini = 910.0*np.ones((mesh.numberofelements,2)) self.Reini = 2.5*np.ones((mesh.numberofelements,2)) self.Gdnini = 0.0*np.ones((mesh.numberofelements,2)) @@ -286,6 +276,6 @@ self.Aini = self.aSnow*np.ones((mesh.numberofelements,2)) self.Tini = 273.15*np.ones((mesh.numberofelements,2)) -# /!\ Default value of Tini must be equal to Tmean but don't know Tmean yet (computed when atmospheric forcings are interpolated on mesh). -# If initialization without restart, this value will be overwritten when snow parameters are retrieved in Element.cpp +# /!\ Default value of Tini must be equal to Tmean but don't know Tmean yet (computed when atmospheric forcings are interpolated on mesh). +# If initialization without restart, this value will be overwritten when snow parameters are retrieved in Element.cpp self.Sizeini = 2*np.ones((mesh.numberofelements,)) #}}} @@ -310,8 +300,8 @@ md = checkfield(md,'fieldname','smb.Tmean','size',[md.mesh.numberofelements],'NaN',1,'Inf',1,'>',273-100,'<',273+100) #-100/100 celsius min/max value - md = checkfield(md,'fieldname','smb.C','size',[md.mesh.numberofelements],'NaN',1,'Inf',1,'> = ',0) - md = checkfield(md,'fieldname','smb.Vmean','size',[md.mesh.numberofelements],'NaN',1,'Inf',1,'> = ',0) - md = checkfield(md,'fieldname','smb.Tz','size',[md.mesh.numberofelements],'NaN',1,'Inf',1,'> = ',0,'< = ',5000) - md = checkfield(md,'fieldname','smb.Vz','size',[md.mesh.numberofelements],'NaN',1,'Inf',1,'> = ',0,'< = ',5000) + md = checkfield(md,'fieldname','smb.C','size',[md.mesh.numberofelements],'NaN',1,'Inf',1,'> = ',0) + md = checkfield(md,'fieldname','smb.Vmean','size',[md.mesh.numberofelements],'NaN',1,'Inf',1,'> = ',0) + md = checkfield(md,'fieldname','smb.Tz','size',[md.mesh.numberofelements],'NaN',1,'Inf',1,'> = ',0,'< = ',5000) + md = checkfield(md,'fieldname','smb.Vz','size',[md.mesh.numberofelements],'NaN',1,'Inf',1,'> = ',0,'< = ',5000) md = checkfield(md,'fieldname','smb.teValue','timeseries',1,'NaN',1,'Inf',1,'>=',0,'<=',1); @@ -320,5 +310,5 @@ md = checkfield(md,'fieldname','smb.swIdx','NaN',1,'Inf',1,'values',[0,1]) md = checkfield(md,'fieldname','smb.denIdx','NaN',1,'Inf',1,'values',[1,2,3,4,5,6,7]) - md = checkfield(md,'fieldname','smb.dsnowIdx','NaN',1,'Inf',1,'values',[0,1,2,3,4]) + md = checkfield(md,'fieldname','smb.dsnowIdx','NaN',1,'Inf',1,'values',[0,1,2,3,4]) md = checkfield(md,'fieldname','smb.zTop','NaN',1,'Inf',1,'> = ',0) @@ -326,5 +316,5 @@ md = checkfield(md,'fieldname','smb.dzMin','NaN',1,'Inf',1,'>',0) md = checkfield(md,'fieldname','smb.zY','NaN',1,'Inf',1,'> = ',1) - md = checkfield(md,'fieldname','smb.outputFreq','NaN',1,'Inf',1,'>',0,'<',10*365) #10 years max + md = checkfield(md,'fieldname','smb.outputFreq','NaN',1,'Inf',1,'>',0,'<',10*365) #10 years max md = checkfield(md,'fieldname','smb.InitDensityScaling','NaN',1,'Inf',1,'> = ',0,'< = ',1) md = checkfield(md,'fieldname','smb.ThermoDeltaTScaling','NaN',1,'Inf',1,'> = ',0,'< = ',1) @@ -342,5 +332,5 @@ md = checkfield(md,'fieldname','smb.t0dry','NaN',1,'Inf',1,'> = ',30,'< = ',30) md = checkfield(md,'fieldname','smb.K','NaN',1,'Inf',1,'> = ',7,'< = ',7) - + #check zTop is < local thickness: @@ -348,5 +338,5 @@ if np.any(he= 1e-10): - error('smb_dt/dt = #f. The number of SMB time steps in one transient core time step has to be an an integer',md.timestepping.time_step/dt) - + error('smb_dt/dt = #f. The number of SMB time steps in one transient core time step has to be an an integer',md.timestepping.time_step/dt) + #process requested outputs outputs = self.requested_outputs @@ -435,6 +425,5 @@ outputscopy=outputs[0:max(0,indices[0]-1)]+self.defaultoutputs(md)+outputs[indices[0]+1:] outputs =outputscopy - + WriteData(fid,prefix,'data',outputs,'name','md.smb.requested_outputs','format','StringArray') # }}} - Index: /issm/trunk-jpl/src/py3/classes/bamggeom.py =================================================================== --- /issm/trunk-jpl/src/py3/classes/bamggeom.py (revision 23708) +++ /issm/trunk-jpl/src/py3/classes/bamggeom.py (revision 23709) @@ -25,5 +25,5 @@ elif len(args) == 1: object=args[0] - for field in object.keys(): + for field in list(object.keys()): if field in vars(self): setattr(self,field,object[field]) Index: /issm/trunk-jpl/src/py3/classes/bamgmesh.py =================================================================== --- /issm/trunk-jpl/src/py3/classes/bamgmesh.py (revision 23708) +++ /issm/trunk-jpl/src/py3/classes/bamgmesh.py (revision 23709) @@ -32,5 +32,5 @@ elif len(args) == 1: object=args[0] - for field in object.keys(): + for field in list(object.keys()): if field in vars(self): setattr(self,field,object[field]) Index: /issm/trunk-jpl/src/py3/classes/frictioncoulomb.py =================================================================== --- /issm/trunk-jpl/src/py3/classes/frictioncoulomb.py (revision 23708) +++ /issm/trunk-jpl/src/py3/classes/frictioncoulomb.py (revision 23709) @@ -5,80 +5,80 @@ class frictioncoulomb(object): - """ - FRICTIONCOULOMB class definition + """ + FRICTIONCOULOMB class definition - Usage: - frictioncoulomb=frictioncoulomb() - """ + Usage: + frictioncoulomb=frictioncoulomb() + """ - def __init__(self): # {{{ - self.coefficient = float('NaN') - self.coefficientcoulomb = float('NaN') - self.p = float('NaN') - self.q = float('NaN') - self.coupling = 0 - self.effective_pressure = float('NaN') - #set defaults - self.setdefaultparameters() + def __init__(self): # {{{ + self.coefficient = float('NaN') + self.coefficientcoulomb = float('NaN') + self.p = float('NaN') + self.q = float('NaN') + self.coupling = 0 + self.effective_pressure = float('NaN') + #set defaults + self.setdefaultparameters() + #}}} - #}}} - def __repr__(self): # {{{ - string="Basal shear stress parameters: Sigma_b = min(coefficient^2 * Neff ^r * |u_b|^(s-1) * u_b,\n coefficientcoulomb^2 * rho_i * g * (h-h_f)), (effective stress Neff=rho_ice*g*thickness+rho_water*g*bed, r=q/p and s=1/p)." + def __repr__(self): # {{{ + string="Basal shear stress parameters: Sigma_b = min(coefficient^2 * Neff ^r * |u_b|^(s-1) * u_b,\n coefficientcoulomb^2 * rho_i * g * (h-h_f)), (effective stress Neff=rho_ice*g*thickness+rho_water*g*bed, r=q/p and s=1/p)." + string="%s\n%s"%(string,fielddisplay(self,"coefficient","power law (Weertman) friction coefficient [SI]")) + string="%s\n%s"%(string,fielddisplay(self,"coefficientcoulomb","Coulomb friction coefficient [SI]")) + string="%s\n%s"%(string,fielddisplay(self,"p","p exponent")) + string="%s\n%s"%(string,fielddisplay(self,"q","q exponent")) + string="%s\n%s"%(string,fielddisplay(self,'coupling','Coupling flag: 0 for default, 1 for forcing(provide md.friction.effective_pressure) and 2 for coupled(not implemented yet)')) + string="%s\n%s"%(string,fielddisplay(self,'effective_pressure','Effective Pressure for the forcing if not coupled [Pa]')) + return string + #}}} + def extrude(self,md): # {{{ + self.coefficient=project3d(md,'vector',self.coefficient,'type','node','layer',1) + self.coefficientcoulomb=project3d(md,'vector',self.coefficientcoulomb,'type','node','layer',1) + self.p=project3d(md,'vector',self.p,'type','element') + self.q=project3d(md,'vector',self.q,'type','element') + if self.coupling==1: + self.effective_pressure=project3d(md,'vector',self.effective_pressure,'type','node','layer',1) + elif self.coupling==2: + raise ValueError('coupling not supported yet') + elif self.coupling > 2: + raise ValueError('md.friction.coupling larger than 2, not supported yet') + return self + #}}} - string="%s\n%s"%(string,fielddisplay(self,"coefficient","power law (Weertman) friction coefficient [SI]")) - string="%s\n%s"%(string,fielddisplay(self,"coefficientcoulomb","Coulomb friction coefficient [SI]")) - string="%s\n%s"%(string,fielddisplay(self,"p","p exponent")) - string="%s\n%s"%(string,fielddisplay(self,"q","q exponent")) - string="%s\n%s"%(string,fielddisplay(self,'coupling','Coupling flag: 0 for default, 1 for forcing(provide md.friction.effective_pressure) and 2 for coupled(not implemented yet)')) - string="%s\n%s"%(string,fielddisplay(self,'effective_pressure','Effective Pressure for the forcing if not coupled [Pa]')) - return string - #}}} - def extrude(self,md): # {{{ - self.coefficient=project3d(md,'vector',self.coefficient,'type','node','layer',1) - self.coefficientcoulomb=project3d(md,'vector',self.coefficientcoulomb,'type','node','layer',1) - self.p=project3d(md,'vector',self.p,'type','element') - self.q=project3d(md,'vector',self.q,'type','element') - if self.coupling==1: - self.effective_pressure=project3d(md,'vector',self.effective_pressure,'type','node','layer',1) - elif self.coupling==2: - raise ValueError('coupling not supported yet') - elif self.coupling > 2: - raise ValueError('md.friction.coupling larger than 2, not supported yet') - return self - #}}} - def setdefaultparameters(self): # {{{ - return self - #}}} - def checkconsistency(self,md,solution,analyses): # {{{ + def setdefaultparameters(self): # {{{ + return self + #}}} - #Early return - if 'StressbalanceAnalysis' not in analyses and 'ThermalAnalysis' not in analyses: - return md + def checkconsistency(self,md,solution,analyses): # {{{ + #Early return + if 'StressbalanceAnalysis' not in analyses and 'ThermalAnalysis' not in analyses: + return md - md = checkfield(md,'fieldname','friction.coefficient','timeseries',1,'NaN',1,'Inf',1) - md = checkfield(md,'fieldname','friction.coefficientcoulomb','timeseries',1,'NaN',1,'Inf',1) - md = checkfield(md,'fieldname','friction.q','NaN',1,'Inf',1,'size',[md.mesh.numberofelements]) - md = checkfield(md,'fieldname','friction.p','NaN',1,'Inf',1,'size',[md.mesh.numberofelements]) - if self.coupling==1: - md = checkfield(md,'fieldname','friction.effective_pressure','NaN',1,'Inf',1,'timeseries',1) - elif self.coupling==2: - raise ValueError('coupling not supported yet') - elif self.coupling > 2: - raise ValueError('md.friction.coupling larger than 2, not supported yet') - return md + md = checkfield(md,'fieldname','friction.coefficient','timeseries',1,'NaN',1,'Inf',1) + md = checkfield(md,'fieldname','friction.coefficientcoulomb','timeseries',1,'NaN',1,'Inf',1) + md = checkfield(md,'fieldname','friction.q','NaN',1,'Inf',1,'size',[md.mesh.numberofelements]) + md = checkfield(md,'fieldname','friction.p','NaN',1,'Inf',1,'size',[md.mesh.numberofelements]) + if self.coupling==1: + md = checkfield(md,'fieldname','friction.effective_pressure','NaN',1,'Inf',1,'timeseries',1) + elif self.coupling==2: + raise ValueError('coupling not supported yet') + elif self.coupling > 2: + raise ValueError('md.friction.coupling larger than 2, not supported yet') + return md + # }}} + def marshall(self,prefix,md,fid): # {{{ + WriteData(fid,prefix,'name','md.friction.law','data',7,'format','Integer') + WriteData(fid,prefix,'object',self,'fieldname','coefficient','format','DoubleMat','mattype',1,'timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts) + WriteData(fid,prefix,'object',self,'fieldname','coefficientcoulomb','format','DoubleMat','mattype',1) + WriteData(fid,prefix,'object',self,'fieldname','p','format','DoubleMat','mattype',2) + WriteData(fid,prefix,'object',self,'fieldname','q','format','DoubleMat','mattype',2) + WriteData(fid,prefix,'class','friction','object',self,'fieldname','coupling','format','Integer') + if self.coupling==1: + WriteData(fid,prefix,'class','friction','object',self,'fieldname','effective_pressure','format','DoubleMat','mattype',1,'timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts) + elif self.coupling==2: + raise ValueError('coupling not supported yet') + elif self.coupling > 2: + raise ValueError('md.friction.coupling larger than 2, not supported yet') # }}} - def marshall(self,prefix,md,fid): # {{{ - WriteData(fid,prefix,'name','md.friction.law','data',7,'format','Integer') - WriteData(fid,prefix,'object',self,'fieldname','coefficient','format','DoubleMat','mattype',1,'timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts) - WriteData(fid,prefix,'object',self,'fieldname','coefficientcoulomb','format','DoubleMat','mattype',1) - WriteData(fid,prefix,'object',self,'fieldname','p','format','DoubleMat','mattype',2) - WriteData(fid,prefix,'object',self,'fieldname','q','format','DoubleMat','mattype',2) - WriteData(fid,prefix,'class','friction','object',self,'fieldname','coupling','format','Integer') - if self.coupling==1: - WriteData(fid,prefix,'class','friction','object',self,'fieldname','effective_pressure','format','DoubleMat','mattype',1,'timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts) - elif self.coupling==2: - raise ValueError('coupling not supported yet') - elif self.coupling > 2: - raise ValueError('md.friction.coupling larger than 2, not supported yet') - # }}} Index: /issm/trunk-jpl/src/py3/classes/frictionshakti.py =================================================================== --- /issm/trunk-jpl/src/py3/classes/frictionshakti.py (revision 23708) +++ /issm/trunk-jpl/src/py3/classes/frictionshakti.py (revision 23709) @@ -5,43 +5,44 @@ class frictionshakti(object): - """ - FRICTIONSHAKTI class definition + """ + FRICTIONSHAKTI class definition - Usage: - friction=frictionshakti() - """ + Usage: + friction=frictionshakti() + """ - def __init__(self,md): # {{{ - self.coefficient = md.friction.coefficient - #set defaults - self.setdefaultparameters() + def __init__(self,md): # {{{ + self.coefficient = md.friction.coefficient + #set defaults + self.setdefaultparameters() + #}}} - #}}} - def __repr__(self): # {{{ - string="Basal shear stress parameters: Sigma_b = coefficient^2 * Neff * u_b\n(effective stress Neff=rho_ice*g*thickness+rho_water*g*(head-b))" + def __repr__(self): # {{{ + string="Basal shear stress parameters: Sigma_b = coefficient^2 * Neff * u_b\n(effective stress Neff=rho_ice*g*thickness+rho_water*g*(head-b))" + string="%s\n%s"%(string,fielddisplay(self,"coefficient","friction coefficient [SI]")) + return string + #}}} - string="%s\n%s"%(string,fielddisplay(self,"coefficient","friction coefficient [SI]")) - return string - #}}} - def extrude(self,md): # {{{ - self.coefficient=project3d(md,'vector',self.coefficient,'type','node','layer',1) - return self - #}}} - def setdefaultparameters(self): # {{{ - return self - #}}} - def checkconsistency(self,md,solution,analyses): # {{{ + def extrude(self,md): # {{{ + self.coefficient=project3d(md,'vector',self.coefficient,'type','node','layer',1) + return self + #}}} - #Early return - if 'StressbalanceAnalysis' not in analyses and 'ThermalAnalysis' not in analyses: - return md + def setdefaultparameters(self): # {{{ + return self + #}}} - md = checkfield(md,'fieldname','friction.coefficient','timeseries',1,'NaN',1,'Inf',1) - return md + def checkconsistency(self,md,solution,analyses): # {{{ + #Early return + if 'StressbalanceAnalysis' not in analyses and 'ThermalAnalysis' not in analyses: + return md + md = checkfield(md,'fieldname','friction.coefficient','timeseries',1,'NaN',1,'Inf',1) + return md + # }}} + + def marshall(self,prefix,md,fid): # {{{ + yts=md.constants.yts + WriteData(fid,prefix,'name','md.friction.law','data',8,'format','Integer') + WriteData(fid,prefix,'object',self,'fieldname','coefficient','format','DoubleMat','mattype',1,'timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts) # }}} - def marshall(self,prefix,md,fid): # {{{ - yts=md.constants.yts - WriteData(fid,prefix,'name','md.friction.law','data',8,'format','Integer') - WriteData(fid,prefix,'object',self,'fieldname','coefficient','format','DoubleMat','mattype',1,'timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts) - # }}} Index: /issm/trunk-jpl/src/py3/classes/hydrologyshakti.py =================================================================== --- /issm/trunk-jpl/src/py3/classes/hydrologyshakti.py (revision 23708) +++ /issm/trunk-jpl/src/py3/classes/hydrologyshakti.py (revision 23709) @@ -30,5 +30,4 @@ #}}} def __repr__(self): # {{{ - string=' hydrologyshakti solution parameters:' string="%s\n%s"%(string,fielddisplay(self,'head','subglacial hydrology water head (m)')) @@ -45,10 +44,12 @@ string="%s\n%s"%(string,fielddisplay(self,'requested_outputs','additional outputs requested')) return string - #}}} + #}}} + def extrude(self,md): # {{{ return self #}}} + def setdefaultparameters(self): # {{{ - # Set under-relaxation parameter to be 1 (no under-relaxation of nonlinear iteration) + # Set under-relaxation parameter to be 1 (no under-relaxation of nonlinear iteration) self.relaxation=1 self.storage=0 @@ -56,4 +57,5 @@ return self #}}} + def defaultoutputs(self,md): # {{{ list = ['HydrologyHead','HydrologyGapHeight','EffectivePressure','HydrologyBasalFlux','DegreeOfChannelization'] @@ -62,5 +64,4 @@ def checkconsistency(self,md,solution,analyses): # {{{ - #Early return if 'HydrologyShaktiAnalysis' not in analyses: @@ -76,10 +77,10 @@ md = checkfield(md,'fieldname','hydrology.neumannflux','timeseries',1,'NaN',1,'Inf',1) md = checkfield(md,'fieldname','hydrology.spchead','size',[md.mesh.numberofvertices]) - md = checkfield(md,'fieldname','hydrology.relaxation','>=',0) + md = checkfield(md,'fieldname','hydrology.relaxation','>=',0) md = checkfield(md,'fieldname','hydrology.storage','>=',0) md = checkfield(md,'fieldname','hydrology.requested_outputs','stringrow',1) - return md # }}} + def marshall(self,prefix,md,fid): # {{{ yts=md.constants.yts Index: /issm/trunk-jpl/src/py3/classes/massfluxatgate.py =================================================================== --- /issm/trunk-jpl/src/py3/classes/massfluxatgate.py (revision 23708) +++ /issm/trunk-jpl/src/py3/classes/massfluxatgate.py (revision 23709) @@ -46,19 +46,19 @@ #}}} def checkconsistency(self,md,solution,analyses): # {{{ - + if not isinstance(self.name, str): raise RuntimeError("massfluxatgate error message: 'name' field should be a string!") - + if not isinstance(self.profilename, str): - raise RuntimeError("massfluxatgate error message: 'profilename' field should be a string!") - + raise RuntimeError("massfluxatgate error message: 'profilename' field should be a string!") + OutputdefinitionStringArray=[] for i in range(1,100): x='Outputdefinition'+str(i) OutputdefinitionStringArray.append(x) - + md = checkfield(md,'field',self.definitionstring,'values',OutputdefinitionStringArray) - - #check the profilename points to a file!: + + #check the profilename points to a file!: if not os.path.isfile(self.profilename): raise RuntimeError("massfluxatgate error message: file name for profile corresponding to gate does not point to a legitimate file on disk!") @@ -67,9 +67,9 @@ # }}} def marshall(self,prefix,md,fid): # {{{ - - #before marshalling, we need to create the segments out of the profilename: + + #before marshalling, we need to create the segments out of the profilename: self.segments=MeshProfileIntersection(md.mesh.elements,md.mesh.x,md.mesh.y,self.profilename)[0] - #ok, marshall name and segments: + #ok, marshall name and segments: WriteData(fid,prefix,'data',self.name,'name','md.massfluxatgate.name','format','String'); WriteData(fid,prefix,'data',self.definitionstring,'name','md.massfluxatgate.definitionstring','format','String'); Index: /issm/trunk-jpl/src/py3/classes/matestar.py =================================================================== --- /issm/trunk-jpl/src/py3/classes/matestar.py (revision 23708) +++ /issm/trunk-jpl/src/py3/classes/matestar.py (revision 23709) @@ -14,23 +14,22 @@ def __init__(self): # {{{ - - rho_ice = 0. - rho_water = 0. - rho_freshwater = 0. - mu_water = 0. - heatcapacity = 0. - latentheat = 0. - thermalconductivity = 0. - temperateiceconductivity = 0. - meltingpoint = 0. - beta = 0. - mixed_layer_capacity = 0. - thermal_exchange_velocity = 0. - rheology_B = float('NaN') - rheology_Ec = float('NaN') - rheology_Es = float('NaN') - rheology_law = '' + rho_ice = 0. + rho_water = 0. + rho_freshwater = 0. + mu_water = 0. + heatcapacity = 0. + latentheat = 0. + thermalconductivity = 0. + temperateiceconductivity = 0. + meltingpoint = 0. + beta = 0. + mixed_layer_capacity = 0. + thermal_exchange_velocity = 0. + rheology_B = float('NaN') + rheology_Ec = float('NaN') + rheology_Es = float('NaN') + rheology_law = '' - #giaivins: + #giaivins: lithosphere_shear_modulus = 0. lithosphere_density = 0. @@ -41,10 +40,10 @@ earth_density = 0 - #set default parameters: + #set default parameters: self.setdefaultparameters() #}}} + def __repr__(self): # {{{ string=" Materials:" - string="%s\n%s"%(string,fielddisplay(self,'rho_ice','ice density [kg/m^3]')) string="%s\n%s"%(string,fielddisplay(self,'rho_water','ocean water density [kg/m^3]')) @@ -68,54 +67,42 @@ string="%s\n%s"%(string,fielddisplay(self,'mantle_density','Mantle density [g/cm^-3]')) string="%s\n%s"%(string,fielddisplay(self,'earth_density','Mantle density [kg/m^-3]')) - return string #}}} + def extrude(self,md): # {{{ self.rheology_B=project3d(md,'vector',self.rheology_B,'type','node') self.rheology_Ec=project3d(md,'vector',self.rheology_Ec,'type','node') self.rheology_Es=project3d(md,'vector',self.rheology_Es,'type','node') - return self + return self #}}} + def setdefaultparameters(self): # {{{ #ice density (kg/m^3) self.rho_ice=917. - #ocean water density (kg/m^3) self.rho_water=1023. - #fresh water density (kg/m^3) self.rho_freshwater=1000. - #water viscosity (N.s/m^2) - self.mu_water=0.001787 - + self.mu_water=0.001787 #ice heat capacity cp (J/kg/K) self.heatcapacity=2093. - #ice latent heat of fusion L (J/kg) self.latentheat=3.34*10**5 - #ice thermal conductivity (W/m/K) self.thermalconductivity=2.4 - #wet ice thermal conductivity (W/m/K) self.temperateiceconductivity=.24 - #the melting point of ice at 1 atmosphere of pressure in K self.meltingpoint=273.15 - #rate of change of melting point with pressure (K/Pa) self.beta=9.8*10**-8 - #mixed layer (ice-water interface) heat capacity (J/kg/K) self.mixed_layer_capacity=3974. - #thermal exchange velocity (ice-water interface) (m/s) self.thermal_exchange_velocity=1.00*10**-4 - #Rheology law: what is the temperature dependence of B with T #available: none, paterson and arrhenius self.rheology_law='Paterson' - # GIA: self.lithosphere_shear_modulus = 6.7*10**10 # (Pa) @@ -123,5 +110,4 @@ self.mantle_shear_modulus = 1.45*10**11 # (Pa) self.mantle_density = 3.34 # (g/cm^-3) - #SLR self.earth_density= 5512 # average density of the Earth, (kg/m^3) @@ -129,4 +115,5 @@ return self #}}} + def checkconsistency(self,md,solution,analyses): # {{{ md = checkfield(md,'fieldname','materials.rho_ice','>',0) @@ -144,4 +131,5 @@ md = checkfield(md,'fieldname','materials.mantle_shear_modulus','>',0,'numel',1) md = checkfield(md,'fieldname','materials.mantle_density','>',0,'numel',1) + if 'SealevelriseAnalysis' in analyses: md = checkfield(md,'fieldname','materials.earth_density','>',0,'numel',1) @@ -149,4 +137,5 @@ return md # }}} + def marshall(self,prefix,md,fid): # {{{ WriteData(fid,prefix,'name','md.materials.type','data',2,'format','Integer') @@ -167,5 +156,4 @@ WriteData(fid,prefix,'object',self,'class','materials','fieldname','rheology_Es','format','DoubleMat','mattype',1) WriteData(fid,prefix,'data',self.rheology_law,'name','md.materials.rheology_law','format','String') - WriteData(fid,prefix,'object',self,'class','materials','fieldname','lithosphere_shear_modulus','format','Double') WriteData(fid,prefix,'object',self,'class','materials','fieldname','lithosphere_density','format','Double','scale',10^3) Index: /issm/trunk-jpl/src/py3/classes/mismipbasalforcings.py =================================================================== --- /issm/trunk-jpl/src/py3/classes/mismipbasalforcings.py (revision 23708) +++ /issm/trunk-jpl/src/py3/classes/mismipbasalforcings.py (revision 23709) @@ -6,90 +6,82 @@ class mismipbasalforcings(object): - """ - MISMIP Basal Forcings class definition + """ + MISMIP Basal Forcings class definition - Usage: - mismipbasalforcings=mismipbasalforcings() - """ + Usage: + mismipbasalforcings=mismipbasalforcings() + """ - def __init__(self): # {{{ + def __init__(self): # {{{ + self.groundedice_melting_rate = float('NaN') + self.meltrate_factor = float('NaN') + self.threshold_thickness = float('NaN') + self.upperdepth_melt = float('NaN') + self.geothermalflux = float('NaN') + self.setdefaultparameters() - self.groundedice_melting_rate = float('NaN') - self.meltrate_factor = float('NaN') - self.threshold_thickness = float('NaN') - self.upperdepth_melt = float('NaN') - self.geothermalflux = float('NaN') - - self.setdefaultparameters() - - #}}} - def __repr__(self): # {{{ - string=" MISMIP+ basal melt parameterization\n" - string="%s\n%s"%(string,fielddisplay(self,"groundedice_melting_rate","basal melting rate (positive if melting) [m/yr]")) - string="%s\n%s"%(string,fielddisplay(self,"meltrate_factor","Melt-rate rate factor [1/yr] (sign is opposite to MISMIP+ benchmark to remain consistent with ISSM convention of positive values for melting)")) - string="%s\n%s"%(string,fielddisplay(self,"threshold_thickness","Threshold thickness for saturation of basal melting [m]")) - string="%s\n%s"%(string,fielddisplay(self,"upperdepth_melt","Depth above which melt rate is zero [m]")) - string="%s\n%s"%(string,fielddisplay(self,"geothermalflux","Geothermal heat flux [W/m^2]")) - return string - #}}} - def extrude(self,md): # {{{ - self.groundedice_melting_rate=project3d(md,'vector',self.groundedice_melting_rate,'type','node','layer',1) - self.geothermalflux=project3d(md,'vector',self.geothermalflux,'type','node','layer',1) #bedrock only gets geothermal flux - return self - #}}} - def initialize(self,md): # {{{ - if np.all(np.isnan(self.groundedice_melting_rate)): - self.groundedice_melting_rate=np.zeros((md.mesh.numberofvertices)) - print(' no basalforcings.groundedice_melting_rate specified: values set as zero') - if np.all(np.isnan(self.geothermalflux)): + #}}} + def __repr__(self): # {{{ + string=" MISMIP+ basal melt parameterization\n" + string="%s\n%s"%(string,fielddisplay(self,"groundedice_melting_rate","basal melting rate (positive if melting) [m/yr]")) + string="%s\n%s"%(string,fielddisplay(self,"meltrate_factor","Melt-rate rate factor [1/yr] (sign is opposite to MISMIP+ benchmark to remain consistent with ISSM convention of positive values for melting)")) + string="%s\n%s"%(string,fielddisplay(self,"threshold_thickness","Threshold thickness for saturation of basal melting [m]")) + string="%s\n%s"%(string,fielddisplay(self,"upperdepth_melt","Depth above which melt rate is zero [m]")) + string="%s\n%s"%(string,fielddisplay(self,"geothermalflux","Geothermal heat flux [W/m^2]")) + return string + #}}} + def extrude(self,md): # {{{ + self.groundedice_melting_rate=project3d(md,'vector',self.groundedice_melting_rate,'type','node','layer',1) + self.geothermalflux=project3d(md,'vector',self.geothermalflux,'type','node','layer',1) #bedrock only gets geothermal flux + return self + #}}} + def initialize(self,md): # {{{ + if np.all(np.isnan(self.groundedice_melting_rate)): + self.groundedice_melting_rate=np.zeros((md.mesh.numberofvertices)) + print(' no basalforcings.groundedice_melting_rate specified: values set as zero') + if np.all(np.isnan(self.geothermalflux)): self.geothermalflux=np.zeros((md.mesh.numberofvertices)) print(" no basalforcings.geothermalflux specified: values set as zero") - return self - #}}} - def setdefaultparameters(self): # {{{ - # default values for melting parameterization - self.meltrate_factor = 0.2 - self.threshold_thickness = 75. - self.upperdepth_melt = -100. - return self - #}}} - def checkconsistency(self,md,solution,analyses): # {{{ + return self + #}}} + def setdefaultparameters(self): # {{{ + # default values for melting parameterization + self.meltrate_factor = 0.2 + self.threshold_thickness = 75. + self.upperdepth_melt = -100. + return self + #}}} + def checkconsistency(self,md,solution,analyses): # {{{ + #Early return + if 'MasstransportAnalysis' in analyses and not (solution=='TransientSolution' and md.transient.ismasstransport==0): + md = checkfield(md,'fieldname','basalforcings.groundedice_melting_rate','NaN',1,'Inf',1,'timeseries',1) + md = checkfield(md,'fieldname','basalforcings.meltrate_factor','>=',0,'numel',[1]) + md = checkfield(md,'fieldname','basalforcings.threshold_thickness','>=',0,'numel',[1]) + md = checkfield(md,'fieldname','basalforcings.upperdepth_melt','<=',0,'numel',[1]) - #Early return - if 'MasstransportAnalysis' in analyses and not (solution=='TransientSolution' and md.transient.ismasstransport==0): + if 'BalancethicknessAnalysis' in analyses: + md = checkfield(md,'fieldname','basalforcings.groundedice_melting_rate','NaN',1,'Inf',1,'size',[md.mesh.numberofvertices]) + md = checkfield(md,'fieldname','basalforcings.meltrate_factor','>=',0,'numel',[1]) + md = checkfield(md,'fieldname','basalforcings.threshold_thickness','>=',0,'numel',[1]) + md = checkfield(md,'fieldname','basalforcings.upperdepth_melt','<=',0,'numel',[1]) - md = checkfield(md,'fieldname','basalforcings.groundedice_melting_rate','NaN',1,'Inf',1,'timeseries',1) - md = checkfield(md,'fieldname','basalforcings.meltrate_factor','>=',0,'numel',[1]) - md = checkfield(md,'fieldname','basalforcings.threshold_thickness','>=',0,'numel',[1]) - md = checkfield(md,'fieldname','basalforcings.upperdepth_melt','<=',0,'numel',[1]) + if 'ThermalAnalysis' in analyses and not (solution=='TransientSolution' and md.transient.isthermal==0): + md = checkfield(md,'fieldname','basalforcings.groundedice_melting_rate','NaN',1,'Inf',1,'timeseries',1) + md = checkfield(md,'fieldname','basalforcings.meltrate_factor','>=',0,'numel',[1]) + md = checkfield(md,'fieldname','basalforcings.threshold_thickness','>=',0,'numel',[1]) + md = checkfield(md,'fieldname','basalforcings.upperdepth_melt','<=',0,'numel',[1]) + md = checkfield(md,'fieldname','basalforcings.geothermalflux','NaN',1,'Inf',1,'timeseries',1,'>=',0) + return md + # }}} + def marshall(self,prefix,md,fid): # {{{ + yts=md.constants.yts + if yts!=365.2422*24.*3600.: + print('WARNING: value of yts for MISMIP+ runs different from ISSM default!') - if 'BalancethicknessAnalysis' in analyses: - - md = checkfield(md,'fieldname','basalforcings.groundedice_melting_rate','NaN',1,'Inf',1,'size',[md.mesh.numberofvertices]) - md = checkfield(md,'fieldname','basalforcings.meltrate_factor','>=',0,'numel',[1]) - md = checkfield(md,'fieldname','basalforcings.threshold_thickness','>=',0,'numel',[1]) - md = checkfield(md,'fieldname','basalforcings.upperdepth_melt','<=',0,'numel',[1]) - - if 'ThermalAnalysis' in analyses and not (solution=='TransientSolution' and md.transient.isthermal==0): - - md = checkfield(md,'fieldname','basalforcings.groundedice_melting_rate','NaN',1,'Inf',1,'timeseries',1) - md = checkfield(md,'fieldname','basalforcings.meltrate_factor','>=',0,'numel',[1]) - md = checkfield(md,'fieldname','basalforcings.threshold_thickness','>=',0,'numel',[1]) - md = checkfield(md,'fieldname','basalforcings.upperdepth_melt','<=',0,'numel',[1]) - md = checkfield(md,'fieldname','basalforcings.geothermalflux','NaN',1,'Inf',1,'timeseries',1,'>=',0) - return md + WriteData(fid,prefix,'name','md.basalforcings.model','data',3,'format','Integer') + WriteData(fid,prefix,'object',self,'fieldname','groundedice_melting_rate','format','DoubleMat','name','md.basalforcings.groundedice_melting_rate','mattype',1,'scale',1./yts,'timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts) + WriteData(fid,prefix,'object',self,'fieldname','geothermalflux','name','md.basalforcings.geothermalflux','format','DoubleMat','mattype',1,'timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts) + WriteData(fid,prefix,'object',self,'fieldname','meltrate_factor','format','Double','scale',1./yts) + WriteData(fid,prefix,'object',self,'fieldname','threshold_thickness','format','Double') + WriteData(fid,prefix,'object',self,'fieldname','upperdepth_melt','format','Double') # }}} - def marshall(self,prefix,md,fid): # {{{ - - yts=md.constants.yts - if yts!=365.2422*24.*3600.: - print('WARNING: value of yts for MISMIP+ runs different from ISSM default!') - - WriteData(fid,prefix,'name','md.basalforcings.model','data',3,'format','Integer') - WriteData(fid,prefix,'object',self,'fieldname','groundedice_melting_rate','format','DoubleMat','name','md.basalforcings.groundedice_melting_rate','mattype',1,'scale',1./yts,'timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts) - WriteData(fid,prefix,'object',self,'fieldname','geothermalflux','name','md.basalforcings.geothermalflux','format','DoubleMat','mattype',1,'timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts) - WriteData(fid,prefix,'object',self,'fieldname','meltrate_factor','format','Double','scale',1./yts) - WriteData(fid,prefix,'object',self,'fieldname','threshold_thickness','format','Double') - WriteData(fid,prefix,'object',self,'fieldname','upperdepth_melt','format','Double') - - # }}} Index: /issm/trunk-jpl/src/py3/classes/model.py =================================================================== --- /issm/trunk-jpl/src/py3/classes/model.py (revision 23708) +++ /issm/trunk-jpl/src/py3/classes/model.py (revision 23709) @@ -218,5 +218,5 @@ # }}} def checkmessage(self,string): # {{{ - print("model not consistent: ", string) + print(("model not consistent: ", string)) self.private.isconsistent=False return self @@ -419,5 +419,5 @@ md2.mesh.elementconnectivity=ElementConnectivity(md2.mesh.elements2d,md2.mesh.vertexconnectivity)[0] segments=contourenvelope(md2) - md2.mesh.vertexonboundary=np.zeros(numberofvertices2/md2.mesh.numberoflayers,bool) + md2.mesh.vertexonboundary=np.zeros(int(numberofvertices2/md2.mesh.numberoflayers),bool) md2.mesh.vertexonboundary[segments[:,0:2]-1]=True md2.mesh.vertexonboundary=np.tile(md2.mesh.vertexonboundary,md2.mesh.numberoflayers) @@ -449,5 +449,5 @@ if md1.results: md2.results=results() - for solutionfield,field in md1.results.__dict__.items(): + for solutionfield,field in list(md1.results.__dict__.items()): if isinstance(field,list): setattr(md2.results,solutionfield,[]) @@ -458,5 +458,5 @@ fieldr=getattr(md2.results,solutionfield)[i] #get subfields - for solutionsubfield,subfield in fieldi.__dict__.items(): + for solutionsubfield,subfield in list(fieldi.__dict__.items()): if np.size(subfield)==numberofvertices1: setattr(fieldr,solutionsubfield,subfield[pos_node]) @@ -472,5 +472,5 @@ fieldr=getattr(md2.results,solutionfield) #get subfields - for solutionsubfield,subfield in field.__dict__.items(): + for solutionsubfield,subfield in list(field.__dict__.items()): if np.size(subfield)==numberofvertices1: setattr(fieldr,solutionsubfield,subfield[pos_node]) @@ -482,5 +482,5 @@ #OutputDefinitions fields if md1.outputdefinition.definitions: - for solutionfield,field in md1.outputdefinition.__dict__.items(): + for solutionfield,field in list(md1.outputdefinition.__dict__.items()): if isinstance(field,list): #get each definition @@ -489,5 +489,5 @@ fieldr=getattr(md2.outputdefinition,solutionfield)[i] #get subfields - for solutionsubfield,subfield in fieldi.__dict__.items(): + for solutionsubfield,subfield in list(fieldi.__dict__.items()): if np.size(subfield)==numberofvertices1: setattr(fieldr,solutionsubfield,subfield[pos_node]) @@ -834,5 +834,5 @@ #OutputDefinitions if md.outputdefinition.definitions: - for solutionfield,field in md.outputdefinition.__dict__.items(): + for solutionfield,field in list(md.outputdefinition.__dict__.items()): if isinstance(field,list): #get each definition @@ -841,5 +841,5 @@ fieldr=getattr(md.outputdefinition,solutionfield)[i] #get subfields - for solutionsubfield,subfield in fieldi.__dict__.items(): + for solutionsubfield,subfield in list(fieldi.__dict__.items()): if np.size(subfield)==md.mesh.numberofvertices: setattr(fieldr,solutionsubfield,project2d(md,subfield,1)) Index: /issm/trunk-jpl/src/py3/classes/organizer.py =================================================================== --- /issm/trunk-jpl/src/py3/classes/organizer.py (revision 23708) +++ /issm/trunk-jpl/src/py3/classes/organizer.py (revision 23709) @@ -6,5 +6,5 @@ from savevars import savevars from model import model -from dbm import whichdb +from dbm.ndbm import whichdb import MatlabFuncs as m @@ -86,5 +86,5 @@ if os.path.exists(path): struc=loadvars(path) - name=name=[key for key in struc.keys()] + name=name=[key for key in list(struc.keys())] md=struc.name[0] else: @@ -115,5 +115,5 @@ raise IOError("Could find neither '%s' nor '%s'" % (path,path2)) else: - print("--> Branching '%s' from trunk '%s'" % (self.prefix,self.trunkprefix)) + print(("--> Branching '%s' from trunk '%s'" % (self.prefix,self.trunkprefix))) md=loadmodel(path2) return md @@ -142,10 +142,10 @@ if 0 in self.requestedsteps: if self._currentstep==1: - print(" prefix: %s" % self.prefix) - print(" step #%i : %s" % (self.steps[self._currentstep-1]['id'],self.steps[self._currentstep-1]['string'])) + print((" prefix: %s" % self.prefix)) + print((" step #%i : %s" % (self.steps[self._currentstep-1]['id'],self.steps[self._currentstep-1]['string']))) #Ok, now if _currentstep is a member of steps, return true if self._currentstep in self.requestedsteps: - print("\n step #%i : %s\n" % (self.steps[self._currentstep-1]['id'],self.steps[self._currentstep-1]['string'])) + print(("\n step #%i : %s\n" % (self.steps[self._currentstep-1]['id'],self.steps[self._currentstep-1]['string']))) bool=True @@ -167,5 +167,5 @@ else: name=os.path.join(self.repository,name) - print("saving model as: '%s'" % name) + print(("saving model as: '%s'" % name)) #check that md is a model Index: /issm/trunk-jpl/src/py3/classes/pairoptions.py =================================================================== --- /issm/trunk-jpl/src/py3/classes/pairoptions.py (revision 23708) +++ /issm/trunk-jpl/src/py3/classes/pairoptions.py (revision 23709) @@ -30,5 +30,5 @@ if self.list: s+=" list: ({}x{}) \n\n".format(len(self.list),2) - for item in self.list.items(): + for item in list(self.list.items()): #if isinstance(item[1],str): s+=" field: {} value: '{}'\n".format((item[0],item[1])) @@ -55,5 +55,5 @@ else: #option is not a string, ignore it - print("WARNING: option number {} is not a string and will be ignored.".format(i+1)) + print(("WARNING: option number {} is not a string and will be ignored.".format(i+1))) # }}} def addfield(self,field,value): # {{{ @@ -61,5 +61,5 @@ if isinstance(field,str): if field in self.list: - print("WARNING: field '{}' with value={} exists and will be overwritten with value={}.".format(field,str(self.list[field]),str(value))) + print(("WARNING: field '{}' with value={} exists and will be overwritten with value={}.".format(field,str(self.list[field]),str(value)))) self.list[field] = value # }}} @@ -72,9 +72,9 @@ def AssignObjectFields(self,obj2): # {{{ """ASSIGNOBJECTFIELDS - assign object fields from options""" - for item in self.list.items(): + for item in list(self.list.items()): if item[0] in dir(obj2): setattr(obj2,item[0],item[1]) else: - print("WARNING: field '%s' is not a property of '%s'." % (item[0],type(obj2))) + print(("WARNING: field '%s' is not a property of '%s'." % (item[0],type(obj2)))) return obj2 # }}} @@ -150,5 +150,5 @@ #warn user if requested if warn: - print("removefield info: option '%s' has been removed from the list of options." % field) + print(("removefield info: option '%s' has been removed from the list of options." % field)) # }}} def marshall(self,md,fid,firstindex): # {{{ Index: /issm/trunk-jpl/src/py3/classes/plotoptions.py =================================================================== --- /issm/trunk-jpl/src/py3/classes/plotoptions.py (revision 23708) +++ /issm/trunk-jpl/src/py3/classes/plotoptions.py (revision 23709) @@ -23,5 +23,5 @@ if self.list: s+=" list: (%ix%i)\n" % (len(self.list),2) - for item in self.list.items(): + for item in list(self.list.items()): #s+=" options of plot number %i\n" % item if isinstance(item[1],str): @@ -50,5 +50,5 @@ else: #option is not a string, ignore it - print("WARNING: option number %d is not a string and will be ignored." % (i+1)) + print(("WARNING: option number %d is not a string and will be ignored." % (i+1))) #get figure number @@ -125,4 +125,4 @@ j=j+1 if j+1>numberofplots: - print("WARNING: too many instances of '%s' in options" % rawlist[i][0]) + print(("WARNING: too many instances of '%s' in options" % rawlist[i][0])) #}}} Index: /issm/trunk-jpl/src/py3/classes/qmu/@dakota_method/dakota_method.py =================================================================== --- /issm/trunk-jpl/src/py3/classes/qmu/@dakota_method/dakota_method.py (revision 23708) +++ /issm/trunk-jpl/src/py3/classes/qmu/@dakota_method/dakota_method.py (revision 23709) @@ -1,4 +1,6 @@ #move this later from helpers import * + +from MatlabFuncs import * import numpy as np @@ -44,12 +46,12 @@ def __init__(self,*args): - self.method = '' - self.type = '' - self.variables = [] - self.lcspec = [] - self.responses = [] - self.ghspec = [] + self.method ='' + self.type ='' + self.variables=[] + self.lcspec =[] + self.responses=[] + self.ghspec =[] #properites - self.params = struct() + self.params =struct() @staticmethod @@ -75,46 +77,45 @@ #given argument was a way of constructing a method else: - mlist=[ - 'dot_bfgs', - 'dot_frcg', - 'dot_mmfd', - 'dot_slp', - 'dot_sqp', - 'npsol_sqp', - 'conmin_frcg', - 'conmin_mfd', - 'optpp_cg', - 'optpp_q_newton', - 'optpp_fd_newton', - 'optpp_newton', - 'optpp_pds', - 'asynch_pattern_search', - 'coliny_cobyla', - 'coliny_direct', - 'coliny_ea', - 'coliny_pattern_search', - 'coliny_solis_wets', - 'ncsu_direct', - 'surrogate_based_local', - 'surrogate_based_global', - 'moga', - 'soga', - 'nl2sol', - 'nlssol_sqp', - 'optpp_g_newton', - 'nond_sampling', - 'nond_local_reliability', - 'nond_global_reliability', - 'nond_polynomial_chaos', - 'nond_stoch_collocation', - 'nond_evidence', - 'dace', - 'fsu_quasi_mc', - 'fsu_cvt', - 'vector_parameter_study', - 'list_parameter_study', - 'centered_parameter_study', - 'multidim_parameter_study', - 'bayes_calibration'] + mlist=['dot_bfgs', + 'dot_frcg', + 'dot_mmfd', + 'dot_slp', + 'dot_sqp', + 'npsol_sqp', + 'conmin_frcg', + 'conmin_mfd', + 'optpp_cg', + 'optpp_q_newton', + 'optpp_fd_newton', + 'optpp_newton', + 'optpp_pds', + 'asynch_pattern_search', + 'coliny_cobyla', + 'coliny_direct', + 'coliny_ea', + 'coliny_pattern_search', + 'coliny_solis_wets', + 'ncsu_direct', + 'surrogate_based_local', + 'surrogate_based_global', + 'moga', + 'soga', + 'nl2sol', + 'nlssol_sqp', + 'optpp_g_newton', + 'nond_sampling', + 'nond_local_reliability', + 'nond_global_reliability', + 'nond_polynomial_chaos', + 'nond_stoch_collocation', + 'nond_evidence', + 'dace', + 'fsu_quasi_mc', + 'fsu_cvt', + 'vector_parameter_study', + 'list_parameter_study', + 'centered_parameter_study', + 'multidim_parameter_study', + 'bayes_calibration'] mlist2=[] @@ -122,7 +123,5 @@ if strncmpi(method,mlist[i],len(method)): mlist2.append(mlist[i]) - - # check for a unique match in the list of methods - + # check for a unique match in the list of methods l = len(mlist2) if l == 0: @@ -134,7 +133,9 @@ # assign the default values for the method + # switch dm.method if dm.method in ['dot_bfgs','dot_frcg']: dm.type ='dot' - dm.variables=['continuous_design','continuous_state'] + dm.variables=['continuous_design', + 'continuous_state'] dm.lcspec =[] dm.responses=['objective_function'] @@ -151,5 +152,6 @@ elif dm.method in ['dot_mmfd','dot_slp','dot_sqp']: dm.type ='dot' - dm.variables=['continuous_design','continuous_state'] + dm.variables=['continuous_design', + 'continuous_state'] dm.lcspec =['linear_inequality_constraint', 'linear_equality_constraint'] @@ -236,5 +238,8 @@ dm.params.max_step=1000. dm.params.gradient_tolerance=1.0e-4 - elif dm.method in ['optpp_q_newton','optpp_fd_newton','optpp_newton']: + + elif dm.method in ['optpp_q_newton', + 'optpp_fd_newton', + 'optpp_newton']: dm.type ='optpp' dm.variables=['continuous_design', @@ -455,6 +460,7 @@ dm.params.vol_boxsize_limit=1.0e-8 - # elif dm.method in ['surrogate_based_local', - # 'surrogate_based_global']: + #if dm.method in ['surrogate_based_local', + #'surrogate_based_global']: + elif dm.method == 'moga': dm.type ='jega' @@ -501,4 +507,5 @@ dm.params.postprocessor_type=False dm.params.orthogonal_distance=[0.01] + elif dm.method == 'soga': dm.type ='jega' @@ -563,4 +570,5 @@ dm.params.covariance=0 dm.params.regression_stressbalances=False + elif dm.method == 'nlssol_sqp': dm.type ='lsq' @@ -619,5 +627,5 @@ dm.responses=['response_function'] dm.ghspec =[] - #not documented, but apparently works + # not documented, but apparently works dm.params.output=False dm.params.seed=False @@ -658,5 +666,5 @@ dm.responses=['response_function'] dm.ghspec =['grad'] - #not documented, but may work + # not documented, but may work dm.params.output=False dm.params.x_gaussian_process=False @@ -860,5 +868,5 @@ else: - raise RuntimeError('Unimplemented method: '+str(dm.method)+'.') + raise RuntimeError('Unimplemented method: {}.'.format(dm.method)) # if more than one argument, issue warning Index: /issm/trunk-jpl/src/py3/classes/qmu/@dakota_method/dmeth_params_set.py =================================================================== --- /issm/trunk-jpl/src/py3/classes/qmu/@dakota_method/dmeth_params_set.py (revision 23708) +++ /issm/trunk-jpl/src/py3/classes/qmu/@dakota_method/dmeth_params_set.py (revision 23709) @@ -22,6 +22,2 @@ return dm - - - - Index: /issm/trunk-jpl/src/py3/classes/qmu/@dakota_method/dmeth_params_write.py =================================================================== --- /issm/trunk-jpl/src/py3/classes/qmu/@dakota_method/dmeth_params_write.py (revision 23708) +++ /issm/trunk-jpl/src/py3/classes/qmu/@dakota_method/dmeth_params_write.py (revision 23709) @@ -1,3 +1,4 @@ from dakota_method import * +from MatlabFuncs import * from IssmConfig import * @@ -21,9 +22,13 @@ # unfortunately this prevents merely looping through the fields # of the parameters structure. + # write method-indepent controls + # param_write(fid,sbeg,'id_method',' = ','\n',dm.params) # param_write(fid,sbeg,'model_pointer',' = ','\n',dm.params) # write method-depent controls + + #switch dm.type if dm.type == 'dot': param_write(fid,sbeg,'max_iterations',' = ','\n',dm.params) @@ -95,17 +100,17 @@ raise RuntimeError('#s'' method must have only one algorithm.', dm.method) - param_write(fid,sbeg,'value_based_line_search','','\n',dm.params) - param_write(fid,sbeg,'gradient_based_line_search','','\n',dm.params) - param_write(fid,sbeg,'trust_region','','\n',dm.params) - param_write(fid,sbeg,'tr_pds','','\n',dm.params) - param_write(fid,sbeg,'max_step',' = ','\n',dm.params) - param_write(fid,sbeg,'gradient_tolerance',' = ','\n',dm.params) - param_write(fid,sbeg,'merit_function',' = ','\n',dm.params) - param_write(fid,sbeg,'central_path',' = ','\n',dm.params) - param_write(fid,sbeg,'steplength_to_boundary',' = ','\n',dm.params) - param_write(fid,sbeg,'centering_parameter',' = ','\n',dm.params) + param_write(fid,sbeg,'value_based_line_search','','\n',dm.params) + param_write(fid,sbeg,'gradient_based_line_search','','\n',dm.params) + param_write(fid,sbeg,'trust_region','','\n',dm.params) + param_write(fid,sbeg,'tr_pds','','\n',dm.params) + param_write(fid,sbeg,'max_step',' = ','\n',dm.params) + param_write(fid,sbeg,'gradient_tolerance',' = ','\n',dm.params) + param_write(fid,sbeg,'merit_function',' = ','\n',dm.params) + param_write(fid,sbeg,'central_path',' = ','\n',dm.params) + param_write(fid,sbeg,'steplength_to_boundary',' = ','\n',dm.params) + param_write(fid,sbeg,'centering_parameter',' = ','\n',dm.params) elif dm.method == 'optpp_pds': - param_write(fid,sbeg,'search_scheme_size',' = ','\n',dm.params) + param_write(fid,sbeg,'search_scheme_size',' = ','\n',dm.params) else: @@ -137,5 +142,4 @@ param_write(fid,sbeg,'output',' ','\n',dm.params) param_write(fid,sbeg,'scaling','','\n',dm.params) - param_write(fid,sbeg,'show_misc_options','','\n',dm.params) param_write(fid,sbeg,'misc_options',' = ','\n',dm.params) @@ -246,6 +250,6 @@ param_write(fid,sbeg,'niching_type',' = ','\n',dm.params) if not isempty(dm.params.radial) and not isempty(dm.params.distance): - raise RuntimeError('#s'' method must have only one niching distance.',dm.method) - + raise RuntimeError('#s'' method must have only one niching distance.', + dm.method) param_write(fid,sbeg,'radial',' = ','\n',dm.params) param_write(fid,sbeg,'distance',' = ','\n',dm.params) @@ -266,5 +270,4 @@ else: raise RuntimeError('Unrecognized '+dm.type+' method: '+dm.method+'.') - elif dm.type == 'lsq': @@ -310,5 +313,6 @@ dm.params.trust_region + dm.params.tr_pds > 1): - raise RuntimeError('#s'' method must have only one algorithm.',dm.method) + raise RuntimeError('#s'' method must have only one algorithm.', + dm.method) param_write(fid,sbeg,'value_based_line_search','','\n',dm.params) @@ -348,5 +352,7 @@ if type(dm.params.mpp_search) == str: if (dm.params.sqp +dm.params.nip > 1): - raise RuntimeError('#s'' method must have only one algorithm.',dm.method) + raise RuntimeError('#s'' method must have only one algorithm.', + dm.method) + param_write(fid,sbeg,'sqp','','\n',dm.params) param_write(fid,sbeg,'nip','','\n',dm.params) @@ -360,5 +366,7 @@ elif dm.method == 'nond_global_reliability': if (dm.params.x_gaussian_process + dm.params.u_gaussian_process != 1): - raise RuntimeError('#s'' method must have one and only one algorithm.',dm.method) + raise RuntimeError('#s'' method must have one and only one algorithm.', + dm.method) + param_write(fid,sbeg,'x_gaussian_process','','\n',dm.params) param_write(fid,sbeg,'u_gaussian_process','','\n',dm.params) @@ -399,9 +407,12 @@ raise RuntimeError('Unrecognized '+dm.type+' method: '+dm.method+'.') + elif dm.type == 'dace': #switch dm.method if dm.method == 'dace': if (dm.params.grid + dm.params.random + dm.params.oas + dm.params.lhs + dm.params.oa_lhs + dm.params.box_behnken + dm.params.central_composite != 1): - raise RuntimeError('#s'' method must have one and only one algorithm.',dm.method) + raise RuntimeError('#s'' method must have one and only one algorithm.', + dm.method) + param_write(fid,sbeg,'grid','','\n',dm.params) param_write(fid,sbeg,'random','','\n',dm.params) @@ -421,4 +432,5 @@ if (dm.params.halton + dm.params.hammersley != 1): raise RuntimeError('#s'' method must have one and only one sequence type.',dm.method) + param_write(fid,sbeg,'halton','','\n',dm.params) param_write(fid,sbeg,'hammersley','','\n',dm.params) @@ -445,5 +457,4 @@ raise RuntimeError('Unrecognized '+dm.type+' method: '+dm.method+'.') - elif dm.type == 'param': param_write(fid,sbeg,'output',' ','\n',dm.params) @@ -457,4 +468,5 @@ param_write(fid,sbeg,'step_length',' = ','\n',dm.params) param_write(fid,sbeg,'num_steps',' = ','\n',dm.params) + elif not isempty(dm.params.step_vector): param_write(fid,sbeg,'step_vector',' = ','\n',dm.params) @@ -474,13 +486,14 @@ raise RuntimeError('Unrecognized '+dm.type+' method: '+dm.method+'.') + elif dm.type == 'bayes': #switch dm.method if dm.method == 'bayes_calibration': - # if (dm.params.queso + - # dm.params.dream + - # dm.params.gpmsa ~= 1) - # raise RuntimeError('''#s'' method must have one and only one bayes type. YOU SUCK', - # dm.method) - # + # if (dm.params.queso + + # dm.params.dream + + # dm.params.gpmsa ~= 1) + # raise RuntimeError('''#s'' method must have one and only one bayes type. YOU SUCK', + # dm.method) + # param_write(fid,sbeg,'queso','','\n',dm.params) param_write(fid,sbeg,'dream','','\n',dm.params) @@ -501,5 +514,4 @@ # loop through each parameter field in the structure fnames=fieldnames(params) - for i in range(np.size(fnames)): param_write(fidi,sbeg,fnames[i],smid,s,params) @@ -517,5 +529,5 @@ return elif isempty(vars(params)[pname]): - print('Warning: param_write:param_empty: Parameter '+pname+' requires input of type '+type(vars(params)[pname])+'.') + print('Warning: param_write:param_empty: Parameter {} requires input of type {}.'.format(pname,type(vars(params)[pname]))) return @@ -523,6 +535,8 @@ if type(vars(params)[pname]) == bool: fidi.write(sbeg+str(pname)+s) + elif type(vars(params)[pname]) in [int,float]: fidi.write(sbeg+str(pname)+smid+str(vars(params)[pname])+s) + elif type(vars(params)[pname]) == list: fidi.write(sbeg+str(pname)+smid+str(vars(params)[pname][0])) @@ -531,7 +545,9 @@ fidi.write(s) + elif type(vars(params)[pname]) == str: fidi.write(sbeg+str(pname)+smid+str(vars(params)[pname])+s) + else: - print('Warning: param_write:param_unrecog: Parameter '+pname+' is of unrecognized type '+type(vars(params)[pname])+'.') + print('Warning: param_write:param_unrecog: Parameter {} is of unrecognized type {}.'.format(pname,type(vars(params)[pname]))) return Index: /issm/trunk-jpl/src/py3/classes/qmu/calibration_function.py =================================================================== --- /issm/trunk-jpl/src/py3/classes/qmu/calibration_function.py (revision 23708) +++ /issm/trunk-jpl/src/py3/classes/qmu/calibration_function.py (revision 23709) @@ -1,4 +1,5 @@ import numpy as np -from rlist_write import rlist_write +from rlist_write import * +from MatlabArray import * class calibration_function(object): @@ -70,9 +71,7 @@ string += ' scale: '+str(self.scale) + '\n' string += ' weight: '+str(self.weight) + '\n' - return string # from here on, cf is either a single, or a 1d vector of, calibration_function - @staticmethod def prop_desc(cf,dstr): @@ -96,5 +95,4 @@ desc = allempty(desc) - return desc Index: /issm/trunk-jpl/src/py3/classes/qmu/least_squares_term.py =================================================================== --- /issm/trunk-jpl/src/py3/classes/qmu/least_squares_term.py (revision 23708) +++ /issm/trunk-jpl/src/py3/classes/qmu/least_squares_term.py (revision 23709) @@ -1,4 +1,5 @@ import numpy as np -from rlist_write import rlist_write +from rlist_write import * +from MatlabArray import * class least_squares_term(object): @@ -58,14 +59,14 @@ for i in range(np.size(lst)): if (np.size(args[2]) > 1): - lst[i].scale = args[2][i] + lst[i].scale = args[2][i] else: - lst[i].scale = args[2] + lst[i].scale = args[2] if (nargin >= 4): for i in range(np.size(lst)): if (np.size(args[3]) > 1): - lst[i].weight = args[3][i] + lst[i].weight = args[3][i] else: - lst[i].weight = args[3] + lst[i].weight = args[3] if (nargin > 4): @@ -82,5 +83,4 @@ string += ' scale: '+str(self.scale) + '\n' string += ' weight: '+str(self.weight) + '\n' - return string @@ -100,5 +100,4 @@ desc = allempty(desc) - return desc @@ -113,5 +112,4 @@ stype = allequal(stype,'none') - return stype @@ -126,5 +124,4 @@ scale = allequal(scale,1.) - return scale @@ -139,5 +136,4 @@ weight = allequal(weight,1.) - return weight Index: /issm/trunk-jpl/src/py3/classes/qmu/normal_uncertain.py =================================================================== --- /issm/trunk-jpl/src/py3/classes/qmu/normal_uncertain.py (revision 23708) +++ /issm/trunk-jpl/src/py3/classes/qmu/normal_uncertain.py (revision 23709) @@ -1,3 +1,5 @@ import numpy as np +#from vlist_write import * +from MatlabArray import * class normal_uncertain(object): @@ -128,5 +130,4 @@ upper = allequal(upper,-np.inf) - return upper @@ -173,4 +174,6 @@ nuv = [struc_class(i,'normal_uncertain','nuv') for i in dvar] + # possible namespace pollution, the above import seems not to work + from vlist_write import vlist_write # write variables vlist_write(fidi,'normal_uncertain','nuv',nuv) Index: /issm/trunk-jpl/src/py3/classes/qmu/objective_function.py =================================================================== --- /issm/trunk-jpl/src/py3/classes/qmu/objective_function.py (revision 23708) +++ /issm/trunk-jpl/src/py3/classes/qmu/objective_function.py (revision 23709) @@ -1,4 +1,5 @@ import numpy as np -from rlist_write import rlist_write +from rlist_write import * +from MatlabArray import * class objective_function(object): @@ -29,4 +30,5 @@ def objective_function(*args): nargin = len(args) + # create a default object if nargin == 0: @@ -40,4 +42,5 @@ shapec = array_size(*args[0:min(nargin,4)]) of = [objective_function() for i in range(shapec[0]) for j in range(shapec[1])] + for i in range(np.size(of)): if (np.size(args[0]) > 1): @@ -72,4 +75,5 @@ return of + def __repr__(self): # display the object @@ -80,5 +84,4 @@ string += ' scale: ' +str(self.scale) + '\n' string += ' weight: ' +str(self.weight) + '\n' - return string @@ -104,5 +107,4 @@ desc = allempty(desc) - return desc @@ -132,5 +134,4 @@ weight = allequal(weight,1.) - return weight @@ -145,5 +146,4 @@ stype = allequal(stype,'none') - return stype @@ -158,5 +158,4 @@ scale = allequal(scale,1.) - return scale Index: /issm/trunk-jpl/src/py3/classes/qmu/response_function.py =================================================================== --- /issm/trunk-jpl/src/py3/classes/qmu/response_function.py (revision 23708) +++ /issm/trunk-jpl/src/py3/classes/qmu/response_function.py (revision 23709) @@ -1,5 +1,6 @@ import numpy as np -from rlist_write import * -from rlev_write import rlev_write +#from rlist_write import * +from rlev_write import * +from MatlabArray import * #move this later @@ -35,4 +36,5 @@ @staticmethod def response_function(*args): + nargin = len(args) # create a default object @@ -75,5 +77,4 @@ return rf - def __repr__(self): # display the object @@ -109,5 +110,4 @@ desc = allempty(desc) - return desc @@ -149,9 +149,6 @@ respl = empty_nd_list(np.size(rf)) - probl = empty_nd_list(np.size(rf)) - rell = empty_nd_list(np.size(rf)) - grell = empty_nd_list(np.size(rf)) @@ -166,5 +163,4 @@ rell = allempty(rell) grell = allempty(grell) - return [respl,probl,rell,grell] Index: /issm/trunk-jpl/src/py3/classes/qmu/uniform_uncertain.py =================================================================== --- /issm/trunk-jpl/src/py3/classes/qmu/uniform_uncertain.py (revision 23708) +++ /issm/trunk-jpl/src/py3/classes/qmu/uniform_uncertain.py (revision 23709) @@ -1,4 +1,4 @@ import numpy as np -from vlist_write import vlist_write +#from vlist_write import * from MatlabArray import * @@ -28,4 +28,5 @@ def uniform_uncertain(*args): nargin = len(args) + # create a default object if nargin == 0: @@ -153,5 +154,6 @@ # collect only the variables of the appropriate class uuv = [struc_class(i,'uniform_uncertain','uuv') for i in dvar] - + # possible namespace pollution, the above import seems not to work + from vlist_write import vlist_write # write variables vlist_write(fidi,'uniform_uncertain','uuv',uuv) Index: /issm/trunk-jpl/src/py3/classes/results.py =================================================================== --- /issm/trunk-jpl/src/py3/classes/results.py (revision 23708) +++ /issm/trunk-jpl/src/py3/classes/results.py (revision 23709) @@ -25,5 +25,5 @@ s+="%s\n" % fielddisplay(self,'SolutionType',"solution type") - for name in self.__dict__.keys(): + for name in list(self.__dict__.keys()): if name not in ['step','time','SolutionType','errlog','outlog']: if isinstance(getattr(self,name),list): Index: /issm/trunk-jpl/src/py3/classes/rifts.py =================================================================== --- /issm/trunk-jpl/src/py3/classes/rifts.py (revision 23708) +++ /issm/trunk-jpl/src/py3/classes/rifts.py (revision 23709) @@ -47,5 +47,5 @@ md.checkmessage("model should be processed for rifts (run meshprocessrifts)!") for i,rift in enumerate(self.riftstruct): - md = checkfield(md,'fieldname',"rifts.riftstruct[%d]['fill']" % i,'values',['Water','Air','Ice','Melange',0,1,2,3]) + md = checkfield(md,'fieldname',"rifts.riftstruct[{}]['fill']".format(i),'values',['Water','Air','Ice','Melange',0,1,2,3]) else: if self.riftstruct and np.any(np.logical_not(isnans(self.riftstruct))): Index: /issm/trunk-jpl/src/py3/classes/toolkits.py =================================================================== --- /issm/trunk-jpl/src/py3/classes/toolkits.py (revision 23708) +++ /issm/trunk-jpl/src/py3/classes/toolkits.py (revision 23709) @@ -35,5 +35,5 @@ def __repr__(self): # {{{ s ="List of toolkits options per analysis:\n\n" - for analysis in vars(self).keys(): + for analysis in list(vars(self).keys()): s+="%s\n" % fielddisplay(self,analysis,'') @@ -56,5 +56,5 @@ # }}} def checkconsistency(self,md,solution,analyses): # {{{ - for analysis in vars(self).keys(): + for analysis in list(vars(self).keys()): if not getattr(self,analysis): md.checkmessage("md.toolkits.%s is empty" % analysis) @@ -83,5 +83,5 @@ #start writing options - for analysis in vars(self).keys(): + for analysis in list(vars(self).keys()): options=getattr(self,analysis) @@ -89,5 +89,5 @@ fid.write("\n+%s\n" % analysis) #append a + to recognize it's an analysis enum #now, write options - for optionname,optionvalue in options.items(): + for optionname,optionvalue in list(options.items()): if not optionvalue: Index: /issm/trunk-jpl/src/py3/classes/verbose.py =================================================================== --- /issm/trunk-jpl/src/py3/classes/verbose.py (revision 23708) +++ /issm/trunk-jpl/src/py3/classes/verbose.py (revision 23709) @@ -67,5 +67,5 @@ #Cast to logicals listproperties=vars(self) - for fieldname,fieldvalue in listproperties.items(): + for fieldname,fieldvalue in list(listproperties.items()): if isinstance(fieldvalue,bool) or isinstance(fieldvalue,(int,float)): setattr(self,fieldname,bool(fieldvalue)) Index: /issm/trunk-jpl/src/py3/consistency/checkfield.py =================================================================== --- /issm/trunk-jpl/src/py3/consistency/checkfield.py (revision 23708) +++ /issm/trunk-jpl/src/py3/consistency/checkfield.py (revision 23709) @@ -1,4 +1,5 @@ import numpy as np import os +from re import findall,split from pairoptions import pairoptions from operator import attrgetter @@ -42,7 +43,17 @@ else: fieldname=options.getfieldvalue('fieldname') - - field=attrgetter(fieldname)(md) + fieldprefix=split(r'\[(.*?)\]',fieldname)[0] + fieldindexes=findall(r'\[(.*?)\]',fieldname) + field=attrgetter(fieldprefix)(md) + for index in fieldindexes: + try: + field=field[index.strip("\'")] + except TypeError: + field=field[int(index)] #looking for an index and not a key + +# that works for py2 +# exec("field=md.{}".format(fieldname)) # exec("field=md.{}".format(fieldname),namespace) + if isinstance(field,(bool,int,float)): @@ -135,4 +146,6 @@ if options.exist('>='): lowerbound = options.getfieldvalue('>=') + if type(lowerbound) is str: + lowerbound=attrgetter(lowerbound)(md) if np.size(lowerbound)>1: #checking elementwise if any(field'): lowerbound=options.getfieldvalue('>') + if type(lowerbound) is str: + lowerbound=attrgetter(lowerbound)(md) if np.size(lowerbound)>1: #checking elementwise if any(field<=upperbound): @@ -172,4 +187,6 @@ if options.exist('<='): upperbound=options.getfieldvalue('<=') + if type(upperbound) is str: + upperbound=attrgetter(upperbound)(md) if np.size(upperbound)>1: #checking elementwise if any(field>upperbound): @@ -189,4 +206,6 @@ if options.exist('<'): upperbound=options.getfieldvalue('<') + if type(upperbound) is str: + upperbound=attrgetter(upperbound)(md) if np.size(upperbound)>1: #checking elementwise if any(field>=upperbound): Index: /issm/trunk-jpl/src/py3/contrib/defleurian/netCDF/export_netCDF.py =================================================================== --- /issm/trunk-jpl/src/py3/contrib/defleurian/netCDF/export_netCDF.py (revision 23708) +++ /issm/trunk-jpl/src/py3/contrib/defleurian/netCDF/export_netCDF.py (revision 23709) @@ -12,5 +12,5 @@ if path.exists(filename): print(('File {} allready exist'.format(filename))) - newname=input('Give a new name or "delete" to replace: ') + newname=eval(input('Give a new name or "delete" to replace: ')) if newname=='delete': remove(filename) Index: /issm/trunk-jpl/src/py3/contrib/defleurian/paraview/enveloppeVTK.py =================================================================== --- /issm/trunk-jpl/src/py3/contrib/defleurian/paraview/enveloppeVTK.py (revision 23708) +++ /issm/trunk-jpl/src/py3/contrib/defleurian/paraview/enveloppeVTK.py (revision 23709) @@ -28,5 +28,5 @@ if os.path.exists(filename): print(('File {} allready exist'.format(filename))) - newname=input('Give a new name or "delete" to replace: ') + newname=eval(input('Give a new name or "delete" to replace: ')) if newname=='delete': filelist = glob.glob(filename+'/*') Index: /issm/trunk-jpl/src/py3/contrib/defleurian/paraview/enveloppeVTK_bleeding.py =================================================================== --- /issm/trunk-jpl/src/py3/contrib/defleurian/paraview/enveloppeVTK_bleeding.py (revision 23708) +++ /issm/trunk-jpl/src/py3/contrib/defleurian/paraview/enveloppeVTK_bleeding.py (revision 23709) @@ -32,5 +32,5 @@ if os.path.exists(filename): print(('File {} allready exist'.format(filename))) - newname=input('Give a new name or "delete" to replace: ') + newname=eval(input('Give a new name or "delete" to replace: ')) if newname=='delete': filelist = glob.glob(filename+'/*') Index: /issm/trunk-jpl/src/py3/contrib/defleurian/paraview/exportVTK.py =================================================================== --- /issm/trunk-jpl/src/py3/contrib/defleurian/paraview/exportVTK.py (revision 23708) +++ /issm/trunk-jpl/src/py3/contrib/defleurian/paraview/exportVTK.py (revision 23709) @@ -28,5 +28,5 @@ if os.path.exists(filename): print(('File {} allready exist'.format(filename))) - newname=input('Give a new name or "delete" to replace: ') + newname=eval(input('Give a new name or "delete" to replace: ')) if newname=='delete': filelist = glob.glob(filename+'/*') @@ -193,5 +193,5 @@ for node in range(0,num_of_points): #paraview does not like NaN, replacing - print(other_struct.__dict__[field][node]) + print((other_struct.__dict__[field][node])) if np.isnan(other_struct.__dict__[field][node]): fid.write('%e\n' % -9999.9999) Index: /issm/trunk-jpl/src/py3/contrib/morlighem/bamg/YamsCall.py =================================================================== --- /issm/trunk-jpl/src/py3/contrib/morlighem/bamg/YamsCall.py (revision 23708) +++ /issm/trunk-jpl/src/py3/contrib/morlighem/bamg/YamsCall.py (revision 23709) @@ -29,19 +29,19 @@ #Compute Hessian t1=time.time() - print("%s" % ' computing Hessian...') + print(("%s" % ' computing Hessian...')) hessian=ComputeHessian(md.mesh.elements,md.mesh.x,md.mesh.y,field,'node') t2=time.time() - print("%s%d%s\n" % (' done (',t2-t1,' seconds)')) + print(("%s%d%s\n" % (' done (',t2-t1,' seconds)'))) #Compute metric t1=time.time() - print("%s" % ' computing metric...') + print(("%s" % ' computing metric...')) metric=ComputeMetric(hessian,scale,epsilon,hmin,hmax,np.empty(0,int)) t2=time.time() - print("%s%d%s\n" % (' done (',t2-t1,' seconds)')) + print(("%s%d%s\n" % (' done (',t2-t1,' seconds)'))) #write files t1=time.time() - print("%s" % ' writing initial mesh files...') + print(("%s" % ' writing initial mesh files...')) np.savetxt('carre0.met',metric) @@ -80,8 +80,8 @@ f.close() t2=time.time() - print("%s%d%s\n" % (' done (',t2-t1,' seconds)')) + print(("%s%d%s\n" % (' done (',t2-t1,' seconds)'))) #call yams - print("%s\n" % ' call Yams...') + print(("%s\n" % ' call Yams...')) if m.ispc(): #windows @@ -96,5 +96,5 @@ #plug new mesh t1=time.time() - print("\n%s" % ' reading final mesh files...') + print(("\n%s" % ' reading final mesh files...')) Tria=np.loadtxt('carre1.tria',int) Coor=np.loadtxt('carre1.coor',float) @@ -107,9 +107,9 @@ numberofelements2=md.mesh.numberofelements t2=time.time() - print("%s%d%s\n\n" % (' done (',t2-t1,' seconds)')) + print(("%s%d%s\n\n" % (' done (',t2-t1,' seconds)'))) #display number of elements - print("\n%s %i" % (' inital number of elements:',numberofelements1)) - print("\n%s %i\n\n" % (' new number of elements:',numberofelements2)) + print(("\n%s %i" % (' inital number of elements:',numberofelements1))) + print(("\n%s %i\n\n" % (' new number of elements:',numberofelements2))) #clean up: Index: /issm/trunk-jpl/src/py3/coordsystems/gmtmask.py =================================================================== --- /issm/trunk-jpl/src/py3/coordsystems/gmtmask.py (revision 23708) +++ /issm/trunk-jpl/src/py3/coordsystems/gmtmask.py (revision 23709) @@ -21,7 +21,7 @@ if recursive: - print(' recursing: num vertices #'+str(lenlat)) + print((' recursing: num vertices #'+str(lenlat))) else: - print('gmtmask: num vertices #'+str(lenlat)) + print(('gmtmask: num vertices #'+str(lenlat))) #Check lat and long size is not more than 50,000 If so, recursively call gmtmask: Index: /issm/trunk-jpl/src/py3/dev/issmversion.py =================================================================== --- /issm/trunk-jpl/src/py3/dev/issmversion.py (revision 23708) +++ /issm/trunk-jpl/src/py3/dev/issmversion.py (revision 23709) @@ -11,8 +11,8 @@ print(' ') -print(IssmConfig('PACKAGE_NAME')[0]+' Version '+IssmConfig('PACKAGE_VERSION')[0]) -print('(website: '+IssmConfig('PACKAGE_URL')[0]+' contact: '+IssmConfig('PACKAGE_BUGREPORT')[0]+')') +print((IssmConfig('PACKAGE_NAME')[0]+' Version '+IssmConfig('PACKAGE_VERSION')[0])) +print(('(website: '+IssmConfig('PACKAGE_URL')[0]+' contact: '+IssmConfig('PACKAGE_BUGREPORT')[0]+')')) print(' ') -print('Build date: '+IssmConfig('PACKAGE_BUILD_DATE')[0]) +print(('Build date: '+IssmConfig('PACKAGE_BUILD_DATE')[0])) print('Copyright (c) 2009-2018 California Institute of Technology') print(' ') Index: /issm/trunk-jpl/src/py3/exp/expcoarsen.py =================================================================== --- /issm/trunk-jpl/src/py3/exp/expcoarsen.py (revision 23708) +++ /issm/trunk-jpl/src/py3/exp/expcoarsen.py (revision 23709) @@ -23,5 +23,5 @@ raise OSError("expcoarsen error message: file '%s' not found!" % oldfile) if os.path.exists(newfile): - choice=input('A file ' + newfile + ' already exists, do you want to modify it? (y/n)') + choice=eval(input('A file ' + newfile + ' already exists, do you want to modify it? (y/n)')) if choice not in 'y': print('no modification done ... exiting') Index: /issm/trunk-jpl/src/py3/geometry/NowickiProfile.py =================================================================== --- /issm/trunk-jpl/src/py3/geometry/NowickiProfile.py (revision 23708) +++ /issm/trunk-jpl/src/py3/geometry/NowickiProfile.py (revision 23709) @@ -28,5 +28,5 @@ #upstream of the GL - for i in range(np.size(x) / 2): + for i in range(int(np.size(x)/2)): ss = np.roots([1, 4 * lamda * beta, 0, 0, 6 * lamda * ms * x[i]**2 + 12 * lamda * q * x[i] - hg**4 - 4 * lamda * beta * hg**3]) @@ -38,5 +38,5 @@ #downstream of the GL - for i in range(np.size(x) / 2, np.size(x)): + for i in range(int(np.size(x)/2), int(np.size(x))): h[i] = (x[i] / (4. * (delta+1) * q) + hg**(-2))**(-0.5) # ice thickness for ice shelf from (3.1) b[i] = sea - h[i] * (1. / (1+delta)) Index: /issm/trunk-jpl/src/py3/io/loadmodel.py =================================================================== --- /issm/trunk-jpl/src/py3/io/loadmodel.py (revision 23708) +++ /issm/trunk-jpl/src/py3/io/loadmodel.py (revision 23709) @@ -1,4 +1,4 @@ from loadvars import loadvars -from dbm import whichdb +from dbm.ndbm import whichdb from netCDF4 import Dataset @@ -27,5 +27,5 @@ #recover model on file and name it md struc=loadvars(path) - name=[key for key in struc.keys()] + name=[key for key in list(struc.keys())] if len(name)>1: raise IOError("loadmodel error message: file '%s' contains several variables. Only one model should be present." % path) Index: /issm/trunk-jpl/src/py3/io/loadvars.py =================================================================== --- /issm/trunk-jpl/src/py3/io/loadvars.py (revision 23708) +++ /issm/trunk-jpl/src/py3/io/loadvars.py (revision 23709) @@ -7,5 +7,5 @@ from os import path from collections import OrderedDict -from dbm import whichdb +from dbm.ndbm import whichdb from model import * @@ -57,20 +57,20 @@ if whichdb(filename): - print("Loading variables from file '%s'." % filename) + print(("Loading variables from file '%s'." % filename)) my_shelf = shelve.open(filename,'r') # 'r' for read-only if nvdict: - for name in nvdict.keys(): + for name in list(nvdict.keys()): try: nvdict[name] = my_shelf[name] - print("Variable '%s' loaded." % name) + print(("Variable '%s' loaded." % name)) except KeyError: value = None - print("Variable '%s' not found." % name) + print(("Variable '%s' not found." % name)) else: - for name in my_shelf.keys(): + for name in list(my_shelf.keys()): nvdict[name] = my_shelf[name] - print("Variable '%s' loaded." % name) + print(("Variable '%s' loaded." % name)) my_shelf.close() Index: /issm/trunk-jpl/src/py3/io/savevars.py =================================================================== --- /issm/trunk-jpl/src/py3/io/savevars.py (revision 23708) +++ /issm/trunk-jpl/src/py3/io/savevars.py (revision 23709) @@ -46,16 +46,16 @@ if os.path.exists(filename): - print("Shelving variables to existing file '%s'." % filename) + print(("Shelving variables to existing file '%s'." % filename)) else: - print("Shelving variables to new file '%s'." % filename) + print(("Shelving variables to new file '%s'." % filename)) my_shelf = shelve.open(filename,'c') # 'c' for create if not exist, else 'n' for new - for name,value in nvdict.items(): + for name,value in list(nvdict.items()): try: my_shelf[name] = value - print("Variable '%s' shelved." % name) + print(("Variable '%s' shelved." % name)) except TypeError: - print("Variable '%s' not shelved." % name) + print(("Variable '%s' not shelved." % name)) my_shelf.close() Index: /issm/trunk-jpl/src/py3/mech/thomasparams.py =================================================================== --- /issm/trunk-jpl/src/py3/mech/thomasparams.py (revision 23708) +++ /issm/trunk-jpl/src/py3/mech/thomasparams.py (revision 23709) @@ -126,5 +126,5 @@ pos=np.nonzero(np.abs((np.abs(a)-2.))<1.e-3) if len(pos)>0: - print('Warning: ', len(pos), ' vertices have alpha within 1e-3 of -2') + print(('Warning: ', len(pos), ' vertices have alpha within 1e-3 of -2')) a[pos]=-2+1e-3 Index: /issm/trunk-jpl/src/py3/mesh/bamg.py =================================================================== --- /issm/trunk-jpl/src/py3/mesh/bamg.py (revision 23708) +++ /issm/trunk-jpl/src/py3/mesh/bamg.py (revision 23709) @@ -572,5 +572,5 @@ if num: - print("WARNING: %d points outside the domain outline have been removed" % num) + print(("WARNING: %d points outside the domain outline have been removed" % num)) """ Index: /issm/trunk-jpl/src/py3/mesh/triangle.py =================================================================== --- /issm/trunk-jpl/src/py3/mesh/triangle.py (revision 23708) +++ /issm/trunk-jpl/src/py3/mesh/triangle.py (revision 23709) @@ -37,5 +37,5 @@ #Check that mesh was not already run, and warn user: if md.mesh.numberofelements: - choice = input('This model already has a mesh. Are you sure you want to go ahead? (y/n)') + choice = eval(input('This model already has a mesh. Are you sure you want to go ahead? (y/n)')) if not m.strcmp(choice,'y'): print('no meshing done ... exiting') Index: /issm/trunk-jpl/src/py3/miscellaneous/fielddisplay.py =================================================================== --- /issm/trunk-jpl/src/py3/miscellaneous/fielddisplay.py (revision 23708) +++ /issm/trunk-jpl/src/py3/miscellaneous/fielddisplay.py (revision 23709) @@ -67,5 +67,5 @@ offset+=' ' - for structure_field,sfield in field.items(): + for structure_field,sfield in list(field.items()): string+=parsedisplay(offset,str(structure_field),sfield,'')+'\n' Index: /issm/trunk-jpl/src/py3/modules/MeshPartition.py =================================================================== --- /issm/trunk-jpl/src/py3/modules/MeshPartition.py (revision 23708) +++ /issm/trunk-jpl/src/py3/modules/MeshPartition.py (revision 23709) @@ -15,5 +15,5 @@ ''' if md == None or numpartitions == None: - print(MeshPartition.__doc__) + print((MeshPartition.__doc__)) raise RuntimeError('Wrong usage (see above)') Index: /issm/trunk-jpl/src/py3/os/issmscpin.py =================================================================== --- /issm/trunk-jpl/src/py3/os/issmscpin.py (revision 23708) +++ /issm/trunk-jpl/src/py3/os/issmscpin.py (revision 23709) @@ -43,6 +43,6 @@ raise OSError("issmscpin error message: could not find ISSM_DIR_WIN environment variable.") - username=input('Username: (quoted string) ') - key=input('Key: (quoted string) ') + username=eval(input('Username: (quoted string) ')) + key=eval(input('Key: (quoted string) ')) for package in packages: Index: /issm/trunk-jpl/src/py3/os/issmscpout.py =================================================================== --- /issm/trunk-jpl/src/py3/os/issmscpout.py (revision 23708) +++ /issm/trunk-jpl/src/py3/os/issmscpout.py (revision 23709) @@ -36,6 +36,6 @@ raise OSError("issmscpout error message: could not find ISSM_DIR_WIN environment variable.") - username=input('Username: (quoted string) ') - key=input('Key: (quoted string) ') + username=eval(input('Username: (quoted string) ')) + key=eval(input('Key: (quoted string) ')) for package in packages: Index: /issm/trunk-jpl/src/py3/os/issmssh.py =================================================================== --- /issm/trunk-jpl/src/py3/os/issmssh.py (revision 23708) +++ /issm/trunk-jpl/src/py3/os/issmssh.py (revision 23709) @@ -29,6 +29,6 @@ raise OSError("issmssh error message: could not find ISSM_DIR_WIN environment variable.") - username=input('Username: (quoted string) ') - key=input('Key: (quoted string) ') + username=eval(input('Username: (quoted string) ')) + key=eval(input('Key: (quoted string) ')) subprocess.call('%s/externalpackages/ssh/plink.exe -ssh -l "%s" -pw "%s" %s "%s"' % (ISSM_DIR,username,key,host,command),shell=True); Index: /issm/trunk-jpl/src/py3/parameterization/setflowequation.py =================================================================== --- /issm/trunk-jpl/src/py3/parameterization/setflowequation.py (revision 23708) +++ /issm/trunk-jpl/src/py3/parameterization/setflowequation.py (revision 23709) @@ -11,6 +11,6 @@ 'SIA','SSA','HO','L1L2','FS' and 'fill' are the possible options that must be followed by the corresponding exp file or flags list - It can either be a domain file (argus type, .exp extension), or an array of element flags. - If user wants every element outside the domain to be + It can either be a domain file (argus type, .exp extension), or an array of element flags. + If user wants every element outside the domain to be setflowequationd, add '~' to the name of the domain file (ex: '~HO.exp'); an empty string '' will be considered as an empty domain @@ -96,5 +96,5 @@ #Then complete with NoneApproximation or the other model used if there is no FS - if any(FSflag): + if any(FSflag): if any(HOflag): #fill with HO HOflag[~FSflag]=True @@ -103,5 +103,5 @@ SSAflag[~FSflag]=True nodeonSSA[md.mesh.elements[np.where(SSAflag),:]-1]=True - else: #fill with none + else: #fill with none noneflag[np.where(~FSflag)]=True @@ -285,3 +285,2 @@ return md - Index: /issm/trunk-jpl/src/py3/parameterization/setmask.py =================================================================== --- /issm/trunk-jpl/src/py3/parameterization/setmask.py (revision 23708) +++ /issm/trunk-jpl/src/py3/parameterization/setmask.py (revision 23709) @@ -10,6 +10,6 @@ SETMASK - establish boundaries between grounded and floating ice. - By default, ice is considered grounded. The contour floatingicename defines nodes - for which ice is floating. The contour groundedicename defines nodes inside an floatingice, + By default, ice is considered grounded. The contour floatingicename defines nodes + for which ice is floating. The contour groundedicename defines nodes inside an floatingice, that are grounded (ie: ice rises, islands, etc ...) All input files are in the Argus format (extension .exp). @@ -39,8 +39,8 @@ #Assign elementonfloatingice, elementongroundedice, vertexongroundedice and vertexonfloatingice. Only change at your own peril! This is synchronized heavily with the GroundingLineMigration module. {{{ elementonfloatingice = FlagElements(md, floatingicename) - elementongroundedice = FlagElements(md, groundedicename) + elementongroundedice = FlagElements(md, groundedicename) - #Because groundedice nodes and elements can be included into an floatingice, we need to update. Remember, all the previous - #arrays come from domain outlines that can intersect one another: + #Because groundedice nodes and elements can be included into an floatingice, we need to update. Remember, all the previous + #arrays come from domain outlines that can intersect one another: elementonfloatingice = np.logical_and(elementonfloatingice,np.logical_not(elementongroundedice)) Index: /issm/trunk-jpl/src/py3/plot/plot_manager.py =================================================================== --- /issm/trunk-jpl/src/py3/plot/plot_manager.py (revision 23708) +++ /issm/trunk-jpl/src/py3/plot/plot_manager.py (revision 23709) @@ -80,5 +80,5 @@ return else: - print("WARNING: '%s' is not implemented or is not a valid string for option 'data'" % data) + print(("WARNING: '%s' is not implemented or is not a valid string for option 'data'" % data)) # }}} # {{{ Gridded plot TODO Index: /issm/trunk-jpl/src/py3/plot/plot_overlay.py =================================================================== --- /issm/trunk-jpl/src/py3/plot/plot_overlay.py (revision 23708) +++ /issm/trunk-jpl/src/py3/plot/plot_overlay.py (revision 23709) @@ -111,5 +111,5 @@ lon_0=0 else: - hemisphere=input('epsg code {} is not supported chose your hemisphere (1 for North, -1 for south)'.format(mesh.epsg)) + hemisphere=eval(input('epsg code {} is not supported chose your hemisphere (1 for North, -1 for south)'.format(mesh.epsg))) lat,lon=xy2ll(xlim,ylim,hemisphere,lon_0,st_lat) Index: /issm/trunk-jpl/src/py3/plot/processdata.py =================================================================== --- /issm/trunk-jpl/src/py3/plot/processdata.py (revision 23708) +++ /issm/trunk-jpl/src/py3/plot/processdata.py (revision 23709) @@ -57,5 +57,5 @@ options.addfielddefault('clim',[lb,ub]) options.addfielddefault('cmap_set_under','1') - print("WARNING: nan's treated as", nanfill, "by default. Change using pairoption 'nan',nan_fill_value in plotmodel call") + print(("WARNING: nan's treated as", nanfill, "by default. Change using pairoption 'nan',nan_fill_value in plotmodel call")) # }}} # {{{ log @@ -118,5 +118,5 @@ spccol=options.getfieldvalue('spccol',0) print('multiple-column spc field; specify column to plot using option "spccol"') - print('column ', spccol, ' plotted for time: ', procdata[-1,spccol]) + print(('column ', spccol, ' plotted for time: ', procdata[-1,spccol])) procdata=procdata[0:-1,spccol] Index: /issm/trunk-jpl/src/py3/qmu/dakota_in_data.py =================================================================== --- /issm/trunk-jpl/src/py3/qmu/dakota_in_data.py (revision 23708) +++ /issm/trunk-jpl/src/py3/qmu/dakota_in_data.py (revision 23709) @@ -1,6 +1,8 @@ #move this stuff elsewhere from helpers import * -from dakota_in_write import dakota_in_write -from dakota_in_params import dakota_in_params + +from dakota_in_write import * +from dakota_in_params import * +from MatlabFuncs import * def dakota_in_data(dmeth,variables,responses,dparams,filei,*args): @@ -43,13 +45,12 @@ # get default set of parameters params=dakota_in_params(struct()) - # merge specified parameters into default set, whether or not # they already exist fnames=fieldnames(dparams) - for i in range(np.size(fnames)): - if not isfield(params,fnames[i]): - print('WARNING: dakota_in_data:unknown_param: No parameter '+str(fnames[i])+' in default parameter set.') - exec(('params.%s = vars(dparams)[fnames[i]]')%(fnames[i])) + for fieldname in fnames: + if not isfield(params,fieldname): + print('WARNING: dakota_in_data:unknown_param: No parameter {} in default parameter set.'.format(str(fieldname))) + exec('params.{} = vars(dparams)[fieldname]'.format(fieldname)) # use matlab even though we are running python Index: /issm/trunk-jpl/src/py3/qmu/dakota_in_params.py =================================================================== --- /issm/trunk-jpl/src/py3/qmu/dakota_in_params.py (revision 23708) +++ /issm/trunk-jpl/src/py3/qmu/dakota_in_params.py (revision 23709) @@ -182,3 +182,2 @@ return params - Index: /issm/trunk-jpl/src/py3/qmu/dakota_in_write.py =================================================================== --- /issm/trunk-jpl/src/py3/qmu/dakota_in_write.py (revision 23708) +++ /issm/trunk-jpl/src/py3/qmu/dakota_in_write.py (revision 23709) @@ -174,6 +174,6 @@ str_name = dmeth.variables[j] - # organize so that multiple instances of the same qmu class - # (2 different variable instances of "normal_uncertain" for example) + # organize so that multiple instances of the same qmu class + # (2 different variable instances of "normal_uncertain" for example) # are in the same dakota_write call regardless of individual size; # but that each class has its own dakota_write call @@ -213,5 +213,4 @@ elif params.system+params.fork+params.direct > 1: raise RuntimeError('Too many interfaces selected.') - if params.system or params.fork: param_write(fidi,'\t','asynchronous','','\n',params) @@ -230,12 +229,12 @@ if len(params.input_filter) != 0: param_write(fidi,'\t ','input_filter',' = ','\n',params) - + if len(params.output_filter) != 0: param_write(fidi,'\t ','output_filter',' = ','\n',params) - + param_write(fidi,'\t ','failure_capture',' ','\n',params) param_write(fidi,'\t ','deactivate',' ','\n',params) - param_write(fidi,'\t ','parameters_file',' = ','\n',params) - param_write(fidi,'\t ','results_file',' = ','\n',params) + param_write(fidi,'\t ','parameters_file',' = \'','\'\n',params) + param_write(fidi,'\t ','results_file',' = \'','\'\n',params) param_write(fidi,'\t ','verbatim', '','\n',params) param_write(fidi,'\t ','aprepro', '','\n',params) @@ -257,14 +256,14 @@ if ext != '': ext='.py' - + params.analysis_components=fullfile(pathstr,name+ext) param_write(fidi,'\t ','analysis_components',' = \'','\'\n',params) - + if len(params.input_filter) != 0: param_write(fidi,'\t ','input_filter',' = ','\n',params) - + if len(params.output_filter) != 0: param_write(fidi,'\t ','output_filter',' = ','\n',params) - + param_write(fidi,'\t ','failure_capture',' ','\n',params) param_write(fidi,'\t ','deactivate',' ','\n',params) @@ -312,5 +311,5 @@ elif (params.numerical_gradients+params.analytic_gradients > 1): raise RuntimeError('Too many gradients selected.') - + if params.numerical_gradients: param_write(fidi,'\t','numerical_gradients','','\n',params) @@ -329,5 +328,2 @@ else: fidi.write('\tno_hessians\n') - - - Index: /issm/trunk-jpl/src/py3/qmu/param_write.py =================================================================== --- /issm/trunk-jpl/src/py3/qmu/param_write.py (revision 23708) +++ /issm/trunk-jpl/src/py3/qmu/param_write.py (revision 23709) @@ -7,5 +7,5 @@ ''' if not isfield(params,pname): - print('WARNING: param_write:param_not_found: Parameter '+str(pname)+' not found in structure.') + print('WARNING: param_write:param_not_found: Parameter {} not found in structure.'.format(pname)) return @@ -19,9 +19,6 @@ elif type(params_pname) in [str]: - fidi.write(sbeg+str(pname)+smid+params_pname+s) + fidi.write(sbeg+pname+smid+params_pname+s) elif type(params_pname) in [int, float]: fidi.write(sbeg+str(pname)+smid+str(params_pname)+s) - - - Index: /issm/trunk-jpl/src/py3/qmu/rlev_write.py =================================================================== --- /issm/trunk-jpl/src/py3/qmu/rlev_write.py (revision 23708) +++ /issm/trunk-jpl/src/py3/qmu/rlev_write.py (revision 23709) @@ -1,7 +1,9 @@ import numpy as np +#move this later from helpers import * -from vector_write import vector_write -from param_write import param_write -from response_function import * +from vector_write import * +from param_write import * +#import relevent qmu classes +from MatlabArray import * def rlevi_write(fidi,ltype,levels): @@ -20,5 +22,4 @@ fidi.write('\n') - fidi.write('\t '+str(ltype)+' =\n') @@ -37,4 +38,5 @@ function to write response levels ''' + from response_function import * if len(dresp) == 0 or len(fieldnames(dresp[0])) == 0: @@ -77,5 +79,7 @@ grell.extend(grelli) + # write response levels + respl=allempty(respl) probl=allempty(probl) Index: /issm/trunk-jpl/src/py3/qmu/setupdesign/QmuSetupVariables.py =================================================================== --- /issm/trunk-jpl/src/py3/qmu/setupdesign/QmuSetupVariables.py (revision 23708) +++ /issm/trunk-jpl/src/py3/qmu/setupdesign/QmuSetupVariables.py (revision 23709) @@ -1,4 +1,5 @@ -from uniform_uncertain import uniform_uncertain -from normal_uncertain import normal_uncertain +from MatlabFuncs import * +from uniform_uncertain import* +from normal_uncertain import * from copy import deepcopy @@ -9,5 +10,5 @@ #decide whether this is a distributed variable, which will drive whether we expand it into npart values, - #or if we just carry it forward as is. + #or if we just carry it forward as is. #ok, key off according to type of descriptor: @@ -18,10 +19,10 @@ if ((type(variables.lower) in [list,np.ndarray] and len(variables.lower) > md.qmu.numberofpartitions) or (type(variables.upper) in [list,np.ndarray] and len(variables.upper) > md.qmu.numberofpartitions)): raise RuntimeError('QmuSetupDesign error message: upper and lower should be either a scalar or a "npart" length vector') - + elif isinstance(variables,normal_uncertain): if type(variables.stddev) in [list,np.ndarray] and len(variables.stddev) > md.qmu.numberofpartitions: raise RuntimeError('QmuSetupDesign error message: stddev should be either a scalar or a "npart" length vector') - #ok, dealing with semi-discrete distributed variable. Distribute according to how many + #ok, dealing with semi-discrete distributed variable. Distribute according to how many #partitions we want @@ -76,2 +77,3 @@ return dvar + Index: /issm/trunk-jpl/src/py3/qmu/vlist_write.py =================================================================== --- /issm/trunk-jpl/src/py3/qmu/vlist_write.py (revision 23708) +++ /issm/trunk-jpl/src/py3/qmu/vlist_write.py (revision 23709) @@ -2,6 +2,9 @@ #move this later from helpers import * -from vector_write import vector_write -from normal_uncertain import normal_uncertain + +from vector_write import * + +from uniform_uncertain import * +from normal_uncertain import * def check(a,l,p): @@ -38,4 +41,5 @@ if dvar == None: return + #from uniform_uncertain import * func = eval(cstring) Index: /issm/trunk-jpl/src/py3/solve/loadresultsfromcluster.py =================================================================== --- /issm/trunk-jpl/src/py3/solve/loadresultsfromcluster.py (revision 23708) +++ /issm/trunk-jpl/src/py3/solve/loadresultsfromcluster.py (revision 23709) @@ -37,7 +37,7 @@ md=loadresultsfromdisk(md,md.miscellaneous.name+'.outbin') else: - print('WARNING, outbin file is empty for run '+md.miscellaneous.name) + print(('WARNING, outbin file is empty for run '+md.miscellaneous.name)) elif not md.qmu.isdakota: - print('WARNING, outbin file does not exist '+md.miscellaneous.name) + print(('WARNING, outbin file does not exist '+md.miscellaneous.name)) #erase the log and output files @@ -88,3 +88,3 @@ os.remove(filename+extension) except OSError: - print('WARNING, no '+extension+' is present for run '+filename) + print(('WARNING, no '+extension+' is present for run '+filename)) Index: /issm/trunk-jpl/src/py3/solve/marshall.py =================================================================== --- /issm/trunk-jpl/src/py3/solve/marshall.py (revision 23708) +++ /issm/trunk-jpl/src/py3/solve/marshall.py (revision 23709) @@ -12,5 +12,5 @@ """ if md.verbose.solution: - print("marshalling file '%s.bin'." % md.miscellaneous.name) + print(("marshalling file '%s.bin'." % md.miscellaneous.name)) #open file for binary writing Index: /issm/trunk-jpl/src/py3/solve/waitonlock.py =================================================================== --- /issm/trunk-jpl/src/py3/solve/waitonlock.py (revision 23708) +++ /issm/trunk-jpl/src/py3/solve/waitonlock.py (revision 23709) @@ -27,5 +27,5 @@ print('solution launched on remote cluster. log in to detect job completion.') - choice=input('Is the job successfully completed? (y/n) ') + choice=eval(input('Is the job successfully completed? (y/n) ')) if not m.strcmp(choice,'y'): print('Results not loaded... exiting') @@ -44,5 +44,5 @@ etime=0 ispresent=0 - print("waiting for '%s' hold on... (Ctrl+C to exit)" % filename) + print(("waiting for '%s' hold on... (Ctrl+C to exit)" % filename)) #loop till file .lock exist or time is up