![(please configure the [header_logo] section in trac.ini)](/trac/issm/chrome/common/trac_banner.png)
Changeset 26660
- Timestamp:
- 11/23/21 12:54:53 (3 years ago)
- Location:
- issm/trunk-jpl/src/m/classes
- Files:
-
- 2 edited
-
stochasticforcing.m (modified) (10 diffs)
-
stochasticforcing.py (modified) (7 diffs)
Legend:
- Unmodified
- Added
- Removed
-
issm/trunk-jpl/src/m/classes/stochasticforcing.m
r26640 r26660 5 5 6 6 classdef stochasticforcing 7 properties (SetAccess=public) 8 isstochasticforcing 9 fields 10 defaultdimension 11 default_id= NaN;12 covariance 13 randomflag 7 properties (SetAccess=public) 8 isstochasticforcing = 0; 9 fields = NaN; 10 defaultdimension = 0; 11 default_id = NaN; 12 covariance = NaN; 13 randomflag = 1; 14 14 end 15 15 methods … … 26 26 end % }}} 27 27 function self = setdefaultparameters(self) % {{{ 28 self.isstochasticforcing 29 self.randomflag 28 self.isstochasticforcing = 0; %stochasticforcing is turned off by default 29 self.randomflag = 1; %true randomness is implemented by default 30 30 end % }}} 31 31 function md = checkconsistency(self,md,solution,analyses) % {{{ … … 35 35 36 36 num_fields = numel(self.fields); 37 37 38 38 %Check that covariance matrix is positive definite 39 39 try … … 43 43 end 44 44 45 %Check that all fields agree with the corresponding md class and if any field needs the default params 46 checkdefaults = false; %need to check defaults only if one of the fielddoes not have its own dimensionality47 structstoch 45 %Check that all fields agree with the corresponding md class and if any field needs the default params 46 checkdefaults = false; %need to check defaults only if one of the fields does not have its own dimensionality 47 structstoch = structstochforcing(); 48 48 for field=self.fields 49 50 49 %Checking agreement of classes 50 if(contains(field,'SMB')) 51 51 mdname = structstoch.mdnames(find(strcmp(structstoch.fields,char(field)))); 52 52 if~(isequal(class(md.smb),char(mdname))) 53 54 55 56 53 error('md.smb does not agree with stochasticforcing field %s', char(field)); 54 end 55 end 56 if(contains(field,'FrontalForcings')) 57 57 mdname = structstoch.mdnames(find(strcmp(structstoch.fields,char(field)))); 58 59 60 61 58 if~(isequal(class(md.frontalforcings),char(mdname))) 59 error('md.frontalforcings does not agree with stochasticforcing field %s', char(field)); 60 end 61 end 62 62 if(contains(field,'Calving')) 63 63 mdname = structstoch.mdnames(find(strcmp(structstoch.fields,char(field)))); 64 65 66 67 64 if~(isequal(class(md.calving),char(mdname))) 65 error('md.calving does not agree with stochasticforcing field %s', char(field)); 66 end 67 end 68 68 if(contains(field,'BasalforcingsFloatingice')) 69 69 mdname = structstoch.mdnames(find(strcmp(structstoch.fields,char(field)))); 70 71 72 73 74 75 76 77 78 70 if~(isequal(class(md.basalforcings),char(mdname))) 71 error('md.basalforcings does not agree with stochasticforcing field %s', char(field)); 72 end 73 end 74 %Checking for specific dimensions 75 if ~(strcmp(field,'SMBautoregression') || strcmp(field,'FrontalForcingsRignotAutoregression')) 76 checkdefaults = true; %field with non-specific dimensionality 77 end 78 end 79 79 80 80 %Retrieve sum of all the field dimensionalities 81 size_tot 82 indSMBar 83 indTFar 81 size_tot = self.defaultdimension*num_fields; 82 indSMBar = -1; %about to check for index of SMBautoregression 83 indTFar = -1; %about to check for index of FrontalForcingsRignotAutoregression 84 84 if any(contains(self.fields,'SMBautoregression')) 85 85 size_tot = size_tot-self.defaultdimension+md.smb.num_basins; 86 86 indSMBar = find(contains(self.fields,'SMBautoregression')); %index of SMBar, now check for consistency with TFar timestep (08Nov2021) 87 88 89 size_tot= size_tot-self.defaultdimension+md.frontalforcings.num_basins;90 indTFar 91 87 end 88 if any(contains(self.fields,'FrontalForcingsRignotAutoregression')) 89 size_tot = size_tot-self.defaultdimension+md.frontalforcings.num_basins; 90 indTFar = find(contains(self.fields,'FrontalForcingsRignotAutoregression')); %index of TFar, now check for consistency with SMBar timestep (08Nov2021) 91 end 92 92 93 93 if(indSMBar~=-1 && indTFar~=-1) %both autoregressive models are used: check autoregressive time step consistency … … 102 102 md = checkfield(md,'fieldname','stochasticforcing.randomflag','numel',[1],'values',[0 1]); 103 103 if(checkdefaults) %need to check the defaults 104 105 104 md = checkfield(md,'fieldname','stochasticforcing.defaultdimension','numel',1,'NaN',1,'Inf',1,'>',0); 105 md = checkfield(md,'fieldname','stochasticforcing.default_id','Inf',1,'>=',0,'<=',self.defaultdimension,'size',[md.mesh.numberofelements,1]); 106 106 end 107 107 end % }}} … … 111 111 fielddisplay(self,'fields','fields with stochasticity applied, ex: {''SMBautoregression''}, or {''FrontalForcingsRignotAutoregression''}'); 112 112 fielddisplay(self,'defaultdimension','dimensionality of the noise terms (does not apply to fields with their specific dimension)'); 113 113 fielddisplay(self,'default_id','id of each element for partitioning of the noise terms (does not apply to fields with their specific partition)'); 114 114 fielddisplay(self,'covariance','covariance matrix for within- and between-fields covariance (units must be squared field units)'); 115 115 fielddisplay(self,'randomflag','whether to apply real randomness (true) or pseudo-randomness with fixed seed (false)'); 116 116 disp('Available fields:'); 117 117 for field=supportedstochforcings() 118 fprintf(' %s 118 fprintf(' %s\n',string(field)); 119 119 end 120 120 end % }}} … … 164 164 WriteData(fid,prefix,'object',self,'fieldname','fields','format','StringArray'); 165 165 WriteData(fid,prefix,'data',dimensions,'name','md.stochasticforcing.dimensions','format','IntMat'); 166 166 WriteData(fid,prefix,'object',self,'fieldname','default_id','data',self.default_id-1,'format','IntMat','mattype',2); %0-indexed 167 167 WriteData(fid,prefix,'object',self,'fieldname','defaultdimension','format','Integer'); 168 168 WriteData(fid,prefix,'data',tempcovariance,'name','md.stochasticforcing.covariance','format','DoubleMat'); … … 173 173 end 174 174 function list = supportedstochforcings() % {{{ 175 176 175 % Defines list of fields supported 176 % by the class md.stochasticforcing 177 177 178 list = structstochforcing(); 178 list = structstochforcing(); 179 179 list = list.fields; 180 180 end % }}} … … 183 183 % supported and corresponding md names 184 184 structure.fields = {... 185 186 187 188 189 185 'DefaultCalving',... 186 'FloatingMeltRate',... 187 'FrontalForcingsRignotAutoregression',... 188 'SMBautoregression' 189 }; 190 190 structure.mdnames = {... 191 191 'calving',... … … 195 195 }; 196 196 end % }}} 197 198 199 200 -
issm/trunk-jpl/src/m/classes/stochasticforcing.py
r26647 r26660 42 42 def setdefaultparameters(self): # {{{ 43 43 # Type of stabilization used 44 self.isstochasticforcing = 0 45 self.fields = [] 46 self.randomflag = 1 44 self.isstochasticforcing = 0 # stochasticforcing is turned off by default 45 self.fields = [] # Need to initialize to list to avoid "RuntimeError: object of type 'float' has no len()" on import of class 46 self.randomflag = 1 # true randomness is implemented by default 47 47 return self 48 48 #}}} … … 55 55 num_fields = len(self.fields) 56 56 57 #Check that covariance matrix is positive definite this is done internaly by linalg 58 np.linalg.cholesky(self.covariance) 57 # Check that covariance matrix is positive definite (this is done internally by linalg) 58 try: 59 np.linalg.cholesky(self.covariance) 60 except LinAlgError: 61 error('md.stochasticforcing.covariance is not positive definite') 59 62 60 # Check that all fields agree with the corresponding md class 61 checkdefaults = False 63 # Check that all fields agree with the corresponding md class and if any field needs the default params 64 checkdefaults = False # Need to check defaults only if one of the fields does not have its own dimensionality 62 65 structstoch = self.structstochforcing() 63 66 for field in self.fields: 64 # Check s are temporarily commented out: need to debug why this does not work (18Nov2021) #67 # Checking agreement of classes 65 68 if 'SMB' in field: 66 69 mdname = structstoch[field] … … 78 81 mdname = structstoch[field] 79 82 if (type(md.basalforcings).__name__ != mdname): 80 raise TypeError('md.basalforcings does not agree with stochasticforcing field {}'.format(mdname)) #Checking for specific dimensions 83 raise TypeError('md.basalforcings does not agree with stochasticforcing field {}'.format(mdname)) 84 # Checking for specific dimensions 81 85 if not (field == 'SMBautoregression' or field == 'FrontalForcingsRignotAutoregression'): 82 checkdefaults = True #field with non-specific dimensionality86 checkdefaults = True # field with non-specific dimensionality 83 87 84 # Retrieve sum of all the field dimensionalities88 # Retrieve sum of all the field dimensionalities 85 89 size_tot = self.defaultdimension * num_fields 86 indSMBar = -1 #about to check for index of SMBautoregression87 indTFar = -1 #about to check for index of FrontalForcingsRignotAutoregression90 indSMBar = -1 # About to check for index of SMBautoregression 91 indTFar = -1 # About to check for index of FrontalForcingsRignotAutoregression 88 92 if ('SMBautoregression' in self.fields): 89 93 size_tot = size_tot - self.defaultdimension + md.smb.num_basins 90 indSMBar = self.fields.index('SMBautoregression') 94 indSMBar = self.fields.index('SMBautoregression') # Index of SMBar, now check for consistency with TFar timestep (08Nov2021) 91 95 if ('FrontalForcingsRignotAutoregression' in self.fields): 92 96 size_tot = size_tot - self.defaultdimension + md.frontalforcings.num_basins 93 indTFar = self.fields.index('FrontalForcingsRignotAutoregression') 94 if (indSMBar != -1 and indTFar != -1): 97 indTFar = self.fields.index('FrontalForcingsRignotAutoregression') # Index of TFar, now check for consistency with SMBar timestep (08Nov2021) 98 if (indSMBar != -1 and indTFar != -1): # Both autoregressive models are used: check autoregressive time step consistency 95 99 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)] 96 100 if((md.smb.ar_timestep != md.frontalforcings.ar_timestep) and np.any(covsum != 0)): … … 99 103 md = checkfield(md, 'fieldname', 'stochasticforcing.isstochasticforcing', 'values', [0, 1]) 100 104 md = checkfield(md, 'fieldname', 'stochasticforcing.fields', 'numel', num_fields, 'cell', 1, 'values', self.supportedstochforcings()) 101 #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 list102 105 md = checkfield(md, 'fieldname', 'stochasticforcing.covariance', 'NaN', 1, 'Inf', 1, 'size', [size_tot, size_tot]) # global covariance matrix 103 106 md = checkfield(md, 'fieldname', 'stochasticforcing.randomflag', 'numel', [1], 'values', [0, 1]) … … 121 124 return md 122 125 else: 123 # Retrieve dimensionality of each field126 # Retrieve dimensionality of each field 124 127 dimensions = self.defaultdimension * np.ones(num_fields) 125 128 for ind, field in enumerate(self.fields): 126 # Checking for specific dimensions129 # Checking for specific dimensions 127 130 if (field == 'SMBautoregression'): 128 131 dimensions[ind] = md.smb.num_basins … … 131 134 132 135 # Scaling covariance matrix (scale column-by-column and row-by-row) 133 scaledfields = ['DefaultCalving', 'SMBautoregression'] 136 scaledfields = ['DefaultCalving', 'SMBautoregression'] # list of fields that need scaling * 1/yts 134 137 tempcovariance = np.copy(self.covariance) 135 138 for i in range(num_fields): … … 140 143 for col in inds: # scale columns corresponding to scaled field 141 144 tempcovariance[:, col] = 1 / yts * tempcovariance[:, col] 142 # Set dummy default_id vector if defaults not used145 # Set dummy default_id vector if defaults not used 143 146 if np.isnan(self.default_id): 144 147 self.default_id = np.zeros(md.mesh.numberofelements)
Note:
See TracChangeset
for help on using the changeset viewer.
