Changeset 26905

Timestamp:

02/25/22 00:19:38 (3 years ago)

Author:

bdef

Message:

CHG: uniformisation of id indexing with matlab, also checking shapes of array that have to be columns and transposing if needed

Location:

issm/trunk-jpl/src/m/classes

Files:

• SMBautoregression.py (modified) (4 diffs)
• autoregressionlinearbasalforcings.py (modified) (2 diffs)
• frontalforcingsrignot.py (modified) (1 diff)
• frontalforcingsrignotautoregression.py (modified) (4 diffs)
• stochasticforcing.py (modified) (1 diff)

issm/trunk-jpl/src/m/classes/SMBautoregression.py

@@ -95 +95 @@
             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])
+            if len(np.shape(self.beta0)) == 1:
+                self.beta0 = np.array([self.beta0])
+                self.beta1 = np.array([self.beta1])
             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; NOTE: As opposed to MATLAB implementation, pass list
@@ -101 +104 @@
             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])
-            if(np.any(np.isnan(md.smb.lapserate_pos)==False) or np.size(md.smb.lapserate_pos)>1):
+
+            if(np.any(np.isnan(self.lapserate_pos) is False) or np.size(self.lapserate_pos) > 1):
+                if len(np.shape(self.lapserate_pos)) == 1:
+                    self.lapserate_pos = np.array([self.lapserate_pos])
                 md = checkfield(md, 'fieldname', 'smb.lapserate_pos', 'NaN', 1, 'Inf', 1, 'size', [1, md.smb.num_basins], 'numel', md.smb.num_basins)
-            if(np.any(np.isnan(md.smb.lapserate_neg)==False) or np.size(md.smb.lapserate_neg)>1):
+
+            if(np.any(np.isnan(self.lapserate_neg) is False) or np.size(self.lapserate_neg) > 1):
+                if len(np.shape(self.lapserate_neg)) == 1:
+                    self.lapserate_neg = np.array([self.lapserate_neg])
                 md = checkfield(md, 'fieldname', 'smb.lapserate_neg', 'NaN', 1, 'Inf', 1, 'size', [1, md.smb.num_basins], 'numel', md.smb.num_basins)
-            if(np.any(np.isnan(md.smb.refelevation)==False) or np.size(md.smb.refelevation)>1):
+
+            if(np.any(np.isnan(self.refelevation) is False) or np.size(self.refelevation) > 1):
+                if len(np.shape(self.refelevation)) == 1:
+                    self.refelevation = np.array([self.refelevation])
                 md = checkfield(md, 'fieldname', 'smb.refelevation', 'NaN', 1, 'Inf', 1, '>=', 0, 'size', [1, md.smb.num_basins], 'numel', md.smb.num_basins)
 
@@ -129 +141 @@
             temprefelevation = np.zeros((md.smb.num_basins)).reshape(1,md.smb.num_basins)
             areas = GetAreas(md.mesh.elements, md.mesh.x, md.mesh.y)
-            for ii,bid in enumerate(np.unique(md.smb.basin_id)):
-                indices = np.where(md.smb.basin_id==bid)[0]
-                elemsh  = np.zeros((len(indices)))
+            for ii, bid in enumerate(np.unique(md.smb.basin_id)):
+                indices = np.where(md.smb.basin_id == bid)[0]
+                elemsh = np.zeros((len(indices)))
                 for jj in range(len(indices)):
-                    elemsh[jj] = np.mean(md.geometry.surface[md.mesh.elements[indices[jj],:]-1])
-                temprefelevation[0,ii] = np.sum(areas[indices]*elemsh)/np.sum(areas[indices])
-            if(np.any(templapserate_pos!=0) or np.any(templapserate_neg!=0)):
+                    elemsh[jj] = np.mean(md.geometry.surface[md.mesh.elements[indices[jj], :] - 1])
+                temprefelevation[0, ii] = np.sum(areas[indices] * elemsh) / np.sum(areas[indices])
+            if(np.any(templapserate_pos != 0) or np.any(templapserate_neg != 0)):
                 print('      smb.refelevation not specified: Reference elevations set to mean surface elevation of basins')
 
@@ -143 +155 @@
         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', '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', 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, 'data', templapserate_pos, 'name', 'md.smb.lapserate_pos', 'format', 'DoubleMat','scale',1/yts,'yts',yts)
-        WriteData(fid, prefix, 'data', templapserate_neg, 'name', 'md.smb.lapserate_neg', 'format', 'DoubleMat','scale',1/yts,'yts',yts)
+        WriteData(fid, prefix, 'data', templapserate_pos, 'name', 'md.smb.lapserate_pos', 'format', 'DoubleMat', 'scale', 1 / yts, 'yts', yts)
+        WriteData(fid, prefix, 'data', templapserate_neg, 'name', 'md.smb.lapserate_neg', 'format', 'DoubleMat', 'scale', 1 / yts, 'yts', yts)
         WriteData(fid, prefix, 'data', temprefelevation, 'name', 'md.smb.refelevation', 'format', 'DoubleMat')
         WriteData(fid, prefix, 'object', self, 'fieldname', 'steps_per_step', 'format', 'Integer')

issm/trunk-jpl/src/m/classes/autoregressionlinearbasalforcings.py

@@ -73 +73 @@
         if 'MasstransportAnalysis' in analyses:
             md = checkfield(md, 'fieldname', 'basalforcings.num_basins', 'numel', 1, 'NaN', 1, 'Inf', 1, '>', 0)
-            md = checkfield(md,'fieldname','basalforcings.groundedice_melting_rate','NaN',1,'Inf',1,'timeseries',1)
-            md = checkfield(md, 'fieldname', 'basalforcings.deepwater_elevation', 'NaN', 1, 'Inf', 1, 'size', [1, md.basalforcings.num_basins], 'numel', md.basalforcings.num_basins) 
-            md = checkfield(md, 'fieldname', 'basalforcings.upperwater_elevation', 'NaN', 1, 'Inf', 1,'<=',0, 'size', [1, md.basalforcings.num_basins], 'numel', md.basalforcings.num_basins) 
-            md = checkfield(md, 'fieldname', 'basalforcings.upperwater_melting_rate', 'NaN', 1, 'Inf', 1,'>=',0, 'size', [1, md.basalforcings.num_basins], 'numel', md.basalforcings.num_basins) 
+            md = checkfield(md, 'fieldname', 'basalforcings.groundedice_melting_rate', 'NaN', 1, 'Inf', 1, 'timeseries', 1)
+
+            if len(np.shape(self.deepwater_elevation)) == 1:
+                self.deepwater_elevation = np.array([self.deepwater_elevation])
+                self.upperwater_elevation = np.array([self.upperwater_elevation])
+                self.upperwater_melting_rate = np.array([self.upperwater_melting_rate])
+            md = checkfield(md, 'fieldname', 'basalforcings.deepwater_elevation', 'NaN', 1, 'Inf', 1, 'size', [1, md.basalforcings.num_basins], 'numel', md.basalforcings.num_basins)
+            md = checkfield(md, 'fieldname', 'basalforcings.upperwater_elevation', 'NaN', 1, 'Inf', 1, '<=', 0, 'size', [1, md.basalforcings.num_basins], 'numel', md.basalforcings.num_basins)
+            md = checkfield(md, 'fieldname', 'basalforcings.upperwater_melting_rate', 'NaN', 1, 'Inf', 1,'>=', 0, 'size', [1, md.basalforcings.num_basins], 'numel', md.basalforcings.num_basins)
             md = checkfield(md, 'fieldname', 'basalforcings.basin_id', 'Inf', 1, '>=', 0, '<=', md.basalforcings.num_basins, 'size', [md.mesh.numberofelements])
+
+            if len(np.shape(self.beta0)) == 1:
+                self.beta0 = np.array([self.beta0])
+                self.beta1 = np.array([self.beta1])
+
             md = checkfield(md, 'fieldname', 'basalforcings.beta0', 'NaN', 1, 'Inf', 1, 'size', [1, md.basalforcings.num_basins], 'numel', md.basalforcings.num_basins) # Scheme fails if passed as column vector
             md = checkfield(md, 'fieldname', 'basalforcings.beta1', 'NaN', 1, 'Inf', 1, 'size', [1, md.basalforcings.num_basins], 'numel', md.basalforcings.num_basins) # Scheme fails if passed as column vector; NOTE: As opposed to MATLAB implementation, pass list
@@ -102 +112 @@
         WriteData(fid, prefix, 'object', self, 'fieldname', 'ar_initialtime', 'format', 'Double', 'scale', yts)
         WriteData(fid, prefix, 'object', self, 'fieldname', 'ar_timestep', 'format', 'Double', 'scale', yts)
-        WriteData(fid, prefix, 'object', self, 'fieldname', 'basin_id', 'data', self.basin_id, 'name', 'md.basalforcings.basin_id', 'format', 'IntMat', 'mattype', 2) # 0-indexed
+        WriteData(fid, prefix, 'object', self, 'fieldname', 'basin_id', 'data', self.basin_id - 1, 'name', 'md.basalforcings.basin_id', 'format', 'IntMat', 'mattype', 2)  # 0-indexed
         WriteData(fid, prefix, 'object', self, 'fieldname', 'beta0', 'format', 'DoubleMat', 'name', 'md.basalforcings.beta0', 'scale', 1 / yts, 'yts', yts)
         WriteData(fid, prefix, 'object', self, 'fieldname', 'beta1', 'format', 'DoubleMat', 'name', 'md.basalforcings.beta1', 'scale', 1 / (yts ** 2), 'yts', yts)

issm/trunk-jpl/src/m/classes/frontalforcingsrignot.py

@@ -64 +64 @@
     def marshall(self, prefix, md, fid):  # {{{
         WriteData(fid, prefix, 'name', 'md.frontalforcings.parameterization', 'data', 2, 'format', 'Integer')
-        WriteData(fid, prefix, 'object', self, 'fieldname', 'basin_id', 'data', self.basin_id, 'name', 'md.frontalforcings.basin_id', 'format', 'IntMat', 'mattype', 2) # 0-indexed
+        WriteData(fid, prefix, 'object', self, 'fieldname', 'basin_id', 'data', self.basin_id - 1, 'name', 'md.frontalforcings.basin_id', 'format', 'IntMat', 'mattype', 2)  # 0-indexed
         WriteData(fid, prefix, 'object', self, 'fieldname', 'num_basins', 'format', 'Integer')
         WriteData(fid, prefix, 'object', self, 'fieldname', 'subglacial_discharge', 'format', 'DoubleMat', 'mattype', 1, 'timeserieslength', md.mesh.numberofvertices + 1, 'yts', md.constants.yts)

issm/trunk-jpl/src/m/classes/frontalforcingsrignotautoregression.py

@@ -1 +1 @@
 # -*- coding: utf-8 -*-
-
 import numpy as np
-
 from checkfield import checkfield
 from fielddisplay import fielddisplay
 from MatlabFuncs import *
 from WriteData import WriteData
+
 
 class frontalforcingsrignotautoregression(object):
@@ -49 +48 @@
         self.num_basins = 0
         self.subglacial_discharge = np.nan
-        self.ar_order = 0.0 # Autoregression model of order 0
+        self.ar_order = 0.0  # Autoregression model of order 0
         return self
     #}}}
@@ -61 +60 @@
         md = checkfield(md, 'fieldname', 'frontalforcings.basin_id', 'Inf', 1, '>=', 0, '<=', md.frontalforcings.num_basins, 'size', [md.mesh.numberofelements])
         md = checkfield(md, 'fieldname', 'frontalforcings.subglacial_discharge', '>=', 0, 'NaN', 1, 'Inf', 1, 'timeseries', 1)
+        if len(np.shape(self.beta0)) == 1:
+            self.beta0 = np.array([self.beta0])
+            self.beta1 = np.array([self.beta1])
         md = checkfield(md, 'fieldname', 'frontalforcings.beta0', 'NaN', 1, 'Inf', 1, 'size', [1, md.frontalforcings.num_basins], 'numel', md.frontalforcings.num_basins)
         md = checkfield(md, 'fieldname', 'frontalforcings.beta1', 'NaN', 1, 'Inf', 1, 'size', [1, md.frontalforcings.num_basins], 'numel', md.frontalforcings.num_basins)
@@ -83 +85 @@
         WriteData(fid, prefix, 'object', self, 'class', 'frontalforcings', 'fieldname', 'ar_initialtime', 'format', 'Double', 'scale', yts)
         WriteData(fid, prefix, 'object', self, 'class', 'frontalforcings', 'fieldname', 'ar_timestep', 'format', 'Double', 'scale', yts)
-        WriteData(fid, prefix, 'object', self, 'class', 'frontalforcings', 'fieldname', 'basin_id', 'data', self.basin_id, 'name', 'md.frontalforcings.basin_id', 'format', 'IntMat', 'mattype', 2) # 0-indexed
+        WriteData(fid, prefix, 'object', self, 'class', 'frontalforcings', 'fieldname', 'basin_id', 'data', self.basin_id - 1, 'name', 'md.frontalforcings.basin_id', 'format', 'IntMat', 'mattype', 2)  # 0-indexed
         WriteData(fid, prefix, 'object', self, 'class', 'frontalforcings', 'fieldname', 'beta0', 'format', 'DoubleMat', 'name', 'md.frontalforcings.beta0')
         WriteData(fid, prefix, 'object', self, 'class', 'frontalforcings', 'fieldname', 'beta1', 'format', 'DoubleMat', 'name', 'md.frontalforcings.beta1', 'scale', 1 / yts, 'yts', yts)

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

@@ -208 +208 @@
             WriteData(fid, prefix, 'object', self, 'fieldname', 'fields', 'format', 'StringArray')
             WriteData(fid, prefix, 'data', dimensions, 'name', 'md.stochasticforcing.dimensions', 'format', 'IntMat', 'mattype', 2)
-            WriteData(fid, prefix, 'object', self, 'fieldname', 'default_id', 'format', 'IntMat', 'mattype', 2)  #12Nov2021 make sure this is zero-indexed!
+            WriteData(fid, prefix, 'object', self, 'fieldname', 'default_id', 'data', self.default_id - 1, 'format', 'IntMat', 'mattype', 2)  #12Nov2021 make sure this is zero-indexed!
             WriteData(fid, prefix, 'object', self, 'fieldname', 'defaultdimension', 'format', 'Integer')
             WriteData(fid, prefix, 'data', tempcovariance, 'name', 'md.stochasticforcing.covariance', 'format', 'DoubleMat')