Changeset 26837


Ignore:
Timestamp:
01/29/22 12:15:51 (3 years ago)
Author:
vverjans
Message:

NEW: added stochasticforcing for autoregressionlinearbasalforcings

Location:
issm/trunk-jpl
Files:
2 added
8 edited

Legend:

Unmodified
Added
Removed

  • issm/trunk-jpl/src/c/analyses/MasstransportAnalysis.cpp

    r26836 r26837  
    221221                case AutoregressionLinearFloatingMeltRateEnum:
    222222                        iomodel->FetchDataToInput(inputs,elements,"md.basalforcings.basin_id",BasalforcingsLinearBasinIdEnum);
     223                        if(isstochastic) iomodel->FetchDataToInput(inputs,elements,"md.stochasticforcing.default_id",StochasticForcingDefaultIdEnum);
    223224                        break;
    224225                default:

  • issm/trunk-jpl/src/c/modules/StochasticForcingx/StochasticForcingx.cpp

    r26825 r26837  
    1111void StochasticForcingx(FemModel* femmodel){/*{{{*/
    1212
    13 
    14         //VV testing (12Nov)
    15         /*
    16         IssmDouble timeVV,dtVV,starttimeVV;
    17         femmodel->parameters->FindParam(&timeVV,TimeEnum);
    18    femmodel->parameters->FindParam(&dtVV,TimesteppingTimeStepEnum);
    19    femmodel->parameters->FindParam(&starttimeVV,TimesteppingStartTimeEnum);
    20         IssmDouble valMean = 0;
    21         IssmDouble valSdev = 0.5;
    22         int seed;
    23         //seed = reCast<int,IssmDouble>((timeVV-starttimeVV)/dtVV);
    24         seed = -1;
    25         IssmDouble rdmVV;
    26         univariateNormal_test0(&rdmVV,valMean,valSdev,seed);
    27         _printf_("VV rdmVV: "<<rdmVV<<'\n');
    28         */
    2913
    3014   /*Retrieve parameters*/
    … …  
    7357
    7458                /*Deal with the autoregressive models*/
    75                 if(fields[j]==SMBautoregressionEnum || fields[j]==FrontalForcingsRignotAutoregressionEnum){
     59                if(fields[j]==SMBautoregressionEnum || fields[j]==FrontalForcingsRignotAutoregressionEnum || fields[j]==BasalforcingsDeepwaterMeltingRateAutoregressionEnum){
    7660                        switch(fields[j]){
    7761                                case SMBautoregressionEnum:
    … …  
    8266                                        dimenum_type   = FrontalForcingsBasinIdEnum;
    8367                                        noiseenum_type = ThermalforcingAutoregressionNoiseEnum;
     68                                        break;
     69                                case BasalforcingsDeepwaterMeltingRateAutoregressionEnum:
     70                                        dimenum_type   = BasalforcingsLinearBasinIdEnum;
     71                                        noiseenum_type = BasalforcingsDeepwaterMeltingRateNoiseEnum;
    8472                                        break;
    8573                        }
    … …  
    9886                                case SMBautoregressionEnum:
    9987                                case FrontalForcingsRignotAutoregressionEnum:
     88                                case BasalforcingsDeepwaterMeltingRateAutoregressionEnum:
    10089                                        /*Already done above*/
    10190                                        break;
    … …  
    185174                                        }
    186175                                        break;
    187                                 //VV working(29Nov2021)
    188176                                case FrictionWaterPressureEnum:
    189177                                        /*Specify that WaterPressure is stochastic*/
    … …  
    208196                                        }
    209197                                        break;
    210 
    211                                 /*
    212                                 case FrictionEffectivePressureEnum:
    213                                         femmodel->parameters->SetParam(true,StochasticForcingIsEffectivePressureEnum);
    214                                         for(Object* &object:femmodel->elements->objects){
    215                   Element* element = xDynamicCast<Element*>(object);
    216                   int numvertices  = element->GetNumberOfVertices();
    217                   IssmDouble Neff[numvertices];
    218                                                 element->GetInputValue(&dimensionid,StochasticForcingDefaultIdEnum);
    219                                                 Gauss* gauss=element->NewGauss();
    220                                                 Friction* friction = new Friction(element);
    221                                                 for(int i=0;i<numvertices;i++){
    222                                                         gauss->GaussVertex(i);
    223                                                         Neff[i] = friction->EffectivePressure(gauss);
    224                                                         Neff[i] = Neff[i]+noisefield[dimensionid];
    225                                                 }
    226                                                 element->AddInput(FrictionEffectivePressureEnum,Neff,P1DGEnum);
    227                                                 delete gauss;
    228                                                 delete friction;
    229                                         }
    230                                         break;
    231                                 */
    232 
    233198                                default:
    234199                                        _error_("Field "<<EnumToStringx(fields[j])<<" does not support stochasticity yet.");

  • issm/trunk-jpl/src/c/shared/Enum/Enum.vim

    r26836 r26837  
    7373syn keyword cConstant BasalforcingsDeepwaterElevationEnum
    7474syn keyword cConstant BasalforcingsDeepwaterMeltingRateEnum
    75 syn keyword cConstant BasalforcingsDeepwaterMeltingRateNoiseEnum
    7675syn keyword cConstant BasalforcingsDtbgEnum
    7776syn keyword cConstant BasalforcingsEnum
    … …  
    603602syn keyword cConstant BasalCrevasseEnum
    604603syn keyword cConstant BasalforcingsDeepwaterMeltingRateAutoregressionEnum
     604syn keyword cConstant BasalforcingsDeepwaterMeltingRateNoiseEnum
    605605syn keyword cConstant BasalforcingsDeepwaterMeltingRateValuesAutoregressionEnum
    606606syn keyword cConstant BasalforcingsFloatingiceMeltingRateEnum

  • issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h

    r26836 r26837  
    6767        BasalforcingsDeepwaterElevationEnum,
    6868        BasalforcingsDeepwaterMeltingRateEnum,
    69         BasalforcingsDeepwaterMeltingRateNoiseEnum,
    7069        BasalforcingsDtbgEnum,
    7170        BasalforcingsEnum,
    … …  
    599598        BasalCrevasseEnum,
    600599        BasalforcingsDeepwaterMeltingRateAutoregressionEnum,
     600        BasalforcingsDeepwaterMeltingRateNoiseEnum,
    601601        BasalforcingsDeepwaterMeltingRateValuesAutoregressionEnum,
    602602        BasalforcingsFloatingiceMeltingRateEnum,

  • issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp

    r26836 r26837  
    7575                case BasalforcingsDeepwaterElevationEnum : return "BasalforcingsDeepwaterElevation";
    7676                case BasalforcingsDeepwaterMeltingRateEnum : return "BasalforcingsDeepwaterMeltingRate";
    77                 case BasalforcingsDeepwaterMeltingRateNoiseEnum : return "BasalforcingsDeepwaterMeltingRateNoise";
    7877                case BasalforcingsDtbgEnum : return "BasalforcingsDtbg";
    7978                case BasalforcingsEnum : return "Basalforcings";
    … …  
    605604                case BasalCrevasseEnum : return "BasalCrevasse";
    606605                case BasalforcingsDeepwaterMeltingRateAutoregressionEnum : return "BasalforcingsDeepwaterMeltingRateAutoregression";
     606                case BasalforcingsDeepwaterMeltingRateNoiseEnum : return "BasalforcingsDeepwaterMeltingRateNoise";
    607607                case BasalforcingsDeepwaterMeltingRateValuesAutoregressionEnum : return "BasalforcingsDeepwaterMeltingRateValuesAutoregression";
    608608                case BasalforcingsFloatingiceMeltingRateEnum : return "BasalforcingsFloatingiceMeltingRate";

  • issm/trunk-jpl/src/c/shared/Enum/Enumjl.vim

    r26836 r26837  
    6666syn keyword juliaConstC BasalforcingsDeepwaterElevationEnum
    6767syn keyword juliaConstC BasalforcingsDeepwaterMeltingRateEnum
    68 syn keyword juliaConstC BasalforcingsDeepwaterMeltingRateNoiseEnum
    6968syn keyword juliaConstC BasalforcingsDtbgEnum
    7069syn keyword juliaConstC BasalforcingsEnum
    … …  
    596595syn keyword juliaConstC BasalCrevasseEnum
    597596syn keyword juliaConstC BasalforcingsDeepwaterMeltingRateAutoregressionEnum
     597syn keyword juliaConstC BasalforcingsDeepwaterMeltingRateNoiseEnum
    598598syn keyword juliaConstC BasalforcingsDeepwaterMeltingRateValuesAutoregressionEnum
    599599syn keyword juliaConstC BasalforcingsFloatingiceMeltingRateEnum

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

    r26825 r26837  
    8484                                        end
    8585                                end
     86                                if(contains(field,'BasalforcingsDeepwaterMeltingRateAutoregression'))
     87                                        mdname = structstoch.mdnames(find(strcmp(structstoch.fields,char(field))));
     88                                        if~(isequal(class(md.basalforcings),char(mdname)))
     89                                                error('md.basalforcings does not agree with stochasticforcing field %s', char(field));
     90                                        end
     91                                end
    8692                                if(contains(field,'WaterPressure'))
    8793               mdname = structstoch.mdnames(find(strcmp(structstoch.fields,char(field))));
    … …  
    97103            end
    98104                                %Checking for specific dimensions
    99                                 if ~(strcmp(field,'SMBautoregression') || strcmp(field,'FrontalForcingsRignotAutoregression'))
     105                                if ~(strcmp(field,'SMBautoregression') || strcmp(field,'FrontalForcingsRignotAutoregression') || strcmp(field,'BasalforcingsDeepwaterMeltingRateAutoregression'))
    100106                                        checkdefaults = true; %field with non-specific dimensionality
    101107                                end
    102108                        end
    103109
    104                         %Retrieve sum of all the field dimensionalities
    105                         size_tot        = self.defaultdimension*num_fields;
    106                         indSMBar        = -1; %about to check for index of SMBautoregression
    107                         indTFar = -1; %about to check for index of FrontalForcingsRignotAutoregression
     110                        %Retrieve all the field dimensionalities
     111                        dimensions = self.defaultdimension*ones(1,num_fields);
     112                        indSMBar   = -1; %about to check for index of SMBautoregression
     113                        indTFar   = -1; %about to check for index of FrontalForcingsRignotAutoregression
     114                        indBDWar          = -1; %about to check for index of BasalforcingsDeepwaterMeltingRateAutoregression
    108115                        if any(contains(self.fields,'SMBautoregression'))
    109                                 size_tot = size_tot-self.defaultdimension+md.smb.num_basins;
    110                                 indSMBar = find(contains(self.fields,'SMBautoregression')); %index of SMBar, now check for consistency with TFar timestep (08Nov2021)
     116                                indSMBar = find(contains(self.fields,'SMBautoregression')); %index of SMBar, now check for consistency with other ar timesteps
     117                                dimensions(indSMBar) = md.smb.num_basins;
    111118                        end
    112119                        if any(contains(self.fields,'FrontalForcingsRignotAutoregression'))
    113                                 size_tot        = size_tot-self.defaultdimension+md.frontalforcings.num_basins;
    114                                 indTFar = find(contains(self.fields,'FrontalForcingsRignotAutoregression')); %index of TFar, now check for consistency with SMBar timestep (08Nov2021)
    115                         end
     120                                indTFar = find(contains(self.fields,'FrontalForcingsRignotAutoregression')); %index of TFar, now check for consistency with other ar timesteps
     121                                dimensions(indTFar) = md.frontalforcings.num_basins;
     122                        end
     123                        if any(contains(self.fields,'BasalforcingsDeepwaterMeltingRateAutoregression'))
     124                                indBDWar        = find(contains(self.fields,'BasalforcingsDeepwaterMeltingRateAutoregression')); %index of BDWar, now check for consistency with other ar timesteps
     125                                dimensions(indBDWar) = md.basalforcings.num_basins;
     126                        end
     127                        size_tot = sum(dimensions);
    116128
    117129                        if(indSMBar~=-1 && indTFar~=-1) %both autoregressive models are used: check autoregressive time step consistency
    118                                 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))))~=0)
     130                                if((md.smb.ar_timestep~=md.frontalforcings.ar_timestep) && any(self.covariance(1+sum(dimensions(1:indSMBar-1)):sum(dimensions(1:indSMBar)),1+sum(dimensions(1:indTFar-1)):sum(dimensions(1:indTFar))),'all')~=0)
    119131                                        error('SMBautoregression and FrontalForcingsRignotAutoregression have different ar_timestep and non-zero covariance');
    120132                                end
    121133                        end
    122 
     134                        if(indSMBar~=-1 && indBDWar~=-1) %both autoregressive models are used: check autoregressive time step consistency
     135                                if((md.smb.ar_timestep~=md.basalforcings.ar_timestep) && any(self.covariance(1+sum(dimensions(1:indSMBar-1)):sum(dimensions(1:indSMBar)),1+sum(dimensions(1:indBDWar-1)):sum(dimensions(1:indBDWar))),'all')~=0)
     136                                        error('SMBautoregression and BasalforcingsDeepwaterMeltingRateAutoregression have different ar_timestep and non-zero covariance');
     137                                end
     138                        end
     139                        if(indTFar~=-1 && indBDWar~=-1) %both autoregressive models are used: check autoregressive time step consistency
     140                                if((md.frontalforcings.ar_timestep~=md.basalforcings.ar_timestep) && any(self.covariance(1+sum(dimensions(1:indTFar-1)):sum(dimensions(1:indTFar)),1+sum(dimensions(1:indBDWar-1)):sum(dimensions(1:indBDWar))),'all')~=0)
     141                                        error('FrontalForcingsRignotAutoregression and BasalforcingsDeepwaterMeltingRateAutoregression have different ar_timestep and non-zero covariance');
     142                                end
     143                        end
    123144                        md = checkfield(md,'fieldname','stochasticforcing.isstochasticforcing','values',[0 1]);
    124145                        md = checkfield(md,'fieldname','stochasticforcing.fields','numel',num_fields,'cell',1,'values',supportedstochforcings());
    … …  
    164185                                                dimensions(ind) = md.frontalforcings.num_basins;
    165186                                        end
     187                                        if(strcmp(field,'BasalforcingsDeepwaterMeltingRateAutoregression'))
     188                                                dimensions(ind) = md.basalforcings.num_basins;
     189                                        end
    166190                                        ind = ind+1;
    167191                                end
    168192
    169193                                %Scaling covariance matrix (scale column-by-column and row-by-row)
    170                                 scaledfields = {'BasalforcingsSpatialDeepwaterMeltingRate','DefaultCalving','FloatingMeltRate','SMBautoregression','SMBforcing'}; %list of fields that need scaling *1/yts
     194                                scaledfields = {'BasalforcingsDeepwaterMeltingRateAutoregression','BasalforcingsSpatialDeepwaterMeltingRate','DefaultCalving','FloatingMeltRate','SMBautoregression','SMBforcing'}; %list of fields that need scaling *1/yts
    171195                                tempcovariance = self.covariance; %copy of covariance to avoid writing back in member variable
    172196                                for i=1:num_fields
    … …  
    185209                                        self.default_id = zeros(md.mesh.numberofelements,1);
    186210                                end
     211
    187212                                WriteData(fid,prefix,'data',num_fields,'name','md.stochasticforcing.num_fields','format','Integer');
    188213                                WriteData(fid,prefix,'object',self,'fieldname','fields','format','StringArray');
    … …  
    207232        % supported and corresponding md names
    208233        structure.fields = {...
     234                'BasalforcingsDeepwaterMeltingRateAutoregression',...
    209235                'BasalforcingsSpatialDeepwaterMeltingRate',...
    210236                'DefaultCalving',...
    … …  
    216242                };
    217243        structure.mdnames = {...
     244                'autoregressionlinearbasalforcings',...
    218245                'spatiallinearbasalforcings',...
    219246                'calving',...

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

    r26825 r26837  
    9090                if (type(md.basalforcings).__name__ != mdname):
    9191                    raise TypeError('md.basalforcings does not agree with stochasticforcing field {}'.format(field))
    92             if 'BasalforcingsDpatialDeepwaterMeltingRate' in field:
     92            if 'BasalforcingsSpatialDeepwaterMeltingRate' in field:
     93                mdname = structstoch[field]
     94                if (type(md.basalforcings).__name__ != mdname):
     95                    raise TypeError('md.basalforcings does not agree with stochasticforcing field {}'.format(field))
     96            if 'BasalforcingsDeepwaterMeltingRateAutoregression' in field:
    9397                mdname = structstoch[field]
    9498                if (type(md.basalforcings).__name__ != mdname):
    … …  
    104108
    105109            # Checking for specific dimensions
    106             if field not in['SMBautoregression', 'FrontalForcingsRignotAutoregression']:
     110            if field not in['SMBautoregression', 'FrontalForcingsRignotAutoregression','BasalforcingsDeepwaterMeltingRateAutoregression']:
    107111                checkdefaults = True  # field with non-specific dimensionality
    108112
    109113        # Retrieve sum of all the field dimensionalities
    110         size_tot = self.defaultdimension * num_fields
    111         indSMBar = -1  # About to check for index of SMBautoregression
    112         indTFar = -1  # About to check for index of FrontalForcingsRignotAutoregression
     114        dimensions = self.defaultdimension*np.ones((num_fields))
     115        indSMBar   = -1  # About to check for index of SMBautoregression
     116        indTFar    = -1  # About to check for index of FrontalForcingsRignotAutoregression
     117        indBDWar   = -1  # About to check for index of BasalforcingsDeepwaterMeltingRateAutoregression
    113118        if ('SMBautoregression' in self.fields):
    114             size_tot = size_tot - self.defaultdimension + md.smb.num_basins
    115             indSMBar = self.fields.index('SMBautoregression')  # Index of SMBar, now check for consistency with TFar timestep (08Nov2021)
     119            indSMBar = self.fields.index('SMBautoregression')  # Index of SMBar, now check for consistency with other timesteps
     120            dimensions[indSMBar] = md.smb.num_basins
    116121        if ('FrontalForcingsRignotAutoregression' in self.fields):
    117             size_tot = size_tot - self.defaultdimension + md.frontalforcings.num_basins
    118             indTFar = self.fields.index('FrontalForcingsRignotAutoregression')  # Index of TFar, now check for consistency with SMBar timestep (08Nov2021)
     122            indTFar = self.fields.index('FrontalForcingsRignotAutoregression')  # Index of TFar, now check for consistency with other timesteps
     123            dimensions[indTFar] = md.frontalforcings.num_basins
     124        if ('BasalforcingsDeepwaterMeltingRateAutoregression' in self.fields):
     125            indBDWar = self.fields.index('BasalforcingsDeepwaterMeltingRateAutoregression')  # Index of BDWar, now check for consistency with other timesteps
     126            dimensions[indTFar] = md.basalforcings.num_basins
     127        size_tot = np.sum(dimensions)
     128
    119129        if (indSMBar != -1 and indTFar != -1):  # Both autoregressive models are used: check autoregressive time step consistency
    120             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)]
     130            covsum = self.covariance[np.sum(dimensions[0:indSMBar]).astype(int):np.sum(dimensions[0:indSMBar + 1]).astype(int), np.sum(dimensions[0:indTFar]).astype(int):np.sum(dimensions[0:indTFar + 1]).astype(int)]
    121131            if((md.smb.ar_timestep != md.frontalforcings.ar_timestep) and np.any(covsum != 0)):
    122132                raise IOError('SMBautoregression and FrontalForcingsRignotAutoregression have different ar_timestep and non-zero covariance')
     133        if (indSMBar != -1 and indBDWar != -1):  # Both autoregressive models are used: check autoregressive time step consistency
     134            covsum = self.covariance[np.sum(dimensions[0:indSMBar]).astype(int):np.sum(dimensions[0:indSMBar + 1]).astype(int), np.sum(dimensions[0:indBDWar]).astype(int):np.sum(dimensions[0:indBDWar + 1]).astype(int)]
     135            if((md.smb.ar_timestep != md.basalforcings.ar_timestep) and np.any(covsum != 0)):
     136                raise IOError('SMBautoregression and BasalforcingsDeepwaterMeltingRateAutoregression have different ar_timestep and non-zero covariance')
     137        if (indTFar != -1 and indBDWar != -1):  # Both autoregressive models are used: check autoregressive time step consistency
     138            covsum = self.covariance[np.sum(dimensions[0:indTFar]).astype(int):np.sum(dimensions[0:indTFar + 1]).astype(int), np.sum(dimensions[0:indBDWar]).astype(int):np.sum(dimensions[0:indBDWar + 1]).astype(int)]
     139            if((md.frontalforcings.ar_timestep != md.basalforcings.ar_timestep) and np.any(covsum != 0)):
     140                raise IOError('FrontalForcingsRignotAutoregression and BasalforcingsDeepwaterMeltingRateAutoregression have different ar_timestep and non-zero covariance')
    123141
    124142        md = checkfield(md, 'fieldname', 'stochasticforcing.isstochasticforcing', 'values', [0, 1])
    … …  
    153171                if (field == 'FrontalForcingsRignotAutoregression'):
    154172                    dimensions[ind] = md.frontalforcings.num_basins
     173                if (field == 'BasalforcingsDeepwaterMeltingRateAutoregression'):
     174                    dimensions[ind] = md.basalforcings.num_basins
    155175
    156176            # Scaling covariance matrix (scale column-by-column and row-by-row)
    157             scaledfields = ['BasalforcingsSpatialDeepwaterMeltingRate','DefaultCalving', 'FloatingMeltRate', 'SMBautoregression', 'SMBforcing']  # list of fields that need scaling * 1/yts
     177            scaledfields = ['BasalforcingsDeepwaterMeltingRateAutoregression','BasalforcingsSpatialDeepwaterMeltingRate','DefaultCalving', 'FloatingMeltRate', 'SMBautoregression', 'SMBforcing']  # list of fields that need scaling * 1/yts
    158178            tempcovariance = np.copy(self.covariance)
    159179            for i in range(num_fields):
    … …  
    191211           supported and corresponding md names
    192212        """
    193         structure = {'BasalforcingsSpatialDeepwaterMeltingRate': 'spatiallinearbasalforcings',
     213        structure = {'BasalforcingsDeepwaterMeltingRateAutoregression': 'autoregressionlinearbasalforcings',
     214                     'BasalforcingsSpatialDeepwaterMeltingRate': 'spatiallinearbasalforcings',
    194215                     'DefaultCalving': 'calving',
    195216                     'FloatingMeltRate': 'basalforcings',
Note: See TracChangeset for help on using the changeset viewer.