Index: /issm/trunk-jpl/src/m/classes/SMBautoregression.m =================================================================== --- /issm/trunk-jpl/src/m/classes/SMBautoregression.m (revision 26485) +++ /issm/trunk-jpl/src/m/classes/SMBautoregression.m (revision 26486) @@ -13,5 +13,5 @@ ar_timestep = 0; phi = NaN; - covmat = NaN; + covmat = NaN; basin_id = NaN; randomflag = 1; @@ -43,13 +43,13 @@ self.ar_order = 1; %dummy 1 value for autoregression self.phi = zeros(self.num_basins,self.ar_order); %autoregression coefficients all set to 0 - disp(' smb.ar_order (order of autoregressive model) not specified: order of autoregressive model set to 0'); + disp(' smb.ar_order (order of autoregressive model) not specified: order of autoregressive model set to 0'); end if (self.ar_initialtime==0) self.ar_initialtime = md.timestepping.start_time; %autoregression model has no prescribed initial time - disp(' smb.ar_initialtime (initial time in the autoregressive model parameterization) not specified: set to md.timestepping.start_time'); + disp(' smb.ar_initialtime (initial time in the autoregressive model parameterization) not specified: set to md.timestepping.start_time'); end if (self.ar_timestep==0) self.ar_timestep = md.timestepping.time_step; %autoregression model has no prescribed time step - disp(' smb.ar_timestep (timestep of autoregressive model) not specified: set to md.timestepping.time_step'); + disp(' smb.ar_timestep (timestep of autoregressive model) not specified: set to md.timestepping.time_step'); end if isnan(self.phi) @@ -114,11 +114,11 @@ WriteData(fid,prefix,'object',self,'class','smb','fieldname','ar_timestep','format','Double','scale',yts); WriteData(fid,prefix,'object',self,'class','smb','fieldname','basin_id','data',self.basin_id-1,'name','md.smb.basin_id','format','IntMat','mattype',2); %0-indexed - WriteData(fid,prefix,'object',self,'class','smb','fieldname','beta0','format','DoubleMat','name','md.smb.beta0','scale',1./yts,'yts',md.constants.yts); - WriteData(fid,prefix,'object',self,'class','smb','fieldname','beta1','format','DoubleMat','name','md.smb.beta1','scale',1./(yts^2),'yts',md.constants.yts); - WriteData(fid,prefix,'object',self,'class','smb','fieldname','phi','format','DoubleMat','name','md.smb.phi','yts',md.constants.yts); - WriteData(fid,prefix,'object',self,'class','smb','fieldname','covmat','format','DoubleMat','name','md.smb.covmat','scale',1./(yts^2),'yts',md.constants.yts); + WriteData(fid,prefix,'object',self,'class','smb','fieldname','beta0','format','DoubleMat','name','md.smb.beta0','scale',1./yts,'yts',yts); + WriteData(fid,prefix,'object',self,'class','smb','fieldname','beta1','format','DoubleMat','name','md.smb.beta1','scale',1./(yts^2),'yts',yts); + WriteData(fid,prefix,'object',self,'class','smb','fieldname','phi','format','DoubleMat','name','md.smb.phi','yts',yts); + WriteData(fid,prefix,'object',self,'class','smb','fieldname','covmat','format','DoubleMat','name','md.smb.covmat','scale',1./(yts^2),'yts',yts); WriteData(fid,prefix,'object',self,'class','smb','fieldname','randomflag','format','Boolean'); - WriteData(fid, prefix, 'object', self, 'fieldname', 'steps_per_step', 'format', 'Integer'); - WriteData(fid, prefix, 'object', self, 'fieldname', 'averaging', 'format', 'Integer'); + WriteData(fid,prefix,'object',self,'fieldname','steps_per_step','format','Integer'); + WriteData(fid,prefix,'object',self,'fieldname','averaging','format','Integer'); %process requested outputs Index: /issm/trunk-jpl/src/m/classes/SMBautoregression.py =================================================================== --- /issm/trunk-jpl/src/m/classes/SMBautoregression.py (revision 26486) +++ /issm/trunk-jpl/src/m/classes/SMBautoregression.py (revision 26486) @@ -0,0 +1,138 @@ +import numpy as np + +from checkfield import * +from fielddisplay import fielddisplay +from project3d import * +from WriteData import * + + +class SMBcomponents(object): + """SMBAUTOREGRESSION class definition + + Usage: + SMBautoregression = SMBautoregression() + """ + + def __init__(self, *args): # {{{ + self.num_basins = 0 + self.beta0 = np.nan + self.beta1 = np.nan + self.ar_order = 0 + self.ar_initialtime = 0 + self.ar_timestep = 0 + self.phi = np.nan + self.covmat = np.nan + self.basin_id = np.nan + self.randomflag = 1 + self.steps_per_step = 1 + self.averaging = 0 + self.requested_outputs = [] + + nargs = len(args) + if nargs == 0: + self.setdefaultparameters() + else: + raise Exception('constructor not supported') + # }}} + + def __repr__(self): # {{{ + s = ' surface forcings parameters:\n' + s += '{}\n'.format(fielddisplay(self, 'num_basins', 'number of different basins [unitless]')) + s += '{}\n'.format(fielddisplay(self, 'basin_id', 'basin number assigned to each element [unitless]')) + s += '{}\n'.format(fielddisplay(self, 'beta0', 'basin-specific intercept values [m ice eq./yr]')) + s += '{}\n'.format(fielddisplay(self, 'beta1', 'basin-specific trend values [m ice eq. yr^(-2)]')) + s += '{}\n'.format(fielddisplay(self, 'ar_order', 'order of the autoregressive model [unitless]')) + s += '{}\n'.format(fielddisplay(self, 'ar_initialtime', 'initial time assumed in the autoregressive model parameterization [yr]')) + s += '{}\n'.format(fielddisplay(self, 'ar_timestep', 'time resolution of the autoregressive model [yr]')) + s += '{}\n'.format(fielddisplay(self, 'phi', 'basin-specific vectors of lag coefficients [unitless]')) + s += '{}\n'.format(fielddisplay(self, 'covmat', 'inter-basin covariance matrix for multivariate normal noise at each time step [m^2 ice eq. yr^(-2)]')) + s += '{}\n'.format(fielddisplay(self, 'randomflag', 'whether to apply real randomness (true) or pseudo-randomness with fixed seed (false)')) + s += '{}\n'.format(fielddisplay(self, 'steps_per_step', 'number of smb steps per time step')) + s += '{}\n'.format(fielddisplay(self, 'averaging', 'averaging methods from short to long steps')) + s += '\t\t{}\n'.format('0: Arithmetic (default)') + s += '\t\t{}\n'.format('1: Geometric') + s += '\t\t{}\n'.format('2: Harmonic') + s += '{}\n'.format(fielddisplay(self, 'requested_outputs', 'additional outputs requested')) + return s + # }}} + + def setdefaultparameters(self): #{{{ + self.ar_order = 0.0 # Autoregression model of order 0 + self.randomflag = 1 + # }}} + + def extrude(self, md): # {{{ + return self # Nothing for now + # }}} + + def defaultoutputs(self, md): # {{{ + return [] + # }}} + + def initialize(self, md): # {{{ + if np.all(np.isnan(self.beta1)): + self.beta1 = np.zeros((1, self.num_basins)) # No trend in SMB + print(' smb.beta1 (trend) not specified: value set to 0') + if self.ar_order == 0: + self.ar_order = 1 # Dummy 1 value for autoregression + self.phi = np.zeros((self.num_basins, self.ar_order)) # Autorgression coefficients all set to 0 + print(' smb.ar_order (order of autoregressive model) not specified: order of autoregressive model set to 0') + if self.ar_initialtime == 0: + self.ar_initialtime = md.timestepping.start_time # Autoregression model has no prescribed initial time + print(' smb.ar_initialtime (initial time in the autoregressive model parameterization) not specified: set to md.timestepping.start_time') + if self.ar_timestep == 0: + self.ar_timestep = md.timestepping.time_step # Autoregression model has no prescribed time step + print(' smb.ar_timestep (timestep of autoregressive model) not specified: set to md.timestepping.time_step') + if np.all(np.isnan(self.phi)): + self.phi = np.zeros((self.num_basins, self.ar_order)) # Autoregression model of order 0 + print(' smb.phi (lag coefficients) not specified: order of autoregressive model set to 0') + if np.all(np.isnan(self.covmat)): + self.covmat = 1e-21 * np.eye(self.num_basins) # No stochasticity and no covariance + print(' smb.covmat not specified: stochasticity set to 0') + return self + # }}} + + def checkconsistency(self, md, solution, analyses): # {{{ + if 'MasstransportAnalysis' in analyses: + md = checkfield(md, 'fieldname', 'smb.num_basins', 'numel', 1, 'NaN', 1, 'Inf', 1, '>', 0) + md = checkfield(md, 'fieldname', 'smb.basin_id', 'Inf', 1, '>=', 0, '<=', md.smb.num_basins, 'size', [md.mesh.numberofelements]) + md = checkfield(md, 'fieldname', 'smb.beta0', 'NaN', 1, 'Inf', 1, 'size', [1, md.smb.num_basins], 'numel', md.smb.num_basins) # Scheme fails if passed as column vector + md = checkfield(md, 'fieldname', 'smb.beta1', 'NaN', 1, 'Inf', 1, 'size', [1, md.smb.num_basins], 'numel', md.smb.num_basins) # Scheme fails if passed as column vector + md = checkfield(md, 'fieldname', 'smb.ar_order', 'numel', 1, 'NaN', 1, 'Inf', 1, '>=', 0) + md = checkfield(md, 'fieldname', 'smb.ar_initialtime', 'numel', 1, 'NaN', 1, 'Inf', 1) + md = checkfield(md, 'fieldname', 'smb.ar_timestep', 'numel', 1, 'NaN', 1, 'Inf', 1, '>=', md.timestepping.time_step) # Autoregression time step cannot be finer than ISSM timestep + md = checkfield(md, 'fieldname', 'smb.phi', 'NaN', 1, 'Inf', 1, 'size', [md.smb.num_basins, md.smb.ar_order]) + md = checkfield(md, 'fieldname', 'smb.covmat', 'NaN', 1, 'Inf', 1, 'size', [md.smb.num_basins, md.smb.num_basins]) + md = checkfield(md, 'fieldname', 'smb.randomflag', 'numel', [1], 'values', [0, 1]) + md = checkfield(md, 'fieldname', 'smb.steps_per_step', '>=', 1, 'numel', [1]) + md = checkfield(md, 'fieldname', 'smb.averaging', 'numel', [1], 'values', [0, 1, 2]) + md = checkfield(md, 'fieldname', 'smb.requested_outputs', 'stringrow', 1) + return md + # }}} + + def marshall(self, prefix, md, fid): # {{{ + yts = md.constants.yts + + WriteData(fid, prefix, 'name', 'md.smb.model', 'data', 55, 'format', 'Integer') + WriteData(fid, prefix, 'object', self, 'class', 'smb', 'fieldname', 'num_basins', 'format', 'Integer') + WriteData(fid, prefix, 'object', self, 'class', 'smb', 'fieldname', 'ar_order', 'format', 'Integer') + WriteData(fid, prefix, 'object', self, 'class', 'smb', 'fieldname', 'ar_initialtime', 'format', 'Double', 'scale', yts) + WriteData(fid, prefix, 'object', self, 'class', 'smb', 'fieldname', 'ar_timestep', 'format', 'Double', 'scale', yts) + WriteData(fid, prefix, 'object', self, 'class', 'smb', 'fieldname', 'basin_id', 'data', self.basin_id, 'name', 'md.smb.basin_id', 'format', 'IntMat', 'mattype', 2) # 0-indexed + WriteData(fid, prefix, 'object', self, 'class', 'smb', 'fieldname', 'beta0', 'format', 'DoubleMat', 'name', 'md.smb.beta0', 'scale', 1 / yts, 'yts', yts) + WriteData(fid, prefix, 'object', self, 'class', 'smb', 'fieldname', 'beta1', 'format', 'DoubleMat', 'name', 'md.smb.beta1', 'scale', 1 / (yts ** 2), 'yts', yts) + WriteData(fid, prefix, 'object', self, 'class', 'smb', 'fieldname', 'phi', 'format', 'DoubleMat', 'name', 'md.smb.phi', 'yts', yts) + WriteData(fid, prefix, 'object', self, 'class', 'smb', 'fieldname', 'covmat', 'format', 'DoubleMat', 'name', 'md.smb.covmat', 'scale', 1 / (yts ** 2), 'yts', yts) + WriteData(fid, prefix, 'object', self, 'class', 'smb', 'fieldname', 'randomflag', 'format', 'Boolean') + WriteData(fid, prefix, 'object', self, 'fieldname', 'steps_per_step', 'format', 'Integer') + WriteData(fid, prefix, 'object', self, 'fieldname', 'averaging', 'format', 'Integer') + + # Process requested outputs + outputs = self.requested_outputs + indices = [i for i, x in enumerate(outputs) if x == 'default'] + if len(indices) > 0: + 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/m/classes/SMBcomponents.m =================================================================== --- /issm/trunk-jpl/src/m/classes/SMBcomponents.m (revision 26485) +++ /issm/trunk-jpl/src/m/classes/SMBcomponents.m (revision 26486) @@ -83,7 +83,7 @@ WriteData(fid,prefix,'name','md.smb.model','data',2,'format','Integer'); - WriteData(fid,prefix,'object',self,'class','smb','fieldname','accumulation','format','DoubleMat','mattype',1,'scale',1./yts,'timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts); - WriteData(fid,prefix,'object',self,'class','smb','fieldname','runoff','format','DoubleMat','mattype',1,'scale',1./yts,'timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts); - WriteData(fid,prefix,'object',self,'class','smb','fieldname','evaporation','format','DoubleMat','mattype',1,'scale',1./yts,'timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts); + WriteData(fid,prefix,'object',self,'class','smb','fieldname','accumulation','format','DoubleMat','mattype',1,'scale',1./yts,'timeserieslength',md.mesh.numberofvertices+1,'yts',yts); + WriteData(fid,prefix,'object',self,'class','smb','fieldname','runoff','format','DoubleMat','mattype',1,'scale',1./yts,'timeserieslength',md.mesh.numberofvertices+1,'yts',yts); + WriteData(fid,prefix,'object',self,'class','smb','fieldname','evaporation','format','DoubleMat','mattype',1,'scale',1./yts,'timeserieslength',md.mesh.numberofvertices+1,'yts',yts); WriteData(fid, prefix, 'object', self, 'fieldname', 'steps_per_step', 'format', 'Integer'); WriteData(fid, prefix, 'object', self, 'fieldname', 'averaging', 'format', 'Integer'); Index: /issm/trunk-jpl/src/m/classes/SMBcomponents.py =================================================================== --- /issm/trunk-jpl/src/m/classes/SMBcomponents.py (revision 26485) +++ /issm/trunk-jpl/src/m/classes/SMBcomponents.py (revision 26486) @@ -61,8 +61,8 @@ if np.all(np.isnan(self.evaporation)): self.evaporation = np.zeros((md.mesh.numberofvertices)) - print(" no SMB.evaporation specified: values set as zero") + print(' no SMB.evaporation specified: values set as zero') if np.all(np.isnan(self.runoff)): self.runoff = np.zeros((md.mesh.numberofvertices)) - print(" no SMB.runoff specified: values set as zero") + print(' no SMB.runoff specified: values set as zero') return self # }}}