Index: /issm/trunk-jpl/src/m/classes/SMBautoregression.py =================================================================== --- /issm/trunk-jpl/src/m/classes/SMBautoregression.py (revision 26552) +++ /issm/trunk-jpl/src/m/classes/SMBautoregression.py (revision 26553) @@ -7,5 +7,5 @@ -class SMBcomponents(object): +class SMBautoregression(object): """SMBAUTOREGRESSION class definition @@ -22,7 +22,5 @@ self.ar_timestep = 0 self.phi = np.nan - self.covmat = np.nan self.basin_id = np.nan - self.randomflag = 1 self.steps_per_step = 1 self.averaging = 0 @@ -59,5 +57,4 @@ def setdefaultparameters(self): #{{{ self.ar_order = 0.0 # Autoregression model of order 0 - self.randomflag = 1 # }}} @@ -87,7 +84,4 @@ self.phi = np.zeros((self.num_basins, self.ar_order)) # Autoregression model of order 0 print(' smb.phi (lag coefficients) not specified: order of autoregressive model set to 0') - if np.all(np.isnan(self.covmat)): - self.covmat = 1e-21 * np.eye(self.num_basins) # No stochasticity and no covariance - print(' smb.covmat not specified: stochasticity set to 0') return self # }}} @@ -98,11 +92,9 @@ md = checkfield(md, 'fieldname', 'smb.basin_id', 'Inf', 1, '>=', 0, '<=', md.smb.num_basins, 'size', [md.mesh.numberofelements]) 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 + 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 md = checkfield(md, 'fieldname', 'smb.ar_order', 'numel', 1, 'NaN', 1, 'Inf', 1, '>=', 0) md = checkfield(md, 'fieldname', 'smb.ar_initialtime', 'numel', 1, 'NaN', 1, 'Inf', 1) 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]) - md = checkfield(md, 'fieldname', 'smb.covmat', 'NaN', 1, 'Inf', 1, 'size', [md.smb.num_basins, md.smb.num_basins]) - md = checkfield(md, 'fieldname', 'smb.randomflag', 'numel', [1], 'values', [0, 1]) md = checkfield(md, 'fieldname', 'smb.steps_per_step', '>=', 1, 'numel', [1]) md = checkfield(md, 'fieldname', 'smb.averaging', 'numel', [1], 'values', [0, 1, 2]) @@ -123,6 +115,4 @@ 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, 'object', self, 'class', 'smb', 'fieldname', 'covmat', 'format', 'DoubleMat', 'name', 'md.smb.covmat', 'scale', 1 / (yts ** 2), 'yts', yts) - WriteData(fid, prefix, 'object', self, 'class', 'smb', 'fieldname', 'randomflag', 'format', 'Boolean') WriteData(fid, prefix, 'object', self, 'fieldname', 'steps_per_step', 'format', 'Integer') WriteData(fid, prefix, 'object', self, 'fieldname', 'averaging', 'format', 'Integer') Index: /issm/trunk-jpl/src/m/classes/calvingvonmises.py =================================================================== --- /issm/trunk-jpl/src/m/classes/calvingvonmises.py (revision 26552) +++ /issm/trunk-jpl/src/m/classes/calvingvonmises.py (revision 26553) @@ -1,20 +1,18 @@ +from checkfield import checkfield from fielddisplay import fielddisplay -from checkfield import checkfield from WriteData import WriteData class calvingvonmises(object): - """ - CALVINGVONMISES class definition + """CALVINGVONMISES class definition - Usage: - calvingvonmises = calvingvonmises() + Usage: + calvingvonmises = calvingvonmises() """ def __init__(self): # {{{ - - self.stress_threshold_groundedice = 0. - self.stress_threshold_floatingice = 0. - self.min_thickness = 0. + self.stress_threshold_groundedice = 0 + self.stress_threshold_floatingice = 0 + self.min_thickness = 0 #set defaults @@ -24,10 +22,9 @@ def __repr__(self): # {{{ - string = ' Calving VonMises parameters:' - string = "%s\n%s" % (string, fielddisplay(self, 'stress_threshold_groundedice', 'sigma_max applied to grounded ice only [Pa]')) - string = "%s\n%s" % (string, fielddisplay(self, 'stress_threshold_floatingice', 'sigma_max applied to floating ice only [Pa]')) - string = "%s\n%s" % (string, fielddisplay(self, 'min_thickness', 'minimum thickness below which no ice is allowed [m]')) - - return string + s = ' Calving VonMises parameters:\n' + s += '{}\n'.format(fielddisplay(self, 'stress_threshold_groundedice', 'sigma_max applied to grounded ice only [Pa]')) + s += '{}\n'.format(fielddisplay(self, 'stress_threshold_floatingice', 'sigma_max applied to floating ice only [Pa]')) + s += '{}\n'.format(fielddisplay(self, 'min_thickness', 'minimum thickness below which no ice is allowed [m]')) + return s #}}} @@ -37,9 +34,9 @@ def setdefaultparameters(self): # {{{ - #Default sigma max + # Default sigma max self.stress_threshold_groundedice = 1e6 self.stress_threshold_floatingice = 150e3 - #turn off min_thickness by default. + # Turn off min_thickness by default self.min_thickness = 0. return self @@ -47,6 +44,6 @@ def checkconsistency(self, md, solution, analyses): # {{{ - #Early return - if solution == 'TransientSolution' or md.transient.ismovingfront == 0: + # Early return + if solution == 'TransientSolution' or not md.transient.ismovingfront: return Index: /issm/trunk-jpl/src/m/classes/frontalforcingsrignot.py =================================================================== --- /issm/trunk-jpl/src/m/classes/frontalforcingsrignot.py (revision 26552) +++ /issm/trunk-jpl/src/m/classes/frontalforcingsrignot.py (revision 26553) @@ -1,2 +1,4 @@ +# -*- coding: utf-8 -*- + import numpy as np @@ -7,5 +9,5 @@ class frontalforcingsrignot(object): - """FRONTAL FORCINGS Rignot class definition + """FRONTALFORCINGSRIGNOT class definition Usage: Index: /issm/trunk-jpl/src/m/classes/frontalforcingsrignotautoregression.m =================================================================== --- /issm/trunk-jpl/src/m/classes/frontalforcingsrignotautoregression.m (revision 26552) +++ /issm/trunk-jpl/src/m/classes/frontalforcingsrignotautoregression.m (revision 26553) @@ -7,11 +7,11 @@ properties (SetAccess=public) num_basins = 0; - beta0 = NaN; - beta1 = NaN; - ar_order = 0; - ar_initialtime = 0; - ar_timestep = 0; - phi = NaN; - basin_id = NaN; + beta0 = NaN; + beta1 = NaN; + ar_order = 0; + ar_initialtime = 0; + ar_timestep = 0; + phi = NaN; + basin_id = NaN; subglacial_discharge = NaN; end @@ -40,8 +40,8 @@ function self = setdefaultparameters(self) % {{{ - basin_id = NaN; - num_basins = 0; - subglacial_discharge = NaN; - self.ar_order = 0.0; %autoregression model of order 0 + self.basin_id = NaN; + self.num_basins = 0; + self.subglacial_discharge = NaN; + self.ar_order = 0.0; %autoregression model of order 0 end % }}} @@ -52,5 +52,5 @@ md = checkfield(md,'fieldname','frontalforcings.num_basins','numel',1,'NaN',1,'Inf',1,'>',0); md = checkfield(md,'fieldname','frontalforcings.basin_id','Inf',1,'>=',0,'<=',md.frontalforcings.num_basins,'size',[md.mesh.numberofelements 1]); - md = checkfield(md,'fieldname','frontalforcings.subglacial_discharge','>=',0,'nan',1,'Inf',1,'timeseries',1); + md = checkfield(md,'fieldname','frontalforcings.subglacial_discharge','>=',0,'NaN',1,'Inf',1,'timeseries',1); 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); @@ -77,5 +77,5 @@ WriteData(fid,prefix,'name','md.frontalforcings.parameterization','data',55,'format','Integer'); WriteData(fid,prefix,'object',self,'class','frontalforcings','fieldname','num_basins','format','Integer'); - WriteData(fid,prefix,'object',self,'class','frontalforcings','fieldname','subglacial_discharge','format','DoubleMat','mattype',1,'timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts); + WriteData(fid,prefix,'object',self,'class','frontalforcings','fieldname','subglacial_discharge','format','DoubleMat','mattype',1,'timeserieslength',md.mesh.numberofvertices+1,'yts',yts); WriteData(fid,prefix,'object',self,'class','frontalforcings','fieldname','ar_order','format','Integer'); WriteData(fid,prefix,'object',self,'class','frontalforcings','fieldname','ar_initialtime','format','Double','scale',yts); @@ -83,6 +83,6 @@ 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',md.constants.yts); - WriteData(fid,prefix,'object',self,'class','frontalforcings','fieldname','phi','format','DoubleMat','name','md.frontalforcings.phi','yts',md.constants.yts); + WriteData(fid,prefix,'object',self,'class','frontalforcings','fieldname','beta1','format','DoubleMat','name','md.frontalforcings.beta1','scale',1./yts,'yts',yts); + WriteData(fid,prefix,'object',self,'class','frontalforcings','fieldname','phi','format','DoubleMat','name','md.frontalforcings.phi','yts',yts); end % }}} end Index: /issm/trunk-jpl/src/m/classes/frontalforcingsrignotautoregression.py =================================================================== --- /issm/trunk-jpl/src/m/classes/frontalforcingsrignotautoregression.py (revision 26553) +++ /issm/trunk-jpl/src/m/classes/frontalforcingsrignotautoregression.py (revision 26553) @@ -0,0 +1,89 @@ +# -*- 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): + """FRONTALFORCINGSRIGNOTAUTOREGRESSION class definition + + Usage: + frontalforcingsrignotautoregression = frontalforcingsrignotautoregression() + """ + + def __init__(self, *args): # {{{ + self.num_basins = 0 + self.beta0 = np.nan + self.beta1 = np.nan + self.ar_order = 0 + self.ar_initialtime = 0 + self.ar_timestep = 0 + self.phi = np.nan + self.basin_id = np.nan + self.subglacial_discharge = np.nan + + if len(args) == 0: + self.setdefaultparameters() + else: + error('constructor not supported') + + def __repr__(self): # {{{ + s = ' Frontalforcings parameters:\n' + s += '{}\n'.format(fielddisplay(self, 'num_basins', 'number of different basins [unitless]')) + s += '{}\n'.format(fielddisplay(self, 'basin_id', 'basin number assigned to each element [unitless]')) + s += '{}\n'.format(fielddisplay(self, 'subglacial_discharge', 'sum of subglacial discharge for each basin [m/d]')) + s += '{}\n'.format(fielddisplay(self, 'beta0', 'basin-specific intercept values [°C]')) + s += '{}\n'.format(fielddisplay(self, 'beta1', 'basin-specific trend values [°C yr^(-1)]')) + s += '{}\n'.format(fielddisplay(self, 'ar_order', 'order of the autoregressive model [unitless]')) + s += '{}\n'.format(fielddisplay(self, 'ar_initialtime', 'initial time assumed in the autoregressive model parameterization [yr]')) + s += '{}\n'.format(fielddisplay(self, 'ar_timestep', 'time resolution of the autoregressive model [yr]')) + s += '{}\n'.format(fielddisplay(self, 'phi', 'basin-specific vectors of lag coefficients [unitless]')) + return s + #}}} + + def setdefaultparameters(self): # {{{ + self.basin_id = np.nan + self.num_basins = 0 + self.subglacial_discharge = np.nan + self.ar_order = 0.0 # Autoregression model of order 0 + return self + #}}} + + def checkconsistency(self, md, solution, analyses): # {{{ + # Early return + if not (solution == 'TransientSolution') or not md.transient.ismovingfront: + return md + + md = checkfield(md, 'fieldname', 'frontalforcings.num_basins', 'numel', 1, 'NaN', 1, 'Inf', 1, '>', 0) + 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) + 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) + md = checkfield(md, 'fieldname', 'frontalforcings.ar_order', 'numel', 1, 'NaN', 1, 'Inf', 1, '>=', 0) + md = checkfield(md, 'fieldname', 'frontalforcings.ar_initialtime', 'numel', 1, 'NaN', 1, 'Inf', 1) + md = checkfield(md, 'fieldname', 'frontalforcings.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', 'frontalforcings.phi', 'NaN', 1, 'Inf', 1, 'size', [md.frontalforcings.num_basins, md.frontalforcings.ar_order]) + return md + # }}} + + def extrude(self, md): # {{{ + # Nothing for now + return self + # }}} + + def marshall(self, prefix, md, fid): # {{{ + yts = md.constants.yts + WriteData(fid, prefix, 'name', 'md.frontalforcings.parameterization', 'data', 55, 'format', 'Integer') + WriteData(fid, prefix, 'object', self, 'class', 'frontalforcings', 'fieldname', 'num_basins', 'format', 'Integer') + WriteData(fid, prefix, 'object', self, 'class', 'frontalforcings', 'fieldname', 'subglacial_discharge', 'format', 'DoubleMat', 'mattype', 1, 'timeserieslength', md.mesh.numberofvertices + 1, 'yts', yts) + WriteData(fid, prefix, 'object', self, 'class', 'frontalforcings', 'fieldname', 'ar_order', 'format', 'Integer') + 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', '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) + WriteData(fid, prefix, 'object', self, 'class', 'frontalforcings', 'fieldname', 'phi', 'format', 'DoubleMat', 'name', 'md.frontalforcings.phi', 'yts', yts) + # }}} Index: /issm/trunk-jpl/src/m/classes/stochasticforcing.m =================================================================== --- /issm/trunk-jpl/src/m/classes/stochasticforcing.m (revision 26552) +++ /issm/trunk-jpl/src/m/classes/stochasticforcing.m (revision 26553) @@ -42,5 +42,5 @@ md = checkfield(md,'fieldname','stochasticforcing.randomflag','numel',[1],'values',[0 1]); - %Check that all fields agree with the corresponding md class + %Check that all fields agree with the corresponding md class for field=self.fields if(contains(field,'SMB')) Index: /issm/trunk-jpl/src/m/classes/stochasticforcing.py =================================================================== --- /issm/trunk-jpl/src/m/classes/stochasticforcing.py (revision 26552) +++ /issm/trunk-jpl/src/m/classes/stochasticforcing.py (revision 26553) @@ -40,4 +40,5 @@ # 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 return self @@ -76,15 +77,5 @@ def marshall(self, prefix, md, fid): # {{{ yts = md.constants.yts - if (type(self.fields) is list): - num_fields = len(self.fields) - # Scaling covariance matrix (scale column-by-column and row-by-row) - scaledfields = ['SMBautoregression'] # list of fields that need scaling * 1/yts - for i in range(num_fields): - if self.fields[i] in scaledfields: - inds = range(1 + np.sum(self.dimensions[0:i]), np.sum(self.dimensions[0:i])) - for row in inds: # scale rows corresponding to scaled field - self.covariance[row, :] = 1 / yts * self.covariance[row, :] - for col in inds: # scale columns corresponding to scaled field - self.covariance[:, col] = 1 / yts * self.covariance[:, col] + num_fields = len(self.fields) WriteData(fid, prefix, 'object', self, 'fieldname', 'isstochasticforcing', 'format', 'Boolean') @@ -92,8 +83,18 @@ return md else: + # Scaling covariance matrix (scale column-by-column and row-by-row) + scaledfields = ['SMBautoregression'] # list of fields that need scaling * 1/yts + tempcovariance = self.covariance + for i in range(num_fields): + if self.fields[i] in scaledfields: + inds = range(int(np.sum(self.dimensions[0:i])), int(np.sum(self.dimensions[0:i + 1]))) + for row in inds: # scale rows corresponding to scaled field + tempcovariance[row, :] = 1 / yts * self.covariance[row, :] + for col in inds: # scale columns corresponding to scaled field + tempcovariance[:, col] = 1 / yts * self.covariance[:, col] 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, 'object', self, 'fieldname', 'covariance', 'format', 'DoubleMat') + WriteData(fid, prefix, 'data', tempcovariance, 'name', 'md.stochasticforcing.covariance', 'format', 'DoubleMat') WriteData(fid, prefix, 'object', self, 'fieldname', 'randomflag', 'format', 'Boolean') # }}} Index: /issm/trunk-jpl/src/m/miscellaneous/MatlabFuncs.py =================================================================== --- /issm/trunk-jpl/src/m/miscellaneous/MatlabFuncs.py (revision 26552) +++ /issm/trunk-jpl/src/m/miscellaneous/MatlabFuncs.py (revision 26553) @@ -50,4 +50,16 @@ return np.degrees(np.arctan2(Y, X)) +# }}} + +def contains(str, pat): #{{{ + """function contains - Determine if pattern is in strings + + Usage: + TF = contains(str, pat) + + TODO: + - Implement 'IgnoreCase' option + """ + # }}} @@ -142,6 +154,23 @@ # }}} +def intersect(A, B): # {{{ + """function intersect - Set intersection of two arrays + + Usage: + C = intersect(A, B) + + NOTE: + - Only the following functionality is currently implemented: + - C = intersect(A,B) returns the data common to both A and B, with no + repetitions. C is in sorted order. + + """ + import numpy as np + + return np.intersect1d(A, B) +#}}} + def isa(A, dataType): # {{{ - """FUNCTION ISA + """function isa NOTE: @@ -216,4 +245,15 @@ def mod(a, m): # {{{ return a % m +# }}} + +def numel(A): # {{{ + """function numel - Number of array elements + + Usage: + n = numel(A)) + """ + import numpy as np + + return np.size(A) # }}} Index: /issm/trunk-jpl/src/m/solve/WriteData.py =================================================================== --- /issm/trunk-jpl/src/m/solve/WriteData.py (revision 26552) +++ /issm/trunk-jpl/src/m/solve/WriteData.py (revision 26553) @@ -67,5 +67,10 @@ #We scale everything but the last line that holds time data[:-1, :] = scale * data[:-1, :] - else: + elif type(data) is list: # Deal with "TypeError: can't multiply sequence by non-int of type 'float'" for type list + scaleddata = [] + for i in range(len(data)): + scaleddata.append(scale * data[i]) + data = scaleddata + else: # data = scale * data if np.ndim(data) > 1 and data.shape[0] == timeserieslength: Index: /issm/trunk-jpl/test/NightlyRun/runme.py =================================================================== --- /issm/trunk-jpl/test/NightlyRun/runme.py (revision 26552) +++ /issm/trunk-jpl/test/NightlyRun/runme.py (revision 26553) @@ -1,3 +1,10 @@ #!/usr/bin/env python +# +# NOTE: Switch to, +# +# #!/usr/bin/python3 +# +# when making Python 3 default +# import argparse from glob import glob Index: /issm/trunk-jpl/test/NightlyRun/test257.m =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test257.m (revision 26552) +++ /issm/trunk-jpl/test/NightlyRun/test257.m (revision 26553) @@ -5,5 +5,5 @@ md=setflowequation(md,'SSA','all'); md.cluster=generic('name',oshostname(),'np',3); -md.transient.requested_outputs={'default','IceVolume','SmbMassBalance',}; +md.transient.requested_outputs={'default','IceVolume','SmbMassBalance'}; ymax = max(md.mesh.y); @@ -66,5 +66,5 @@ (md.results.TransientSolution(1).Thickness),... (md.results.TransientSolution(1).IceVolume),... - (md.results.TransientSolution(1).SmbMassBalance),... + (md.results.TransientSolution(1).SmbMassBalance),... (md.results.TransientSolution(2).Vx),... (md.results.TransientSolution(2).Vy),... @@ -72,5 +72,5 @@ (md.results.TransientSolution(2).Thickness),... (md.results.TransientSolution(2).IceVolume),... - (md.results.TransientSolution(2).SmbMassBalance),... + (md.results.TransientSolution(2).SmbMassBalance),... (md.results.TransientSolution(3).Vx),... (md.results.TransientSolution(3).Vy),... @@ -78,4 +78,4 @@ (md.results.TransientSolution(3).Thickness),... (md.results.TransientSolution(3).IceVolume),... - (md.results.TransientSolution(3).SmbMassBalance),... + (md.results.TransientSolution(3).SmbMassBalance),... }; Index: /issm/trunk-jpl/test/NightlyRun/test257.py =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test257.py (revision 26553) +++ /issm/trunk-jpl/test/NightlyRun/test257.py (revision 26553) @@ -0,0 +1,94 @@ +#Test Name: SquareShelfSMBautoregression +from __future__ import division # TODO: Remove this import after ISSM is updated to use Python 3 by default + +import numpy as np + +from MatlabFuncs import * +from model import * +from parameterize import * +from setflowequation import * +from setmask import * +from SMBautoregression import * +from solve import * +from triangle import * + + +md = triangle(model(), '../Exp/Square.exp', 80000) +md = setmask(md, 'all', '') +md = parameterize(md, '../Par/SquareShelfConstrained.py') +md = setflowequation(md, 'SSA', 'all') +md.cluster = generic('name', oshostname(), 'np', 3) +md.transient.requested_outputs = ['default', 'IceVolume', 'SmbMassBalance'] + +ymax = max(md.mesh.y) +xmax = max(md.mesh.x) +# Generate basin IDs for 3 basins +idbasin = np.zeros((md.mesh.numberofelements,)) +iid1 = np.where(md.mesh.y >= 2/3 * ymax)[0] +iid2 = intersect(np.where(md.mesh.y < 2/3 * ymax)[0], np.where(md.mesh.x >= 1/3 * xmax)[0]) +iid3 = intersect(np.where(md.mesh.y < 2/3 * ymax)[0], np.where(md.mesh.x < 1/3 * xmax)[0]) +for ii in range(md.mesh.numberofelements): + for vertex in range(3): + if md.mesh.elements[ii][vertex] - 1 in iid1: # one vertex in basin 1; NOTE: offset because of 1-based vertex indexing + idbasin[ii] = 1 + if idbasin[ii] == 0: # no vertex was found in basin 1 + for vertex in range(3): + if md.mesh.elements[ii][vertex] - 1 in iid2: # one vertex in basin 2; NOTE: offset because of 1-based vertex indexing + idbasin[ii] = 2 + if idbasin[ii] == 0: # no vertex was found in basin 1 and 2 + idbasin[ii] = 3 + +# SMB parameters +md.timestepping.start_time = 0 +md.timestepping.time_step = 1 +md.timestepping.final_time = 5 +md.smb = SMBautoregression() +md.smb.num_basins = 3 # number of basins +md.smb.basin_id = idbasin - 1 # prescribe basin ID number to elements; # NOTE: offset because of 1-based vertex indexing +md.smb.beta0 = np.array([[0.5, 1.2, 1.5]]) # intercept values of SMB in basins [m ice eq./yr] +md.smb.beta1 = np.array([[0.0, 0.01, -0.01]]) # trend values of SMB in basins [m ice eq./yr^2] +md.smb.ar_initialtime = md.timestepping.start_time +md.smb.ar_order = 4 +md.smb.ar_timestep = 2.0 #timestep of the autoregressive model [yr] +md.smb.phi = np.array([[0.2, 0.1, 0.05, 0.01], [0.4, 0.2, -0.2, 0.1], [0.4, -0.4, 0.1, -0.1]]) + +# Stochastic forcing +md.stochasticforcing.isstochasticforcing = 1 +md.stochasticforcing.fields = ['SMBautoregression'] +md.stochasticforcing.dimensions = [md.smb.num_basins] # dimension of each field +md.stochasticforcing.covariance = np.array([[0.15, 0.08, -0.02], [0.08, 0.12, -0.05], [-0.02, -0.05, 0.1]]) # global covariance among- and between-fields +md.stochasticforcing.randomflag = 0 # fixed random seeds + +md = solve(md, 'Transient') + +# Fields and tolerances to track changes +field_names = [ + 'Vx1', 'Vy1', 'Vel1', 'Thickness1', 'Volume1', 'SmbMassBalance1', + 'Vx2', 'Vy2', 'Vel2', 'Thickness2', 'Volume2', 'SmbMassBalance2', + 'Vx3', 'Vy3', 'Vel3', 'Thickness3', 'Volume3', 'SmbMassBalance3' +] +field_tolerances = [ + 1e-13, 1e-13, 1e-13, 1e-13, 1e-13, 1e-13, + 1e-13, 1e-13, 1e-13, 1e-13, 1e-13, 1e-13, + 1e-13, 1e-13, 1e-13, 1e-13, 1e-13, 1e-13 +] +field_values = [ + md.results.TransientSolution[0].Vx, + md.results.TransientSolution[0].Vy, + md.results.TransientSolution[0].Vel, + md.results.TransientSolution[0].Thickness, + md.results.TransientSolution[0].IceVolume, + md.results.TransientSolution[0].SmbMassBalance, + md.results.TransientSolution[1].Vx, + md.results.TransientSolution[1].Vy, + md.results.TransientSolution[1].Vel, + md.results.TransientSolution[1].Thickness, + md.results.TransientSolution[1].IceVolume, + md.results.TransientSolution[1].SmbMassBalance, + md.results.TransientSolution[2].Vx, + md.results.TransientSolution[2].Vy, + md.results.TransientSolution[2].Vel, + md.results.TransientSolution[2].Thickness, + md.results.TransientSolution[2].IceVolume, + md.results.TransientSolution[2].SmbMassBalance +] Index: /issm/trunk-jpl/test/NightlyRun/test542.m =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test542.m (revision 26552) +++ /issm/trunk-jpl/test/NightlyRun/test542.m (revision 26553) @@ -34,4 +34,5 @@ end end + md.transient.ismovingfront = 1; md.levelset.spclevelset = NaN(md.mesh.numberofvertices,1); Index: /issm/trunk-jpl/test/NightlyRun/test542.py =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test542.py (revision 26553) +++ /issm/trunk-jpl/test/NightlyRun/test542.py (revision 26553) @@ -0,0 +1,91 @@ +#Test Name: PigTranFrontalforcingsrignot +import numpy as np + +from calvingvonmises import * +from frontalforcingsrignot import * +from MatlabFuncs import * +from model import * +from parameterize import * +from setflowequation import * +from setmask import * +from solve import * +from triangle import * + + +md = triangle(model(), '../Exp/Pig.exp', 10000) +md = setmask(md, '../Exp/PigShelves.exp', '../Exp/PigIslands.exp') +md = parameterize(md, '../Par/Pig.py') +md = setflowequation(md, 'SSA', 'all') +md.timestepping.time_step = 2 +md.timestepping.final_time = 50 + +# Separate domain in 2 basins +idbasin = np.zeros((md.mesh.numberofelements,)) +iid1 = np.where(md.mesh.x <= -1.6e6)[0] +for ii in range(md.mesh.numberofelements): + for vertex in range(3): + if md.mesh.elements[ii][vertex] - 1 in iid1: # one vertex in basin 1; NOTE: offset because of 1-based vertex indexing + idbasin[ii] = 1 + if idbasin[ii] == 0: # no vertex was found in basin 1 + idbasin[ii] = 2 + +# Calving and frontalforcings parameters +md.mask.ice_levelset = 1e4 * (md.mask.ice_levelset + 0.5) +md.calving = calvingvonmises() +md.frontalforcings = frontalforcingsrignot() +md.frontalforcings.num_basins = 2 +md.frontalforcings.basin_id = idbasin - 1 # NOTE: offset because of 1-based vertex indexing +md.frontalforcings.subglacial_discharge = 0.1 * np.ones((md.mesh.numberofvertices,)) +md.frontalforcings.thermalforcing = 0.5 * np.ones((md.mesh.numberofvertices,)) +for elem in range(md.mesh.numberofelements): + if idbasin[elem] == 2: + md.frontalforcings.thermalforcing[md.mesh.elements[elem, 0:3] - 1] = 1.5 # NOTE: offset because of 1-based vertex indexing + +md.transient.ismovingfront = 1 +md.levelset.spclevelset = np.full((md.mesh.numberofvertices,), np.nan) +md.levelset.migration_max = 1e10 + +md.transient.requested_outputs = ['default', 'CalvingMeltingrate'] +md.cluster = generic('name', oshostname(), 'np', 2) +md = solve(md, 'Transient') + +# Fields and tolerances to track changes +field_names = [ + 'Vx1', 'Vy1', 'Vel1', 'Pressure1', 'Bed1', 'Surface1', 'Thickness1', 'MaskIceLevelset1', 'CalvingMetlingRate1', + 'Vx2', 'Vy2', 'Vel2', 'Pressure2', 'Bed2', 'Surface2', 'Thickness2', 'MaskIceLevelset2', 'CalvingMetlingRate2', + 'Vx10', 'Vy10', 'Vel10', 'Pressure10', 'Bed10', 'Surface10', 'Thickness10', 'MaskIceLevelset10', 'CalvingMetlingRate10' +] +field_tolerances = [ + 1e-11, 2e-11, 2e-11, 1e-12, 2e-11, 6e-12, 9e-12, 1e-12, 1e-9, + 2e-11, 1e-11, 1e-11, 9e-12, 2e-1, 2e-11, 2e-11, 1e-11, 1e-9, + 2e-10, 1e-10, 1e-10, 1e-10, 1e-10, 1e-10, 1e-10, 1e-9, 1e-9 +] +field_values = [ + md.results.TransientSolution[0].Vx, + md.results.TransientSolution[0].Vy, + md.results.TransientSolution[0].Vel, + md.results.TransientSolution[0].Pressure, + md.results.TransientSolution[0].Base, + md.results.TransientSolution[0].Surface, + md.results.TransientSolution[0].Thickness, + md.results.TransientSolution[0].MaskIceLevelset, + md.results.TransientSolution[0].CalvingMeltingrate, + md.results.TransientSolution[1].Vx, + md.results.TransientSolution[1].Vy, + md.results.TransientSolution[1].Vel, + md.results.TransientSolution[1].Pressure, + md.results.TransientSolution[1].Base, + md.results.TransientSolution[1].Surface, + md.results.TransientSolution[1].Thickness, + md.results.TransientSolution[1].MaskIceLevelset, + md.results.TransientSolution[1].CalvingMeltingrate, + md.results.TransientSolution[9].Vx, + md.results.TransientSolution[9].Vy, + md.results.TransientSolution[9].Vel, + md.results.TransientSolution[9].Pressure, + md.results.TransientSolution[9].Base, + md.results.TransientSolution[9].Surface, + md.results.TransientSolution[9].Thickness, + md.results.TransientSolution[9].MaskIceLevelset, + md.results.TransientSolution[9].CalvingMeltingrate +] Index: /issm/trunk-jpl/test/NightlyRun/test543.m =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test543.m (revision 26552) +++ /issm/trunk-jpl/test/NightlyRun/test543.m (revision 26553) @@ -40,5 +40,5 @@ md.stochasticforcing.fields = [{'FrontalForcingsRignotAutoregression'}]; md.stochasticforcing.dimensions = [md.frontalforcings.num_basins]; %dimension of each field -md.stochasticforcing.covariance = 1e-4*[[1.5,0.5];[0.5,0.4]];; %global covariance among- and between-fields +md.stochasticforcing.covariance = 1e-4*[[1.5,0.5];[0.5,0.4]]; %global covariance among- and between-fields md.stochasticforcing.randomflag = 0; %determines true/false randomness Index: /issm/trunk-jpl/test/NightlyRun/test543.py =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test543.py (revision 26553) +++ /issm/trunk-jpl/test/NightlyRun/test543.py (revision 26553) @@ -0,0 +1,102 @@ +#Test Name: PigTranFrontalforcingsrignotautoregression +import numpy as np + +from calvingvonmises import * +from frontalforcingsrignotautoregression import * +from MatlabFuncs import * +from model import * +from parameterize import * +from setflowequation import * +from setmask import * +from solve import * +from triangle import * + + +md = triangle(model(), '../Exp/Pig.exp', 10000) +md = setmask(md, '../Exp/PigShelves.exp', '../Exp/PigIslands.exp') +md = parameterize(md, '../Par/Pig.py') +md = setflowequation(md, 'SSA', 'all') +md.timestepping.start_time = 0 +md.timestepping.time_step = 1 +md.timestepping.final_time = 25 + +# Separate domain in 2 basins +idbasin = np.zeros((md.mesh.numberofelements,)) +iid1 = np.where(md.mesh.x <= -1.6e6)[0] +for ii in range(md.mesh.numberofelements): + for vertex in range(3): + if md.mesh.elements[ii][vertex] - 1 in iid1: # one vertex in basin 1; NOTE: offset because of 1-based vertex indexing + idbasin[ii] = 1 + if idbasin[ii] == 0: # no vertex was found in basin 1 + for vertex in range(3): + idbasin[ii] = 2 + +# Calving and frontalforcings parameters +md.mask.ice_levelset = 1e4 * (md.mask.ice_levelset + 0.5) +md.calving = calvingvonmises() +md.frontalforcings = frontalforcingsrignotautoregression() +md.frontalforcings.num_basins = 2 +md.frontalforcings.basin_id = idbasin - 1 # NOTE: offset because of 1-based vertex indexing +md.frontalforcings.subglacial_discharge = 0.1 * np.ones((md.mesh.numberofvertices,)) +md.frontalforcings.beta0 = np.array([[0.05, 0.01]]) # intercept values of TF in basins [C] +md.frontalforcings.beta1 = np.array([[0.001, 0.0001]]) # trend values of TF in basins [C/yr] +md.frontalforcings.ar_initialtime = md.timestepping.start_time # initial time in the AR model parameterization [yr] +md.frontalforcings.ar_order = 4 +md.frontalforcings.ar_timestep = 2 # timestep of the autoregressive model [yr] +md.frontalforcings.phi = np.array([[0.1, -0.1, 0.01, -0.01], [0.2, -0.2, 0.1, 0.0]]) # autoregressive parameters + +# Stochastic forcing +md.stochasticforcing.isstochasticforcing = 1 +md.stochasticforcing.fields = ['FrontalForcingsRignotAutoregression'] +md.stochasticforcing.dimensions = [md.frontalforcings.num_basins] # dimension of each field +md.stochasticforcing.covariance = 1e-4 * np.array([[1.5, 0.5], [0.5, 0.4]]) # global covariance among- and between-fields +md.stochasticforcing.randomflag = 0 # determines true/false randomness + +md.transient.ismovingfront = 1 +md.levelset.spclevelset = np.full((md.mesh.numberofvertices,), np.nan) +md.levelset.migration_max = 1e10 + +md.transient.requested_outputs = ['default', 'CalvingMeltingrate'] +md.cluster = generic('name', oshostname(), 'np', 2) +md = solve(md, 'Transient') + +# Fields and tolerances to track changes +field_names = [ + 'Vx1', 'Vy1', 'Vel1', 'Pressure1', 'Bed1', 'Surface1', 'Thickness1', 'MaskIceLevelset1', 'CalvingMetlingRate1', + 'Vx2', 'Vy2', 'Vel2', 'Pressure2', 'Bed2', 'Surface2', 'Thickness2', 'MaskIceLevelset2', 'CalvingMetlingRate2', + 'Vx10', 'Vy10', 'Vel10', 'Pressure10', 'Bed10', 'Surface10', 'Thickness10', 'MaskIceLevelset10', 'CalvingMetlingRate10' +] +field_tolerances = [ + 1e-11, 2e-11, 2e-11, 1e-12, 2e-11, 6e-12, 9e-12, 1e-12, 1e-9, + 2e-11, 1e-11, 1e-11, 9e-12, 2e-1, 2e-11, 2e-11, 1e-11, 1e-9, + 2e-10, 1e-10, 1e-10, 1e-10, 1e-10, 1e-10, 1e-10, 1e-9, 1e-9 +] +field_values = [ + md.results.TransientSolution[0].Vx, + md.results.TransientSolution[0].Vy, + md.results.TransientSolution[0].Vel, + md.results.TransientSolution[0].Pressure, + md.results.TransientSolution[0].Base, + md.results.TransientSolution[0].Surface, + md.results.TransientSolution[0].Thickness, + md.results.TransientSolution[0].MaskIceLevelset, + md.results.TransientSolution[0].CalvingMeltingrate, + md.results.TransientSolution[1].Vx, + md.results.TransientSolution[1].Vy, + md.results.TransientSolution[1].Vel, + md.results.TransientSolution[1].Pressure, + md.results.TransientSolution[1].Base, + md.results.TransientSolution[1].Surface, + md.results.TransientSolution[1].Thickness, + md.results.TransientSolution[1].MaskIceLevelset, + md.results.TransientSolution[1].CalvingMeltingrate, + md.results.TransientSolution[9].Vx, + md.results.TransientSolution[9].Vy, + md.results.TransientSolution[9].Vel, + md.results.TransientSolution[9].Pressure, + md.results.TransientSolution[9].Base, + md.results.TransientSolution[9].Surface, + md.results.TransientSolution[9].Thickness, + md.results.TransientSolution[9].MaskIceLevelset, + md.results.TransientSolution[9].CalvingMeltingrate +]