Index: /issm/trunk-jpl/src/m/boundaryconditions/getlovenumbers.py =================================================================== --- /issm/trunk-jpl/src/m/boundaryconditions/getlovenumbers.py (revision 26300) +++ /issm/trunk-jpl/src/m/boundaryconditions/getlovenumbers.py (revision 26301) @@ -5,21 +5,25 @@ def getlovenumbers(*args): #{{{ - ''' - GETLOVENUMBERS - provide love numbers retrieved from: + """GETLOVENUMBERS - provide love numbers retrieved from: http://www.srosat.com/iag-jsg/loveNb.php in a chosen reference frame - Usage: - series = love_numbers('type', 'loadingverticaldisplacement', 'referenceframe', 'CM', 'maxdeg', 1001) + Usage: + series = love_numbers('type', 'loadingverticaldisplacement', 'referenceframe', 'CM', 'maxdeg', 1001) - - type = one of 'loadingverticaldisplacement', - 'loadinggravitationalpotential', 'loadinghorizontaldisplacement', - 'tidalverticaldisplacement', 'tidalgravitationalpotential', - 'tidalhorizontaldisplacement' - - reference_frame = one of 'CM' (default) and 'CF' - - maxdeg = default 1000 + - type = one of 'loadingverticaldisplacement', + 'loadinggravitationalpotential', 'loadinghorizontaldisplacement', + 'tidalverticaldisplacement', 'tidalgravitationalpotential', + 'tidalhorizontaldisplacement' + - reference_frame = one of 'CM' (default) and 'CF' + - maxdeg = default 1000 - Example: - h = love_number - ''' + Example: + h = getlovenumbers('type', 'loadingverticaldisplacement', 'referenceframe', 'CM', 'maxdeg', maxdeg) + k = getlovenumbers('type', 'loadinggravitationalpotential', 'referenceframe', 'CM', 'maxdeg', maxdeg) + l = getlovenumbers('type', 'loadinghorizontaldisplacement', 'referenceframe', 'CM', 'maxdeg', maxdeg) + th = getlovenumbers('type', 'tidalverticaldisplacement', 'referenceframe', 'CM', 'maxdeg', maxdeg) + tk = getlovenumbers('type', 'tidalgravitationalpotential', 'referenceframe', 'CM', 'maxdeg', maxdeg) + tl = getlovenumbers('type', 'tidalhorizontaldisplacement', 'referenceframe', 'CM', 'maxdeg', maxdeg) + """ TYPES = [ @@ -32,9 +36,12 @@ ] - #recover options + # Recover options options = pairoptions(*args) type = options.getfieldvalue('type') frame = options.getfieldvalue('referenceframe', 'CM') maxdeg = options.getfieldvalue('maxdeg', 1000) + + if (maxdeg > 10000): + raise Exception('PREM love numbers computed only for deg < 10,000. Request lower maxdeg') if type not in TYPES: @@ -10045,8 +10052,8 @@ ]) - #cut off series at maxdeg + # Cut off series at maxdeg love_numbers = np.delete(love_numbers, range(maxdeg + 1, len(love_numbers)), axis=0) - #retrieve right type + # Retrieve right type if type == 'loadingverticaldisplacement': series = love_numbers[:, 0] @@ -10062,5 +10069,5 @@ series = love_numbers[:, 5] else: - raise Exception("love_numbers error message: unknown type: {}".format(type)) + raise Exception('love_numbers error message: unknown type: {}'.format(type)) #choose degree 1 term for CF reference system @@ -10075,5 +10082,5 @@ series[1] = 0.134 else: - raise Exception("love_numbers error message: unknown reference frame: {}".format(frame)) + raise Exception('love_numbers error message: unknown reference frame: {}'.format(frame)) return series Index: /issm/trunk-jpl/src/m/classes/amr.m =================================================================== --- /issm/trunk-jpl/src/m/classes/amr.m (revision 26300) +++ /issm/trunk-jpl/src/m/classes/amr.m (revision 26301) @@ -84,5 +84,5 @@ self.gradation=1.5; - %other criterias + %other criteria self.groundingline_resolution=500.; self.groundingline_distance=0.; Index: /issm/trunk-jpl/src/m/classes/amr.py =================================================================== --- /issm/trunk-jpl/src/m/classes/amr.py (revision 26300) +++ /issm/trunk-jpl/src/m/classes/amr.py (revision 26301) @@ -5,6 +5,5 @@ class amr(object): - """ - AMR Class definition + """AMR Class definition Usage: @@ -12,75 +11,85 @@ """ - def __init__(self): # {{{ - self.hmin = 0. - self.hmax = 0. + def __init__(self): #{{{ + self.hmin = 0 + self.hmax = 0 self.fieldname = '' - self.err = 0. - self.keepmetric = 0. - self.gradation = 0. - self.groundingline_resolution = 0. - self.groundingline_distance = 0. - self.icefront_resolution = 0. - self.icefront_distance = 0. - self.thicknesserror_resolution = 0. - self.thicknesserror_threshold = 0. - self.thicknesserror_groupthreshold = 0. - self.thicknesserror_maximum = 0. - self.deviatoricerror_resolution = 0. - self.deviatoricerror_threshold = 0. - self.deviatoricerror_groupthreshold = 0. - self.deviatoricerror_maximum = 0. - self.restart = 0. - #set defaults + self.err = 0 + self.keepmetric = 0 + self.gradation = 0 + self.groundingline_resolution = 0 + self.groundingline_distance = 0 + self.icefront_resolution = 0 + self.icefront_distance = 0 + self.thicknesserror_resolution = 0 + self.thicknesserror_threshold = 0 + self.thicknesserror_groupthreshold = 0 + self.thicknesserror_maximum = 0 + self.deviatoricerror_resolution = 0 + self.deviatoricerror_threshold = 0 + self.deviatoricerror_groupthreshold = 0 + self.deviatoricerror_maximum = 0 + self.restart = 0 + self.setdefaultparameters() #}}} - def __repr__(self): # {{{ - string = " amr parameters:" - string = "%s\n%s" % (string, fielddisplay(self, "hmin", "minimum element length")) - string = "%s\n%s" % (string, fielddisplay(self, "hmax", "maximum element length")) - string = "%s\n%s" % (string, fielddisplay(self, "fieldname", "name of input that will be used to compute the metric (should be an input of FemModel)")) - string = "%s\n%s" % (string, fielddisplay(self, "keepmetric", "indicates whether the metric should be kept every remeshing time")) - string = "%s\n%s" % (string, fielddisplay(self, "gradation", "maximum ratio between two adjacent edges")) - string = "%s\n%s" % (string, fielddisplay(self, "groundingline_resolution", "element length near the grounding line")) - string = "%s\n%s" % (string, fielddisplay(self, "groundingline_distance", "distance around the grounding line which elements will be refined")) - string = "%s\n%s" % (string, fielddisplay(self, "icefront_resolution", "element length near the ice front")) - string = "%s\n%s" % (string, fielddisplay(self, "icefront_distance", "distance around the ice front which elements will be refined")) - string = "%s\n%s" % (string, fielddisplay(self, "thicknesserror_resolution", "element length when thickness error estimator is used")) - string = "%s\n%s" % (string, fielddisplay(self, "thicknesserror_threshold", "maximum threshold thickness error permitted")) - string = "%s\n%s" % (string, fielddisplay(self, "thicknesserror_groupthreshold", "maximum group threshold thickness error permitted")) - string = "%s\n%s" % (string, fielddisplay(self, "thicknesserror_maximum", "maximum thickness error permitted")) - string = "%s\n%s" % (string, fielddisplay(self, "deviatoricerror_resolution", "element length when deviatoric stress error estimator is used")) - string = "%s\n%s" % (string, fielddisplay(self, "deviatoricerror_threshold", "maximum threshold deviatoricstress error permitted")) - string = "%s\n%s" % (string, fielddisplay(self, "deviatoricerror_groupthreshold", "maximum group threshold deviatoric stress error permitted")) - string = "%s\n%s" % (string, fielddisplay(self, "deviatoricerror_maximum", "maximum deviatoricstress error permitted")) - string = "%s\n%s" % (string, fielddisplay(self, "restart", "indicates if ReMesh() will call before first time step")) - return string + def __repr__(self): #{{{ + s = ' amr parameters:\n' + s += '{}\n'.format(fielddisplay(self, 'hmin', 'minimum element length')) + s += '{}\n'.format(fielddisplay(self, 'hmax', 'maximum element length')) + s += '{}\n'.format(fielddisplay(self, 'fieldname', 'name of input that will be used to compute the metric (should be an input of FemModel)')) + s += '{}\n'.format(fielddisplay(self, 'keepmetric', 'indicates whether the metric should be kept every remeshing time')) + s += '{}\n'.format(fielddisplay(self, 'gradation', 'maximum ratio between two adjacent edges')) + s += '{}\n'.format(fielddisplay(self, 'groundingline_resolution', 'element length near the grounding line')) + s += '{}\n'.format(fielddisplay(self, 'groundingline_distance', 'distance around the grounding line which elements will be refined')) + s += '{}\n'.format(fielddisplay(self, 'icefront_resolution', 'element length near the ice front')) + s += '{}\n'.format(fielddisplay(self, 'icefront_distance', 'distance around the ice front which elements will be refined')) + s += '{}\n'.format(fielddisplay(self, 'thicknesserror_resolution', 'element length when thickness error estimator is used')) + s += '{}\n'.format(fielddisplay(self, 'thicknesserror_threshold', 'maximum threshold thickness error permitted')) + s += '{}\n'.format(fielddisplay(self, 'thicknesserror_groupthreshold', 'maximum group threshold thickness error permitted')) + s += '{}\n'.format(fielddisplay(self, 'thicknesserror_maximum', 'maximum thickness error permitted')) + s += '{}\n'.format(fielddisplay(self, 'deviatoricerror_resolution', 'element length when deviatoric stress error estimator is used')) + s += '{}\n'.format(fielddisplay(self, 'deviatoricerror_threshold', 'maximum threshold deviatoricstress error permitted')) + s += '{}\n'.format(fielddisplay(self, 'deviatoricerror_groupthreshold', 'maximum group threshold deviatoric stress error permitted')) + s += '{}\n'.format(fielddisplay(self, 'deviatoricerror_maximum', 'maximum deviatoricstress error permitted')) + s += '{}\n'.format(fielddisplay(self, 'restart', 'indicates if ReMesh() will call before first time step')) + return s #}}} - def setdefaultparameters(self): # {{{ - self.hmin = 100. - self.hmax = 100.e3 + def setdefaultparameters(self): #{{{ + self.hmin = 100 + self.hmax = 100e3 + + # Fields self.fieldname = 'Vel' - self.err = 3. + self.err = 3 + + # Keep metric? self.keepmetric = 1 + + # Control of element lengths self.gradation = 1.5 - self.groundingline_resolution = 500. + + # Other criteria + self.groundingline_resolution = 500 self.groundingline_distance = 0 - self.icefront_resolution = 500. + self.icefront_resolution = 500 self.icefront_distance = 0 - self.thicknesserror_resolution = 500. + self.thicknesserror_resolution = 500 self.thicknesserror_threshold = 0 self.thicknesserror_groupthreshold = 0 self.thicknesserror_maximum = 0 - self.deviatoricerror_resolution = 500. + self.deviatoricerror_resolution = 500 self.deviatoricerror_threshold = 0 self.deviatoricerror_groupthreshold = 0 self.deviatoricerror_maximum = 0 - self.restart = 0. + + # Is restart? This calls femmodel->ReMesh() before first time step. + self.restart = 0 return self #}}} - def checkconsistency(self, md, solution, analyses): # {{{ + def checkconsistency(self, md, solution, analyses): #{{{ md = checkfield(md, 'fieldname', 'amr.hmax', 'numel', [1], '>', 0, 'NaN', 1) md = checkfield(md, 'fieldname', 'amr.hmin', 'numel', [1], '>', 0, '<', self.hmax, 'NaN', 1) @@ -101,7 +110,7 @@ md = checkfield(md, 'fieldname', 'amr.restart', 'numel', [1], '>=', 0, '<=', 1, 'NaN', 1) return md - # }}} + # }}} - def marshall(self, prefix, md, fid): # {{{ + def marshall(self, prefix, md, fid): #{{{ WriteData(fid, prefix, 'name', 'md.amr.type', 'data', 1, 'format', 'Integer') WriteData(fid, prefix, 'object', self, 'fieldname', 'hmin', 'format', 'Double') @@ -124,3 +133,3 @@ WriteData(fid, prefix, 'object', self, 'fieldname', 'deviatoricerror_maximum', 'format', 'Double') WriteData(fid, prefix, 'object', self, 'class', 'amr', 'fieldname', 'restart', 'format', 'Integer') - # }}} + #}}} Index: /issm/trunk-jpl/src/m/classes/dsl.m =================================================================== --- /issm/trunk-jpl/src/m/classes/dsl.m (revision 26300) +++ /issm/trunk-jpl/src/m/classes/dsl.m (revision 26301) @@ -13,8 +13,4 @@ end methods - function self = extrude(self,md) % {{{ - self.sea_surface_height_above_geoid=project3d(md,'vector',self.sea_surface_height_above_geoid,'type','node','layer',1); - self.sea_water_pressure_at_sea_floor=project3d(md,'vector',self.sea_water_pressure_at_sea_floor,'type','node','layer',1); - end % }}} function self = dsl(varargin) % {{{ switch nargin @@ -34,4 +30,12 @@ end % }}} + function disp(self) % {{{ + + disp(sprintf(' dsl parameters:')); + fielddisplay(self,'global_average_thermosteric_sea_level','Corresponds to zostoga field in CMIP5 archives. Specified as a temporally variable quantity (in m).'); + fielddisplay(self,'sea_surface_height_above_geoid','Corresponds to zos field in CMIP5 archives. Spatial average is 0. Specified as a spatio-temporally variable quantity (in m).'); + fielddisplay(self,'sea_water_pressure_at_sea_floor','Corresponds to bpo field in CMIP5 archives. Specified as a spatio-temporally variable quantity (in m equivalent, not in Pa!).'); + + end % }}} function md = checkconsistency(self,md,solution,analyses) % {{{ @@ -49,26 +53,15 @@ end % }}} - function disp(self) % {{{ - - disp(sprintf(' dsl parameters:')); - fielddisplay(self,'global_average_thermosteric_sea_level','Corresponds to zostoga field in CMIP5 archives. Specified as a temporally variable quantity (in m).'); - fielddisplay(self,'sea_surface_height_above_geoid','Corresponds to zos field in CMIP5 archives. Spatial average is 0. Specified as a spatio-temporally variable quantity (in m).'); - fielddisplay(self,'sea_water_pressure_at_sea_floor','Corresponds to bpo field in CMIP5 archives. Specified as a spatio-temporally variable quantity (in m equivalent, not in Pa!).'); + function marshall(self,prefix,md,fid) % {{{ + yts=md.constants.yts; + WriteData(fid,prefix,'name','md.dsl.model','data',1,'format','Integer'); + WriteData(fid,prefix,'object',self,'fieldname','global_average_thermosteric_sea_level','format','DoubleMat','mattype',2,'timeseries',1,'timeserieslength',2,'yts',yts); %mattype 2, because we are sending a GMSL value identical everywhere on each element. + WriteData(fid,prefix,'object',self,'fieldname','sea_surface_height_above_geoid','format','DoubleMat','mattype',1,'timeserieslength',md.mesh.numberofvertices+1,'yts',yts); %mattype 1 because we specify DSL at vertex locations. + WriteData(fid,prefix,'object',self,'fieldname','sea_water_pressure_at_sea_floor','format','DoubleMat','mattype',1,'timeserieslength',md.mesh.numberofvertices+1,'yts',yts); %mattype 1 because we specify bottom pressure at vertex locations. end % }}} - function marshall(self,prefix,md,fid) % {{{ - - WriteData(fid,prefix,'name','md.dsl.model','data',1,'format','Integer'); - WriteData(fid,prefix,'object',self,'fieldname','global_average_thermosteric_sea_level','format','DoubleMat','mattype',2,'timeseries',1,'timeserieslength',2,'yts',md.constants.yts); %mattype 2, because we are sending a GMSL value identical everywhere on each element. - WriteData(fid,prefix,'object',self,'fieldname','sea_surface_height_above_geoid','format','DoubleMat','mattype',1,'timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts); %mattype 1 because we specify DSL at vertex locations. - WriteData(fid,prefix,'object',self,'fieldname','sea_water_pressure_at_sea_floor','format','DoubleMat','mattype',1,'timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts); %mattype 1 because we specify bottom pressure at vertex locations. - - end % }}} - function savemodeljs(self,fid,modelname) % {{{ - - writejs1Darray(fid,[modelname '.dsl.global_average_thermosteric_sea_level'],self.global_average_thermosteric_sea_level); - writejs1Darray(fid,[modelname '.dsl.sea_surface_height_above_geoid'],self.sea_surface_height_above_geoid); - writejs1Darray(fid,[modelname '.dsl.sea_water_pressure_at_sea_floor'],self.sea_water_pressure_at_sea_floor); - + function self = extrude(self,md) % {{{ + self.sea_surface_height_above_geoid=project3d(md,'vector',self.sea_surface_height_above_geoid,'type','node','layer',1); + self.sea_water_pressure_at_sea_floor=project3d(md,'vector',self.sea_water_pressure_at_sea_floor,'type','node','layer',1); end % }}} function self = initialize(self,md) % {{{ @@ -87,5 +80,11 @@ end end % }}} - + function savemodeljs(self,fid,modelname) % {{{ + + writejs1Darray(fid,[modelname '.dsl.global_average_thermosteric_sea_level'],self.global_average_thermosteric_sea_level); + writejs1Darray(fid,[modelname '.dsl.sea_surface_height_above_geoid'],self.sea_surface_height_above_geoid); + writejs1Darray(fid,[modelname '.dsl.sea_water_pressure_at_sea_floor'],self.sea_water_pressure_at_sea_floor); + + end % }}} end end Index: /issm/trunk-jpl/src/m/classes/dsl.py =================================================================== --- /issm/trunk-jpl/src/m/classes/dsl.py (revision 26300) +++ /issm/trunk-jpl/src/m/classes/dsl.py (revision 26301) @@ -8,14 +8,19 @@ class dsl(object): - """DSL - slass definition + """DSL - class definition Usage: - dsl = dsl() + dsl = dsl() # dynamic sea level class, based on CMIP5 outputs """ - def __init__(self): #{{{ + def __init__(self, *args): #{{{ self.global_average_thermosteric_sea_level = np.nan # Corresponds to zostoga field in CMIP5 archives. Specified as a temporally variable quantity (in m). self.sea_surface_height_above_geoid = np.nan # Corresponds to zos field in CMIP5 archives. Spatial average is 0. Specified as a spatio-temporally variable quantity (in m). self.sea_water_pressure_at_sea_floor = np.nan # Corresponds to bpo field in CMIP5 archives. Specified as a spatio-temporally variable quantity (in m equivalent, not in Pa!). + + if len(args) == 0: + self.setdefaultparameters() + else: + raise Exception('constructor not supported') #}}} @@ -28,25 +33,13 @@ #}}} - def defaultoutputs(self, md): #{{{ - return [] - #}}} - - def initialize(self, md): #{{{ - if np.isnan(self.global_average_thermosteric_sea_level): - self.global_average_thermosteric_sea_level = np.array([0, 0]).reshape(-1, 1) - print(' no dsl.global_average_thermosteric_sea_level specified: transient values set at zero') - - if np.isnan(self.sea_surface_height_above_geoid): - self.sea_surface_height_above_geoid = np.array(np.zeros(md.mesh.numberofvertices)).reshape(-1, 1) - print(' no dsl.sea_surface_height_above_geoid specified: transient values set at zero') - - if np.isnan(self.sea_water_pressure_at_sea_floor): - self.sea_water_pressure_at_sea_floor = np.array(np.zeros(md.mesh.numberofvertices)).reshape(-1, 1) - print(' no dsl.sea_water_pressure_at_sea_floor specified: transient values set at zero') + def setdefaultparameters(self): #{{{ + self.global_average_thermosteric_sea_level = np.nan + self.sea_surface_height_above_geoid = np.nan + self.sea_water_pressure_at_sea_floor = np.nan #}}} def checkconsistency(self, md, solution, analyses): #{{{ # Early return - if ('SealevelriseAnalysis' not in analyses) or (solution == 'TransientSolution' and not md.transient.isslc): + if ('SealevelchangeAnalysis' not in analyses) or (solution == 'TransientSolution' and not md.transient.isslc): return md @@ -61,15 +54,30 @@ # }}} + def marshall(self, prefix, md, fid): #{{{ + yts = md.constants.yts + WriteData(fid, prefix, 'name', 'md.dsl.model', 'data', 1, 'format', 'Integer') + WriteData(fid, prefix, 'object', self, 'fieldname', 'global_average_thermosteric_sea_level', 'format', 'DoubleMat', 'mattype', 2, 'timeseries', 1, 'yts', yts) # mattype 2, because we are sending a GMSL value identical everywhere on each element. + WriteData(fid, prefix, 'object', self, 'fieldname', 'sea_surface_height_above_geoid', 'format', 'DoubleMat', 'mattype', 1, 'timeserieslength', md.mesh.numberofvertices + 1, 'yts', yts) # mattype 1 because we specify DSL at vertex locations. + WriteData(fid, prefix, 'object', self, 'fieldname', 'sea_water_pressure_at_sea_floor', 'format', 'DoubleMat', 'mattype', 1, 'timeserieslength', md.mesh.numberofvertices + 1, 'yts', yts) # mattype 1 because we specify bottom pressure at vertex locations. + # }}} + def extrude(self, md): #{{{ - self.sea_surface_height_above_geoid = project3d(md, 'vector', self.sea_surface_height_above_geoid, 'type', 'node') - self.sea_water_pressure_at_sea_floor = project3d(md, 'vector', self.sea_water_pressure_at_sea_floor, 'type', 'node') + self.sea_surface_height_above_geoid = project3d(md, 'vector', self.sea_surface_height_above_geoid, 'type', 'node', 'layer', 1) + self.sea_water_pressure_at_sea_floor = project3d(md, 'vector', self.sea_water_pressure_at_sea_floor, 'type', 'node', 'layer', 1) return self #}}} - def marshall(self, prefix, md, fid): #{{{ - yts = md.constants.yts - WriteData(fid, prefix, 'name', 'md.dsl.model', 'data', 1, 'format', 'Integer') - WriteData(fid, prefix, 'object', self, 'class', 'dsl', 'fieldname', 'global_average_thermosteric_sea_level', 'format', 'DoubleMat', 'mattype', 2, 'timeserieslength', 1+1, 'yts', md.constants.yts, 'scale', 1e-3/md.constants.yts) # mattype 2, because we are sending a GMSL value identical everywhere on each element. - WriteData(fid, prefix, 'object', self, 'class', 'dsl', 'fieldname', 'sea_water_pressure_at_sea_floor', 'format', 'DoubleMat', 'mattype', 1, 'timeserieslength', md.mesh.numberofvertices+1, 'yts', md.constants.yts, 'scale', 1e-3/md.constants.yts) # mattype 1 because we specify DSL at vertex locations. - WriteData(fid, prefix, 'object', self, 'class', 'dsl', 'fieldname', 'sea_surface_height_above_geoid', 'format', 'DoubleMat', 'mattype', 1, 'timeserieslength', md.mesh.numberofvertices+1, 'yts', md.constants.yts) # mattype 1 because we specify bottom pressure at vertex locations. - # }}} + + def initialize(self, md): #{{{ + if np.isnan(self.global_average_thermosteric_sea_level): + self.global_average_thermosteric_sea_level = np.array([0, 0]).reshape(-1, 1) + print(' no dsl.global_average_thermosteric_sea_level specified: transient values set at zero') + + if np.isnan(self.sea_surface_height_above_geoid): + self.sea_surface_height_above_geoid = np.append(np.zeros((md.mesh.numberofvertices, 1)), 0).reshape(-1, 1) + print(' no dsl.sea_surface_height_above_geoid specified: transient values set at zero') + + if np.isnan(self.sea_water_pressure_at_sea_floor): + self.sea_water_pressure_at_sea_floor = np.append(np.zeros((md.mesh.numberofvertices, 1)), 0).reshape(-1, 1) + print(' no dsl.sea_water_pressure_at_sea_floor specified: transient values set at zero') + #}}} Index: /issm/trunk-jpl/src/m/classes/fourierlove.m =================================================================== --- /issm/trunk-jpl/src/m/classes/fourierlove.m (revision 26300) +++ /issm/trunk-jpl/src/m/classes/fourierlove.m (revision 26301) @@ -6,23 +6,23 @@ classdef fourierlove properties (SetAccess=public) - nfreq = 0; - frequencies = 0; - sh_nmax = 0; - sh_nmin = 0; - g0 = 0; - r0 = 0; - mu0 = 0; - Gravitational_Constant = 0; - allow_layer_deletion = 0; - underflow_tol = 0; - integration_steps_per_layer= 0; - istemporal = 0; - n_temporal_iterations = 0; - time = 0; - love_kernels = 0; - forcing_type = 0; - inner_core_boundary = 0; - core_mantle_boundary = 0; - complex_computation = 0; + nfreq = 0; + frequencies = 0; + sh_nmax = 0; + sh_nmin = 0; + g0 = 0; + r0 = 0; + mu0 = 0; + Gravitational_Constant = 0; + allow_layer_deletion = 0; + underflow_tol = 0; + integration_steps_per_layer = 0; + istemporal = 0; + n_temporal_iterations = 0; + time = 0; + love_kernels = 0; + forcing_type = 0; + inner_core_boundary = 0; + core_mantle_boundary = 0; + complex_computation = 0; end @@ -34,6 +34,4 @@ end methods - function self = extrude(self,md) % {{{ - end % }}} function self = fourierlove(varargin) % {{{ switch nargin @@ -53,5 +51,5 @@ self.g0=9.81; % m/s^2; self.r0=6371*1e3; %m; - self.mu0=10^11; % Pa + self.mu0=1e11; % Pa self.Gravitational_Constant=6.67259e-11; % m^3 kg^-1 s^-2 self.allow_layer_deletion=1; @@ -68,22 +66,22 @@ end % }}} function disp(self) % {{{ - fielddisplay(self,'nfreq','number of frequencies sampled (default 1, elastic) [Hz]'); + fielddisplay(self,'nfreq','number of frequencies sampled (default: 1, elastic) [Hz]'); fielddisplay(self,'frequencies','frequencies sampled (convention defaults to 0 for the elastic case) [Hz]'); - fielddisplay(self,'sh_nmax','maximum spherical harmonic degree (default 256, .35 deg, or 40 km at equator)'); - fielddisplay(self,'sh_nmin','minimum spherical harmonic degree (default 1)'); - fielddisplay(self,'g0','adimensioning constant for gravity (default 10) [m/s^2]'); - fielddisplay(self,'r0','adimensioning constant for radius (default 6371*10^3) [m]'); - fielddisplay(self,'mu0','adimensioning constant for stress (default 10^11) [Pa]'); - fielddisplay(self,'Gravitational_Constant','Newtonian constant of gravitation (default 6.67259e-11 [m^3 kg^-1 s^-2])'); - fielddisplay(self,'allow_layer_deletion','allow for migration of the integration boundary with increasing spherical harmonics degree (default 1)'); - fielddisplay(self,'underflow_tol','threshold of deep to surface love number ratio to trigger the deletion of layers (default 2.2204460492503131E-016)'); - fielddisplay(self,'integration_steps_per_layer','number of radial steps to propagate the yi system from the bottom to the top of each layer (default 100)'); - fielddisplay(self,'istemporal',{'1 for time-dependent love numbers, 0 for frequency-dependent or elastic love numbers (default 0)', 'If 1: use fourierlove function build_frequencies_from_time to meet consistency'}); - fielddisplay(self,'n_temporal_iterations','max number of iterations in the inverse Laplace transform. Also the number of spectral samples per time step requested (default 8)'); - fielddisplay(self,'time','time vector for deformation if istemporal (default 0) [s]'); - fielddisplay(self,'love_kernels','compute love numbers at depth? (default 0)'); - fielddisplay(self,'forcing_type',{'integer indicating the nature and depth of the forcing for the Love number calculation (default 11) :','1: Inner core boundary -- Volumic Potential','2: Inner core boundary -- Pressure','3: Inner core boundary -- Loading','4: Inner core boundary -- Tangential traction','5: Core mantle boundary -- Volumic Potential','6: Core mantle boundary -- Pressure','7: Core mantle boundary -- Loading','8: Core mantle boundary -- Tangential traction','9: Surface -- Volumic Potential','10: Surface -- Pressure','11: Surface -- Loading','12: Surface -- Tangential traction '}); - fielddisplay(self,'inner_core_boundary','interface index in materials.radius locating forcing. Only used for forcing_type 1--4 (default 1)'); - fielddisplay(self,'core_mantle_boundary','interface index in materials.radius locating forcing. Only used for forcing_type 5--8 (default 2)'); + fielddisplay(self,'sh_nmax','maximum spherical harmonic degree (default: 256, .35 deg, or 40 km at equator)'); + fielddisplay(self,'sh_nmin','minimum spherical harmonic degree (default: 1)'); + fielddisplay(self,'g0','adimensioning constant for gravity (default: 10) [m/s^2]'); + fielddisplay(self,'r0','adimensioning constant for radius (default: 6371*10^3) [m]'); + fielddisplay(self,'mu0','adimensioning constant for stress (default: 10^11) [Pa]'); + fielddisplay(self,'Gravitational_Constant','Newtonian constant of gravitation (default: 6.67259e-11 [m^3 kg^-1 s^-2])'); + fielddisplay(self,'allow_layer_deletion','allow for migration of the integration boundary with increasing spherical harmonics degree (default: 1)'); + fielddisplay(self,'underflow_tol','threshold of deep to surface love number ratio to trigger the deletion of layers (default: 1e-16)'); + fielddisplay(self,'integration_steps_per_layer','number of radial steps to propagate the yi system from the bottom to the top of each layer (default: 100)'); + fielddisplay(self,'istemporal',{'1 for time-dependent love numbers, 0 for frequency-dependent or elastic love numbers (default: 0)', 'If 1: use fourierlove function build_frequencies_from_time to meet consistency'}); + fielddisplay(self,'n_temporal_iterations','max number of iterations in the inverse Laplace transform. Also the number of spectral samples per time step requested (default: 8)'); + fielddisplay(self,'time','time vector for deformation if istemporal (default: 0) [s]'); + fielddisplay(self,'love_kernels','compute love numbers at depth? (default: 0)'); + fielddisplay(self,'forcing_type',{'integer indicating the nature and depth of the forcing for the Love number calculation (default: 11):','1: Inner core boundary -- Volumic Potential','2: Inner core boundary -- Pressure','3: Inner core boundary -- Loading','4: Inner core boundary -- Tangential traction','5: Core mantle boundary -- Volumic Potential','6: Core mantle boundary -- Pressure','7: Core mantle boundary -- Loading','8: Core mantle boundary -- Tangential traction','9: Surface -- Volumic Potential','10: Surface -- Pressure','11: Surface -- Loading','12: Surface -- Tangential traction '}); + fielddisplay(self,'inner_core_boundary','interface index in materials.radius locating forcing. Only used for forcing_type 1--4 (default: 1)'); + fielddisplay(self,'core_mantle_boundary','interface index in materials.radius locating forcing. Only used for forcing_type 5--8 (default: 2)'); end % }}} @@ -112,6 +110,6 @@ md = checkfield(md,'fieldname','love.time','NaN',1,'Inf',1,'numel',md.love.nfreq/2/md.love.n_temporal_iterations); end - if md.love.sh_nmin<=1 & (md.love.forcing_type==9 || md.love.forcing_type==5 || md.love.forcing_type==1) - error('Degree 1 not supported for Volumetric Potential forcing. Use sh_min>=2 for this kind of calculation.') + if md.love.sh_nmin<=1 & (md.love.forcing_type==1 || md.love.forcing_type==5 || md.love.forcing_type==9) + error(['Degree 1 not supported for forcing type ' num2str(md.love.forcing_type) '. Use sh_min>=2 for this kind of calculation.']) end @@ -130,5 +128,5 @@ end % }}} function marshall(self,prefix,md,fid) % {{{ - + WriteData(fid,prefix,'object',self,'fieldname','nfreq','format','Integer'); WriteData(fid,prefix,'object',self,'fieldname','frequencies','format','DoubleMat','mattype',3); @@ -152,4 +150,6 @@ end % }}} + function self = extrude(self,md) % {{{ + end % }}} function savemodeljs(self,fid,modelname) % {{{ error('not implemented yet!'); Index: /issm/trunk-jpl/src/m/classes/fourierlove.py =================================================================== --- /issm/trunk-jpl/src/m/classes/fourierlove.py (revision 26300) +++ /issm/trunk-jpl/src/m/classes/fourierlove.py (revision 26301) @@ -1,2 +1,4 @@ +import numpy as np + from fielddisplay import fielddisplay from checkfield import checkfield @@ -5,10 +7,11 @@ class fourierlove(object): - """FOURIERLOVE - Fourier Love Number class definition + """FOURIERLOVE - class definition Usage: - fourierlove = fourierlove() + md.love = fourierlove() """ - def __init__(self): # {{{ + + def __init__(self): #{{{ self.nfreq = 0 self.frequencies = 0 @@ -18,88 +21,123 @@ self.r0 = 0 self.mu0 = 0 + self.Gravitational_Constant = 0 self.allow_layer_deletion = 0 + self.underflow_tol = 0 + self.integration_steps_per_layer = 0 + self.istemporal = 0 + self.n_temporal_iterations = 0 + self.time = 0 self.love_kernels = 0 self.forcing_type = 0 + self.inner_core_boundary = 0 + self.core_mantle_boundary = 0 + self.complex_computation = 0 self.setdefaultparameters() #}}} - def __repr__(self): # {{{ - # TODO: - # - Convert all formatting to calls to .format (see any - # already converted .__repr__ method for examples) - # - string = ' Fourier Love class:' - string = "%s\n%s" % (string, fielddisplay(self, 'nfreq', 'number of frequencies sampled (default 1, elastic) [Hz]')) - string = "%s\n%s" % (string, fielddisplay(self, 'frequencies', 'frequencies sampled (convention defaults to 0 for the elastic case) [Hz]')) - string = "%s\n%s" % (string, fielddisplay(self, 'sh_nmax', 'maximum spherical harmonic degree (default 256, .35 deg, or 40 km at equator)')) - string = "%s\n%s" % (string, fielddisplay(self, 'sh_nmin', 'minimum spherical harmonic degree (default 1)')) - string = "%s\n%s" % (string, fielddisplay(self, 'g0', 'adimensioning constant for gravity (default 10) [m / s^2]')) - string = "%s\n%s" % (string, fielddisplay(self, 'r0', 'adimensioning constant for radius (default 6378e3) [m]')) - string = "%s\n%s" % (string, fielddisplay(self, 'mu0', 'adimensioning constant for stress (default 1.0e11) [Pa]')) - string = "%s\n%s" % (string, fielddisplay(self, 'allow_layer_deletion', 'allow for migration of the integration boundary with increasing spherical harmonics degree (default 1)')) - string = "%s\n%s" % (string, fielddisplay(self, 'love_kernels', 'compute love numbers at depth? (default 0)')) - string = "%s\n%s" % (string, fielddisplay(self, 'forcing_type', 'integer indicating the nature and depth of the forcing for the Love number calculation (default 11) :')) - string = "%s\n%s" % (string, ' 1: Inner core boundary -- Volumic Potential') - string = "%s\n%s" % (string, ' 2: Inner core boundary -- Pressure') - string = "%s\n%s" % (string, ' 3: Inner core boundary -- Loading') - string = "%s\n%s" % (string, ' 4: Inner core boundary -- Tangential traction') - string = "%s\n%s" % (string, ' 5: Core mantle boundary -- Volumic Potential') - string = "%s\n%s" % (string, ' 6: Core mantle boundary -- Pressure') - string = "%s\n%s" % (string, ' 7: Core mantle boundary -- Loading') - string = "%s\n%s" % (string, ' 8: Core mantle boundary -- Tangential traction') - string = "%s\n%s" % (string, ' 9: Surface-- Volumic Potential') - string = "%s\n%s" % (string, ' 10: Surface-- Pressure') - string = "%s\n%s" % (string, ' 11: Surface-- Loading') - string = "%s\n%s" % (string, ' 12: Surface-- Tangential traction ') + def __repr__(self): #{{{ + s = ' Fourier Love class:\n' + s += '{}\n'.format(fielddisplay(self, 'nfreq', 'number of frequencies sampled (default: 1, elastic) [Hz]')) + s += '{}\n'.format(fielddisplay(self, 'frequencies', 'frequencies sampled (convention defaults to 0 for the elastic case) [Hz]')) + s += '{}\n'.format(fielddisplay(self, 'sh_nmax', 'maximum spherical harmonic degree (default: 256, .35 deg, or 40 km at equator)')) + s += '{}\n'.format(fielddisplay(self, 'sh_nmin', 'minimum spherical harmonic degree (default: 1)')) + s += '{}\n'.format(fielddisplay(self, 'g0', 'adimensioning constant for gravity (default: 10) [m / s^2]')) + s += '{}\n'.format(fielddisplay(self, 'r0', 'adimensioning constant for radius (default: 6378e3) [m]')) + s += '{}\n'.format(fielddisplay(self, 'mu0', 'adimensioning constant for stress (default: 1.0e11) [Pa]')) + s += '{}\n'.format(fielddisplay(self, 'allow_layer_deletion', 'allow for migration of the integration boundary with increasing spherical harmonics degree (default: 1)')) + s += '{}\n'.format(fielddisplay(self, 'Gravitational_Constant', 'Newtonian constant of gravitation (default: 6.67259e-11 [m^3 kg^-1 s^-2])')) + s += '{}\n'.format(fielddisplay(self, 'allow_layer_deletion', 'allow for migration of the integration boundary with increasing spherical harmonics degree (default: 1)')) + s += '{}\n'.format(fielddisplay(self, 'underflow_tol', 'threshold of deep to surface love number ratio to trigger the deletion of layers (default: 1e-16)')) + s += '{}\n'.format(fielddisplay(self, 'integration_steps_per_layer', 'number of radial steps to propagate the yi system from the bottom to the top of each layer (default: 100)')) + s += '{}\n'.format(fielddisplay(self, 'istemporal', {'1 for time-dependent love numbers, 0 for frequency-dependent or elastic love numbers (default: 0)', 'If 1: use fourierlove function build_frequencies_from_time to meet consistency'})) + s += '{}\n'.format(fielddisplay(self, 'n_temporal_iterations', 'max number of iterations in the inverse Laplace transform. Also the number of spectral samples per time step requested (default: 8)')) + s += '{}\n'.format(fielddisplay(self, 'time', 'time vector for deformation if istemporal (default: 0) [s]')) + s += '{}\n'.format(fielddisplay(self, 'love_kernels', 'compute love numbers at depth? (default: 0)')) + s += '{}\n'.format(fielddisplay(self, 'forcing_type', 'integer indicating the nature and depth of the forcing for the Love number calculation (default: 11):')) + s += '{}\n'.format(' 1: Inner core boundary -- Volumic Potential') + s += '{}\n'.format(' 2: Inner core boundary -- Pressure') + s += '{}\n'.format(' 3: Inner core boundary -- Loading') + s += '{}\n'.format(' 4: Inner core boundary -- Tangential traction') + s += '{}\n'.format(' 5: Core mantle boundary -- Volumic Potential') + s += '{}\n'.format(' 6: Core mantle boundary -- Pressure') + s += '{}\n'.format(' 7: Core mantle boundary -- Loading') + s += '{}\n'.format(' 8: Core mantle boundary -- Tangential traction') + s += '{}\n'.format(' 9: Surface-- Volumic Potential') + s += '{}\n'.format(' 10: Surface-- Pressure') + s += '{}\n'.format(' 11: Surface-- Loading') + s += '{}\n'.format(' 12: Surface-- Tangential traction ') + s += '{}\n'.format(fielddisplay(self, 'inner_core_boundary', 'interface index in materials.radius locating forcing. Only used for forcing_type 1--4 (default: 1)')) + s += '{}\n'.format(fielddisplay(self, 'core_mantle_boundary', 'interface index in materials.radius locating forcing. Only used for forcing_type 5--8 (default: 2)')) - return string + return s #}}} - def extrude(self, md): # {{{ - return self + def setdefaultparameters(self): #{{{ + # We setup an elastic love number computation by default + self.nfreq = 1 + self.frequencies = [0] # Hz + self.sh_nmax = 256 # .35 degree, 40 km at the equator + self.sh_nmin = 1 + # Work on matlab script for computing g0 for given Earth's structure + self.g0 = 9.81 # m/s^2 + self.r0 = 6371 * 1e3 # m + self.mu0 = 1e11 # Pa + self.Gravitational_Constant = 6.67259e-11 # m^3 kg^-1 s^-2 + self.allow_layer_deletion = 1 + self.underflow_tol = 1e-16 # Threshold of deep to surface love number ratio to trigger the deletion of layer + self.integration_steps_per_layer = 100 + self.istemporal = 0 + self.n_temporal_iterations = 8 + self.time = [0] # s + self.love_kernels = 0 + self.forcing_type = 11 # Surface loading + self.inner_core_boundary = 1 + self.core_mantle_boundary = 2 + self.complex_computation = 0 #}}} - def setdefaultparameters(self): # {{{ - #we setup an elastic love number computation by default. - self.nfreq = 1 - self.frequencies = [0] #Hz - self.sh_nmax = 256 # .35 degree, 40 km at the equator. - self.sh_nmin = 1 - self.g0 = 9.81 # m/s^2 - self.r0 = 6371e3 #m - self.mu0 = 1e11 # Pa - self.allow_layer_deletion = 1 - self.love_kernels = 0 - self.forcing_type = 11 + def checkconsistency(self, md, solution, analyses): #{{{ + if 'LoveAnalysis' not in analyses: + return md - return self + md = checkfield(md, 'fieldname', 'love.nfreq', 'NaN', 1, 'Inf', 1, 'numel', 1, '>', 0) + md = checkfield(md, 'fieldname', 'love.frequencies', 'NaN', 1, 'Inf', 1, 'numel', md.love.nfreq) + md = checkfield(md, 'fieldname', 'love.sh_nmax', 'NaN', 1, 'Inf', 1, 'numel', 1, '>', 0) + md = checkfield(md, 'fieldname', 'love.sh_nmin', 'NaN', 1, 'Inf', 1, 'numel', 1, '>', 0) + md = checkfield(md, 'fieldname', 'love.g0', 'NaN', 1, 'Inf', 1, 'numel', 1, '>', 0) + md = checkfield(md, 'fieldname', 'love.r0', 'NaN', 1, 'Inf', 1, 'numel', 1, '>', 0) + md = checkfield(md, 'fieldname', 'love.mu0', 'NaN', 1, 'Inf', 1, 'numel', 1, '>', 0) + md = checkfield(md, 'fieldname', 'love.Gravitational_Constant', 'NaN', 1, 'Inf', 1, 'numel', 1, '>', 0) + md = checkfield(md, 'fieldname', 'love.allow_layer_deletion', 'values', [0, 1]) + md = checkfield(md, 'fieldname', 'love.underflow_tol', 'NaN', 1, 'Inf', 1, 'numel', 1, '>', 0) + md = checkfield(md, 'fieldname', 'love.integration_steps_per_layer', 'NaN', 1, 'Inf', 1, 'numel', 1, '>', 0) + md = checkfield(md, 'fieldname', 'love.love_kernels', 'values', [0, 1]) + md = checkfield(md, 'fieldname', 'love.forcing_type', 'NaN', 1, 'Inf', 1, 'numel', 1, '>', 0, '<=', 12) + md = checkfield(md, 'fieldname', 'love.complex_computation', 'NaN', 1, 'Inf', 1, 'numel', 1, 'values', [0, 1]) + + md = checkfield(md, 'fieldname', 'love.istemporal', 'values', [0, 1]) + if md.love.istemporal: + md = checkfield(md, 'fieldname', 'love.n_temporal_iterations', 'NaN', 1, 'Inf', 1, 'numel', 1, '>', 0) + md = checkfield(md, 'fieldname', 'love.time', 'NaN', 1, 'Inf', 1, 'numel', md.love.nfreq / 2 / md.love.n_temporal_iterations) + if md.love.sh_nmin <= 1 and (md.love.forcing_type == 1 or md.love.forcing_type == 5 or md.love.forcing_type == 9): + raise RuntimeError('Degree 1 not supported for forcing type {}. Use sh_min >= 2 for this kind of calculation.'.format(md.love.forcing_type)) + + # Need 'litho' material + if not hasattr(md.materials, 'materials') or 'litho' not in md.materials.nature: + raise RuntimeError('Need a \'litho\' material to run a Fourier Love number analysis') + + mat = np.where(np.array(nature) == 'litho')[0] + if md.love.forcing_type <= 4: + md = checkfield(md, 'fieldname', 'love.inner_core_boundary', 'NaN', 1, 'Inf', 1, 'numel', 1, '>', 0, '<=', md.materials[mat].numlayers) + elif md.love.forcing_type <= 8: + md = checkfield(md, 'fieldname', 'love.core_mantle_boundary', 'NaN', 1, 'Inf', 1, 'numel', 1, '>', 0, '<=', md.materials[mat].numlayers) + + return md #}}} - def checkconsistency(self, md, solution, analyses): # {{{ - if 'LoveAnalysis' not in analyses: - return md - md = checkfield(md, 'fieldname', 'love.nfreq', 'NaN', 1, 'Inf', 1, 'numel', [1], '>', 0) - md = checkfield(md, 'fieldname', 'love.frequencies', 'NaN', 1, 'Inf', 1, 'numel', [md.love.nfreq]) - md = checkfield(md, 'fieldname', 'love.sh_nmax', 'NaN', 1, 'Inf', 1, 'numel', [1], '>', 0) - md = checkfield(md, 'fieldname', 'love.sh_nmin', 'NaN', 1, 'Inf', 1, 'numel', [1], '>', 0) - md = checkfield(md, 'fieldname', 'love.g0', 'NaN', 1, 'Inf', 1, 'numel', [1], '>', 0) - md = checkfield(md, 'fieldname', 'love.r0', 'NaN', 1, 'Inf', 1, 'numel', [1], '>', 0) - md = checkfield(md, 'fieldname', 'love.mu0', 'NaN', 1, 'Inf', 1, 'numel', [1], '>', 0) - md = checkfield(md, 'fieldname', 'love.allow_layer_deletion', 'values', [0, 1]) - md = checkfield(md, 'fieldname', 'love.love_kernels', 'values', [0, 1]) - md = checkfield(md, 'fieldname', 'love.forcing_type', 'NaN', 1, 'Inf', 1, 'numel', [1], '>', 0, '<=', 12) - if md.love.sh_nmin <= 1 and md.love.forcing_type == 9: - raise RuntimeError("Degree 1 not supported for Volumetric Potential forcing. Use sh_min >= 2 for this kind of calculation.") - - # need 'litho' material - if not hasattr(md.materials, 'materials') or 'litho' not in md.materials.nature: - raise RuntimeError('Need a \'litho\' material to run a Fourier Love number analysis') - return md - # }}} - - def marshall(self, prefix, md, fid): # {{{ + def marshall(self, prefix, md, fid): #{{{ WriteData(fid, prefix, 'object', self, 'fieldname', 'nfreq', 'format', 'Integer') - WriteData(fid, prefix, 'object', self, 'fieldname', 'frequencies', 'format', 'DoubleMat', 'mattype', 3) + WriteData(fid, prefix, 'object', self, 'fieldname', 'frequencies', 'format', 'DoubleMat', 'mattype',3) WriteData(fid, prefix, 'object', self, 'fieldname', 'sh_nmax', 'format', 'Integer') WriteData(fid, prefix, 'object', self, 'fieldname', 'sh_nmin', 'format', 'Integer') @@ -107,6 +145,30 @@ WriteData(fid, prefix, 'object', self, 'fieldname', 'r0', 'format', 'Double') WriteData(fid, prefix, 'object', self, 'fieldname', 'mu0', 'format', 'Double') + WriteData(fid, prefix, 'object', self, 'fieldname', 'Gravitational_Constant', 'format', 'Double') WriteData(fid, prefix, 'object', self, 'fieldname', 'allow_layer_deletion', 'format', 'Boolean') + WriteData(fid, prefix, 'object', self, 'fieldname', 'underflow_tol', 'format', 'Double') + WriteData(fid, prefix, 'object', self, 'fieldname', 'integration_steps_per_layer', 'format', 'Integer') + WriteData(fid, prefix, 'object', self, 'fieldname', 'istemporal', 'format', 'Boolean') + WriteData(fid, prefix, 'object', self, 'fieldname', 'n_temporal_iterations', 'format', 'Integer') + WriteData(fid, prefix, 'object', self, 'fieldname', 'complex_computation', 'format', 'Boolean') + # Note: no need to marshall the time vector, we have frequencies WriteData(fid, prefix, 'object', self, 'fieldname', 'love_kernels', 'format', 'Boolean') WriteData(fid, prefix, 'object', self, 'fieldname', 'forcing_type', 'format', 'Integer') - # }}} + WriteData(fid, prefix, 'object', self, 'fieldname', 'inner_core_boundary', 'format', 'Integer') + WriteData(fid, prefix, 'object', self, 'fieldname', 'core_mantle_boundary', 'format', 'Integer') + #}}} + + def extrude(self, md): #{{{ + return self + #}}} + + def build_frequencies_from_time(self): #{{{ + if not self.istemporal: + raise RuntimeError('cannot build frequencies for temporal love numbers if love.istemporal==0') + print('Temporal love numbers: Overriding md.love.nfreq and md.love.frequencies') + self.nfreq = len(self.time) * 2 * self.n_temporal_iterations + self.frequencies = np.zeros((self.nfreq, 1)) + for i in range(len(self.time)): + for j in range(2 * self.n_temporal_iterations): + self.frequencies[(i - 1) * 2 * self.n_temporal_iterations + j] = j * np.log(2) / self.time[i] / 2 / np.pi + #}}} Index: /issm/trunk-jpl/src/m/classes/geometry.py =================================================================== --- /issm/trunk-jpl/src/m/classes/geometry.py (revision 26300) +++ /issm/trunk-jpl/src/m/classes/geometry.py (revision 26301) @@ -37,5 +37,5 @@ def setdefaultparameters(self): #{{{ - return self + return #}}} Index: /issm/trunk-jpl/src/m/classes/hydrologyshreve.m =================================================================== --- /issm/trunk-jpl/src/m/classes/hydrologyshreve.m (revision 26300) +++ /issm/trunk-jpl/src/m/classes/hydrologyshreve.m (revision 26301) @@ -11,6 +11,4 @@ end methods - function self = extrude(self,md) % {{{ - end % }}} function self = hydrologyshreve(varargin) % {{{ switch nargin @@ -23,8 +21,11 @@ end end % }}} - function list = defaultoutputs(self,md) % {{{ - list = {'Watercolumn','HydrologyWaterVx','HydrologyWaterVy'}; - end % }}} + function disp(self) % {{{ + disp(sprintf(' hydrologyshreve solution parameters:')); + fielddisplay(self,'spcwatercolumn','water thickness constraints (NaN means no constraint) [m]'); + fielddisplay(self,'stabilization','artificial diffusivity (default: 1). can be more than 1 to increase diffusivity.'); + fielddisplay(self,'requested_outputs','additional outputs requested'); + end % }}} function self = setdefaultparameters(self) % {{{ @@ -34,5 +35,5 @@ end % }}} function md = checkconsistency(self,md,solution,analyses) % {{{ - + %Early return if ~ismember('HydrologyShreveAnalysis',analyses) | (strcmp(solution,'TransientSolution') & md.transient.ishydrology==0), return; end @@ -41,10 +42,6 @@ md = checkfield(md,'fieldname','hydrology.stabilization','>=',0); end % }}} - function disp(self) % {{{ - disp(sprintf(' hydrologyshreve solution parameters:')); - fielddisplay(self,'spcwatercolumn','water thickness constraints (NaN means no constraint) [m]'); - fielddisplay(self,'stabilization','artificial diffusivity (default is 1). can be more than 1 to increase diffusivity.'); - fielddisplay(self,'requested_outputs','additional outputs requested'); - + function list = defaultoutputs(self,md) % {{{ + list = {'Watercolumn','HydrologyWaterVx','HydrologyWaterVy'}; end % }}} function marshall(self,prefix,md,fid) % {{{ @@ -52,11 +49,13 @@ WriteData(fid,prefix,'object',self,'fieldname','spcwatercolumn','format','DoubleMat','mattype',1,'timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts); WriteData(fid,prefix,'object',self,'fieldname','stabilization','format','Double'); - outputs = self.requested_outputs; - pos = find(ismember(outputs,'default')); - if ~isempty(pos), - outputs(pos) = []; %remove 'default' from outputs - outputs = [outputs defaultoutputs(self,md)]; %add defaults - end - WriteData(fid,prefix,'data',outputs,'name','md.hydrology.requested_outputs','format','StringArray'); + outputs = self.requested_outputs; + pos = find(ismember(outputs,'default')); + if ~isempty(pos), + outputs(pos) = []; %remove 'default' from outputs + outputs = [outputs defaultoutputs(self,md)]; %add defaults + end + WriteData(fid,prefix,'data',outputs,'name','md.hydrology.requested_outputs','format','StringArray'); + end % }}} + function self = extrude(self,md) % {{{ end % }}} function savemodeljs(self,fid,modelname) % {{{ Index: /issm/trunk-jpl/src/m/classes/hydrologyshreve.py =================================================================== --- /issm/trunk-jpl/src/m/classes/hydrologyshreve.py (revision 26300) +++ /issm/trunk-jpl/src/m/classes/hydrologyshreve.py (revision 26301) @@ -1,2 +1,4 @@ +import numpy as np + from fielddisplay import fielddisplay from checkfield import checkfield @@ -11,41 +13,33 @@ """ - def __init__(self): # {{{ - self.spcwatercolumn = float('NaN') + def __init__(self, *args): #{{{ + self.spcwatercolumn = np.nan self.stabilization = 0 self.requested_outputs = [] - self.setdefaultparameters() + nargs = len(args) + if nargs == 0: + self.setdefaultparameters() + elif nargs == 1: + self.setdefaultparameters(args) + else: + raise RuntimeError('constructor not supported') #}}} - def __repr__(self): # {{{ - # TODO: - # - Convert all formatting to calls to .format (see any - # already converted .__repr__ method for examples) - # - string = ' hydrologyshreve solution parameters:' - string = "%s\n%s" % (string, fielddisplay(self, 'spcwatercolumn', 'water thickness constraints (NaN means no constraint) [m]')) - string = "%s\n%s" % (string, fielddisplay(self, 'stabilization', 'artificial diffusivity (default is 1). can be more than 1 to increase diffusivity.')) - string = "%s\n%s" % (string, fielddisplay(self, 'requested_outputs', 'additional outputs requested')) - return string + def __repr__(self): #{{{ + s = ' hydrologyshreve solution parameters:\n' + s += '{}\n'.format(fielddisplay(self, 'spcwatercolumn', 'water thickness constraints (NaN means no constraint) [m]')) + s += '{}\n'.format(fielddisplay(self, 'stabilization', 'artificial diffusivity (default: 1). can be more than 1 to increase diffusivity.')) + s += '{}\n'.format(fielddisplay(self, 'requested_outputs', 'additional outputs requested')) + return s #}}} - def extrude(self, md): # {{{ - return self + def setdefaultparameters(self): #{{{ + # Type of stabilization to use 0:nothing 1:artificial_diffusivity + self.stabilization = 1 + self.requested_outputs = ['default'] #}}} - def setdefaultparameters(self): # {{{ - #Type of stabilization to use 0:nothing 1:artificial_diffusivity - self.stabilization = 1 - self.requested_outputs = ['default'] - return self - #}}} - - def defaultoutputs(self, md): # {{{ - list = ['Watercolumn', 'HydrologyWaterVx', 'HydrologyWaterVy'] - return list - #}}} - - def checkconsistency(self, md, solution, analyses): # {{{ + def checkconsistency(self, md, solution, analyses): #{{{ #Early return if 'HydrologyShreveAnalysis' not in analyses or (solution == 'TransientSolution' and not md.transient.ishydrology): @@ -54,19 +48,24 @@ md = checkfield(md, 'fieldname', 'hydrology.spcwatercolumn', 'Inf', 1, 'timeseries', 1) md = checkfield(md, 'fieldname', 'hydrology.stabilization', '>=', 0) - md = checkfield(md, 'fieldname', 'hydrology.requested_outputs', 'stringrow', 1) return md # }}} - def marshall(self, prefix, md, fid): # {{{ + def defaultoutputs(self, md): #{{{ + return ['Watercolumn', 'HydrologyWaterVx', 'HydrologyWaterVy'] + #}}} + + def marshall(self, prefix, md, fid): #{{{ WriteData(fid, prefix, 'name', 'md.hydrology.model', 'data', 2, 'format', 'Integer') WriteData(fid, prefix, 'object', self, 'fieldname', 'spcwatercolumn', 'format', 'DoubleMat', 'mattype', 1, 'timeserieslength', md.mesh.numberofvertices + 1, 'yts', md.constants.yts) WriteData(fid, prefix, 'object', self, 'fieldname', 'stabilization', 'format', 'Double') - #process requested outputs outputs = self.requested_outputs indices = [i for i, x in enumerate(outputs) if x == 'default'] - if len(indices) > 0: + if len(indices): outputscopy = outputs[0:max(0, indices[0] - 1)] + self.defaultoutputs(md) + outputs[indices[0] + 1:] outputs = outputscopy WriteData(fid, prefix, 'data', outputs, 'name', 'md.hydrology.requested_outputs', 'format', 'StringArray') + # }}} - # }}} + def extrude(self, md): #{{{ + return self + #}}} Index: /issm/trunk-jpl/src/m/classes/initialization.m =================================================================== --- /issm/trunk-jpl/src/m/classes/initialization.m (revision 26300) +++ /issm/trunk-jpl/src/m/classes/initialization.m (revision 26301) @@ -27,24 +27,4 @@ end methods - function self = extrude(self,md) % {{{ - self.vx=project3d(md,'vector',self.vx,'type','node'); - self.vy=project3d(md,'vector',self.vy,'type','node'); - self.vz=project3d(md,'vector',self.vz,'type','node'); - self.vel=project3d(md,'vector',self.vel,'type','node'); - self.temperature=project3d(md,'vector',self.temperature,'type','node'); - self.enthalpy=project3d(md,'vector',self.enthalpy,'type','node'); - self.waterfraction=project3d(md,'vector',self.waterfraction,'type','node'); - self.watercolumn=project3d(md,'vector',self.watercolumn,'type','node','layer',1); - self.sediment_head=project3d(md,'vector',self.sediment_head,'type','node','layer',1); - self.epl_head=project3d(md,'vector',self.epl_head,'type','node','layer',1); - self.epl_thickness=project3d(md,'vector',self.epl_thickness,'type','node','layer',1); - self.sealevel=project3d(md,'vector',self.sealevel,'type','node','layer',1); - self.bottompressure=project3d(md,'vector',self.bottompressure,'type','node','layer',1); - self.dsl=project3d(md,'vector',self.dsl,'type','node','layer',1); - self.str=project3d(md,'vector',self.str,'type','node','layer',1); - - %Lithostatic pressure by default - self.pressure=md.constants.g*md.materials.rho_ice*(md.geometry.surface-md.mesh.z); - end % }}} function self = initialization(varargin) % {{{ switch nargin @@ -201,4 +181,24 @@ end end % }}} + function self = extrude(self,md) % {{{ + self.vx=project3d(md,'vector',self.vx,'type','node'); + self.vy=project3d(md,'vector',self.vy,'type','node'); + self.vz=project3d(md,'vector',self.vz,'type','node'); + self.vel=project3d(md,'vector',self.vel,'type','node'); + self.temperature=project3d(md,'vector',self.temperature,'type','node'); + self.enthalpy=project3d(md,'vector',self.enthalpy,'type','node'); + self.waterfraction=project3d(md,'vector',self.waterfraction,'type','node'); + self.watercolumn=project3d(md,'vector',self.watercolumn,'type','node','layer',1); + self.sediment_head=project3d(md,'vector',self.sediment_head,'type','node','layer',1); + self.epl_head=project3d(md,'vector',self.epl_head,'type','node','layer',1); + self.epl_thickness=project3d(md,'vector',self.epl_thickness,'type','node','layer',1); + self.sealevel=project3d(md,'vector',self.sealevel,'type','node','layer',1); + self.bottompressure=project3d(md,'vector',self.bottompressure,'type','node','layer',1); + self.dsl=project3d(md,'vector',self.dsl,'type','node','layer',1); + self.str=project3d(md,'vector',self.str,'type','node','layer',1); + + %Lithostatic pressure by default + self.pressure=md.constants.g*md.materials.rho_ice*(md.geometry.surface-md.mesh.z); + end % }}} function savemodeljs(self,fid,modelname) % {{{ Index: /issm/trunk-jpl/src/m/classes/initialization.py =================================================================== --- /issm/trunk-jpl/src/m/classes/initialization.py (revision 26300) +++ /issm/trunk-jpl/src/m/classes/initialization.py (revision 26301) @@ -11,30 +11,31 @@ Usage: - initialization = initialization() + initialization = initialization() """ - def __init__(self): # {{{ + def __init__(self): #{{{ self.vx = np.nan self.vy = np.nan self.vz = np.nan self.vel = np.nan - self.enthalpy = np.nan self.pressure = np.nan self.temperature = np.nan + self.enthalpy = np.nan self.waterfraction = np.nan - self.watercolumn = np.nan self.sediment_head = np.nan self.epl_head = np.nan self.epl_thickness = np.nan + self.watercolumn = np.nan self.hydraulic_potential = np.nan self.channelarea = np.nan self.sealevel = np.nan self.bottompressure = np.nan + self.dsl = np.nan + self.str = np.nan self.sample = np.nan - #set defaults self.setdefaultparameters() #}}} - def __repr__(self): # {{{ + def __repr__(self): #{{{ s = ' initial field values:\n' s += '{}\n'.format(fielddisplay(self, 'vx', 'x component of velocity [m / yr]')) @@ -52,35 +53,11 @@ s += '{}\n'.format(fielddisplay(self, 'hydraulic_potential', 'Hydraulic potential (for GlaDS) [Pa]')) s += '{}\n'.format(fielddisplay(self, 'channelarea', 'subglaciale water channel area (for GlaDS) [m2]')) - s += '{}\n'.format(fielddisplay(self,'sample','Realization of a Gaussian random field')) + s += '{}\n'.format(fielddisplay(self,'sample', 'Realization of a Gaussian random field')) return s #}}} - def extrude(self, md): # {{{ - self.vx = project3d(md, 'vector', self.vx, 'type', 'node') - self.vy = project3d(md, 'vector', self.vy, 'type', 'node') - self.vz = project3d(md, 'vector', self.vz, 'type', 'node') - self.vel = project3d(md, 'vector', self.vel, 'type', 'node') - self.temperature = project3d(md, 'vector', self.temperature, 'type', 'node') - self.enthalpy = project3d(md, 'vector', self.enthalpy, 'type', 'node') - self.waterfraction = project3d(md, 'vector', self.waterfraction, 'type', 'node') - self.watercolumn = project3d(md, 'vector', self.watercolumn, 'type', 'node') - self.sediment_head = project3d(md, 'vector', self.sediment_head, 'type', 'node', 'layer', 1) - self.epl_head = project3d(md, 'vector', self.epl_head, 'type', 'node', 'layer', 1) - self.epl_thickness = project3d(md, 'vector', self.epl_thickness, 'type', 'node', 'layer', 1) - self.sealevel = project3d(md, 'vector', self.sealevel, 'type', 'node', 'layer', 1) - self.bottompressure = project3d(md, 'vector', self.bottompressure, 'type', 'node', 'layer', 1) - - #Lithostatic pressure by default - if np.ndim(md.geometry.surface) == 2: - print('Reshaping md.geometry.surface for your convenience but you should fix it in your model set up') - self.pressure = md.constants.g * md.materials.rho_ice * (md.geometry.surface.reshape(-1, ) - md.mesh.z) - else: - self.pressure = md.constants.g * md.materials.rho_ice * (md.geometry.surface - md.mesh.z) - - return self + def setdefaultparameters(self): #{{{ + return #}}} - def setdefaultparameters(self): # {{{ - return self - #}}} - def checkconsistency(self, md, solution, analyses): # {{{ + def checkconsistency(self, md, solution, analyses): #{{{ if 'StressbalanceAnalysis' in analyses and not solution == 'TransientSolution' and not md.transient.isstressbalance: if not np.any(np.logical_or(np.isnan(md.initialization.vx), np.isnan(md.initialization.vy))): @@ -130,15 +107,18 @@ if 'SamplingAnalysis' in analyses and not solution == 'TransientSolution' and not md.transient.issampling: if np.any(np.isnan(md.initialization.sample)): - md = checkfield(md,'fieldname','initialization.sample','NaN',1,'Inf',1,'size',[md.mesh.numberofvertices]) + md = checkfield(md, 'fieldname', 'initialization.sample', 'NaN', 1,'Inf', 1, 'size', [md.mesh.numberofvertices]) return md - # }}} - def marshall(self, prefix, md, fid): # {{{ + #}}} + def marshall(self, prefix, md, fid): #{{{ yts = md.constants.yts - WriteData(fid, prefix, 'object', self, 'fieldname', 'vx', 'format', 'DoubleMat', 'mattype', 1, 'scale', 1. / yts) - WriteData(fid, prefix, 'object', self, 'fieldname', 'vy', 'format', 'DoubleMat', 'mattype', 1, 'scale', 1. / yts) - WriteData(fid, prefix, 'object', self, 'fieldname', 'vz', 'format', 'DoubleMat', 'mattype', 1, 'scale', 1. / yts) + + WriteData(fid, prefix, 'object', self, 'fieldname', 'vx', 'format', 'DoubleMat', 'mattype', 1, 'scale', 1 / yts) + WriteData(fid, prefix, 'object', self, 'fieldname', 'vy', 'format', 'DoubleMat', 'mattype', 1, 'scale', 1 / yts) + WriteData(fid, prefix, 'object', self, 'fieldname', 'vz', 'format', 'DoubleMat', 'mattype', 1, 'scale', 1 / yts) WriteData(fid, prefix, 'object', self, 'fieldname', 'pressure', 'format', 'DoubleMat', 'mattype', 1) - WriteData(fid, prefix, 'object', self, 'fieldname', 'sealevel', 'format', 'DoubleMat', 'mattype', 1, 'timeserieslength', md.mesh.numberofvertices + 1, 'yts', md.constants.yts) + WriteData(fid, prefix, 'object', self, 'fieldname', 'sealevel', 'format', 'DoubleMat', 'mattype', 1,'timeserieslength', md.mesh.numberofvertices + 1, 'yts', yts) WriteData(fid, prefix, 'object', self, 'fieldname', 'bottompressure', 'format', 'DoubleMat', 'mattype', 1) + WriteData(fid, prefix, 'object', self, 'fieldname', 'str', 'format', 'DoubleMat', 'mattype', 1) + WriteData(fid, prefix, 'object', self, 'fieldname', 'dsl', 'format', 'DoubleMat', 'mattype', 1) WriteData(fid, prefix, 'object', self, 'fieldname', 'temperature', 'format', 'DoubleMat', 'mattype', 1) WriteData(fid, prefix, 'object', self, 'fieldname', 'waterfraction', 'format', 'DoubleMat', 'mattype', 1) @@ -149,12 +129,38 @@ WriteData(fid, prefix, 'object', self, 'fieldname', 'channelarea', 'format', 'DoubleMat', 'mattype', 1) WriteData(fid, prefix, 'object', self, 'fieldname', 'hydraulic_potential', 'format', 'DoubleMat', 'mattype', 1) - WriteData(fid,prefix,'object',self,'fieldname','sample','format','DoubleMat','mattype',1) + WriteData(fid, prefix, 'object', self, 'fieldname', 'sample', 'format', 'DoubleMat', 'mattype', 1) + if md.thermal.isenthalpy: if (np.size(self.enthalpy) <= 1): + # Reconstruct enthalpy tpmp = md.materials.meltingpoint - md.materials.beta * md.initialization.pressure - pos = np.nonzero(md.initialization.waterfraction > 0.)[0] + pos = np.where(md.initialization.waterfraction > 0)[0] self.enthalpy = md.materials.heatcapacity * (md.initialization.temperature - md.constants.referencetemperature) - self.enthalpy[pos] = md.materials.heatcapacity * (tpmp[pos].reshape(-1, ) - md.constants.referencetemperature) + md.materials.latentheat * md.initialization.waterfraction[pos].reshape(-1, ) + self.enthalpy[pos] = md.materials.heatcapacity * (tpmp[pos].reshape(-1, 1) - md.constants.referencetemperature) + md.materials.latentheat * md.initialization.waterfraction[pos].reshape(-1, 1) WriteData(fid, prefix, 'data', self.enthalpy, 'format', 'DoubleMat', 'mattype', 1, 'name', 'md.initialization.enthalpy') - # }}} + #}}} + def extrude(self, md): #{{{ + self.vx = project3d(md, 'vector', self.vx, 'type', 'node') + self.vy = project3d(md, 'vector', self.vy, 'type', 'node') + self.vz = project3d(md, 'vector', self.vz, 'type', 'node') + self.vel = project3d(md, 'vector', self.vel, 'type', 'node') + self.temperature = project3d(md, 'vector', self.temperature, 'type', 'node') + self.enthalpy = project3d(md, 'vector', self.enthalpy, 'type', 'node') + self.waterfraction = project3d(md, 'vector', self.waterfraction, 'type', 'node') + self.watercolumn = project3d(md, 'vector', self.watercolumn, 'type', 'node') + self.sediment_head = project3d(md, 'vector', self.sediment_head, 'type', 'node', 'layer', 1) + self.epl_head = project3d(md, 'vector', self.epl_head, 'type', 'node', 'layer', 1) + self.epl_thickness = project3d(md, 'vector', self.epl_thickness, 'type', 'node', 'layer', 1) + self.sealevel = project3d(md, 'vector', self.sealevel, 'type', 'node', 'layer', 1) + self.bottompressure = project3d(md, 'vector', self.bottompressure, 'type', 'node', 'layer', 1) + + # Lithostatic pressure by default + if np.ndim(md.geometry.surface) == 2: + print('Reshaping md.geometry.surface for your convenience but you should fix it in your model set up') + self.pressure = md.constants.g * md.materials.rho_ice * (md.geometry.surface.reshape(-1, 1) - md.mesh.z) + else: + self.pressure = md.constants.g * md.materials.rho_ice * (md.geometry.surface - md.mesh.z) + + return self + #}}} Index: /issm/trunk-jpl/src/m/classes/lovenumbers.m =================================================================== --- /issm/trunk-jpl/src/m/classes/lovenumbers.m (revision 26300) +++ /issm/trunk-jpl/src/m/classes/lovenumbers.m (revision 26301) @@ -16,12 +16,12 @@ %tidal love numbers for computing rotational feedback: - th = []; - tk = []; - tl = []; - tk2secular = 0; %deg 2 secular number. + th = []; + tk = []; + tl = []; + tk2secular = 0; %deg 2 secular number. %time/frequency for visco-elastic love numbers timefreq = []; - istime = 1; + istime = 1; end @@ -33,6 +33,22 @@ self=setdefaultparameters(self,maxdeg,referenceframe); end % }}} + function disp(self) % {{{ + disp(sprintf(' lovenumbers parameters:')); + + fielddisplay(self,'h','load Love number for radial displacement'); + fielddisplay(self,'k','load Love number for gravitational potential perturbation'); + fielddisplay(self,'l','load Love number for horizontal displacements'); + + fielddisplay(self,'th','tidal load Love number (deg 2)'); + fielddisplay(self,'tk','tidal load Love number (deg 2)'); + fielddisplay(self,'tl','tidal load Love number (deg 2)'); + fielddisplay(self,'tk2secular','secular fluid Love number'); + + fielddisplay(self,'istime','time (default: 1) or frequency love numbers (0)'); + fielddisplay(self,'timefreq','time/frequency vector (yr or 1/yr)'); + + end % }}} function self = setdefaultparameters(self,maxdeg,referenceframe) % {{{ - + %initialize love numbers: self.h=getlovenumbers('type','loadingverticaldisplacement','referenceframe',referenceframe,'maxdeg',maxdeg); @@ -45,5 +61,5 @@ %secular fluid love number: self.tk2secular=0.942; - + %time: self.istime=1; %temporal love numbers by default @@ -73,5 +89,5 @@ error('lovenumbers error message: love numbers should be provided at the same level of accuracy'); end - + ntf=length(self.timefreq); if( size(self.h,2) ~= ntf | size(self.k,2) ~= ntf | size(self.l,2) ~= ntf | size(self.th,2) ~= ntf | size(self.tk,2) ~= ntf | size(self.tl,2) ~= ntf ), @@ -82,20 +98,4 @@ function list=defaultoutputs(self,md) % {{{ list = {}; - end % }}} - function disp(self) % {{{ - disp(sprintf(' lovenumbers parameters:')); - - fielddisplay(self,'h','load Love number for radial displacement'); - fielddisplay(self,'k','load Love number for gravitational potential perturbation'); - fielddisplay(self,'l','load Love number for horizontal displacements'); - - fielddisplay(self,'th','tidal load Love number (deg 2)'); - fielddisplay(self,'tk','tidal load Love number (deg 2)'); - fielddisplay(self,'tl','tidal load Love number (deg 2)'); - fielddisplay(self,'tk2secular','secular fluid Love number'); - - fielddisplay(self,'istime','time (default=1) or frequency love numbers (0)'); - fielddisplay(self,'timefreq','time/frequency vector (yr or 1/yr)'); - end % }}} function marshall(self,prefix,md,fid) % {{{ Index: /issm/trunk-jpl/src/m/classes/lovenumbers.py =================================================================== --- /issm/trunk-jpl/src/m/classes/lovenumbers.py (revision 26300) +++ /issm/trunk-jpl/src/m/classes/lovenumbers.py (revision 26301) @@ -10,6 +10,9 @@ Usage: - lovenumbers = lovenumbers() #will setup love numbers deg 1001 by default - lovenumbers = lovenumbers('maxdeg', 10001); #supply numbers of degrees required (here, 10001) + lovenumbers = lovenumbers() + lovenumbers = lovenumbers('maxdeg', 10000, 'referenceframe', 'CF'); + + Choose numbers of degrees required (1000 by default) and reference frame + (between CF and CM; CM by default) """ @@ -26,8 +29,26 @@ self.tk2secular = 0 # deg 2 secular number + # Time/frequency for visco-elastic love numbers + self.timefreq = [] + self.istime = 1 + options = pairoptions(*args) maxdeg = options.getfieldvalue('maxdeg', 1000) referenceframe = options.getfieldvalue('referenceframe', 'CM') self.setdefaultparameters(maxdeg, referenceframe) + #}}} + + def __repr__(self): #{{{ + s = ' lovenumbers parameters:\n' + s += '{}\n'.format(fielddisplay(self, 'h', 'load Love number for radial displacement')) + s += '{}\n'.format(fielddisplay(self, 'k', 'load Love number for gravitational potential perturbation')) + s += '{}\n'.format(fielddisplay(self, 'l', 'load Love number for horizontal displacements')) + s += '{}\n'.format(fielddisplay(self, 'th', 'tidal load Love number (deg 2)')) + s += '{}\n'.format(fielddisplay(self, 'tk', 'tidal load Love number (deg 2)')) + s += '{}\n'.format(fielddisplay(self, 'tl', 'tidal load Love number (deg 2)')) + s += '{}\n'.format(fielddisplay(self, 'tk2secular', 'secular fluid Love number')) + s += '{}\n'.format(fielddisplay(self, 'istime', 'time (default: 1) or frequency love numbers (0)')) + s += '{}\n'.format(fielddisplay(self, 'timefreq', 'time/frequency vector (yr or 1/yr)')) + return s #}}} @@ -43,9 +64,13 @@ # Secular fluid love number self.tk2secular = 0.942 + + # Time + self.istime = 1 # Temporal love numbers by default + self.timefreq = 0 # Elastic case by default return self #}}} def checkconsistency(self, md, solution, analyses): #{{{ - if ('SealevelriseAnalysis' not in analyses) or (solution == 'TransientSolution' and not md.transient.isslc): + if ('SealevelchangeAnalysis' not in analyses) or (solution == 'TransientSolution' and not md.transient.isslc): return md = checkfield(md, 'fieldname', 'solidearth.lovenumbers.h', 'NaN', 1, 'Inf', 1) @@ -56,7 +81,15 @@ md = checkfield(md, 'fieldname', 'solidearth.lovenumbers.tk', 'NaN', 1, 'Inf', 1) md = checkfield(md, 'fieldname', 'solidearth.lovenumbers.tk2secular', 'NaN', 1, 'Inf', 1) + md = checkfield(md, 'fieldname', 'solidearth.lovenumbers.timefreq', 'NaN', 1, 'Inf', 1) + md = checkfield(md, 'fieldname', 'solidearth.lovenumbers.istime', 'NaN', 1, 'Inf', 1, 'values', [0, 1]) + # Check that love numbers are provided at the same level of accuracy if (self.h.shape[0] != self.k.shape[0]) or (self.h.shape[0] != self.l.shape[0]): raise ValueError('lovenumbers error message: love numbers should be provided at the same level of accuracy') + + ntf = len(self.timefreq) + if (self.h.shape[1] != ntf or self.k.shape[1] != ntf or self.l.shape[1] != ntf or self.th.shape[1] != ntf or self.tk.shape[1] != ntf or self.tl.shape[1] != ntf): + raise ValueError('lovenumbers error message: love numbers should have as many time/frequency steps as the time/frequency vector') + return md #}}} @@ -64,15 +97,4 @@ def defaultoutputs(self, md): #{{{ return[] - #}}} - - def __repr__(self): #{{{ - s = ' lovenumbers parameters:\n' - s += '{}\n'.format(fielddisplay(self, 'h', 'load Love number for radial displacement')) - s += '{}\n'.format(fielddisplay(self, 'k', 'load Love number for gravitational potential perturbation')) - s += '{}\n'.format(fielddisplay(self, 'l', 'load Love number for horizontal displacements')) - s += '{}\n'.format(fielddisplay(self, 'th', 'tidal load Love number (deg 2)')) - s += '{}\n'.format(fielddisplay(self, 'tk', 'tidal load Love number (deg 2)')) - s += '{}\n'.format(fielddisplay(self, 'tk2secular', 'secular fluid Love number')) - return s #}}} @@ -86,4 +108,11 @@ WriteData(fid, prefix, 'object', self, 'fieldname', 'tl', 'name', 'md.solidearth.lovenumbers.tl', 'format', 'DoubleMat', 'mattype', 1) WriteData(fid, prefix, 'object', self, 'data', self.tk2secular, 'fieldname', 'lovenumbers.tk2secular', 'format', 'Double') + + if (self.istime): + scale = md.constants.yts + else: + scale = 1.0 / md.constants.yts + WriteData(fid, prefix, 'object', self, 'fieldname', 'istime', 'name', 'md.solidearth.lovenumbers.istime', 'format', 'Boolean') + WriteData(fid, prefix, 'object', self, 'fieldname', 'timefreq', 'name', 'md.solidearth.lovenumbers.timefreq', 'format', 'DoubleMat', 'mattype', 1, 'scale', scale); #}}} Index: /issm/trunk-jpl/src/m/classes/materials.m =================================================================== --- /issm/trunk-jpl/src/m/classes/materials.m (revision 26300) +++ /issm/trunk-jpl/src/m/classes/materials.m (revision 26301) @@ -204,5 +204,5 @@ fielddisplay(self,'ebm_delta','array describing each layer''s amplitude of the transient relaxation (ratio between elastic rigity to pre-maxwell relaxation rigity), only for EBM rheology (numlayers)'); fielddisplay(self,'ebm_taul','array describing each layer''s starting period for transient relaxation, only for EBM rheology (numlayers) [s]'); - fielddisplay(self,'ebm_tauh','array describing each layer''s array describing each layer''s end period for transient relaxation, only for Burgers rheology (numlayers) [s]'); + fielddisplay(self,'ebm_tauh','array describing each layer''s array describing each layer''s end period for transient relaxation, only for Burgers rheology (numlayers) [s]'); @@ -256,5 +256,5 @@ error('materials checkconsistency error message: Litho burgers_viscosity or burgers_mu has NaN values, inconsistent with rheologymodel choice'); end - if md.materials.rheologymodel(i)==2 & (isnan(md.materials.ebm_alpha(i)) | isnan(md.materials.ebm_delta(i)) | isnan(md.materials.ebm_taul(i)) | isnan(md.materials.ebm_tauh(i)) ), + if md.materials.rheologymodel(i)==2 & (isnan(md.materials.ebm_alpha(i)) | isnan(md.materials.ebm_delta(i)) | isnan(md.materials.ebm_taul(i)) | isnan(md.materials.ebm_tauh(i))), error('materials checkconsistency error message: Litho ebm_alpha, ebm_delta, ebm_taul or ebm_tauh has NaN values, inconsistent with rheologymodel choice'); end @@ -284,5 +284,5 @@ %1: MatdamageiceEnum 2: MatestarEnum 3: MaticeEnum 4: MatenhancediceEnum 5: MaterialsEnum WriteData(fid,prefix,'name','md.materials.nature','data',naturetointeger(self.nature),'format','IntMat','mattype',3); - WriteData(fid,prefix,'name','md.materials.type','data',5,'format','Integer'); %DANGER, this can evolve if you have classes. + WriteData(fid,prefix,'name','md.materials.type','data',5,'format','Integer'); %DANGER: this can evolve if you have classes. for i=1:length(self.nature), nat=self.nature{i}; @@ -326,4 +326,5 @@ end earth_density=earth_density/self.radius(self.numlayers+1)^3; + disp(earth_density) self.earth_density=earth_density; case 'hydro' @@ -541,5 +542,5 @@ vs = vs + t1*(r2^3) + t2*(r2^4) + t3*(r2^5) + t4*(r2^6) - ( t1*(r1^3) + t2*(r1^4) + t3*(r1^5) + t4*(r1^6) ); - end + end ro = ro*3 / (rad(j+1)^3-rad(j)^3); vp = vp*3 /(rad(j+1)^3-rad(j)^3); @@ -557,5 +558,4 @@ end % }}} - end end @@ -584,4 +584,2 @@ end end % }}} - - Index: /issm/trunk-jpl/src/m/classes/materials.py =================================================================== --- /issm/trunk-jpl/src/m/classes/materials.py (revision 26300) +++ /issm/trunk-jpl/src/m/classes/materials.py (revision 26301) @@ -25,4 +25,5 @@ raise RuntimeError("materials constructor error message: nature of the material not supported yet! ('ice' or 'litho' or 'hydro')") + # Start filling in the dynamic fields (not truly dynamic under Python) for i in range(len(self.nature)): nat = self.nature[i] @@ -67,5 +68,4 @@ self.earth_density = 0 - # Set default parameters self.setdefaultparameters() #}}} @@ -105,9 +105,9 @@ s += '{}\n'.format(fielddisplay(self, 'ebm_delta', 'array describing each layer\'s amplitude of the transient relaxation (ratio between elastic rigity to pre-maxwell relaxation rigity), only for EBM rheology (numlayers)')) s += '{}\n'.format(fielddisplay(self, 'ebm_taul', 'array describing each layer\'s starting period for transient relaxation, only for EBM rheology (numlayers) [s]')) - s += '{}\n'.format(fielddisplay(self, 'ebm_tauh', 'array describing each layer''s array describing each layer\'s end period for transient relaxation, only for Burgers rheology (numlayers) [s]')) + s += '{}\n'.format(fielddisplay(self, 'ebm_tauh', 'array describing each layer''s array describing each layer\'s end period for transient relaxation, only for Burgers rheology (numlayers) [s]')) s += '{}\n'.format(fielddisplay(self, 'rheologymodel', 'array describing whether we adopt a Maxwell (0), Burgers (1) or EBM (2) rheology (default: 0)')) s += '{}\n'.format(fielddisplay(self, 'density', 'array describing each layer\'s density (numlayers) [kg/m^3]')) - s += '{}\n'.format(fielddisplay(self, 'issolid', 'array describing whether the layer is solid or liquid (default 1) (numlayers)')) + s += '{}\n'.format(fielddisplay(self, 'issolid', 'array describing whether the layer is solid or liquid (default: 1) (numlayers)')) elif nat == 'hydro': s += 'Hydro:\n' @@ -127,11 +127,11 @@ if nat == 'ice': # Ice density (kg/m^3) - self.rho_ice = 917. + self.rho_ice = 917 # Ocean water density (kg/m^3) - self.rho_water = 1023. + self.rho_water = 1023 # Fresh water density (kg/m^3) - self.rho_freshwater = 1000. + self.rho_freshwater = 1000 # Water viscosity (N.s/m^2) @@ -139,5 +139,5 @@ # Ice heat capacity cp (J/kg/K) - self.heatcapacity = 2093. + self.heatcapacity = 2093 # Ice latent heat of fusion L (J/kg) @@ -160,5 +160,5 @@ # Mixed layer (ice-water interface) heat capacity (J/kg/K) - self.mixed_layer_capacity = 3974. + self.mixed_layer_capacity = 3974 # Thermal exchange velocity (ice-water interface) (m/s) @@ -193,15 +193,15 @@ self.ebm_tauh = [np.nan, np.nan] self.rheologymodel = [0, 0] - self.density = [5.51e3, 5.50e3] # (Pa) # Mantle and lithosphere density [kg/m^3] - self.issolid = [1, 1] # Is layer solid or liquid? + self.density = [5.51e3, 5.50e3] # (Pa) # Mantle and lithosphere density [kg/m^3] + self.issolid = [1, 1] # Is layer solid or liquid? elif nat == 'hydro': # Ice density (kg/m^3) - self.rho_ice = 917. + self.rho_ice = 917 # Ocean water density (kg/m^3) - self.rho_water = 1023. + self.rho_water = 1023 # Fresh water density (kg/m^3) - self.rho_freshwater = 1000. + self.rho_freshwater = 1000 else: raise RuntimeError("materials setdefaultparameters error message: nature of the material not supported yet! ('ice' or 'litho' or 'hydro')") @@ -251,5 +251,5 @@ ind = np.where(md.materials.issolid == 0)[0] if np.sum(np.in1d(np.diff(ind),1) >= 1): # If there are at least two consecutive indices that contain issolid = 0 - raise RuntimeError('Fluid layers detected at layers #{} but having 2 or more adjacent fluid layers is not supported yet. Consider merging them.'.format(i)) + raise RuntimeError('Fluid layers detected at layers #' + indices + ', but having 2 or more adjacent fluid layers is not supported yet. Consider merging them.') elif nat == 'hydro': @@ -267,5 +267,5 @@ #1: MatdamageiceEnum 2: MatestarEnum 3: MaticeEnum 4: MatenhancediceEnum 5: MaterialsEnum WriteData(fid, prefix, 'name', 'md.materials.nature', 'data', naturetointeger(self.nature), 'format', 'IntMat', 'mattype', 3) - WriteData(fid, prefix, 'name', 'md.materials.type', 'data', 5, 'format', 'Integer') #DANGER, this can evolve if you have classes + WriteData(fid, prefix, 'name', 'md.materials.type', 'data', 5, 'format', 'Integer') #DANGER: this can evolve if you have classes for i in range(len(self.nature)): nat = self.nature[i] @@ -307,4 +307,5 @@ earth_density = earth_density + (pow(self.radius[i + 1], 3) - pow(self.radius[i], 3)) * self.density[i] earth_density = earth_density / pow(self.radius[self.numlayers], 3) + print(earth_density) self.earth_density = earth_density elif nat == 'hydro': Index: /issm/trunk-jpl/src/m/classes/model.m =================================================================== --- /issm/trunk-jpl/src/m/classes/model.m (revision 26300) +++ /issm/trunk-jpl/src/m/classes/model.m (revision 26301) @@ -39,8 +39,8 @@ steadystate = 0; transient = 0; - levelset = 0; + levelset = 0; calving = 0; frontalforcings = 0; - love = 0; + love = 0; esa = 0; sampling = 0; @@ -49,5 +49,5 @@ inversion = 0; qmu = 0; - amr = 0; + amr = 0; results = 0; outputdefinition = 0; @@ -207,4 +207,48 @@ end %}}} + function disp(self) % {{{ + disp(sprintf('%19s: %-22s -- %s','mesh' ,['[1x1 ' class(self.mesh) ']'],'mesh properties')); + disp(sprintf('%19s: %-22s -- %s','mask' ,['[1x1 ' class(self.mask) ']'],'defines grounded and floating elements')); + disp(sprintf('%19s: %-22s -- %s','geometry' ,['[1x1 ' class(self.geometry) ']'],'surface elevation, bedrock topography, ice thickness,...')); + disp(sprintf('%19s: %-22s -- %s','constants' ,['[1x1 ' class(self.constants) ']'],'physical constants')); + disp(sprintf('%19s: %-22s -- %s','smb' ,['[1x1 ' class(self.smb) ']'],'surface mass balance')); + disp(sprintf('%19s: %-22s -- %s','basalforcings' ,['[1x1 ' class(self.basalforcings) ']'],'bed forcings')); + disp(sprintf('%19s: %-22s -- %s','materials' ,['[1x1 ' class(self.materials) ']'],'material properties')); + disp(sprintf('%19s: %-22s -- %s','damage' ,['[1x1 ' class(self.damage) ']'],'parameters for damage evolution solution')); + disp(sprintf('%19s: %-22s -- %s','friction' ,['[1x1 ' class(self.friction) ']'],'basal friction/drag properties')); + disp(sprintf('%19s: %-22s -- %s','flowequation' ,['[1x1 ' class(self.flowequation) ']'],'flow equations')); + disp(sprintf('%19s: %-22s -- %s','timestepping' ,['[1x1 ' class(self.timestepping) ']'],'time stepping for transient models')); + disp(sprintf('%19s: %-22s -- %s','initialization' ,['[1x1 ' class(self.initialization) ']'],'initial guess/state')); + disp(sprintf('%19s: %-22s -- %s','rifts' ,['[1x1 ' class(self.rifts) ']'],'rifts properties')); + disp(sprintf('%19s: %-22s -- %s','solidearth' ,['[1x1 ' class(self.solidearth) ']'],'solidearth inputs and settings')); + disp(sprintf('%19s: %-22s -- %s','dsl' ,['[1x1 ' class(self.dsl) ']'],'dynamic sea-level ')); + disp(sprintf('%19s: %-22s -- %s','debug' ,['[1x1 ' class(self.debug) ']'],'debugging tools (valgrind, gprof)')); + disp(sprintf('%19s: %-22s -- %s','verbose' ,['[1x1 ' class(self.verbose) ']'],'verbosity level in solve')); + disp(sprintf('%19s: %-22s -- %s','settings' ,['[1x1 ' class(self.settings) ']'],'settings properties')); + disp(sprintf('%19s: %-22s -- %s','toolkits' ,['[1x1 ' class(self.toolkits) ']'],'PETSc options for each solution')); + disp(sprintf('%19s: %-22s -- %s','cluster' ,['[1x1 ' class(self.cluster) ']'],'cluster parameters (number of CPUs...)')); + disp(sprintf('%19s: %-22s -- %s','balancethickness',['[1x1 ' class(self.balancethickness) ']'],'parameters for balancethickness solution')); + disp(sprintf('%19s: %-22s -- %s','stressbalance' ,['[1x1 ' class(self.stressbalance) ']'],'parameters for stressbalance solution')); + disp(sprintf('%19s: %-22s -- %s','groundingline' ,['[1x1 ' class(self.groundingline) ']'],'parameters for groundingline solution')); + disp(sprintf('%19s: %-22s -- %s','hydrology' ,['[1x1 ' class(self.hydrology) ']'],'parameters for hydrology solution')); + disp(sprintf('%19s: %-22s -- %s','masstransport' ,['[1x1 ' class(self.masstransport) ']'],'parameters for masstransport solution')); + disp(sprintf('%19s: %-22s -- %s','thermal' ,['[1x1 ' class(self.thermal) ']'],'parameters for thermal solution')); + disp(sprintf('%19s: %-22s -- %s','steadystate' ,['[1x1 ' class(self.steadystate) ']'],'parameters for steadystate solution')); + disp(sprintf('%19s: %-22s -- %s','transient' ,['[1x1 ' class(self.transient) ']'],'parHwoameters for transient solution')); + disp(sprintf('%19s: %-22s -- %s','levelset' ,['[1x1 ' class(self.levelset) ']'],'parameters for moving boundaries (level-set method)')); + disp(sprintf('%19s: %-22s -- %s','calving' ,['[1x1 ' class(self.calving) ']'],'parameters for calving')); + disp(sprintf('%19s: %-22s -- %s','frontalforcings' ,['[1x1 ' class(self.frontalforcings) ']'],'parameters for frontalforcings')); + disp(sprintf('%19s: %-22s -- %s','esa' ,['[1x1 ' class(self.esa) ']'],'parameters for elastic adjustment solution')); + disp(sprintf('%19s: %-22s -- %s','love' ,['[1x1 ' class(self.love) ']'],'parameters for love solution')); + disp(sprintf('%19s: %-22s -- %s','sampling' ,['[1x1 ' class(self.sampling) ']'],'parameters for stochastic sampler')); + disp(sprintf('%19s: %-22s -- %s','autodiff' ,['[1x1 ' class(self.autodiff) ']'],'automatic differentiation parameters')); + disp(sprintf('%19s: %-22s -- %s','inversion' ,['[1x1 ' class(self.inversion) ']'],'parameters for inverse methods')); + disp(sprintf('%19s: %-22s -- %s','qmu' ,['[1x1 ' class(self.qmu) ']'],'Dakota properties')); + disp(sprintf('%19s: %-22s -- %s','amr' ,['[1x1 ' class(self.amr) ']'],'adaptive mesh refinement properties')); + disp(sprintf('%19s: %-22s -- %s','outputdefinition',['[1x1 ' class(self.outputdefinition) ']'],'output definition')); + disp(sprintf('%19s: %-22s -- %s','results' ,['[1x1 ' class(self.results) ']'],'model results')); + disp(sprintf('%19s: %-22s -- %s','radaroverlay' ,['[1x1 ' class(self.radaroverlay) ']'],'radar image for plot overlay')); + disp(sprintf('%19s: %-22s -- %s','miscellaneous' ,['[1x1 ' class(self.miscellaneous) ']'],'miscellaneous fields')); + end % }}} function md = setdefaultparameters(md,planet) % {{{ @@ -1133,5 +1177,5 @@ md.mesh.average_vertex_connectivity=100; end - end % }}} + end % }}} function md = structtomodel(md,structmd) % {{{ @@ -1584,95 +1628,52 @@ end end % }}} - function disp(self) % {{{ - disp(sprintf('%19s: %-22s -- %s','mesh' ,['[1x1 ' class(self.mesh) ']'],'mesh properties')); - disp(sprintf('%19s: %-22s -- %s','mask' ,['[1x1 ' class(self.mask) ']'],'defines grounded and floating elements')); - disp(sprintf('%19s: %-22s -- %s','geometry' ,['[1x1 ' class(self.geometry) ']'],'surface elevation, bedrock topography, ice thickness,...')); - disp(sprintf('%19s: %-22s -- %s','constants' ,['[1x1 ' class(self.constants) ']'],'physical constants')); - disp(sprintf('%19s: %-22s -- %s','smb' ,['[1x1 ' class(self.smb) ']'],'surface mass balance')); - disp(sprintf('%19s: %-22s -- %s','basalforcings' ,['[1x1 ' class(self.basalforcings) ']'],'bed forcings')); - disp(sprintf('%19s: %-22s -- %s','materials' ,['[1x1 ' class(self.materials) ']'],'material properties')); - disp(sprintf('%19s: %-22s -- %s','damage' ,['[1x1 ' class(self.damage) ']'],'parameters for damage evolution solution')); - disp(sprintf('%19s: %-22s -- %s','friction' ,['[1x1 ' class(self.friction) ']'],'basal friction/drag properties')); - disp(sprintf('%19s: %-22s -- %s','flowequation' ,['[1x1 ' class(self.flowequation) ']'],'flow equations')); - disp(sprintf('%19s: %-22s -- %s','timestepping' ,['[1x1 ' class(self.timestepping) ']'],'time stepping for transient models')); - disp(sprintf('%19s: %-22s -- %s','initialization' ,['[1x1 ' class(self.initialization) ']'],'initial guess/state')); - disp(sprintf('%19s: %-22s -- %s','rifts' ,['[1x1 ' class(self.rifts) ']'],'rifts properties')); - disp(sprintf('%19s: %-22s -- %s','solidearth' ,['[1x1 ' class(self.solidearth) ']'],'solidearth inputs and settings')); - disp(sprintf('%19s: %-22s -- %s','dsl' ,['[1x1 ' class(self.dsl) ']'],'dynamic sea-level ')); - disp(sprintf('%19s: %-22s -- %s','debug' ,['[1x1 ' class(self.debug) ']'],'debugging tools (valgrind, gprof)')); - disp(sprintf('%19s: %-22s -- %s','verbose' ,['[1x1 ' class(self.verbose) ']'],'verbosity level in solve')); - disp(sprintf('%19s: %-22s -- %s','settings' ,['[1x1 ' class(self.settings) ']'],'settings properties')); - disp(sprintf('%19s: %-22s -- %s','toolkits' ,['[1x1 ' class(self.toolkits) ']'],'PETSc options for each solution')); - disp(sprintf('%19s: %-22s -- %s','cluster' ,['[1x1 ' class(self.cluster) ']'],'cluster parameters (number of CPUs...)')); - disp(sprintf('%19s: %-22s -- %s','balancethickness',['[1x1 ' class(self.balancethickness) ']'],'parameters for balancethickness solution')); - disp(sprintf('%19s: %-22s -- %s','stressbalance' ,['[1x1 ' class(self.stressbalance) ']'],'parameters for stressbalance solution')); - disp(sprintf('%19s: %-22s -- %s','groundingline' ,['[1x1 ' class(self.groundingline) ']'],'parameters for groundingline solution')); - disp(sprintf('%19s: %-22s -- %s','hydrology' ,['[1x1 ' class(self.hydrology) ']'],'parameters for hydrology solution')); - disp(sprintf('%19s: %-22s -- %s','masstransport' ,['[1x1 ' class(self.masstransport) ']'],'parameters for masstransport solution')); - disp(sprintf('%19s: %-22s -- %s','thermal' ,['[1x1 ' class(self.thermal) ']'],'parameters for thermal solution')); - disp(sprintf('%19s: %-22s -- %s','steadystate' ,['[1x1 ' class(self.steadystate) ']'],'parameters for steadystate solution')); - disp(sprintf('%19s: %-22s -- %s','transient' ,['[1x1 ' class(self.transient) ']'],'parHwoameters for transient solution')); - disp(sprintf('%19s: %-22s -- %s','levelset' ,['[1x1 ' class(self.levelset) ']'],'parameters for moving boundaries (level-set method)')); - disp(sprintf('%19s: %-22s -- %s','calving' ,['[1x1 ' class(self.calving) ']'],'parameters for calving')); - disp(sprintf('%19s: %-22s -- %s','frontalforcings' ,['[1x1 ' class(self.frontalforcings) ']'],'parameters for frontalforcings')); - disp(sprintf('%19s: %-22s -- %s','esa' ,['[1x1 ' class(self.esa) ']'],'parameters for elastic adjustment solution')); - disp(sprintf('%19s: %-22s -- %s','love' ,['[1x1 ' class(self.love) ']'],'parameters for love solution')); - disp(sprintf('%19s: %-22s -- %s','sampling' ,['[1x1 ' class(self.sampling) ']'],'parameters for stochastic sampler')); - disp(sprintf('%19s: %-22s -- %s','autodiff' ,['[1x1 ' class(self.autodiff) ']'],'automatic differentiation parameters')); - disp(sprintf('%19s: %-22s -- %s','inversion' ,['[1x1 ' class(self.inversion) ']'],'parameters for inverse methods')); - disp(sprintf('%19s: %-22s -- %s','qmu' ,['[1x1 ' class(self.qmu) ']'],'Dakota properties')); - disp(sprintf('%19s: %-22s -- %s','amr' ,['[1x1 ' class(self.amr) ']'],'adaptive mesh refinement properties')); - disp(sprintf('%19s: %-22s -- %s','outputdefinition',['[1x1 ' class(self.outputdefinition) ']'],'output definition')); - disp(sprintf('%19s: %-22s -- %s','results' ,['[1x1 ' class(self.results) ']'],'model results')); - disp(sprintf('%19s: %-22s -- %s','radaroverlay' ,['[1x1 ' class(self.radaroverlay) ']'],'radar image for plot overlay')); - disp(sprintf('%19s: %-22s -- %s','miscellaneous' ,['[1x1 ' class(self.miscellaneous) ']'],'miscellaneous fields')); - end % }}} function memory(self) % {{{ - disp(sprintf('\nMemory imprint:\n')); - - fields=properties('model'); - mem=0; - - for i=1:length(fields), - field=self.(fields{i}); - s=whos('field'); - mem=mem+s.bytes/1e6; - disp(sprintf('%19s: %6.2f Mb',fields{i},s.bytes/1e6)); + disp(sprintf('\nMemory imprint:\n')); + + fields=properties('model'); + mem=0; + + for i=1:length(fields), + field=self.(fields{i}); + s=whos('field'); + mem=mem+s.bytes/1e6; + disp(sprintf('%19s: %6.2f Mb',fields{i},s.bytes/1e6)); + end + disp(sprintf('%19s--%10s','--------------','--------------')); + disp(sprintf('%19s: %g Mb','Total',mem)); end - disp(sprintf('%19s--%10s','--------------','--------------')); - disp(sprintf('%19s: %g Mb','Total',mem)); - end % }}} + % }}} function netcdf(self,filename) % {{{ - %NETCDF - save model as netcdf - % - % Usage: - % netcdf(md,filename) - % - % Example: - % netcdf(md,'model.nc'); - - disp('Saving model as NetCDF'); - %1. Create NetCDF file - ncid=netcdf.create(filename,'CLOBBER'); - netcdf.putAtt(ncid,netcdf.getConstant('NC_GLOBAL'),'Conventions','CF-1.4'); - netcdf.putAtt(ncid,netcdf.getConstant('NC_GLOBAL'),'Title',['ISSM model (' self.miscellaneous.name ')']); - netcdf.putAtt(ncid,netcdf.getConstant('NC_GLOBAL'),'Author',getenv('USER')); - netcdf.putAtt(ncid,netcdf.getConstant('NC_GLOBAL'),'Date',datestr(now)); - - %Preallocate variable id, needed to write variables in netcdf file - var_id=zeros(1000,1);%preallocate - - for step=1:2, - counter=0; - [var_id,counter]=structtonc(ncid,'md',self,0,var_id,counter,step); - if step==1, netcdf.endDef(ncid); end - end - - if counter>1000, - warning(['preallocation of var_id need to be updated from ' num2str(1000) ' to ' num2str(counter)]); - end - - netcdf.close(ncid) + %NETCDF - save model as netcdf + % + % Usage: + % netcdf(md,filename) + % + % Example: + % netcdf(md,'model.nc'); + + disp('Saving model as NetCDF'); + %1. Create NetCDF file + ncid=netcdf.create(filename,'CLOBBER'); + netcdf.putAtt(ncid,netcdf.getConstant('NC_GLOBAL'),'Conventions','CF-1.4'); + netcdf.putAtt(ncid,netcdf.getConstant('NC_GLOBAL'),'Title',['ISSM model (' self.miscellaneous.name ')']); + netcdf.putAtt(ncid,netcdf.getConstant('NC_GLOBAL'),'Author',getenv('USER')); + netcdf.putAtt(ncid,netcdf.getConstant('NC_GLOBAL'),'Date',datestr(now)); + + %Preallocate variable id, needed to write variables in netcdf file + var_id=zeros(1000,1);%preallocate + + for step=1:2, + counter=0; + [var_id,counter]=structtonc(ncid,'md',self,0,var_id,counter,step); + if step==1, netcdf.endDef(ncid); end + end + + if counter>1000, + warning(['preallocation of var_id need to be updated from ' num2str(1000) ' to ' num2str(counter)]); + end + + netcdf.close(ncid) end % }}} function xylim(self) % {{{ @@ -1682,38 +1683,38 @@ end % }}} function md=upload(md) % {{{ - %the goal of this routine is to upload the model onto a server, and to empty it. - %So first, save the model with a unique name and upload the file to the server: - random_part=fix(rand(1)*10000); - id=[md.miscellaneous.name '-' regexprep(datestr(now),'[^\w'']','') '-' num2str(random_part) '-' getenv('USER') '-' oshostname() '.upload']; - eval(['save ' id ' md']); - - %Now, upload the file: - issmscpout(md.settings.upload_server,md.settings.upload_path,md.settings.upload_login,md.settings.upload_port,{id},1); - - %Now, empty this model of everything except settings, and record name of file we just uploaded! - settings_back=md.settings; - md=model(); - md.settings=settings_back; - md.settings.upload_filename=id; - - %get locally rid of file that was uploaded - eval(['delete ' id]); + %the goal of this routine is to upload the model onto a server, and to empty it. + %So first, save the model with a unique name and upload the file to the server: + random_part=fix(rand(1)*10000); + id=[md.miscellaneous.name '-' regexprep(datestr(now),'[^\w'']','') '-' num2str(random_part) '-' getenv('USER') '-' oshostname() '.upload']; + eval(['save ' id ' md']); + + %Now, upload the file: + issmscpout(md.settings.upload_server,md.settings.upload_path,md.settings.upload_login,md.settings.upload_port,{id},1); + + %Now, empty this model of everything except settings, and record name of file we just uploaded! + settings_back=md.settings; + md=model(); + md.settings=settings_back; + md.settings.upload_filename=id; + + %get locally rid of file that was uploaded + eval(['delete ' id]); end % }}} function md=download(md) % {{{ - %the goal of this routine is to download the internals of the current model from a server, because - %this model is empty, except for the settings which tell us where to go and find this model! - - %Download the file: - issmscpin(md.settings.upload_server, md.settings.upload_login, md.settings.upload_port, md.settings.upload_path, {md.settings.upload_filename}); - - name=md.settings.upload_filename; - - %Now, load this model: - md=loadmodel(md.settings.upload_filename); - - %get locally rid of file that was downloaded - eval(['delete ' name]); + %the goal of this routine is to download the internals of the current model from a server, because + %this model is empty, except for the settings which tell us where to go and find this model! + + %Download the file: + issmscpin(md.settings.upload_server, md.settings.upload_login, md.settings.upload_port, md.settings.upload_path, {md.settings.upload_filename}); + + name=md.settings.upload_filename; + + %Now, load this model: + md=loadmodel(md.settings.upload_filename); + + %get locally rid of file that was downloaded + eval(['delete ' name]); end % }}} Index: /issm/trunk-jpl/src/m/classes/model.py =================================================================== --- /issm/trunk-jpl/src/m/classes/model.py (revision 26300) +++ /issm/trunk-jpl/src/m/classes/model.py (revision 26301) @@ -76,22 +76,28 @@ class model(object): - #properties + """MODEL - class definition + + Usage: + md = model() + """ + def __init__(self, *args): #{{{ self.mesh = None self.mask = None + + self.geometry = None self.constants = None - self.geometry = None - self.initialization = None self.smb = None self.basalforcings = None - self.friction = None - self.rifts = None - self.solidearth = None - self.dsl = None - self.timestepping = None - self.groundingline = None self.materials = None self.damage = None + self.friction = None self.flowequation = None + self.timestepping = None + self.initialization = None + self.rifts = None + self.dsl = None + self.solidearth = None + self.debug = None self.verbose = None @@ -99,6 +105,8 @@ self.toolkits = None self.cluster = None + self.balancethickness = None self.stressbalance = None + self.groundingline = None self.hydrology = None self.masstransport = None @@ -112,71 +120,22 @@ self.esa = None self.sampling = None + self.autodiff = None self.inversion = None self.qmu = None self.amr = None - self.radaroverlay = None self.results = None self.outputdefinition = None + self.radaroverlay = None self.miscellaneous = None self.private = None - nargs = len(args) - - if nargs == 0: + if len(args) == 0: self.setdefaultparameters('earth') else: - self.setdefaultparameters(args[0]) options = pairoptions(*args) planet = options.getfieldvalue('planet', 'earth') self.setdefaultparameters(planet) -#}}} - - def properties(self): # {{{ - # ordered list of properties since vars(self) is random - return ['mesh', - 'mask', - 'constants', - 'geometry', - 'initialization', - 'smb', - 'basalforcings', - 'friction', - 'rifts', - 'solidearth', - 'dsl', - 'timestepping', - 'groundingline', - 'materials', - 'damage', - 'flowequation', - 'debug', - 'verbose', - 'settings', - 'toolkits', - 'cluster', - 'balancethickness', - 'stressbalance', - 'hydrology', - 'masstransport', - 'thermal', - 'steadystate', - 'transient', - 'levelset', - 'calving', - 'frontalforcings', - 'love', - 'esa', - 'sampling', - 'autodiff', - 'inversion', - 'qmu', - 'amr', - 'radaroverlay', - 'results', - 'outputdefinition', - 'miscellaneous', - 'private'] - # }}} + #}}} def __repr__(obj): #{{{ @@ -185,6 +144,4 @@ # already converted .__repr__ method for examples) # - - #print "Here %s the number: %d" % ("is", 37) s = "%19s: %-22s -- %s" % ("mesh", "[%s %s]" % ("1x1", obj.mesh.__class__.__name__), "mesh properties") s = "%s\n%s" % (s, "%19s: %-22s -- %s" % ("mask", "[%s %s]" % ("1x1", obj.mask.__class__.__name__), "defines grounded and floating elements")) @@ -230,5 +187,52 @@ s = "%s\n%s" % (s, "%19s: %-22s -- %s" % ("miscellaneous", "[%s %s]" % ("1x1", obj.miscellaneous.__class__.__name__), "miscellaneous fields")) return s - # }}} + #}}} + + def properties(self): #{{{ + # ordered list of properties since vars(self) is random + return ['mesh', + 'mask', + 'geometry', + 'constants', + 'smb', + 'basalforcings', + 'materials', + 'damage', + 'friction', + 'flowequation', + 'timestepping', + 'initialization', + 'rifts', + 'dsl', + 'solidearth', + 'debug', + 'verbose', + 'settings', + 'toolkits', + 'cluster', + 'balancethickness', + 'stressbalance', + 'groundingline', + 'hydrology', + 'masstransport', + 'thermal', + 'steadystate', + 'transient', + 'levelset', + 'calving', + 'frontalforcings', + 'love', + 'esa', + 'sampling', + 'autodiff', + 'inversion', + 'qmu', + 'amr', + 'results', + 'outputdefinition', + 'radaroverlay', + 'miscellaneous', + 'private'] + #}}} def setdefaultparameters(self, planet): #{{{ @@ -278,12 +282,12 @@ #}}} - def checkmessage(self, string): # {{{ + def checkmessage(self, string): #{{{ print("model not consistent: {}".format(string)) self.private.isconsistent = False return self - # }}} + #}}} #@staticmethod - def extract(self, area): # {{{ + def extract(self, area): #{{{ """EXTRACT - extract a model according to an Argus contour or flag list @@ -562,7 +566,7 @@ return md2 - # }}} - - def extrude(md, *args): # {{{ + #}}} + + def extrude(md, *args): #{{{ """EXTRUDE - vertically extrude a 2d mesh @@ -749,5 +753,5 @@ return md - # }}} + #}}} def collapse(md): #{{{ Index: /issm/trunk-jpl/src/m/classes/sampling.m =================================================================== --- /issm/trunk-jpl/src/m/classes/sampling.m (revision 26300) +++ /issm/trunk-jpl/src/m/classes/sampling.m (revision 26301) @@ -6,105 +6,96 @@ classdef sampling properties (SetAccess=public) - kappa = NaN; - tau = 0; - beta = NaN; - phi = 0; - alpha = 0; - robin = 0; - seed = 0; - requested_outputs = {}; - end + kappa = NaN; + tau = 0; + beta = NaN; + phi = 0; + alpha = 0; + robin = 0; + seed = 0; + requested_outputs = {}; + end methods function self = sampling(varargin) % {{{ - switch nargin + switch nargin case 0 self=setdefaultparameters(self); otherwise error('constructor not supported'); - end + end end % }}} - function list = defaultoutputs(self,md) % {{{ + function disp(self) % {{{ - list = {}; + disp(sprintf(' Sampling parameters:')); - end % }}} + disp(sprintf('\n %s','Parameters of PDE operator (kappa^2 I-Laplacian)^(alpha/2)(tau):')); + fielddisplay(self,'kappa','coefficient of the identity operator'); + fielddisplay(self,'tau','scaling coefficient of the solution (default: 1.0)'); + fielddisplay(self,'alpha','exponent in PDE operator, (default: 2.0, BiLaplacian covariance operator)'); + + disp(sprintf('\n %s','Parameters of Robin boundary conditions nabla () \cdot normvec + beta ():')); + fielddisplay(self,'robin','Apply Robin boundary conditions (1 if applied and 0 for homogenous Neumann boundary conditions) (default: 0)'); + fielddisplay(self,'beta','Coefficient in Robin boundary conditions (to be defined for robin = 1)'); + + disp(sprintf('\n %s','Parameters for first-order autoregressive process (X_t = phi X_{t-1} + noise) (if transient):')); + fielddisplay(self,'phi','Temporal correlation factor (|phi|<1 for stationary process, phi = 1 for random walk process) (default 0)'); + + disp(sprintf('\n %s','Other parameters of stochastic sampler:')); + fielddisplay(self,'seed','Seed for pseudorandom number generator (given seed if >=0 and random seed if <0) (default: -1)'); + fielddisplay(self,'requested_outputs','additional outputs requested (not implemented yet)'); + + end % }}}',' function self = setdefaultparameters(self) % {{{ - - %Scaling coefficient - self.tau=1; - %Apply Robin boundary conditions - self.robin=0; - - %Temporal correlation factor - self.phi=0; - - %Exponent in fraction SPDE (default=2, biLaplacian covariance + %Scaling coefficient + self.tau=1; + + %Apply Robin boundary conditions + self.robin=0; + + %Temporal correlation factor + self.phi=0; + + %Exponent in fraction SPDE (default: 2, biLaplacian covariance %operator) self.alpha=2; % Default - - %Seed for pseudorandom number generator (default -1 for random seed) - self.seed=-1; - + + %Seed for pseudorandom number generator (default: -1, for random seed) + self.seed=-1; + %default output self.requested_outputs={'default'}; end % }}} - function md = setparameters(self,md,lc,sigma) % {{{ - - nu = self.alpha-1; - KAPPA = sqrt(8*nu)/lc; - TAU = sqrt(gamma(nu)/(gamma(self.alpha)*(4*pi)*KAPPA^(2*nu)*sigma^2)); - md.sampling.kappa = KAPPA*ones(md.mesh.numberofvertices,1); - md.sampling.tau = TAU; - + function list = defaultoutputs(self,md) % {{{ + + list = {}; + end % }}} function md = checkconsistency(self,md,solution,analyses) % {{{ - - if ~ismember('SamplingAnalysis',analyses), return; end - - md = checkfield(md,'fieldname','sampling.kappa','NaN',1,'Inf',1,'size',[md.mesh.numberofvertices 1],'>',0); - md = checkfield(md,'fieldname','sampling.tau','NaN',1,'Inf',1,'numel',1,'>',0); - md = checkfield(md,'fieldname','sampling.robin','numel',1,'values',[0 1]); - if(md.sampling.robin) - md = checkfield(md,'fieldname','sampling.beta','NaN',1,'Inf',1,'size',[md.mesh.numberofvertices 1],'>',0); - end - md = checkfield(md,'fieldname','sampling.phi','NaN',1,'Inf',1,'numel',1,'>=',0); - md = checkfield(md,'fieldname','sampling.alpha','NaN',1,'Inf',1,'numel',1,'>',0); - md = checkfield(md,'fieldname','sampling.seed','NaN',1,'Inf',1,'numel',1); - md = checkfield(md,'fieldname','sampling.requested_outputs','stringrow',1); - end % }}} - function disp(self) % {{{ + if ~ismember('SamplingAnalysis',analyses), return; end - disp(sprintf(' Sampling parameters:')); + md = checkfield(md,'fieldname','sampling.kappa','NaN',1,'Inf',1,'size',[md.mesh.numberofvertices 1],'>',0); + md = checkfield(md,'fieldname','sampling.tau','NaN',1,'Inf',1,'numel',1,'>',0); + md = checkfield(md,'fieldname','sampling.robin','numel',1,'values',[0 1]); + if(md.sampling.robin) + md = checkfield(md,'fieldname','sampling.beta','NaN',1,'Inf',1,'size',[md.mesh.numberofvertices 1],'>',0); + end + md = checkfield(md,'fieldname','sampling.phi','NaN',1,'Inf',1,'numel',1,'>=',0); + md = checkfield(md,'fieldname','sampling.alpha','NaN',1,'Inf',1,'numel',1,'>',0); + md = checkfield(md,'fieldname','sampling.seed','NaN',1,'Inf',1,'numel',1); + md = checkfield(md,'fieldname','sampling.requested_outputs','stringrow',1); - disp(sprintf('\n %s','Parameters of PDE operator (kappa^2 I-Laplacian)^(alpha/2)(tau):')); - fielddisplay(self,'kappa','coefficient of the identity operator'); - fielddisplay(self,'tau','scaling coefficient of the solution (default 1.0)'); - fielddisplay(self,'alpha','exponent in PDE operator, (default 2.0, BiLaplacian covariance operator)'); - - disp(sprintf('\n %s','Parameters of Robin boundary conditions nabla () \cdot normvec + beta ():')); - fielddisplay(self,'robin','Apply Robin boundary conditions (1 if applied and 0 for homogenous Neumann boundary conditions) (default 0)'); - fielddisplay(self,'beta','Coefficient in Robin boundary conditions (to be defined for robin = 1)'); - - disp(sprintf('\n %s','Parameters for first-order autoregressive process (X_t = phi X_{t-1} + noise) (if transient):')); - fielddisplay(self,'phi','Temporal correlation factor (|phi|<1 for stationary process, phi = 1 for random walk process) (default 0)'); - - disp(sprintf('\n %s','Other parameters of stochastic sampler:')); - fielddisplay(self,'seed','Seed for pseudorandom number generator (given seed if >=0 and random seed if <0) (default -1)'); - fielddisplay(self,'requested_outputs','additional outputs requested (not implemented yet)'); - end % }}} function marshall(self,prefix,md,fid) % {{{ - WriteData(fid,prefix,'object',self,'fieldname','kappa','format','DoubleMat','mattype',1); - WriteData(fid,prefix,'object',self,'fieldname','tau','format','Double'); - WriteData(fid,prefix,'object',self,'fieldname','beta','format','DoubleMat','mattype',1); - WriteData(fid,prefix,'object',self,'fieldname','phi','format','Double'); - WriteData(fid,prefix,'object',self,'fieldname','alpha','format','Integer'); - WriteData(fid,prefix,'object',self,'fieldname','robin','format','Boolean'); - WriteData(fid,prefix,'object',self,'fieldname','seed','format','Integer'); - + WriteData(fid,prefix,'object',self,'fieldname','kappa','format','DoubleMat','mattype',1); + WriteData(fid,prefix,'object',self,'fieldname','tau','format','Double'); + WriteData(fid,prefix,'object',self,'fieldname','beta','format','DoubleMat','mattype',1); + WriteData(fid,prefix,'object',self,'fieldname','phi','format','Double'); + WriteData(fid,prefix,'object',self,'fieldname','alpha','format','Integer'); + WriteData(fid,prefix,'object',self,'fieldname','robin','format','Boolean'); + WriteData(fid,prefix,'object',self,'fieldname','seed','format','Integer'); + %process requested outputs outputs = self.requested_outputs; @@ -116,13 +107,22 @@ WriteData(fid,prefix,'data',outputs,'name','md.sampling.requested_outputs','format','StringArray'); end % }}} + function md = setparameters(self,md,lc,sigma) % {{{ + + nu = self.alpha-1; + KAPPA = sqrt(8*nu)/lc; + TAU = sqrt(gamma(nu)/(gamma(self.alpha)*(4*pi)*KAPPA^(2*nu)*sigma^2)); + md.sampling.kappa = KAPPA*ones(md.mesh.numberofvertices,1); + md.sampling.tau = TAU; + + end % }}} function savemodeljs(self,fid,modelname) % {{{ - - writejsdouble(fid,[modelname '.sampling.kappa'],self.kappa); - writejsdouble(fid,[modelname '.sampling.tau'],self.tau); - writejsdouble(fid,[modelname '.sampling.beta'],self.beta); - writejsdouble(fid,[modelname '.sampling.phi'],self.beta); - writejsdouble(fid,[modelname '.sampling.alpha'],self.alpha); - writejsdouble(fid,[modelname '.sampling.robin'],self.robin); - writejsdouble(fid,[modelname '.sampling.seed'],self.seed); + + writejsdouble(fid,[modelname '.sampling.kappa'],self.kappa); + writejsdouble(fid,[modelname '.sampling.tau'],self.tau); + writejsdouble(fid,[modelname '.sampling.beta'],self.beta); + writejsdouble(fid,[modelname '.sampling.phi'],self.beta); + writejsdouble(fid,[modelname '.sampling.alpha'],self.alpha); + writejsdouble(fid,[modelname '.sampling.robin'],self.robin); + writejsdouble(fid,[modelname '.sampling.seed'],self.seed); writejscellstring(fid,[modelname '.sampling.requested_outputs'],self.requested_outputs); Index: /issm/trunk-jpl/src/m/classes/sampling.py =================================================================== --- /issm/trunk-jpl/src/m/classes/sampling.py (revision 26300) +++ /issm/trunk-jpl/src/m/classes/sampling.py (revision 26301) @@ -1,3 +1,5 @@ import numpy as np + +from math import * from checkfield import checkfield @@ -14,9 +16,8 @@ """ - def __init__(self): # {{{ - self.kappa = float('NaN') + def __init__(self, *args): #{{{ + self.kappa = np.nan self.tau = 0 - self.beta = float('NaN') - self.hydrostatic_adjustment = 0 + self.beta = np.nan self.phi = 0 self.alpha = 0 @@ -25,70 +26,85 @@ self.requested_outputs = [] - # Set defaults - self.setdefaultparameters() - + if len(args) == 0: + self.setdefaultparameters() + else: + raise RuntimeError('constructor not supported') #}}} - def __repr__(self): # {{{ + def __repr__(self): #{{{ s = ' Sampling parameters::\n' - s += '{}\n'.format(fielddisplay(self,'kappa','coefficient of the identity operator')) - s += '{}\n'.format(fielddisplay(self,'tau','scaling coefficient of the solution (default 1.0)')) - s += '{}\n'.format(fielddisplay(self,'alpha','exponent in PDE operator, (default 2.0, BiLaplacian covariance operator)')) - s += '{}\n'.format(disp(sprintf('\n %s','Parameters of Robin boundary conditions nabla () \cdot normvec + beta ():'))) - s += '{}\n'.format(fielddisplay(self,'beta','Coefficient in Robin boundary conditions (to be defined for robin = 1)')) - s += '{}\n'.format(fielddisplay(self,'phi','Temporal correlation factor (|phi|<1 for stationary process, phi = 1 for random walk process) (default 0)')) - s += '{}\n'.format(fielddisplay(self,'seed','Seed for pseudorandom number generator (given seed if >=0 and random seed if <0) (default -1)')) + s += ' Parameters of PDE operator (kappa^2 I-Laplacian)^(alpha/2)(tau):\n' + s += '{}\n'.format(fielddisplay(self, 'kappa', 'coefficient of the identity operator')) + s += '{}\n'.format(fielddisplay(self, 'tau', 'scaling coefficient of the solution (default: 1.0)')) + s += '{}\n'.format(fielddisplay(self, 'alpha', 'exponent in PDE operator, (default: 2.0, BiLaplacian covariance operator)')) + + s += ' Parameters of Robin boundary conditions nabla () \cdot normvec + beta ():\n' + s += '{}\n'.format(fielddisplay(self, 'robin', 'Apply Robin boundary conditions (1 if applied and 0 for homogenous Neumann boundary conditions) (default: 0)')) + s += '{}\n'.format(fielddisplay(self, 'beta', 'Coefficient in Robin boundary conditions (to be defined for robin = 1)')) + + s += ' Parameters for first-order autoregressive process (X_t = phi X_{t-1} + noise) (if transient):\n' + s += '{}\n'.format(fielddisplay(self, 'phi', 'Temporal correlation factor (|phi|<1 for stationary process, phi = 1 for random walk process) (default 0)')) + + s += ' Other parameters of stochastic sampler:\n' + s += '{}\n'.format(fielddisplay(self, 'seed', 'Seed for pseudorandom number generator (given seed if >=0 and random seed if <0) (default: -1)')) s += '{}\n'.format(fielddisplay(self, 'requested_outputs', 'additional outputs requested (not implemented yet)')) + return s #}}} - def defaultoutputs(self, md): # {{{ - return [] - - #}}} - def setdefaultparameters(self): # {{{ + def setdefaultparameters(self): #{{{ # Scaling coefficient self.tau = 1 + # Apply Robin boundary conditions self.robin = 0 + # Temporal correlation factor self.phi = 0 - # Exponent in fraction SPDE (default=2, biLaplacian covariance operator) + + # Exponent in fraction SPDE (default: 2, biLaplacian covariance operator) self.alpha = 2 - # Seed for pseudorandom number generator (default -1 for random seed) - self.alpha = -1 + + # Seed for pseudorandom number generator (default: -1, for random seed) + self.seed = -1 + # Default output self.requested_outputs = ['default'] + return self #}}} - def checkconsistency(self, md, solution, analyses): # {{{ - # Early return + def defaultoutputs(self, md): #{{{ + return [] + #}}} + + def checkconsistency(self, md, solution, analyses): #{{{ if ('SamplingAnalysis' not in analyses): return md - md = checkfield(md,'fieldname','sampling.kappa','NaN',1,'Inf',1,'size',[md.mesh.numberofvertices],'>',0) - md = checkfield(md,'fieldname','sampling.tau','NaN',1,'Inf',1,'numel',1,'>',0) - md = checkfield(md,'fieldname','sampling.robin','numel',1,'values',[0, 1]) + md = checkfield(md, 'fieldname', 'sampling.kappa', 'NaN', 1, 'Inf', 1, 'size', [md.mesh.numberofvertices], '>', 0) + md = checkfield(md, 'fieldname', 'sampling.tau', 'NaN', 1, 'Inf', 1, 'numel', 1, '>', 0) + md = checkfield(md, 'fieldname', 'sampling.robin', 'numel', 1, 'values', [0, 1]) if md.sampling.robin: - md = checkfield(md,'fieldname','sampling.beta','NaN',1,'Inf',1,'size',[md.mesh.numberofvertices],'>',0) - md = checkfield(md,'fieldname','sampling.phi','NaN',1,'Inf',1,'numel',1,'>=',0) - md = checkfield(md,'fieldname','sampling.alpha','NaN',1,'Inf',1,'numel',1,'>',0) - md = checkfield(md,'fieldname','sampling.seed','NaN',1,'Inf',1,'numel',1) - md = checkfield(md,'fieldname','sampling.requested_outputs','stringrow',1) + md = checkfield(md, 'fieldname', 'sampling.beta', 'NaN', 1, 'Inf', 1, 'size', [md.mesh.numberofvertices], '>', 0) + end + md = checkfield(md, 'fieldname', 'sampling.phi', 'NaN', 1, 'Inf', 1, 'numel', 1, '>=', 0) + md = checkfield(md, 'fieldname', 'sampling.alpha', 'NaN', 1, 'Inf', 1, 'numel', 1, '>', 0) + md = checkfield(md, 'fieldname', 'sampling.seed', 'NaN', 1, 'Inf', 1, 'numel', 1) + md = checkfield(md, 'fieldname', 'sampling.requested_outputs', 'stringrow', 1) return md - # }}} + #}}} - def marshall(self, prefix, md, fid): # {{{ - WriteData(fid,prefix,'object',self,'fieldname','kappa','format','DoubleMat','mattype',1) - WriteData(fid,prefix,'object',self,'fieldname','tau','format','Double') - WriteData(fid,prefix,'object',self,'fieldname','beta','format','DoubleMat','mattype',1) - WriteData(fid,prefix,'object',self,'fieldname','phi','format','Double') - WriteData(fid,prefix,'object',self,'fieldname','alpha','format','Integer') - WriteData(fid,prefix,'object',self,'fieldname','robin','format','Boolean') - WriteData(fid,prefix,'object',self,'fieldname','seed','format','Integer') + def marshall(self, prefix, md, fid): #{{{ + WriteData(fid, prefix, 'object', self, 'fieldname', 'kappa', 'format', 'DoubleMat', 'mattype', 1) + WriteData(fid, prefix, 'object', self, 'fieldname', 'tau', 'format', 'Double') + WriteData(fid, prefix, 'object', self, 'fieldname', 'beta', 'format', 'DoubleMat', 'mattype', 1) + WriteData(fid, prefix, 'object', self, 'fieldname', 'phi', 'format', 'Double') + WriteData(fid, prefix, 'object', self, 'fieldname', 'alpha', 'format', 'Integer') + WriteData(fid, prefix, 'object', self, 'fieldname', 'robin', 'format', 'Boolean') + WriteData(fid, prefix, 'object', self, 'fieldname', 'seed', 'format', 'Integer') - #process requested outputs + # Process requested outputs outputs = self.requested_outputs indices = [i for i, x in enumerate(outputs) if x == 'default'] @@ -97,4 +113,13 @@ outputs = outputscopy WriteData(fid, prefix, 'data', outputs, 'name', 'md.sampling.requested_outputs', 'format', 'StringArray') + #}}} - # }}} + def setparameters(self, md, lc, sigma): #{{{ + nu = self.alpha - 1 + KAPPA = pow((8 * nu), 0.5) / lc + TAU = pow((math.gamma(nu) / math.gamma(self.alpha) * (4 * np.pi) * pow(KAPPA, 2 * nu) * pow(sigma, 2)), 0.5) + md.sampling.kappa = KAPPA * np.ones((md.mesh.numberofvertices, 1)) + md.sampling.tau = TAU + + return md + #}}} Index: /issm/trunk-jpl/src/m/classes/solidearth.m =================================================================== --- /issm/trunk-jpl/src/m/classes/solidearth.m (revision 26300) +++ /issm/trunk-jpl/src/m/classes/solidearth.m (revision 26301) @@ -45,4 +45,22 @@ end end % }}} + function disp(self) % {{{ + disp(sprintf(' solidearth inputs, forcings and settings:')); + + fielddisplay(self,'planetradius','planet radius [m]'); + fielddisplay(self,'transitions','indices into parts of the mesh that will be icecaps'); + fielddisplay(self,'requested_outputs','additional outputs requested'); + fielddisplay(self,'partitionice','ice partition vector for barystatic contribution'); + fielddisplay(self,'partitionhydro','hydro partition vector for barystatic contribution'); + fielddisplay(self,'partitionocean','ocean partition vector for barystatic contribution'); + if isempty(self.external), fielddisplay(self,'external','external solution, of the type solidearthsolution'); end + self.settings.disp(); + self.lovenumbers.disp(); + self.rotational.disp(); + if ~isempty(self.external), + self.external.disp(); + end + + end % }}} function self = setdefaultparameters(self,planet) % {{{ @@ -86,22 +104,4 @@ function list=defaultoutputs(self,md) % {{{ list = {'Sealevel'}; - end % }}} - function disp(self) % {{{ - disp(sprintf(' solidearth inputs, forcings and settings:')); - - fielddisplay(self,'planetradius','planet radius [m]'); - fielddisplay(self,'transitions','indices into parts of the mesh that will be icecaps'); - fielddisplay(self,'requested_outputs','additional outputs requested'); - fielddisplay(self,'partitionice','ice partition vector for barystatic contribution'); - fielddisplay(self,'partitionhydro','hydro partition vector for barystatic contribution'); - fielddisplay(self,'partitionocean','ocean partition vector for barystatic contribution'); - if isempty(self.external), fielddisplay(self,'external','external solution, of the type solidearthsolution'); end - self.settings.disp(); - self.lovenumbers.disp(); - self.rotational.disp(); - if ~isempty(self.external), - self.external.disp(); - end - end % }}} function marshall(self,prefix,md,fid) % {{{ Index: /issm/trunk-jpl/src/m/classes/solidearthsettings.m =================================================================== --- /issm/trunk-jpl/src/m/classes/solidearthsettings.m (revision 26300) +++ /issm/trunk-jpl/src/m/classes/solidearthsettings.m (revision 26301) @@ -92,4 +92,25 @@ end % }}} + function disp(self) % {{{ + disp(sprintf(' solidearth settings:')); + + fielddisplay(self,'reltol','sea level change relative convergence criterion (default, NaN: not applied)'); + fielddisplay(self,'abstol','sea level change absolute convergence criterion(default, NaN: not applied)'); + fielddisplay(self,'maxiter','maximum number of nonlinear iterations'); + fielddisplay(self,'grdocean','does this planet have an ocean, if set to 1: global water mass is conserved in GRD module (default: 1)'); + fielddisplay(self,'ocean_area_scaling','correction for model representation of ocean area (default: No correction)'); + fielddisplay(self,'sealevelloading','enables surface loading from sea-level change (default: 1)'); + fielddisplay(self,'isgrd','compute GRD patterns (default: 1)'); + fielddisplay(self,'compute_bp_grd','compute GRD patterns for bottom pressure loads (default: 1)'); + fielddisplay(self,'runfrequency','how many time steps we skip before we run solidearthsettings solver during transient (default: 1)'); + fielddisplay(self,'selfattraction','enables surface mass load to perturb the gravity field'); + fielddisplay(self,'elastic','enables elastic deformation from surface loading'); + fielddisplay(self,'viscous','enables viscous deformation from surface loading'); + fielddisplay(self,'rotation','enables polar motion to feedback on the GRD fields'); + fielddisplay(self,'degacc','accuracy (default: .01 deg) for numerical discretization of the Green''s functions'); + fielddisplay(self,'timeacc','time accuracy (default: 1 yr) for numerical discretization of the Green''s functions'); + fielddisplay(self,'grdmodel','type of deformation model, 0 for no GRD, 1 for spherical GRD model (SESAW model), 2 for half-space planar GRD (visco-elastic model from Ivins)'); + fielddisplay(self,'cross_section_shape','1: square-edged (default). 2: elliptical. See iedge in GiaDeflectionCore'); + end % }}} function md = checkconsistency(self,md,solution,analyses) % {{{ @@ -135,25 +156,4 @@ end % }}} - function disp(self) % {{{ - disp(sprintf(' solidearth settings:')); - - fielddisplay(self,'reltol','sea level change relative convergence criterion (default, NaN: not applied)'); - fielddisplay(self,'abstol','sea level change absolute convergence criterion(default, NaN: not applied)'); - fielddisplay(self,'maxiter','maximum number of nonlinear iterations'); - fielddisplay(self,'grdocean','does this planet have an ocean, if set to 1: global water mass is conserved in GRD module (default: 1)'); - fielddisplay(self,'ocean_area_scaling','correction for model representation of ocean area (default: No correction)'); - fielddisplay(self,'sealevelloading','enables surface loading from sea-level change (default: 1)'); - fielddisplay(self,'isgrd','compute GRD patterns (default: 1)'); - fielddisplay(self,'compute_bp_grd','compute GRD patterns for bottom pressure loads (default: 1)'); - fielddisplay(self,'runfrequency','how many time steps we skip before we run solidearthsettings solver during transient (default: 1)'); - fielddisplay(self,'selfattraction','enables surface mass load to perturb the gravity field'); - fielddisplay(self,'elastic','enables elastic deformation from surface loading'); - fielddisplay(self,'viscous','enables viscous deformation from surface loading'); - fielddisplay(self,'rotation','enables polar motion to feedback on the GRD fields'); - fielddisplay(self,'degacc','accuracy (default: .01 deg) for numerical discretization of the Green''s functions'); - fielddisplay(self,'timeacc','time accuracy (default: 1 yr) for numerical discretization of the Green''s functions'); - fielddisplay(self,'grdmodel','type of deformation model, 0 for no GRD, 1 for spherical GRD model (SESAW model), 2 for half-space planar GRD (visco-elastic model from Ivins)'); - fielddisplay(self,'cross_section_shape','1: square-edged (default). 2: elliptical. See iedge in GiaDeflectionCore'); - end % }}} function marshall(self,prefix,md,fid) % {{{ WriteData(fid,prefix,'object',self,'fieldname','reltol','name','md.solidearth.settings.reltol','format','Double'); Index: /issm/trunk-jpl/src/m/classes/solidearthsettings.py =================================================================== --- /issm/trunk-jpl/src/m/classes/solidearthsettings.py (revision 26300) +++ /issm/trunk-jpl/src/m/classes/solidearthsettings.py (revision 26301) @@ -121,7 +121,7 @@ if self.grdmodel == 2: raise RuntimeException('model requires a 2D mesh to run gia Ivins computations (change mesh from mesh3dsurface to mesh2d)') - else: - if self.grdmodel == 1: - raise RuntimeException('model requires a 3D surface mesh to run GRD computations (change mesh from mesh2d to mesh3dsurface)') + else: + if self.grdmodel == 1: + raise RuntimeException('model requires a 3D surface mesh to run GRD computations (change mesh from mesh2d to mesh3dsurface)') if self.sealevelloading and not self.grdocean: raise RuntimeException('solidearthsettings checkconsistency error message: need grdocean on if sealevelloading flag is set') @@ -134,22 +134,22 @@ def marshall(self, prefix, md, fid): #{{{ - WriteData(fid, prefix, 'object', self, 'fieldname', 'reltol', 'name', 'md.solidearth.settings.reltol', 'format', 'Double') - WriteData(fid, prefix, 'object', self, 'fieldname', 'abstol', 'name', 'md.solidearth.settings.abstol', 'format', 'Double') - WriteData(fid, prefix, 'object', self, 'fieldname', 'maxiter', 'name', 'md.solidearth.settings.maxiter', 'format', 'Integer') - WriteData(fid, prefix, 'object', self, 'fieldname', 'selfattraction', 'name', 'md.solidearth.settings.selfattraction', 'format', 'Boolean') - WriteData(fid, prefix, 'object', self, 'fieldname', 'elastic', 'name', 'md.solidearth.settings.elastic', 'format', 'Boolean') - WriteData(fid, prefix, 'object', self, 'fieldname', 'viscous', 'name', 'md.solidearth.settings.viscous', 'format', 'Boolean') - WriteData(fid, prefix, 'object', self, 'fieldname', 'rotation', 'name', 'md.solidearth.settings.rotation', 'format', 'Boolean') - WriteData(fid, prefix, 'object', self, 'fieldname', 'grdocean', 'name', 'md.solidearth.settings.grdocean', 'format', 'Boolean') - WriteData(fid, prefix, 'object', self, 'fieldname', 'ocean_area_scaling', 'name', 'md.solidearth.settings.ocean_area_scaling', 'format', 'Integer') - WriteData(fid, prefix, 'object', self, 'fieldname', 'runfrequency', 'name', 'md.solidearth.settings.runfrequency', 'format', 'Integer') - WriteData(fid, prefix, 'object', self, 'fieldname', 'degacc', 'name', 'md.solidearth.settings.degacc', 'format', 'Double') - WriteData(fid, prefix, 'object', self, 'fieldname', 'timeacc', 'name', 'md.solidearth.settings.timeacc', 'format', 'Double', 'scale', md.constants.yts) - WriteData(fid, prefix, 'object', self, 'fieldname', 'horiz', 'name', 'md.solidearth.settings.horiz', 'format', 'Integer') - WriteData(fid, prefix, 'object', self, 'fieldname', 'sealevelloading', 'name', 'md.solidearth.settings.sealevelloading', 'format', 'Integer') - WriteData(fid, prefix, 'object', self, 'fieldname','isgrd', 'name', 'md.solidearth.settings.isgrd', 'format', 'Integer') - WriteData(fid, prefix, 'object', self, 'fieldname', 'compute_bp_grd', 'name', 'md.solidearth.settings.compute_bp_grd', 'format', 'Integer') - WriteData(fid, prefix, 'object', self, 'fieldname', 'grdmodel', 'name', 'md.solidearth.settings.grdmodel', 'format', 'Integer') - WriteData(fid, prefix, 'object', self, 'fieldname', 'cross_section_shape', 'name', 'md.solidearth.settings.cross_section_shape', 'format', 'Integer') + WriteData(fid, prefix, 'object', self, 'fieldname', 'reltol', 'name', 'md.solidearth.settings.reltol', 'format', 'Double'); + WriteData(fid, prefix, 'object', self, 'fieldname', 'abstol', 'name', 'md.solidearth.settings.abstol', 'format', 'Double'); + WriteData(fid, prefix, 'object', self, 'fieldname', 'maxiter', 'name', 'md.solidearth.settings.maxiter', 'format', 'Integer'); + WriteData(fid, prefix, 'object', self, 'fieldname', 'selfattraction', 'name', 'md.solidearth.settings.selfattraction', 'format', 'Boolean'); + WriteData(fid, prefix, 'object', self, 'fieldname', 'elastic', 'name', 'md.solidearth.settings.elastic', 'format', 'Boolean'); + WriteData(fid, prefix, 'object', self, 'fieldname', 'viscous', 'name', 'md.solidearth.settings.viscous', 'format', 'Boolean'); + WriteData(fid, prefix, 'object', self, 'fieldname', 'rotation', 'name', 'md.solidearth.settings.rotation', 'format', 'Boolean'); + WriteData(fid, prefix, 'object', self, 'fieldname', 'grdocean', 'name', 'md.solidearth.settings.grdocean', 'format', 'Boolean'); + WriteData(fid, prefix, 'object', self, 'fieldname', 'ocean_area_scaling', 'name', 'md.solidearth.settings.ocean_area_scaling', 'format', 'Boolean'); + WriteData(fid, prefix, 'object', self, 'fieldname', 'runfrequency', 'name', 'md.solidearth.settings.runfrequency', 'format', 'Integer'); + WriteData(fid, prefix, 'object', self, 'fieldname', 'degacc', 'name', 'md.solidearth.settings.degacc', 'format', 'Double'); + WriteData(fid, prefix, 'object', self, 'fieldname', 'timeacc', 'name', 'md.solidearth.settings.timeacc', 'format', 'Double', 'scale',md.constants.yts); + WriteData(fid, prefix, 'object', self, 'fieldname', 'horiz', 'name', 'md.solidearth.settings.horiz', 'format', 'Integer'); + WriteData(fid, prefix, 'object', self, 'fieldname', 'sealevelloading', 'name', 'md.solidearth.settings.sealevelloading', 'format', 'Integer'); + WriteData(fid, prefix, 'object', self, 'fieldname', 'isgrd', 'name', 'md.solidearth.settings.isgrd', 'format', 'Integer'); + WriteData(fid, prefix, 'object', self, 'fieldname', 'compute_bp_grd', 'name', 'md.solidearth.settings.compute_bp_grd', 'format', 'Integer'); + WriteData(fid, prefix, 'object', self, 'fieldname', 'grdmodel', 'name', 'md.solidearth.settings.grdmodel', 'format', 'Integer'); + WriteData(fid, prefix, 'object', self, 'fieldname', 'cross_section_shape', 'name', 'md.solidearth.settings.cross_section_shape', 'format', 'Integer'); #}}} Index: /issm/trunk-jpl/src/m/classes/solidearthsolution.py =================================================================== --- /issm/trunk-jpl/src/m/classes/solidearthsolution.py (revision 26301) +++ /issm/trunk-jpl/src/m/classes/solidearthsolution.py (revision 26301) @@ -0,0 +1,80 @@ +import numpy as np + +from checkfield import checkfield +from fielddisplay import fielddisplay +from WriteData import WriteData + + +class solidearthsolution(object): + """SOLIDEARTHSOLUTION class definition + + Usage: + solidearthsolution = solidearthsolution() + """ + + def __init__(self, *args): #{{{ + self.displacementeast = None + self.displacementnorth = None + self.displacementup = None + self.geoid = None + + if len(args) == 0: + self.setdefaultparameters() + else: + raise RuntimeError('constructor not supported') + #}}} + + def __repr__(self): #{{{ + s = ' units for time series is (yr)\n' + s += '{}\n'.format(fielddisplay(self, 'displacementeast', 'solid-Earth Eastwards bedrock displacement series (m)')) + s += '{}\n'.format(fielddisplay(self, 'displacementnorth', 'solid-Earth Northwards bedrock displacement time series (m)')) + s += '{}\n'.format(fielddisplay(self, 'displacementup', 'solid-Earth bedrock uplift time series (m)')) + s += '{}\n'.format(fielddisplay(self, 'geoid', 'solid-Earth geoid time series (m)')) + + return s + #}}} + + def setdefaultparameters(self): #{{{ + self.displacementeast = [] + self.displacementnorth = [] + self.displacementup = [] + self.geoid = [] + #}}} + + def checkconsistency(self, md, solution, analyses): #{{{ + md = checkfield(md, 'fieldname', 'solidearth.external.displacementeast', 'Inf', 1, 'timeseries', 1) + md = checkfield(md, 'fieldname', 'solidearth.external.displacementnorth', 'Inf', 1, 'timeseries', 1) + md = checkfield(md, 'fieldname', 'solidearth.external.displacementup', 'Inf', 1, 'timeseries', 1) + md = checkfield(md, 'fieldname', 'solidearth.external.geoid', 'Inf', 1, 'timeseries', 1) + #}}} + + def marshall(self, prefix, md, fid): #{{{ + yts = md.constants.yts + + # Transform our time series into time series rates + if np.shape(self.displacementeast, 1) == 1: + print('External solidearthsolution warning: only one time step provided, assuming the values are rates per year') + displacementeast_rate = np.append(np.array(self.displacementeast).reshape(-1, 1), 0) + displacementnorth_rate = np.append(np.array(self.displacementnorth).reshape(-1, 1), 0) + displacementup_rate = np.append(np.array(self.displacementup).reshape(-1, 1), 0) + geoid_rate = np.append(np.array(self.geoid).reshape(-1, 1), 0) + else: + time = self.displacementeast[-1, :] + dt = np.diff(time, 1, 1) + displacementeast_rate = np.diff(self.displacementeast[0:-2, :], 1, 1) / dt + displacementeast_rate.append(time[0:-2]) + displacementnorth_rate = np.diff(self.displacementnorth[0:-2, :], 1, 1) / dt + displacementnorth_rate.append(time[0:-2]) + displacementup_rate = np.diff(self.displacementup[0:-2, :], 1, 1) / dt + displacementup_rate.append(time[0:-2]) + geoid_rate = np.diff(self.geoid[0:-2, :], 1, 1) / dt + geoid_rate.append(time[0:-2]) + WriteData(fid, prefix, 'object', self, 'fieldname', 'displacementeast', 'data', displacementeast_rate, 'format', 'DoubleMat', 'name', 'md.solidearth.external.displacementeast', 'mattype', 1, 'scale', 1 / yts,'timeserieslength', md.mesh.numberofvertices + 1, 'yts', yts); + WriteData(fid, prefix, 'object', self, 'fieldname', 'displacementup', 'data', displacementup_rate,'format', 'DoubleMat', 'name', 'md.solidearth.external.displacementup', 'mattype', 1, 'scale', 1 / yts,'timeserieslength', md.mesh.numberofvertices + 1, 'yts', yts); + WriteData(fid, prefix, 'object', self, 'fieldname', 'displacementnorth', 'data', displacementnorth_rate,'format', 'DoubleMat', 'name', 'md.solidearth.external.displacementnorth', 'mattype', 1, 'scale', 1 / yts,'timeserieslength', md.mesh.numberofvertices + 1, 'yts', yts); + WriteData(fid, prefix, 'object', self, 'fieldname', 'geoid', 'data', geoid_rate,'format', 'DoubleMat', 'name', 'md.solidearth.external.geoid', 'mattype', 1, 'scale', 1 / yts,'timeserieslength', md.mesh.numberofvertices + 1, 'yts', yts); + #}}} + + def extrude(self, md): #{{{ + return self + #}}} Index: /issm/trunk-jpl/src/m/miscellaneous/MatlabFuncs.py =================================================================== --- /issm/trunk-jpl/src/m/miscellaneous/MatlabFuncs.py (revision 26300) +++ /issm/trunk-jpl/src/m/miscellaneous/MatlabFuncs.py (revision 26301) @@ -1,18 +1,69 @@ -def oshostname(): - import socket +def acosd(X): #{{{ + """ function acosd - Inverse cosine in degrees - return socket.gethostname() + Usage: + Y = acosd(X) + """ + import numpy as np + + return np.degrees(np.arccos(X)) +#}}} + +def asind(X): #{{{ + """ function asind - Inverse sine in degrees + + Usage: + Y = asind(X) + """ + import numpy as np + + return np.degrees(np.arcsin(X)) +#}}} + +def atand(X): #{{{ + """ function atand - Inverse tangent in degrees + + Usage: + Y = atand(X) + """ + import numpy as np + + return np.degrees(np.arctan(X)) +#}}} -def ispc(): - import platform +def atan2d(Y, X): #{{{ + """ function atan2d - Four-quadrant inverse tangent in degrees - if 'Windows' in platform.system(): - return True + Usage: + D = atan2d(Y, X) + """ + import numpy as np + + return np.degrees(np.arctan2(Y, X)) +#}}} + +def det(a): #{{{ + if a.shape == (1, ): + return a[0] + elif a.shape == (1, 1): + return a[0, 0] + elif a.shape == (2, 2): + return a[0, 0] * a[1, 1] - a[0, 1] * a[1, 0] else: - return False + raise TypeError("MatlabFunc.det only implemented for shape (2, 2), not for shape %s." % str(a.shape)) +#}}} +def heaviside(x): #{{{ + import numpy as np -def ismac(): + y = np.zeros_like(x) + y[np.nonzero(x > 0.)] = 1. + y[np.nonzero(x == 0.)] = 0.5 + + return y +#}}} + +def ismac(): #{{{ import platform @@ -21,39 +72,7 @@ else: return False +#}}} - -def strcmp(s1, s2): - - if s1 == s2: - return True - else: - return False - - -def strncmp(s1, s2, n): - - if s1[0:n] == s2[0:n]: - return True - else: - return False - - -def strcmpi(s1, s2): - - if s1.lower() == s2.lower(): - return True - else: - return False - - -def strncmpi(s1, s2, n): - - if s1.lower()[0:n] == s2.lower()[0:n]: - return True - else: - return False - - -def ismember(a, s): +def ismember(a, s): #{{{ import numpy as np @@ -66,5 +85,4 @@ if not isinstance(a, np.ndarray): b = [item in s for item in a] - else: if not isinstance(s, np.ndarray): @@ -72,24 +90,26 @@ for i, item in enumerate(a.flat): b.flat[i] = item in s - else: b = np.in1d(a.flat, s.flat).reshape(a.shape) return b +#}}} +def ispc(): #{{{ + import platform -def det(a): + if 'Windows' in platform.system(): + return True + else: + return False +#}}} - if a.shape == (1, ): - return a[0] - elif a.shape == (1, 1): - return a[0, 0] - elif a.shape == (2, 2): - return a[0, 0] * a[1, 1] - a[0, 1] * a[1, 0] - else: - raise TypeError("MatlabFunc.det only implemented for shape (2, 2), not for shape %s." % str(a.shape)) +def oshostname(): #{{{ + import socket + return socket.gethostname() +#}}} -def sparse(ivec, jvec, svec, m=0, n=0, nzmax=0): +def sparse(ivec, jvec, svec, m=0, n=0, nzmax=0): #{{{ import numpy as np @@ -105,12 +125,31 @@ return a +#}}} +def strcmp(s1, s2): #{{{ + if s1 == s2: + return True + else: + return False +#}}} -def heaviside(x): - import numpy as np +def strcmpi(s1, s2): #{{{ + if s1.lower() == s2.lower(): + return True + else: + return False +#}}} - y = np.zeros_like(x) - y[np.nonzero(x > 0.)] = 1. - y[np.nonzero(x == 0.)] = 0.5 +def strncmp(s1, s2, n): #{{{ + if s1[0:n] == s2[0:n]: + return True + else: + return False +#}}} - return y +def strncmpi(s1, s2, n): #{{{ + if s1.lower()[0:n] == s2.lower()[0:n]: + return True + else: + return False +#}}} Index: /issm/trunk-jpl/src/m/miscellaneous/PythonFuncs.py =================================================================== --- /issm/trunk-jpl/src/m/miscellaneous/PythonFuncs.py (revision 26300) +++ /issm/trunk-jpl/src/m/miscellaneous/PythonFuncs.py (revision 26301) @@ -2,6 +2,5 @@ -def logical_and_n(*arg): - +def logical_and_n(*arg): #{{{ if len(arg): result = arg[0] @@ -9,11 +8,9 @@ result = np.logical_and(result, item) return result - else: return None +#}}} - -def logical_or_n(*arg): - +def logical_or_n(*arg): #{{{ if len(arg): result = arg[0] @@ -21,5 +18,5 @@ result = np.logical_or(result, item) return result - else: return None +#}}} Index: /issm/trunk-jpl/src/m/solve/marshall.m =================================================================== --- /issm/trunk-jpl/src/m/solve/marshall.m (revision 26300) +++ /issm/trunk-jpl/src/m/solve/marshall.m (revision 26301) @@ -45,5 +45,5 @@ % Uncomment the following to make a copy of the binary input file for debugging -% purposes (can be fed into scripts/ReadBin.py). +% purposes (can be fed into scripts/BinRead.py). % copyfile([md.miscellaneous.name '.bin'], [md.miscellaneous.name '.m.bin']) Index: /issm/trunk-jpl/src/m/solve/marshall.py =================================================================== --- /issm/trunk-jpl/src/m/solve/marshall.py (revision 26300) +++ /issm/trunk-jpl/src/m/solve/marshall.py (revision 26301) @@ -46,5 +46,5 @@ # Uncomment the following to make a copy of the binary input file for - # debugging purposes (can be fed into scripts/ReadBin.py). + # debugging purposes (can be fed into scripts/BinRead.py). # copy_cmd = "cp {}.bin {}.py.bin".format(md.miscellaneous.name, md.miscellaneous.name) # subprocess.call(copy_cmd, shell=True) Index: /issm/trunk-jpl/test/NightlyRun/test2002.m =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test2002.m (revision 26300) +++ /issm/trunk-jpl/test/NightlyRun/test2002.m (revision 26301) @@ -15,4 +15,5 @@ md.masstransport.spcthickness=[md.geometry.thickness;0]; md.smb.mass_balance=zeros(md.mesh.numberofvertices,1); + %antarctica xe=md.mesh.x(md.mesh.elements)*[1;1;1]/3; @@ -26,4 +27,5 @@ md.masstransport.spcthickness(md.mesh.elements(pos,:))= md.masstransport.spcthickness(md.mesh.elements(pos,:))-100; posant=pos; + %greenland pos=find(late>60 & late<90 & longe>-75 & longe<-15); @@ -96,4 +98,5 @@ Seustatic=md.results.TransientSolution.Sealevel; Beustatic=md.results.TransientSolution.Bed; +pause %eustatic + selfattraction run: Index: /issm/trunk-jpl/test/NightlyRun/test2002.py =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test2002.py (revision 26300) +++ /issm/trunk-jpl/test/NightlyRun/test2002.py (revision 26301) @@ -1,4 +1,6 @@ #Test Name: EarthSlc import numpy as np + +from MatlabFuncs import * from gmshplanet import * @@ -12,97 +14,136 @@ -#mesh earth: +# Mesh earth +# +# NOTE: In MATLAB, we currently use cached mesh to account for differences in +# mesh generated under Linux versus under macOS +# md = model() -md.mesh = gmshplanet('radius', 6.371012 * 1e3, 'resolution', 700.) #700 km resolution mesh +md.mesh = gmshplanet('radius', 6.371012 * 1e3, 'resolution', 700.) # 700 km resolution mesh -#parameterize solidearth solution: -#solidearth loading: -md.solidearth.surfaceload.icethicknesschange = np.zeros((md.mesh.numberofelements, 1)) -md.solidearth.initialsealevel = np.zeros((md.mesh.numberofvertices, 1)) -md.dsl.global_average_thermosteric_sea_level_change = np.zeros((2, 1)) -md.dsl.sea_surface_height_change_above_geoid = np.zeros((md.mesh.numberofvertices + 1, 1)) -md.dsl.sea_water_pressure_change_at_sea_floor = np.zeros((md.mesh.numberofvertices + 1, 1)) +# Geometry for the bed, arbitrary thickness of 100 +md.geometry.bed = np.zeros((md.mesh.numberofvertices, 1)) +md.geometry.base = md.geometry.bed +md.geometry.thickness = 100 * np.ones((md.mesh.numberofvertices, 1)) +md.geometry.surface = md.geometry.bed + md.geometry.thickness -#antarctica -late = md.mesh.lat[md.mesh.elements - 1].sum(axis=1) / 3 -longe = md.mesh.long[md.mesh.elements - 1].sum(axis=1) / 3 +# Solidearth loading #{{{ +md.masstransport.spcthickness = np.append(md.geometry.thickness, 0) +md.smb.mass_balance = np.zeros((md.mesh.numberofvertices, 1)) + +# Antarctica +xe = md.mesh.x[md.mesh.elements - 1].sum(axis=1) / 3 +ye = md.mesh.y[md.mesh.elements - 1].sum(axis=1) / 3 +ze = md.mesh.z[md.mesh.elements - 1].sum(axis=1) / 3 +re = pow((pow(xe, 2) + pow(ye, 2) + pow(ze, 2)), 0.5) + +late = asind(ze / re) +longe = atan2d(ye, xe) pos = np.where(late < -80)[0] -md.solidearth.surfaceload.icethicknesschange[pos] = -100 -#greenland -pos = np.where(np.logical_and.reduce((late > 70, late < 80, longe > -60, longe < -30)))[0] -md.solidearth.surfaceload.icethicknesschange[pos] = -100 +md.masstransport.spcthickness[md.mesh.elements[pos]] = md.masstransport.spcthickness[md.mesh.elements[pos]] - 100 +posant = pos -#elastic loading from love numbers: -md.solidearth.lovenumbers = lovenumbers('maxdeg', 100) +# Greenland +pos = np.where(np.logical_and.reduce((late > 60, late < 90, longe > -75, longe < -15)))[0] +md.masstransport.spcthickness[md.mesh.elements[pos]] = md.masstransport.spcthickness[md.mesh.elements[pos]] - 100 +posgre = pos + +# Elastic loading from love numbers: +md.solidearth.lovenumbers = lovenumbers('maxdeg', 1000) #}}} -#mask: {{{ +# Mask: {{{ mask = gmtmask(md.mesh.lat, md.mesh.long) +oceanmask = -1 * np.ones((md.mesh.numberofvertices, 1)) +pos = np.where(mask == 0)[0] +oceanmask[pos] = 1 + icemask = np.ones((md.mesh.numberofvertices, 1)) -pos = np.where(mask == 0)[0] -icemask[pos] = -1 -pos = np.where(mask[md.mesh.elements - 1].sum(axis=1) < 3)[0] -icemask[md.mesh.elements[pos, :] - 1] = -1 +icemask[md.mesh.elements[posant]] = -1 +icemask[md.mesh.elements[posgre]] = -1 + md.mask.ice_levelset = icemask -md.mask.ocean_levelset = -icemask +md.mask.ocean_levelset = oceanmask +#}}} -#make sure that the elements that have loads are fully grounded -pos = np.nonzero(md.solidearth.surfaceload.icethicknesschange)[0] -md.mask.ocean_levelset[md.mesh.elements[pos, :] - 1] = 1 +# Time stepping {{{ +md.timestepping.start_time = 0 +md.timestepping.time_step = 1 +md.timestepping.final_time = 1 +#}}} -#make sure wherever there is an ice load, that the mask is set to ice: -#pos = np.nonzero(md.solidearth.surfaceload.icethicknesschange)[0] # TODO: Do we need to do this twice? -md.mask.ice_levelset[md.mesh.elements[pos, :] - 1] = -1 -# }}} +# Masstransport +md.basalforcings.groundedice_melting_rate = np.zeros((md.mesh.numberofvertices, 1)) +md.basalforcings.floatingice_melting_rate = np.zeros((md.mesh.numberofvertices, 1)) +md.initialization.vx = np.zeros((md.mesh.numberofvertices, 1)) +md.initialization.vy = np.zeros((md.mesh.numberofvertices, 1)) +md.initialization.sealevel = np.zeros((md.mesh.numberofvertices, 1)) +md.initialization.str = 0 -md.solidearth.settings.ocean_area_scaling = 0 +# Materials +md.materials = materials('hydro') -#geometry for the bed, arbitrary -md.geometry.bed = -np.ones((md.mesh.numberofvertices, 1)) - -#materials -md.materials=materials('hydro') - -#Miscellaneous +# Miscellaneous md.miscellaneous.name = 'test2002' -#Solution parameters +# Solution parameters +md.cluster.np = 3 md.solidearth.settings.reltol = np.nan md.solidearth.settings.abstol = 1e-3 -md.solidearth.settings.computesealevelchange = 1 +md.solidearth.settings.sealevelloading = 1 +md.solidearth.settings.isgrd = 1 +md.solidearth.settings.ocean_area_scaling = 0 +md.solidearth.settings.grdmodel = 1 -#max number of iterations reverted back to 10 (i.e., the original default value) +# Physics +md.transient.issmb = 0 +md.transient.isstressbalance = 0 +md.transient.isthermal = 0 +md.transient.ismasstransport = 1 +md.transient.isslc = 1 +md.solidearth.requested_outputs = ['Sealevel', 'Bed'] + +# Max number of iterations reverted back to 10 (i.e., the original default value) md.solidearth.settings.maxiter = 10 -#eustatic run: -md.solidearth.settings.rigid = 0 +# Eustatic run +md.solidearth.settings.selfattraction = 0 md.solidearth.settings.elastic = 0 md.solidearth.settings.rotation = 0 -md = solve(md, 'Sealevelrise') -Seustatic = md.results.SealevelriseSolution.Sealevel +md.solidearth.settings.viscous = 0 -#eustatic + rigid run: -md.solidearth.settings.rigid = 1 +md = solve(md, 'Transient') +Seustatic = md.results.TransientSolution.Sealevel +Beustatic = md.results.TransientSolution.Bed + +# Eustatic + selfattraction run +md.solidearth.settings.selfattraction = 1 md.solidearth.settings.elastic = 0 md.solidearth.settings.rotation = 0 -md = solve(md, 'Sealevelrise') -Srigid = md.results.SealevelriseSolution.Sealevel +md.solidearth.settings.viscous = 0 +md = solve(md, 'tr') +Sselfattraction = md.results.TransientSolution.Sealevel +Bselfattraction = md.results.TransientSolution.Bed -#eustatic + rigid + elastic run: -md.solidearth.settings.rigid = 1 +# Eustatic + selfattraction + elastic run +md.solidearth.settings.selfattraction = 1 md.solidearth.settings.elastic = 1 md.solidearth.settings.rotation = 0 -md = solve(md, 'Sealevelrise') -Selastic = md.results.SealevelriseSolution.Sealevel +md.solidearth.settings.viscous = 0 +md = solve(md, 'tr') +Selastic = md.results.TransientSolution.Sealevel +Belastic = md.results.TransientSolution.Bed -#eustatic + rigid + elastic + rotation run: -md.solidearth.settings.rigid = 1 +# Eustatic + selfattraction + elastic + rotation run +md.solidearth.settings.selfattraction = 1 md.solidearth.settings.elastic = 1 md.solidearth.settings.rotation = 1 -md = solve(md, 'Sealevelrise') -Srotation = md.results.SealevelriseSolution.Sealevel +md.solidearth.settings.viscous = 0 +md = solve(md, 'tr') +Srotation = md.results.TransientSolution.Sealevel +Brotation = md.results.TransientSolution.Bed -#Fields and tolerances to track changes -field_names = ['Eustatic', 'Rigid', 'Elastic', 'Rotation'] -field_tolerances = [1e-13, 1e-13, 1e-13, 1e-13] -field_values = [Seustatic, Srigid, Selastic, Srotation] +# Fields and tolerances to track changes +field_names = ['Seustatic', 'Sselfattraction', 'Selastic', 'Srotation', 'Beustatic', 'Bselfattraction', 'Belastic', 'Brotation'] +field_tolerances = [1e-13, 1e-13, 1e-13, 1e-13, 1e-13, 1e-13, 1e-13, 1e-13] +field_values = [Seustatic, Sselfattraction, Selastic, Srotation, Beustatic, Bselfattraction, Belastic, Brotation]