Index: /issm/trunk-jpl/src/m/classes/stochasticforcing.m =================================================================== --- /issm/trunk-jpl/src/m/classes/stochasticforcing.m (revision 26659) +++ /issm/trunk-jpl/src/m/classes/stochasticforcing.m (revision 26660) @@ -5,11 +5,11 @@ classdef stochasticforcing - properties (SetAccess=public) - isstochasticforcing = 0; - fields = NaN; - defaultdimension = 0; - default_id = NaN; - covariance = NaN; - randomflag = 1; + properties (SetAccess=public) + isstochasticforcing = 0; + fields = NaN; + defaultdimension = 0; + default_id = NaN; + covariance = NaN; + randomflag = 1; end methods @@ -26,6 +26,6 @@ end % }}} function self = setdefaultparameters(self) % {{{ - self.isstochasticforcing = 0; %stochasticforcing is turned off by default - self.randomflag = 1; %true randomness is implemented by default + self.isstochasticforcing = 0; %stochasticforcing is turned off by default + self.randomflag = 1; %true randomness is implemented by default end % }}} function md = checkconsistency(self,md,solution,analyses) % {{{ @@ -35,5 +35,5 @@ num_fields = numel(self.fields); - + %Check that covariance matrix is positive definite try @@ -43,51 +43,51 @@ end - %Check that all fields agree with the corresponding md class and if any field needs the default params - checkdefaults = false; %need to check defaults only if one of the field does not have its own dimensionality - structstoch = structstochforcing(); + %Check that all fields agree with the corresponding md class and if any field needs the default params + checkdefaults = false; %need to check defaults only if one of the fields does not have its own dimensionality + structstoch = structstochforcing(); for field=self.fields - %Checking agreement of classes - if(contains(field,'SMB')) + %Checking agreement of classes + if(contains(field,'SMB')) mdname = structstoch.mdnames(find(strcmp(structstoch.fields,char(field)))); if~(isequal(class(md.smb),char(mdname))) - error('md.smb does not agree with stochasticforcing field %s', char(field)); - end - end - if(contains(field,'FrontalForcings')) + error('md.smb does not agree with stochasticforcing field %s', char(field)); + end + end + if(contains(field,'FrontalForcings')) mdname = structstoch.mdnames(find(strcmp(structstoch.fields,char(field)))); - if~(isequal(class(md.frontalforcings),char(mdname))) - error('md.frontalforcings does not agree with stochasticforcing field %s', char(field)); - end - end + if~(isequal(class(md.frontalforcings),char(mdname))) + error('md.frontalforcings does not agree with stochasticforcing field %s', char(field)); + end + end if(contains(field,'Calving')) mdname = structstoch.mdnames(find(strcmp(structstoch.fields,char(field)))); - if~(isequal(class(md.calving),char(mdname))) - error('md.calving does not agree with stochasticforcing field %s', char(field)); - end - end + if~(isequal(class(md.calving),char(mdname))) + error('md.calving does not agree with stochasticforcing field %s', char(field)); + end + end if(contains(field,'BasalforcingsFloatingice')) mdname = structstoch.mdnames(find(strcmp(structstoch.fields,char(field)))); - if~(isequal(class(md.basalforcings),char(mdname))) - error('md.basalforcings does not agree with stochasticforcing field %s', char(field)); - end - end - %Checking for specific dimensions - if ~(strcmp(field,'SMBautoregression') || strcmp(field,'FrontalForcingsRignotAutoregression')) - checkdefaults = true; %field with non-specific dimensionality - end - end + if~(isequal(class(md.basalforcings),char(mdname))) + error('md.basalforcings does not agree with stochasticforcing field %s', char(field)); + end + end + %Checking for specific dimensions + if ~(strcmp(field,'SMBautoregression') || strcmp(field,'FrontalForcingsRignotAutoregression')) + checkdefaults = true; %field with non-specific dimensionality + end + end %Retrieve sum of all the field dimensionalities - size_tot = self.defaultdimension*num_fields; - indSMBar = -1; %about to check for index of SMBautoregression - indTFar = -1; %about to check for index of FrontalForcingsRignotAutoregression + size_tot = self.defaultdimension*num_fields; + indSMBar = -1; %about to check for index of SMBautoregression + indTFar = -1; %about to check for index of FrontalForcingsRignotAutoregression if any(contains(self.fields,'SMBautoregression')) - size_tot = size_tot-self.defaultdimension+md.smb.num_basins; + size_tot = size_tot-self.defaultdimension+md.smb.num_basins; indSMBar = find(contains(self.fields,'SMBautoregression')); %index of SMBar, now check for consistency with TFar timestep (08Nov2021) - end - if any(contains(self.fields,'FrontalForcingsRignotAutoregression')) - size_tot = size_tot-self.defaultdimension+md.frontalforcings.num_basins; - indTFar = find(contains(self.fields,'FrontalForcingsRignotAutoregression')); %index of TFar, now check for consistency with SMBar timestep (08Nov2021) - end + end + if any(contains(self.fields,'FrontalForcingsRignotAutoregression')) + size_tot = size_tot-self.defaultdimension+md.frontalforcings.num_basins; + indTFar = find(contains(self.fields,'FrontalForcingsRignotAutoregression')); %index of TFar, now check for consistency with SMBar timestep (08Nov2021) + end if(indSMBar~=-1 && indTFar~=-1) %both autoregressive models are used: check autoregressive time step consistency @@ -102,6 +102,6 @@ md = checkfield(md,'fieldname','stochasticforcing.randomflag','numel',[1],'values',[0 1]); if(checkdefaults) %need to check the defaults - md = checkfield(md,'fieldname','stochasticforcing.defaultdimension','numel',1,'NaN',1,'Inf',1,'>',0); - md = checkfield(md,'fieldname','stochasticforcing.default_id','Inf',1,'>=',0,'<=',self.defaultdimension,'size',[md.mesh.numberofelements,1]); + md = checkfield(md,'fieldname','stochasticforcing.defaultdimension','numel',1,'NaN',1,'Inf',1,'>',0); + md = checkfield(md,'fieldname','stochasticforcing.default_id','Inf',1,'>=',0,'<=',self.defaultdimension,'size',[md.mesh.numberofelements,1]); end end % }}} @@ -111,10 +111,10 @@ fielddisplay(self,'fields','fields with stochasticity applied, ex: {''SMBautoregression''}, or {''FrontalForcingsRignotAutoregression''}'); fielddisplay(self,'defaultdimension','dimensionality of the noise terms (does not apply to fields with their specific dimension)'); - fielddisplay(self,'default_id','id of each element for partitioning of the noise terms (does not apply to fields with their specific partition)'); + fielddisplay(self,'default_id','id of each element for partitioning of the noise terms (does not apply to fields with their specific partition)'); fielddisplay(self,'covariance','covariance matrix for within- and between-fields covariance (units must be squared field units)'); fielddisplay(self,'randomflag','whether to apply real randomness (true) or pseudo-randomness with fixed seed (false)'); disp('Available fields:'); for field=supportedstochforcings() - fprintf(' %s \n',string(field)); + fprintf(' %s\n',string(field)); end end % }}} @@ -164,5 +164,5 @@ WriteData(fid,prefix,'object',self,'fieldname','fields','format','StringArray'); WriteData(fid,prefix,'data',dimensions,'name','md.stochasticforcing.dimensions','format','IntMat'); - WriteData(fid,prefix,'object',self,'fieldname','default_id','data',self.default_id-1,'format','IntMat','mattype',2); %0-indexed + WriteData(fid,prefix,'object',self,'fieldname','default_id','data',self.default_id-1,'format','IntMat','mattype',2); %0-indexed WriteData(fid,prefix,'object',self,'fieldname','defaultdimension','format','Integer'); WriteData(fid,prefix,'data',tempcovariance,'name','md.stochasticforcing.covariance','format','DoubleMat'); @@ -173,8 +173,8 @@ end function list = supportedstochforcings() % {{{ - % Defines list of fields supported - % by the class md.stochasticforcing + % Defines list of fields supported + % by the class md.stochasticforcing - list = structstochforcing(); + list = structstochforcing(); list = list.fields; end % }}} @@ -183,9 +183,9 @@ % supported and corresponding md names structure.fields = {... - 'DefaultCalving',... - 'FloatingMeltRate',... - 'FrontalForcingsRignotAutoregression',... - 'SMBautoregression' - }; + 'DefaultCalving',... + 'FloatingMeltRate',... + 'FrontalForcingsRignotAutoregression',... + 'SMBautoregression' + }; structure.mdnames = {... 'calving',... @@ -195,6 +195,2 @@ }; end % }}} - - - - Index: /issm/trunk-jpl/src/m/classes/stochasticforcing.py =================================================================== --- /issm/trunk-jpl/src/m/classes/stochasticforcing.py (revision 26659) +++ /issm/trunk-jpl/src/m/classes/stochasticforcing.py (revision 26660) @@ -42,7 +42,7 @@ def setdefaultparameters(self): # {{{ # Type of stabilization used - self.isstochasticforcing = 0 # stochasticforcing is turned off by default - self.fields = [] # Need to initialize to list to avoid "RuntimeError: object of type 'float' has no len()" on import of class - self.randomflag = 1 # true randomness is implemented by default + self.isstochasticforcing = 0 # stochasticforcing is turned off by default + self.fields = [] # Need to initialize to list to avoid "RuntimeError: object of type 'float' has no len()" on import of class + self.randomflag = 1 # true randomness is implemented by default return self #}}} @@ -55,12 +55,15 @@ num_fields = len(self.fields) - #Check that covariance matrix is positive definite this is done internaly by linalg - np.linalg.cholesky(self.covariance) + # Check that covariance matrix is positive definite (this is done internally by linalg) + try: + np.linalg.cholesky(self.covariance) + except LinAlgError: + error('md.stochasticforcing.covariance is not positive definite') - # Check that all fields agree with the corresponding md class - checkdefaults = False + # Check that all fields agree with the corresponding md class and if any field needs the default params + checkdefaults = False # Need to check defaults only if one of the fields does not have its own dimensionality structstoch = self.structstochforcing() for field in self.fields: - # Checks are temporarily commented out: need to debug why this does not work (18Nov2021) # + # Checking agreement of classes if 'SMB' in field: mdname = structstoch[field] @@ -78,19 +81,20 @@ mdname = structstoch[field] if (type(md.basalforcings).__name__ != mdname): - raise TypeError('md.basalforcings does not agree with stochasticforcing field {}'.format(mdname)) #Checking for specific dimensions + raise TypeError('md.basalforcings does not agree with stochasticforcing field {}'.format(mdname)) + # Checking for specific dimensions if not (field == 'SMBautoregression' or field == 'FrontalForcingsRignotAutoregression'): - checkdefaults = True #field with non-specific dimensionality + checkdefaults = True # field with non-specific dimensionality - #Retrieve sum of all the field dimensionalities + # Retrieve sum of all the field dimensionalities size_tot = self.defaultdimension * num_fields - indSMBar = -1 #about to check for index of SMBautoregression - indTFar = -1 #about to check for index of FrontalForcingsRignotAutoregression + indSMBar = -1 # About to check for index of SMBautoregression + indTFar = -1 # About to check for index of FrontalForcingsRignotAutoregression if ('SMBautoregression' in self.fields): size_tot = size_tot - self.defaultdimension + md.smb.num_basins - indSMBar = self.fields.index('SMBautoregression') + indSMBar = self.fields.index('SMBautoregression') # Index of SMBar, now check for consistency with TFar timestep (08Nov2021) if ('FrontalForcingsRignotAutoregression' in self.fields): size_tot = size_tot - self.defaultdimension + md.frontalforcings.num_basins - indTFar = self.fields.index('FrontalForcingsRignotAutoregression') - if (indSMBar != -1 and indTFar != -1): + indTFar = self.fields.index('FrontalForcingsRignotAutoregression') # Index of TFar, now check for consistency with SMBar timestep (08Nov2021) + if (indSMBar != -1 and indTFar != -1): # Both autoregressive models are used: check autoregressive time step consistency covsum = self.covariance[np.sum(self.defaultdimensions[0:indSMBar]).astype(int):np.sum(self.defaultdimensions[0:indSMBar + 1]).astype(int), np.sum(self.defaultdimensions[0:indTFar]).astype(int):np.sum(self.defaultdimensions[0:indTFar + 1]).astype(int)] if((md.smb.ar_timestep != md.frontalforcings.ar_timestep) and np.any(covsum != 0)): @@ -99,5 +103,4 @@ md = checkfield(md, 'fieldname', 'stochasticforcing.isstochasticforcing', 'values', [0, 1]) md = checkfield(md, 'fieldname', 'stochasticforcing.fields', 'numel', num_fields, 'cell', 1, 'values', self.supportedstochforcings()) - #md = checkfield(md, 'fieldname', 'stochasticforcing.dimensions', 'NaN', 1, 'Inf', 1, '>', 0, 'size', [num_fields]) # specific dimension for each field; NOTE: As opposed to MATLAB implementation, pass list md = checkfield(md, 'fieldname', 'stochasticforcing.covariance', 'NaN', 1, 'Inf', 1, 'size', [size_tot, size_tot]) # global covariance matrix md = checkfield(md, 'fieldname', 'stochasticforcing.randomflag', 'numel', [1], 'values', [0, 1]) @@ -121,8 +124,8 @@ return md else: - #Retrieve dimensionality of each field + # Retrieve dimensionality of each field dimensions = self.defaultdimension * np.ones(num_fields) for ind, field in enumerate(self.fields): - #Checking for specific dimensions + # Checking for specific dimensions if (field == 'SMBautoregression'): dimensions[ind] = md.smb.num_basins @@ -131,5 +134,5 @@ # Scaling covariance matrix (scale column-by-column and row-by-row) - scaledfields = ['DefaultCalving', 'SMBautoregression'] # list of fields that need scaling * 1/yts + scaledfields = ['DefaultCalving', 'SMBautoregression'] # list of fields that need scaling * 1/yts tempcovariance = np.copy(self.covariance) for i in range(num_fields): @@ -140,5 +143,5 @@ for col in inds: # scale columns corresponding to scaled field tempcovariance[:, col] = 1 / yts * tempcovariance[:, col] - #Set dummy default_id vector if defaults not used + # Set dummy default_id vector if defaults not used if np.isnan(self.default_id): self.default_id = np.zeros(md.mesh.numberofelements)