Changeset 26604 for issm/trunk-jpl/src/m/classes/stochasticforcing.m

Timestamp:

11/11/21 09:33:47 (3 years ago)

Author:

vverjans

Message:

CHG: making stochasticforcing.m more generic, extending stochasticforcing capability to DefaultCalving

File:

• issm/trunk-jpl/src/m/classes/stochasticforcing.m (modified) (6 diffs)

issm/trunk-jpl/src/m/classes/stochasticforcing.m

@@ -8 +8 @@
                 isstochasticforcing  = 0;
                 fields               = NaN;
-                dimensions           = NaN;
+                defaultdimension     = 0;
+      default_id           = NaN;
                 covariance           = NaN;
                 randomflag           = 1;
@@ -34 +35 @@
 
                         num_fields = numel(self.fields);
-                        size_tot   = sum(self.dimensions);
+                        
+                        %Check that covariance matrix is positive definite
+                        try
+                                chol(self.covariance);
+                        catch
+                                error('md.stochasticforcing.covariance is not positive definite');
+                        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
+                        for field=self.fields
+            %Checking agreement of classes
+            if(contains(field,'SMB'))
+               if~(isequal(class(md.smb),char(field)))
+                  error('md.smb does not agree with stochasticforcing field %s', char(field));
+               end
+            end
+            if(contains(field,'frontalforcings'))
+               if~(isequal(class(md.frontalforcings),char(field)))
+                  error('md.frontalforcings 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
+                        if any(contains(self.fields,'SMBautoregression'))
+            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
+
+                        if(indSMBar~=-1 && indTFar~=-1) %both autoregressive models are used: check autoregressive time step consistency
+                                if((md.smb.ar_timestep~=md.frontalforcings.ar_timestep) && any(self.covariance(1+sum(self.dimensions(1:indSMBar-1)):sum(self.dimensions(1:indSMBar)),1+sum(self.dimensions(1:indTFar-1)):sum(self.dimensions(1:indTFar)))))
+                                        error('SMBautoregression and FrontalForcingsRignotAutoregression have different ar_timestep and non-zero covariance');
+                                end
+                        end
 
                         md = checkfield(md,'fieldname','stochasticforcing.isstochasticforcing','values',[0 1]);
                         md = checkfield(md,'fieldname','stochasticforcing.fields','numel',num_fields,'cell',1,'values',supportedstochforcings()); %VV check here 'cell' (19Oct2021)
-                        md = checkfield(md,'fieldname','stochasticforcing.dimensions','NaN',1,'Inf',1,'>',0,'size',[1,num_fields]); %specific dimension for each field
                         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]);
-
-                        %Check that all fields agree with the corresponding md class
-                        for field=self.fields
-                                if(contains(field,'SMB'))
-                                        if~(isequal(class(md.smb),char(field)))
-                                                error('md.smb does not agree with stochasticforcing field %s', char(field));
-                                        end
-                                end
-                                if(contains(field,'frontalforcings'))
-                                        if~(isequal(class(md.frontalforcings),char(field)))
-                                                error('md.frontalforcings does not agree with stochasticforcing field %s', char(field));
-                                        end
-                                end
-                        end
+                        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]);
+         end
                 end % }}}
                 function disp(self) % {{{
@@ -60 +95 @@
                         fielddisplay(self,'isstochasticforcing','is stochasticity activated?');
                         fielddisplay(self,'fields','fields with stochasticity applied, ex: {''SMBautoregression''}, or {''FrontalForcingsRignotAutoregression''}');
-                        fielddisplay(self,'dimensions','dimensionality of each field');
+                        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,'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)');
@@ -77 +113 @@
                                 return
                         else
+
+                                %Retrieve dimensionality of each field
+                                dimensions = self.defaultdimension*ones(1,num_fields);
+                                ind = 1;
+                                for field=self.fields
+                                        %Checking for specific dimensions
+                                        if(strcmp(field,'SMBautoregression'))
+                                                dimensions(ind) = md.smb.num_basins;
+                                        end
+                                        if(strcmp(field,'FrontalForcingsRignotAutoregression'))
+                                                dimensions(ind) = md.frontalforcings.num_basins;
+                                        end
+                                        ind = ind+1;
+                                end
+
                                 %Scaling covariance matrix (scale column-by-column and row-by-row)
-                                scaledfields = {'SMBautoregression'}; %list of fields that need scaling *1/yts
+                                scaledfields = {'DefaultCalving','SMBautoregression'}; %list of fields that need scaling *1/yts
                                 tempcovariance = self.covariance; %copy of covariance to avoid writing back in member variable
                                 for i=1:num_fields
                                         if any(strcmp(scaledfields,self.fields(i)))
-                                                inds = [1+sum(self.dimensions(1:i-1)):1:sum(self.dimensions(1:i))];
+                                                inds = [1+sum(dimensions(1:i-1)):1:sum(dimensions(1:i))];
                                                 for row=inds %scale rows corresponding to scaled field
                                                         tempcovariance(row,:) = 1./yts*tempcovariance(row,:);
@@ -93 +144 @@
                                 WriteData(fid,prefix,'data',num_fields,'name','md.stochasticforcing.num_fields','format','Integer');
                                 WriteData(fid,prefix,'object',self,'fieldname','fields','format','StringArray');
-                                WriteData(fid,prefix,'object',self,'fieldname','dimensions','format','IntMat');
+                                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','defaultdimension','format','Integer');
                                 WriteData(fid,prefix,'data',tempcovariance,'name','md.stochasticforcing.covariance','format','DoubleMat');
                                 WriteData(fid,prefix,'object',self,'fieldname','randomflag','format','Boolean');
@@ -105 +158 @@
 
    list = {...
-      'SMBautoregression',...
-      'FrontalForcingsRignotAutoregression'
+      'DefaultCalving',...
+                'FrontalForcingsRignotAutoregression',...
+      'SMBautoregression'
       };
 end % }}}