Index: /issm/trunk-jpl/src/m/boundaryconditions/SetIceSheetBC.m =================================================================== --- /issm/trunk-jpl/src/m/boundaryconditions/SetIceSheetBC.m (revision 26058) +++ /issm/trunk-jpl/src/m/boundaryconditions/SetIceSheetBC.m (revision 26059) @@ -31,8 +31,8 @@ %Initialize surface and basal forcings md.smb = initialize(md.smb,md); -md.basalforcings = initialize(md.basalforcings,md); +md.basalforcings = initialize(md.basalforcings,md); %Initialize ocean forcings and sealevel -md.dsl= initialize(md.dsl,md); +md.dsl = initialize(md.dsl,md); %Deal with other boundary conditions Index: /issm/trunk-jpl/src/m/boundaryconditions/SetIceSheetBC.py =================================================================== --- /issm/trunk-jpl/src/m/boundaryconditions/SetIceSheetBC.py (revision 26058) +++ /issm/trunk-jpl/src/m/boundaryconditions/SetIceSheetBC.py (revision 26059) @@ -3,11 +3,10 @@ def SetIceSheetBC(md): - """ - SETICESHEETBC - Create the boundary conditions for stressbalance and thermal models for an IceSheet with no Ice Front + """SETICESHEETBC - Create the boundary conditions for stressbalance and thermal models for an IceSheet with no Ice Front - Usage: - md = SetIceSheetBC(md) + Usage: + md = SetIceSheetBC(md) - See also: SETICESHELFBC, SETMARINEICESHEETBC + See also: SETICESHELFBC, SETMARINEICESHEETBC """ @@ -33,7 +32,10 @@ #No ice front -> do nothing - #Create zeros basalforcings and smb + #Initialize surface and basal forcings md.smb.initialize(md) md.basalforcings.initialize(md) + + #Initialize ocean forcings and sealevel + md.dsl.initialize(md) #Deal with other boundary conditions Index: /issm/trunk-jpl/src/m/classes/dsl.m =================================================================== --- /issm/trunk-jpl/src/m/classes/dsl.m (revision 26058) +++ /issm/trunk-jpl/src/m/classes/dsl.m (revision 26059) @@ -7,7 +7,7 @@ properties (SetAccess=public) - global_average_thermosteric_sea_level; %corresponds to zostoga field in CMIP5 archives. Specified as a temporally variable quantity (in m) - sea_surface_height_above_geoid; %corresponds to zos field in CMIP5 archives. Spatial average is 0. Specified as a spatio-temporally variable (in m) - 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!) + global_average_thermosteric_sea_level; %Corresponds to zostoga field in CMIP5 archives. Specified as a temporally variable quantity (in m). + 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). + 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 @@ -51,7 +51,7 @@ 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 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)'); + 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 % }}} Index: /issm/trunk-jpl/src/m/classes/dsl.py =================================================================== --- /issm/trunk-jpl/src/m/classes/dsl.py (revision 26058) +++ /issm/trunk-jpl/src/m/classes/dsl.py (revision 26059) @@ -15,22 +15,20 @@ def __init__(self): #{{{ - self.global_average_thermosteric_sea_level_change = 0 #corresponds to zostoga field in CMIP5 archives. Specified as a temporally variable global rate (mm/yr) - self.sea_surface_height_change_above_geoid = float('NaN') #corresponds to zos field in CMIP5 archives. Spatial average is 0. Specified as a spatio-temporally variable rate (mm/yr) - self.sea_water_pressure_change_at_sea_floor = float('NaN') #corresponds to bpo field in CMIP5 archives. Specified as a spatio-temporally variable rate (in Pa/yr) - self.compute_fingerprints = 0; #do we use the sea water pressure change to compute fingerprints and correct sea_surface_height_change_above_geoid + 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!). #}}} def __repr__(self): #{{{ s = ' dsl parameters:\n' - s += '{}\n'.format(fielddisplay(self, 'global_average_thermosteric_sea_level_change', 'corresponds to zostoga field in CMIP5 archives. Specified as a temporally variable global rate (mm/yr)')) - s += '{}\n'.format(fielddisplay(self, 'sea_surface_height_change_above_geoid', 'corresponds to zos field in CMIP5 archives. Spatial average is 0. Specified as a spatio-temporally variable rate (mm/yr)')) - s += '{}\n'.format(fielddisplay(self, 'sea_water_pressure_change_at_sea_floor', 'corresponds to bpo field in CMIP5 archives. Specified as a spatio-temporally variable rate (in Pa/yr)')) - s += '{}\n'.format(fielddisplay(self, 'compute_fingerprints', 'do we use the sea water pressure change to compute fingerprints and correct sea_surface_height_change_above_geoid')) + s += '{}\n'.format(fielddisplay(self, 'global_average_thermosteric_sea_level', 'Corresponds to zostoga field in CMIP5 archives. Specified as a temporally variable quantity (in m).')) + s += '{}\n'.format(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).')) + s += '{}\n'.format(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!).')) return s #}}} def extrude(self, md): #{{{ - self.sea_surface_height_change_above_geoid = project3d(md, 'vector', self.sea_surface_height_change_above_geoid, 'type', 'node') - self.sea_water_pressure_change_at_sea_floor = project3d(md, 'vector', self.sea_water_pressure_change_at_sea_floor, 'type', 'node') + 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') return self #}}} @@ -42,16 +40,14 @@ def checkconsistency(self, md, solution, analyses): #{{{ # Early return - if 'SealevelriseAnalysis' not in analyses: + if ('SealevelriseAnalysis' not in analyses) or (solution == 'TransientSolution' and not md.transient.isslc): return md - if solution == 'TransientSolution' and md.transient.isslc == 0: - return md - md = checkfield(md, 'fieldname', 'dsl.global_average_thermosteric_sea_level_change', 'NaN', 1, 'Inf', 1) - md = checkfield(md, 'fieldname', 'dsl.sea_surface_height_change_above_geoid', 'NaN', 1, 'Inf', 1, 'timeseries', 1) - md = checkfield(md, 'fieldname', 'dsl.sea_water_pressure_change_at_sea_floor', 'NaN', 1, 'Inf', 1, 'timeseries', 1) - md = checkfield(md, 'fieldname', 'dsl.compute_fingerprints', 'NaN', 1, 'Inf', 1, 'values', [0, 1]) - if self.compute_fingerprints: - # Check if geodetic flag of slr is on - if not md.slr.geodetic: - raise RuntimeError('DSL checkconsistency error message: if bottom pressure fingerprints computations are requested, slr class should have geodetic flag on') + + md = checkfield(md, 'fieldname', 'dsl.global_average_thermosteric_sea_level', 'NaN', 1, 'Inf', 1) + md = checkfield(md, 'fieldname', 'dsl.sea_surface_height_above_geoid', 'NaN', 1, 'Inf', 1, 'timeseries', 1) + md = checkfield(md, 'fieldname', 'dsl.sea_water_pressure_at_sea_floor', 'NaN', 1, 'Inf', 1, 'timeseries', 1) + + if md.solidearth.settings.compute_bp_grd: + md = checkfield(md, 'fieldname', dsl.sea_water_pressure_at_sea_floor, 'empty', 1) + return md # }}} @@ -60,7 +56,6 @@ yts = md.constants.yts WriteData(fid, prefix, 'name', 'md.dsl.model', 'data', 1, 'format', 'Integer') - WriteData(fid, prefix, 'object', self, 'class', 'dsl', 'fieldname', 'compute_fingerprints', 'format', 'Integer') - WriteData(fid, prefix, 'object', self, 'class', 'dsl', 'fieldname', 'global_average_thermosteric_sea_level_change', 'format', 'DoubleMat', 'mattype', 1, 'timeserieslength', 1+1, 'yts', md.constants.yts, 'scale', 1e-3/md.constants.yts) - WriteData(fid, prefix, 'object', self, 'class', 'dsl', 'fieldname', 'sea_water_pressure_change_at_sea_floor', 'format', 'DoubleMat', 'mattype', 1, 'timeserieslength', md.mesh.numberofvertices+1, 'yts', md.constants.yts, 'scale', 1e-3/md.constants.yts) - WriteData(fid, prefix, 'object', self, 'class', 'dsl', 'fieldname', 'sea_surface_height_change_above_geoid', 'format', 'DoubleMat', 'mattype', 1, 'timeserieslength', md.mesh.numberofvertices+1, 'yts', md.constants.yts) + 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. # }}} Index: /issm/trunk-jpl/src/m/classes/dslmme.m =================================================================== --- /issm/trunk-jpl/src/m/classes/dslmme.m (revision 26058) +++ /issm/trunk-jpl/src/m/classes/dslmme.m (revision 26059) @@ -8,7 +8,7 @@ modelid; %index into the multi-model ensemble, determine which field will be used. - global_average_thermosteric_sea_level; %corresponds to zostoga fields in CMIP5 archives. Specified as a temporally quantity (in m) for each ensemble. - sea_surface_height_above_geoid; %corresponds to zos fields in CMIP5 archives. Spatial average is 0. Specified as a spatio-temporally quantity (in m) for each ensemble - 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!) for each ensemble + global_average_thermosteric_sea_level; %Corresponds to zostoga field in CMIP5 archives. Specified as a temporally variable quantity (in m) for each ensemble. + 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) for each ensemble. + 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!) for each ensemble. end @@ -53,7 +53,7 @@ disp(sprintf(' dsl mme parameters:')); fielddisplay(self,'modelid','index into the multi-model ensemble, determine which field will be used.'); - fielddisplay(self,'global_average_thermosteric_sea_level','corresponds to zostoga fields in CMIP5 archives. Specified as a temporally variable quantity (in m) for each ensemble.'); - fielddisplay(self,'sea_surface_height_above_geoid','corresponds to zos fields in CMIP5 archives. Spatial average is 0. Specified as a spatio-temporally quantity (in m) for each ensemble.'); - fielddisplay(self,'sea_water_pressure_at_sea_floor','corresponds to bpo fields in CMIP5 archives. Specified as a spatio-temporally quantity (in m) for each ensemble.'); + fielddisplay(self,'global_average_thermosteric_sea_level','Corresponds to zostoga field in CMIP5 archives. Specified as a temporally variable quantity (in m) for each ensemble.'); + 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) for each ensemble.'); + 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!) for each ensemble.'); end % }}} function marshall(self,prefix,md,fid) % {{{ Index: /issm/trunk-jpl/src/m/classes/dslmme.py =================================================================== --- /issm/trunk-jpl/src/m/classes/dslmme.py (revision 26058) +++ /issm/trunk-jpl/src/m/classes/dslmme.py (revision 26059) @@ -15,13 +15,12 @@ def __init__(self, *args): #{{{ - self.modelid = 0 #index into the multi-model ensemble - self.global_average_thermosteric_sea_level_change = [] #corresponds to zostoga fields in CMIP5 archives. Specified as a temporally variable global rate (mm/yr) for each ensemble. - self.sea_surface_height_change_above_geoid = [] #corresponds to zos fields in CMIP5 archives. Spatial average is 0. Specified as a spatio-temporally variable rate (mm/yr) for each ensemble - self.sea_water_pressure_change_at_sea_floor = [] #corresponds to bpo fields in CMIP5 archives. Specified as a spatio-temporally variable rate (in mm/yr equivalent, not in Pa/yr!) for each ensemble - self.compute_fingerprints = 0 #corresponds to zos fields in CMIP5 archives. Spatial average is 0. Specified as a spatio-temporally variable rate (mm/yr) for each ensemble - - nargin = len(args) + self.modelid = 0 # Index into the multi-model ensemble + self.global_average_thermosteric_sea_level = [] # Corresponds to zostoga field in CMIP5 archives. Specified as a temporally variable quantity (in m) for each ensemble. + 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) for each ensemble. + 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!) for each ensemble. - if nargin == 0: + nargs = len(args) + + if nargs == 0: self.setdefaultparameters() else: @@ -32,8 +31,7 @@ s = ' dsl mme parameters:\n' s += '{}\n'.format(fielddisplay(self, 'modelid', 'index into the multi-model ensemble, determines which field will be used.')) - s += '{}\n'.format(fielddisplay(self, 'global_average_thermosteric_sea_level_change', 'corresponds to zostoga fields in CMIP5 archives. Specified as a temporally variable global rate (mm/yr) for each ensemble.')) - s += '{}\n'.format(fielddisplay(self, 'sea_surface_height_change_above_geoid', 'corresponds to zos fields in CMIP5 archives. Spatial average is 0. Specified as a spatio-temporally variable rate (mm/yr) for each ensemble.')) - s += '{}\n'.format(fielddisplay(self, 'sea_water_pressure_change_at_sea_floor', 'corresponds to bpo fields in CMIP5 archives. Specified as a spatio-temporally variable rate (in mm/yr) for each ensemble.')) - s += '{}\n'.format(fielddisplay(self, 'compute_fingerprints', 'do we use the sea water pressure change to compute fingerprints and correct sea_surface_height_change_above_geoid')) + s += '{}\n'.format(fielddisplay(self, 'global_average_thermosteric_sea_level', 'Corresponds to zostoga field in CMIP5 archives. Specified as a temporally variable quantity (in m) for each ensemble.')) + s += '{}\n'.format(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) for each ensemble.')) + s += '{}\n'.format(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!) for each ensemble.')) return s #}}} @@ -44,15 +42,17 @@ def checkconsistency(self, md, solution, analyses): # {{{ - if ('SealevelriseAnalysis' not in analyses) or (solution == 'TransientSolution' and not md.transient.isslr): + # Early return + if ('SealevelriseAnalysis' not in analyses) or (solution == 'TransientSolution' and not md.transient.isslc): return md - for i in range(len(self.global_average_thermosteric_sea_level_change)): - md = checkfield(md, 'field', self.global_average_thermosteric_sea_level_change[i], 'NaN', 1, 'Inf', 1) - md = checkfield(md, 'field', self.sea_surface_height_change_above_geoid[i], 'NaN', 1, 'Inf', 1, 'timeseries', 1) - md = checkfield(md, 'field', self.sea_water_pressure_change_at_sea_floor[i], 'NaN', 1, 'Inf', 1, 'timeseries', 1) - md = checkfield(md, 'field', self.modelid, 'NaN', 1, 'Inf', 1, '>=', 1, '<=',len(self.global_average_thermosteric_sea_level_change)) - if self.compute_fingerprints: - #check geodetic flag of slr is on - if not md.solidearth.settings.computesealevelchange: - raise Exception('DSL checkconsistency error message: if bottom pressure fingerprints computations are requested, slr class should have geodetic flag on') + + for i in range(len(self.global_average_thermosteric_sea_level)): + md = checkfield(md, 'field', self.global_average_thermosteric_sea_level[i], 'NaN', 1, 'Inf', 1) + md = checkfield(md, 'field', self.sea_surface_height_above_geoid[i], 'NaN', 1, 'Inf', 1, 'timeseries', 1) + md = checkfield(md, 'field', self.sea_water_pressure_at_sea_floor[i], 'NaN', 1, 'Inf', 1, 'timeseries', 1) + md = checkfield(md, 'field', self.modelid, 'NaN', 1, 'Inf', 1, '>=', 1, '<=',len(self.global_average_thermosteric_sea_level)) + + if self.solidearth.settings.compute_bp_grd: + md = checkfield(md, 'field', self.sea_water_pressure_at_sea_floor, 'empty', 1) + return md #}}} @@ -60,16 +60,15 @@ def marshall(self, prefix, md, fid): #{{{ WriteData(fid, prefix, 'name', 'md.dsl.model', 'data', 2, 'format', 'Integer') - WriteData(fid, prefix, 'object', self, 'fieldname', 'compute_fingerprints', 'format', 'Integer') WriteData(fid, prefix, 'object', self, 'fieldname', 'modelid', 'format', 'Double') - WriteData(fid, prefix, 'name', 'md.dsl.nummodels', 'data', len(self.global_average_thermosteric_sea_level_change), 'format', 'Integer') - WriteData(fid, prefix, 'object', self, 'fieldname', 'global_average_thermosteric_sea_level_change', 'format', 'MatArray', 'timeseries', 1, 'timeserieslength', 2, 'yts', md.constants.yts, 'scale', 1e-3 / md.constants.yts) - WriteData(fid, prefix, 'object', self, 'fieldname', 'sea_water_pressure_change_at_sea_floor', 'format', 'MatArray', 'timeserieslength', md.mesh.numberofvertices + 1, 'yts', md.constants.yts, 'scale', 1e-3) - WriteData(fid, prefix, 'object', self, 'fieldname', 'sea_surface_height_change_above_geoid', 'format', 'MatArray', 'timeserieslength', md.mesh.numberofvertices + 1, 'yts', md.constants.yts, 'scale', 1e-3 / md.constants.yts) + WriteData(fid, prefix, 'name', 'md.dsl.nummodels', 'data', len(self.global_average_thermosteric_sea_level), 'format', 'Integer') + WriteData(fid, prefix, 'object', self, 'fieldname', 'global_average_thermosteric_sea_level', 'format', 'MatArray', 'timeseries', 1, 'timeserieslength', 2, 'yts', md.constants.yts, 'scale', 1e-3 / md.constants.yts) + WriteData(fid, prefix, 'object', self, 'fieldname', 'sea_water_pressure_at_sea_floor', 'format', 'MatArray', 'timeserieslength', md.mesh.numberofvertices + 1, 'yts', md.constants.yts, 'scale', 1e-3) + WriteData(fid, prefix, 'object', self, 'fieldname', 'sea_surface_height_above_geoid', 'format', 'MatArray', 'timeserieslength', md.mesh.numberofvertices + 1, 'yts', md.constants.yts, 'scale', 1e-3 / md.constants.yts) #}}} def extrude(self, md): #{{{ - for i in range(len(self.global_average_thermosteric_sea_level_change)): - self.sea_surface_height_change_above_geoid[i] = project3d(md, 'vector', self.self.sea_surface_height_change_above_geoid[i], 'type', 'node', 'layer', 1) - self.sea_water_pressure_change_at_sea_floor[i] = project3d(md, 'vector', self.sea_water_pressure_change_at_sea_floor[i], 'type', 'node', 'layer', 1) + for i in range(len(self.global_average_thermosteric_sea_level)): + self.sea_surface_height_above_geoid[i] = project3d(md, 'vector', self.self.sea_surface_height_above_geoid[i], 'type', 'node', 'layer', 1) + self.sea_water_pressure_at_sea_floor[i] = project3d(md, 'vector', self.sea_water_pressure_at_sea_floor[i], 'type', 'node', 'layer', 1) return self Index: /issm/trunk-jpl/src/m/classes/fourierlove.m =================================================================== --- /issm/trunk-jpl/src/m/classes/fourierlove.m (revision 26058) +++ /issm/trunk-jpl/src/m/classes/fourierlove.m (revision 26059) @@ -80,10 +80,6 @@ %need 'litho' material: - if ~isa(md.materials,'materials') + if ~isa(md.materials,'materials') | ~sum(strcmpi(md.materials.nature,'litho')) error('Need a ''litho'' material to run a Fourier Love number analysis'); - else - if ~sum(strcmpi(md.materials.nature,'litho')), - error('Need a ''litho'' material to run a Fourier Love number analysis'); - end end Index: /issm/trunk-jpl/src/m/classes/fourierlove.py =================================================================== --- /issm/trunk-jpl/src/m/classes/fourierlove.py (revision 26058) +++ /issm/trunk-jpl/src/m/classes/fourierlove.py (revision 26059) @@ -5,9 +5,8 @@ class fourierlove(object): - """ - Fourier Love Number class definition + """FOURIERLOVE - Fourier Love Number class definition - Usage: - fourierlove = fourierlove() + Usage: + fourierlove = fourierlove() """ def __init__(self): # {{{ @@ -23,9 +22,12 @@ self.forcing_type = 0 - #set defaults 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]')) @@ -76,4 +78,6 @@ 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]) @@ -89,4 +93,7 @@ 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 # }}} Index: /issm/trunk-jpl/src/m/classes/geometry.m =================================================================== --- /issm/trunk-jpl/src/m/classes/geometry.m (revision 26058) +++ /issm/trunk-jpl/src/m/classes/geometry.m (revision 26059) @@ -85,8 +85,8 @@ end % }}} function marshall(self,prefix,md,fid) % {{{ - - if (size(self.thickness==md.mesh.numberofvertices) | (self.thickness==md.mesh.numberofvertices+1)), + + if (size(self.thickness)==md.mesh.numberofvertices) | (size(self.thickness)==md.mesh.numberofvertices+1)), WriteData(fid,prefix,'object',self,'fieldname','thickness','format','DoubleMat','mattype',1,'timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts); - elseif (size(self.thickness==md.mesh.numberofelements) | (self.thickness==md.mesh.numberofelements+1)), + elseif (size(self.thickness)==md.mesh.numberofelements) | (size(self.thickness)==md.mesh.numberofelements+1)), WriteData(fid,prefix,'object',self,'fieldname','thickness','format','DoubleMat','mattype',1,'timeserieslength',md.mesh.numberofelements+1,'yts',md.constants.yts); else Index: /issm/trunk-jpl/src/m/classes/geometry.py =================================================================== --- /issm/trunk-jpl/src/m/classes/geometry.py (revision 26058) +++ /issm/trunk-jpl/src/m/classes/geometry.py (revision 26059) @@ -51,10 +51,6 @@ def checkconsistency(self, md, solution, analyses): #{{{ - if (solution == 'TransientSolution' and md.transient.isgia) or (solution == 'GiaSolution'): - md = checkfield(md, 'fieldname', 'geometry.thickness', 'timeseries', 1, 'NaN', 1, 'Inf', 1) - elif solution == 'SealevelriseSolution': - md = checkfield(md, 'fieldname', 'geometry.bed', 'NaN', 1, 'Inf', 1, 'size', [md.mesh.numberofvertices]) - elif solution == 'LoveSolution': - return + if solution == 'LoveSolution': + return md else: md = checkfield(md, 'fieldname', 'geometry.surface', 'NaN', 1, 'Inf', 1, 'size', [md.mesh.numberofvertices]) @@ -71,11 +67,16 @@ md.checkmessage('equality base = bed on grounded ice violated') md = checkfield(md, 'fieldname', 'geometry.bed', 'NaN', 1, 'Inf', 1, 'size', [md.mesh.numberofvertices]) - return md # }}} def marshall(self, prefix, md, fid): #{{{ + if (len(self.thickness) == md.mesh.numberofvertices) or (len(self.thickness) == md.mesh.numberofvertices + 1): + WriteData(fid, prefix, 'object', self, 'fieldname', 'thickness', 'format', 'DoubleMat', 'mattype', 1, 'timeserieslength', md.mesh.numberofvertices + 1, 'yts', md.constants.yts) + elif (len(self.thickness) == md.mesh.numberofelements) or (len(self.thickness) == md.mesh.numberofelements + 1): + WriteData(fid, prefix, 'object', self, 'fieldname', 'thickness', 'format', 'DoubleMat', 'mattype', 1, 'timeserieslength', md.mesh.numberofvertices + 1, 'yts', md.constants.yts) + else: + raise RuntimeError('geometry thickness time series should be a vertex or element time series') + WriteData(fid, prefix, 'object', self, 'fieldname', 'surface', 'format', 'DoubleMat', 'mattype', 1) - WriteData(fid, prefix, 'object', self, 'fieldname', 'thickness', 'format', 'DoubleMat', 'mattype', 1, 'timeserieslength', md.mesh.numberofvertices + 1, 'yts', md.constants.yts) WriteData(fid, prefix, 'object', self, 'fieldname', 'base', 'format', 'DoubleMat', 'mattype', 1) WriteData(fid, prefix, 'object', self, 'fieldname', 'bed', 'format', 'DoubleMat', 'mattype', 1) Index: /issm/trunk-jpl/src/m/classes/hydrologydc.py =================================================================== --- /issm/trunk-jpl/src/m/classes/hydrologydc.py (revision 26058) +++ /issm/trunk-jpl/src/m/classes/hydrologydc.py (revision 26059) @@ -7,9 +7,8 @@ class hydrologydc(object): - """ - Hydrologydc class definition + """HYDROLOGYDC class definition Usage: - hydrologydc = hydrologydc() + hydrologydc = hydrologydc() """ @@ -49,9 +48,12 @@ self.eplflip_lock = 0 - #set defaults self.setdefaultparameters() # }}} def __repr__(self): # {{{ + # TODO: + # - Convert all formatting to calls to .format (see any + # already converted .__repr__ method for examples) + # string = ' hydrology Dual Porous Continuum Equivalent parameters:' string = ' - general parameters' Index: /issm/trunk-jpl/src/m/classes/hydrologyshreve.m =================================================================== --- /issm/trunk-jpl/src/m/classes/hydrologyshreve.m (revision 26058) +++ /issm/trunk-jpl/src/m/classes/hydrologyshreve.m (revision 26059) @@ -31,5 +31,5 @@ %Type of stabilization to use 0:nothing 1:artificial_diffusivity self.stabilization = 1; - self.requested_outputs = {'default'}; + self.requested_outputs = {'default'}; end % }}} function md = checkconsistency(self,md,solution,analyses) % {{{ Index: /issm/trunk-jpl/src/m/classes/hydrologyshreve.py =================================================================== --- /issm/trunk-jpl/src/m/classes/hydrologyshreve.py (revision 26058) +++ /issm/trunk-jpl/src/m/classes/hydrologyshreve.py (revision 26059) @@ -5,9 +5,8 @@ class hydrologyshreve(object): - """ - HYDROLOGYSHREVE class definition + """HYDROLOGYSHREVE class definition - Usage: - hydrologyshreve = hydrologyshreve() + Usage: + hydrologyshreve = hydrologyshreve() """ @@ -16,11 +15,13 @@ self.stabilization = 0 self.requested_outputs = [] - #set defaults + self.setdefaultparameters() - #}}} 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]')) @@ -48,5 +49,5 @@ def checkconsistency(self, md, solution, analyses): # {{{ #Early return - if 'HydrologyShreveAnalysis' not in analyses: + if 'HydrologyShreveAnalysis' not in analyses or (solution == 'TransientSolution' and not md.transient.ishydrology): return md @@ -61,5 +62,5 @@ 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 + #process requested outputs outputs = self.requested_outputs indices = [i for i, x in enumerate(outputs) if x == 'default'] Index: /issm/trunk-jpl/src/m/classes/hydrologytws.m =================================================================== --- /issm/trunk-jpl/src/m/classes/hydrologytws.m (revision 26058) +++ /issm/trunk-jpl/src/m/classes/hydrologytws.m (revision 26059) @@ -24,5 +24,5 @@ function list = defaultoutputs(self,md) % {{{ list = {''}; - end % }}} + end % }}} function self = setdefaultparameters(self) % {{{ @@ -46,14 +46,14 @@ 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); - 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 savemodeljs(self,fid,modelname) % {{{ - + writejs1Darray(fid,[modelname '.hydrology.spcwatercolumn'],self.spcwatercolumn); Index: /issm/trunk-jpl/src/m/classes/hydrologytws.py =================================================================== --- /issm/trunk-jpl/src/m/classes/hydrologytws.py (revision 26059) +++ /issm/trunk-jpl/src/m/classes/hydrologytws.py (revision 26059) @@ -0,0 +1,64 @@ +import numpy as np + +from structtoobj import structtoobj + +class hydrologytws(object): + """HYDROLOGYTWS class definition + + Usage: + hydrologytws = hydrologytws() + """ + + def __init__(self): # {{{ + self.spcwatercolumn = np.nan + self.requested_outputs = np.nan + + nargs = len(args) + if nargs == 0: + self.setdefaultparameters() + elif nargs == 1: + self = structtoobj(self, args[0]) + else: + raise RuntimeError('constructor not supported') + #}}} + + def __repr__(self): # {{{ + s = ' hydrologytws solution parameters:\n' + s += '{}\n'.format(fielddisplay(self, 'spcwatercolumn', 'water thickness constraints (NaN means no constraint) [m]')) + s += '{}\n'.format(fielddisplay(self, 'requested_outputs', 'additional outputs requested')) + return s + #}}} + + def defaultoutputs(self, md): # {{{ + return [''] + #}}} + + def setdefaultparameters(self): # {{{ + self.requested_outputs = ['defualt'] + return self + #}}} + + def extrude(self, md): # {{{ + return self + #}}} + + def checkconsistency(self, md, solution, analyses): # {{{ + # Early return + if 'HydrologyTwsAnalysis' not in analyses: + return + md = checkfield(md, 'fieldname', 'hydrology.spcwatercolumn', 'Inf', 1, 'timeseries', 1) + #}}} + + 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) + outputs = self.requested_outputs + pos = find(ismember(outputs,'default')) + if not len(pos), + outputs[pos] = []; # remove 'default' from outputs + outputs.extend(defaultoutputs(self, md)) # add defaults + end + WriteData(fid, prefix, 'data', outputs, 'name', 'md.hydrology.requested_outputs', 'format', 'StringArray') + # }}} + + Index: /issm/trunk-jpl/src/m/classes/initialization.m =================================================================== --- /issm/trunk-jpl/src/m/classes/initialization.m (revision 26058) +++ /issm/trunk-jpl/src/m/classes/initialization.m (revision 26059) @@ -22,5 +22,5 @@ sealevel = NaN; bottompressure = NaN; - sample = NaN; + sample = NaN; end methods @@ -52,5 +52,4 @@ end % }}} function self = setdefaultparameters(self) % {{{ - end % }}} function md = checkconsistency(self,md,solution,analyses) % {{{ @@ -98,8 +97,5 @@ if ismember('HydrologyShreveAnalysis',analyses), if isa(md.hydrology,'hydrologyshreve'), - if strcmp(solution,'TransientSolution') & md.transient.ishydrology, - md = checkfield(md,'fieldname','initialization.watercolumn','NaN',1,'Inf',1,'size',[md.mesh.numberofvertices 1]); - end - if strcmp(solution,'HydrologySolution'), + if (strcmp(solution,'TransientSolution') & md.transient.ishydrology) | strcmp(solution,'HydrologySolution'), md = checkfield(md,'fieldname','initialization.watercolumn','NaN',1,'Inf',1,'size',[md.mesh.numberofvertices 1]); end @@ -111,5 +107,5 @@ end end - if ismember('SealevelchangeAnalysis',analyses) & ~(strcmp(solution,'TransientSolution') & md.transient.isslc == 0), + if ismember('SealevelchangeAnalysis',analyses), if strcmp(solution,'TransientSolution') & md.transient.isslc, md = checkfield(md,'fieldname','initialization.sealevel','NaN',1,'Inf',1,'size',[md.mesh.numberofvertices 1]); @@ -134,10 +130,10 @@ md = checkfield(md,'fieldname','initialization.epl_thickness','NaN',1,'Inf',1,'size',[md.mesh.numberofvertices 1]); end - end - end - if ismember('SamplingAnalysis',analyses) & ~(strcmp(solution,'TransientSolution') & md.transient.issampling == 0), - if ~isnan(md.initialization.sample) - md = checkfield(md,'fieldname','initialization.sample','NaN',1,'Inf',1,'size',[md.mesh.numberofvertices 1]); - end + end + end + if ismember('SamplingAnalysis',analyses) & ~(strcmp(solution,'TransientSolution') & md.transient.issampling == 0), + if ~isnan(md.initialization.sample) + md = checkfield(md,'fieldname','initialization.sample','NaN',1,'Inf',1,'size',[md.mesh.numberofvertices 1]); + end end end % }}} @@ -159,5 +155,5 @@ fielddisplay(self,'hydraulic_potential','Hydraulic potential (for GlaDS) [Pa]'); fielddisplay(self,'channelarea','subglacial water channel area (for GlaDS) [m2]'); - fielddisplay(self,'sample','Realization of a Gaussian random field'); + fielddisplay(self,'sample','Realization of a Gaussian random field'); end % }}} function marshall(self,prefix,md,fid) % {{{ @@ -179,6 +175,6 @@ 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 numel(self.enthalpy) <= 1, @@ -195,5 +191,5 @@ end % }}} function savemodeljs(self,fid,modelname) % {{{ - + writejs1Darray(fid,[modelname '.initialization.vx'],self.vx); writejs1Darray(fid,[modelname '.initialization.vy'],self.vy); @@ -210,6 +206,6 @@ writejs1Darray(fid,[modelname '.initialization.hydraulic_potential'],self.hydraulic_potential); writejs1Darray(fid,[modelname '.initialization.channel'],self.channelarea); - writejs1Darray(fid,[modelname '.initialization.sample'],self.sample); - + writejs1Darray(fid,[modelname '.initialization.sample'],self.sample); + end % }}} end Index: /issm/trunk-jpl/src/m/classes/initialization.py =================================================================== --- /issm/trunk-jpl/src/m/classes/initialization.py (revision 26058) +++ /issm/trunk-jpl/src/m/classes/initialization.py (revision 26059) @@ -29,4 +29,6 @@ self.hydraulic_potential = np.nan self.channelarea = np.nan + self.sealevel = np.nan + self.bottompressure = np.nan self.sample = np.nan @@ -67,4 +69,6 @@ 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 @@ -90,4 +94,7 @@ md = checkfield(md, 'fieldname', 'initialization.vx', 'NaN', 1, 'Inf', 1, 'size', [md.mesh.numberofvertices]) md = checkfield(md, 'fieldname', 'initialization.vy', 'NaN', 1, 'Inf', 1, 'size', [md.mesh.numberofvertices]) + if 'OceantransportAnalysis' in analyses: + if solution == 'TransientSolution' and md.transient.isslc and md.transient.isoceantransport: + md = checkfield(md, 'fieldname', 'initialization.bottompressure', 'NaN', 1, 'Inf', 1, 'size', [md.mesh.numberofvertices]) if 'BalancethicknessAnalysis' in analyses and solution == 'BalancethicknessSolution': md = checkfield(md, 'fieldname', 'initialization.vx', 'NaN', 1, 'Inf', 1, 'size', [md.mesh.numberofvertices]) @@ -112,13 +119,12 @@ if 'HydrologyShreveAnalysis' in analyses: if hasattr(md.hydrology, 'hydrologyshreve'): + if (solution == 'TransientSolution' and md.transient.ishydrology) or solution == 'HydrologySolution': + md = checkfield(md, 'fieldname', 'initialization.watercolumn', 'NaN', 1, 'Inf', 1, 'size', [md.mesh.numberofvertices]) + if 'HydrologyTwsAnalysis' in analyses: + if hasattr(md.hydrology, 'hydrologytws'): md = checkfield(md, 'fieldname', 'initialization.watercolumn', 'NaN', 1, 'Inf', 1, 'size', [md.mesh.numberofvertices]) - if 'HydrologyDCInefficientAnalysis' in analyses: - if hasattr(md.hydrology, 'hydrologydc'): - md = checkfield(md, 'fieldname', 'initialization.sediment_head', 'NaN', 1, 'Inf', 1, 'size', [md.mesh.numberofvertices]) - if 'HydrologyDCEfficientAnalysis' in analyses: - if hasattr(md.hydrology, 'hydrologydc'): - if md.hydrology.isefficientlayer == 1: - md = checkfield(md, 'fieldname', 'initialization.epl_head', 'NaN', 1, 'Inf', 1, 'size', [md.mesh.numberofvertices]) - md = checkfield(md, 'fieldname', 'initialization.epl_thickness', 'NaN', 1, 'Inf', 1, 'size', [md.mesh.numberofvertices]) + if 'SealevelchangeAnalysis' in analyses: + if solution == 'TransientSolution' and md.transient.isslc: + md = checkfield(md, 'fieldname', 'initialization.sealevel', 'NaN', 1, 'Inf', 1, 'size', [md.mesh.numberofvertices]) if 'HydrologyGlaDSAnalysis' in analyses: if hasattr(md.hydrology, 'hydrologyglads'): @@ -138,4 +144,6 @@ 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) + WriteData(fid, prefix, 'object', self, 'fieldname', 'bottompressure', '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) Index: /issm/trunk-jpl/src/m/classes/matdamageice.m =================================================================== --- /issm/trunk-jpl/src/m/classes/matdamageice.m (revision 26058) +++ /issm/trunk-jpl/src/m/classes/matdamageice.m (revision 26059) @@ -6,23 +6,23 @@ classdef matdamageice properties (SetAccess=public) - rho_ice = 0.; - rho_water = 0.; - rho_freshwater = 0.; - mu_water = 0.; - heatcapacity = 0.; - latentheat = 0.; - thermalconductivity = 0.; - temperateiceconductivity = 0.; + rho_ice = 0.; + rho_water = 0.; + rho_freshwater = 0.; + mu_water = 0.; + heatcapacity = 0.; + latentheat = 0.; + thermalconductivity = 0.; + temperateiceconductivity = 0.; effectiveconductivity_averaging = 0.; - meltingpoint = 0.; - beta = 0.; - mixed_layer_capacity = 0.; - thermal_exchange_velocity = 0.; - rheology_B = NaN; - rheology_n = NaN; - rheology_law = ''; + meltingpoint = 0.; + beta = 0.; + mixed_layer_capacity = 0.; + thermal_exchange_velocity = 0.; + rheology_B = NaN; + rheology_n = NaN; + rheology_law = ''; %slc - earth_density = 0; + earth_density = 0; end Index: /issm/trunk-jpl/src/m/classes/matdamageice.py =================================================================== --- /issm/trunk-jpl/src/m/classes/matdamageice.py (revision 26058) +++ /issm/trunk-jpl/src/m/classes/matdamageice.py (revision 26059) @@ -31,16 +31,15 @@ self.rheology_law = '' - #giaivins: - self.lithosphere_shear_modulus = 0. - self.lithosphere_density = 0. - self.mantle_shear_modulus = 0. - self.mantle_density = 0. + #SLC + self.earth_density = 5512 # average density of the Earth, (kg / m^3) - #SLR - self.earth_density = 5512 # average density of the Earth, (kg / m^3) self.setdefaultparameters() #}}} def __repr__(self): # {{{ + # TODO: + # - Convert all formatting to calls to .format (see any + # already converted .__repr__ method for examples) + # string = " Materials:" string = "%s\n%s" % (string, fielddisplay(self, "rho_ice", "ice density [kg / m^3]")) @@ -60,8 +59,4 @@ string = "%s\n%s" % (string, fielddisplay(self, "rheology_n", "Glen's flow law exponent")) string = "%s\n%s" % (string, fielddisplay(self, "rheology_law", "law for the temperature dependance of the rheology: 'None', 'BuddJacka', 'Cuffey', 'CuffeyTemperate', 'Paterson', 'Arrhenius' or 'LliboutryDuval'")) - string = "%s\n%s" % (string, fielddisplay(self, "lithosphere_shear_modulus", "Lithosphere shear modulus [Pa]")) - string = "%s\n%s" % (string, fielddisplay(self, "lithosphere_density", "Lithosphere density [g / cm^ - 3]")) - string = "%s\n%s" % (string, fielddisplay(self, "mantle_shear_modulus", "Mantle shear modulus [Pa]")) - string = "%s\n%s" % (string, fielddisplay(self, "mantle_density", "Mantle density [g / cm^ - 3]")) string = "%s\n%s" % (string, fielddisplay(self, "earth_density", "Mantle density [kg / m^ - 3]")) return string @@ -105,11 +100,5 @@ self.rheology_law = 'Paterson' - # GIA: - self.lithosphere_shear_modulus = 6.7e10 # (Pa) - self.lithosphere_density = 3.32 # (g / cm^ - 3) - self.mantle_shear_modulus = 1.45e11 # (Pa) - self.mantle_density = 3.34 # (g / cm^ - 3) - - #SLR + #SLC self.earth_density = 5512 #average density of the Earth, (kg / m^3) return self @@ -130,10 +119,4 @@ md = checkfield(md, 'fieldname', 'materials.mantle_density', '>', 0, 'numel', [1]) md = checkfield(md, 'fieldname', 'materials.earth_density', '>', 0, 'numel', [1]) - - if 'GiaAnalysis' in analyses: - md = checkfield(md, 'fieldname', 'materials.lithosphere_shear_modulus', '>', 0, 'numel', 1) - md = checkfield(md, 'fieldname', 'materials.lithosphere_density', '>', 0, 'numel', 1) - md = checkfield(md,'fieldname', 'materials.mantle_shear_modulus', '>', 0, 'numel', 1) - md = checkfield(md,'fieldname', 'materials.mantle_density', '>', 0, 'numel', 1) if 'SealevelriseAnalysis' in analyses: @@ -161,8 +144,4 @@ WriteData(fid, prefix, 'object', self, 'class', 'materials', 'fieldname', 'rheology_n', 'format', 'DoubleMat', 'mattype', 2) WriteData(fid, prefix, 'data', self.rheology_law, 'name', 'md.materials.rheology_law', 'format', 'String') - WriteData(fid, prefix, 'object', self, 'class', 'materials', 'fieldname', 'lithosphere_shear_modulus', 'format', 'Double') - WriteData(fid, prefix, 'object', self, 'class', 'materials', 'fieldname', 'lithosphere_density', 'format', 'Double', 'scale', 10.**3.) - WriteData(fid, prefix, 'object', self, 'class', 'materials', 'fieldname', 'mantle_shear_modulus', 'format', 'Double') - WriteData(fid, prefix, 'object', self, 'class', 'materials', 'fieldname', 'mantle_density', 'format', 'Double', 'scale', 10.**3.) WriteData(fid, prefix, 'object', self, 'class', 'materials', 'fieldname', 'earth_density', 'format', 'Double') Index: /issm/trunk-jpl/src/m/classes/matenhancedice.m =================================================================== --- /issm/trunk-jpl/src/m/classes/matenhancedice.m (revision 26058) +++ /issm/trunk-jpl/src/m/classes/matenhancedice.m (revision 26059) @@ -6,24 +6,24 @@ classdef matenhancedice properties (SetAccess=public) - rho_ice = 0.; - rho_water = 0.; - rho_freshwater = 0.; - mu_water = 0.; - heatcapacity = 0.; - latentheat = 0.; - thermalconductivity = 0.; - temperateiceconductivity = 0.; + rho_ice = 0.; + rho_water = 0.; + rho_freshwater = 0.; + mu_water = 0.; + heatcapacity = 0.; + latentheat = 0.; + thermalconductivity = 0.; + temperateiceconductivity = 0.; effectiveconductivity_averaging = 0.; - meltingpoint = 0.; - beta = 0.; - mixed_layer_capacity = 0.; - thermal_exchange_velocity = 0.; - rheology_E = NaN; - rheology_B = NaN; - rheology_n = NaN; - rheology_law = ''; + meltingpoint = 0.; + beta = 0.; + mixed_layer_capacity = 0.; + thermal_exchange_velocity = 0.; + rheology_E = NaN; + rheology_B = NaN; + rheology_n = NaN; + rheology_law = ''; - %slr - earth_density = 0; + %SLC + earth_density = 0; end Index: /issm/trunk-jpl/src/m/classes/matenhancedice.py =================================================================== --- /issm/trunk-jpl/src/m/classes/matenhancedice.py (revision 26058) +++ /issm/trunk-jpl/src/m/classes/matenhancedice.py (revision 26059) @@ -1,14 +1,13 @@ +from checkfield import checkfield from fielddisplay import fielddisplay from project3d import project3d -from checkfield import checkfield from WriteData import WriteData class matenhancedice(object): - """ - MATICE class definition + """MATICE class definition - Usage: - matenhancedice = matenhancedice() + Usage: + matenhancedice = matenhancedice() """ @@ -32,16 +31,15 @@ self.rheology_law = '' - #GIA - self.lithosphere_shear_modulus = 0. - self.lithosphere_density = 0. - self.mantle_shear_modulus = 0. - self.mantle_density = 0. + #SLC + self.earth_density = 0 # average density of the Earth, (kg/m^3) - #SLR - self.earth_density = 0 # average density of the Earth, (kg/m^3) self.setdefaultparameters() #}}} def __repr__(self): #{{{ + # TODO: + # - Convert all formatting to calls to .format (see any + # already converted .__repr__ method for examples) + # s = " Materials:" s = "%s\n%s" % (s, fielddisplay(self, "rho_ice", "ice density [kg/m^3]")) @@ -62,8 +60,4 @@ s = "%s\n%s" % (s, fielddisplay(self, "rheology_n", "Glen's flow law exponent")) s = "%s\n%s" % (s, fielddisplay(self, "rheology_law", "law for the temperature dependance of the rheology: 'None', 'BuddJacka', 'Cuffey', 'CuffeyTemperate', 'Paterson', 'Arrhenius' or 'LliboutryDuval'")) - s = "%s\n%s" % (s, fielddisplay(self, "lithosphere_shear_modulus", "Lithosphere shear modulus [Pa]")) - s = "%s\n%s" % (s, fielddisplay(self, "lithosphere_density", "Lithosphere density [g/cm^-3]")) - s = "%s\n%s" % (s, fielddisplay(self, "mantle_shear_modulus", "Mantle shear modulus [Pa]")) - s = "%s\n%s" % (s, fielddisplay(self, "mantle_density", "Mantle density [g/cm^-3]")) s = "%s\n%s" % (s, fielddisplay(self, "earth_density", "Mantle density [kg/m^-3]")) @@ -109,5 +103,5 @@ self.rheology_law = 'Paterson' - # GIA: + #GIA self.lithosphere_shear_modulus = 6.7 * 10**10 # (Pa) self.lithosphere_density = 3.32 # (g / cm^ - 3) @@ -115,5 +109,5 @@ self.mantle_density = 3.34 # (g / cm^ - 3) - #SLR + #SLC self.earth_density = 5512 #average density of the Earth, (kg / m^3) @@ -132,10 +126,4 @@ md = checkfield(md, 'fieldname', 'materials.effectiveconductivity_averaging', 'numel', [1], 'values', [0, 1, 2]) - if 'GiaAnalysis' in analyses: - md = checkfield(md, 'fieldname', 'materials.lithosphere_shear_modulus', '>', 0, 'numel', 1) - md = checkfield(md, 'fieldname', 'materials.lithosphere_density', '>', 0, 'numel', 1) - md = checkfield(md, 'fieldname', 'materials.mantle_shear_modulus', '>', 0, 'numel', 1) - md = checkfield(md, 'fieldname', 'materials.mantle_density', '>', 0, 'numel', 1) - if 'SealevelriseAnalysis' in analyses: md = checkfield(md, 'fieldname', 'materials.earth_density', '>', 0, 'numel', 1) @@ -162,8 +150,4 @@ WriteData(fid, prefix, 'object', self, 'class', 'materials', 'fieldname', 'rheology_n', 'format', 'DoubleMat', 'mattype', 2) WriteData(fid, prefix, 'data', self.rheology_law, 'name', 'md.materials.rheology_law', 'format', 'String') - WriteData(fid, prefix, 'object', self, 'class', 'materials', 'fieldname', 'lithosphere_shear_modulus', 'format', 'Double') - WriteData(fid, prefix, 'object', self, 'class', 'materials', 'fieldname', 'lithosphere_density', 'format', 'Double', 'scale', 10**3) - WriteData(fid, prefix, 'object', self, 'class', 'materials', 'fieldname', 'mantle_shear_modulus', 'format', 'Double') - WriteData(fid, prefix, 'object', self, 'class', 'materials', 'fieldname', 'mantle_density', 'format', 'Double', 'scale', 10**3) WriteData(fid, prefix, 'object', self, 'class', 'materials', 'fieldname', 'earth_density', 'format', 'Double') # }}} Index: /issm/trunk-jpl/src/m/classes/materials.m =================================================================== --- /issm/trunk-jpl/src/m/classes/materials.m (revision 26058) +++ /issm/trunk-jpl/src/m/classes/materials.m (revision 26059) @@ -122,5 +122,4 @@ self.rheology_B = 1*1e8; self.rheology_n = 3; - case 'litho' Index: /issm/trunk-jpl/src/m/classes/materials.py =================================================================== --- /issm/trunk-jpl/src/m/classes/materials.py (revision 26058) +++ /issm/trunk-jpl/src/m/classes/materials.py (revision 26059) @@ -8,10 +8,9 @@ class materials(object): - ''' - MATERIALS class definition - - Usage: - materials = materials() - ''' + """MATERIALS class definition + + Usage: + materials = materials() + """ def __init__(self, *args): #{{{ @@ -38,4 +37,5 @@ setattr(self, 'thermalconductivity', 0) setattr(self, 'temperateiceconductivity', 0) + setattr(self, 'effectiveconductivity_averaging', 0) setattr(self, 'meltingpoint', 0) setattr(self, 'beta', 0) @@ -60,7 +60,7 @@ setattr(self, 'rho_water', 0) setattr(self, 'rho_freshwater', 0) - setattr(self, 'earth_density', 0) else: raise RuntimeError("materials constructor error message: nature of the material not supported yet! ('ice' or 'litho' or 'hydro')") + setattr(self, 'earth_density', 0) #set default parameters: @@ -74,35 +74,56 @@ #ice density (kg/m^3) self.rho_ice = 917. + #ocean water density (kg/m^3) self.rho_water = 1023. + #fresh water density (kg/m^3) self.rho_freshwater = 1000. + #water viscosity (N.s/m^2) self.mu_water = 0.001787 + #ice heat capacity cp (J/kg/K) self.heatcapacity = 2093. + #ice latent heat of fusion L (J / kg) self.latentheat = 3.34e5 + #ice thermal conductivity (W/m/K) self.thermalconductivity = 2.4 + #wet ice thermal conductivity (W/m/K) self.temperateiceconductivity = 0.24 + + #computation of effective conductivity + self.effectiveconductivity_averaging = 1 + #the melting point of ice at 1 atmosphere of pressure in K self.meltingpoint = 273.15 + #rate of change of melting point with pressure (K/Pa) self.beta = 9.8e-8 + #mixed layer (ice-water interface) heat capacity (J/kg/K) self.mixed_layer_capacity = 3974. + #thermal exchange velocity (ice-water interface) (m/s) self.thermal_exchange_velocity = 1.00e-4 + #Rheology law: what is the temperature dependence of B with T #available: none, paterson and arrhenius self.rheology_law = 'Paterson' - elif nat == 'litho': - #we default to a configuration that enables running GIA solutions using giacaron and / or giaivins. + + #Rheology fields default + self.rheology_B = 1e8 + self.rheology_n = 3 + elif nat == 'litho': + #we default to a configuration that enables running GIA solutions using giacaron and/or giaivins. self.numlayers = 2 + #center of the earth (approximation, must not be 0), then the lab (lithosphere / asthenosphere boundary) then the surface #(with 1d3 to avoid numerical singularities) self.radius = [1e3, 6278e3, 6378e3] + self.viscosity = [1e21, 1e40] #mantle and lithosphere viscosity (respectively) [Pa.s] self.lame_mu = [1.45e11, 6.7e10] # (Pa) #lithosphere and mantle shear modulus (respectively) [Pa] @@ -116,12 +137,15 @@ #ice density (kg/m^3) self.rho_ice = 917. + #ocean water density (kg/m^3) self.rho_water = 1023. - #average density of the Earth (kg/m^3) - self.earth_density = 5512 + #fresh water density (kg/m^3) self.rho_freshwater = 1000. else: raise RuntimeError("materials setdefaultparameters error message: nature of the material not supported yet! ('ice' or 'litho' or 'hydro')") + + #average density of the Earth (kg/m^3) + self.earth_density = 5512 #}}} @@ -223,5 +247,4 @@ 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 - for i in range(len(self.nature)): nat = self.nature[i] @@ -236,4 +259,5 @@ WriteData(fid, prefix, 'object', self, 'class', 'materials', 'fieldname', 'thermalconductivity', 'format', 'Double') WriteData(fid, prefix, 'object', self, 'class', 'materials', 'fieldname', 'temperateiceconductivity', 'format', 'Double') + WriteData(fid, prefix, 'object', self, 'class', 'materials', 'fieldname', 'effectiveconductivity_averaging', 'format', 'Integer') WriteData(fid, prefix, 'object', self, 'class', 'materials', 'fieldname', 'meltingpoint', 'format', 'Double') WriteData(fid, prefix, 'object', self, 'class', 'materials', 'fieldname', 'beta', 'format', 'Double') @@ -254,11 +278,17 @@ WriteData(fid, prefix, 'object', self, 'class', 'materials', 'fieldname', 'burgers_viscosity', 'format', 'DoubleMat', 'mattype', 3) WriteData(fid, prefix, 'object', self, 'class', 'materials', 'fieldname', 'burgers_mu', 'format', 'DoubleMat', 'mattype', 3) + # Compute earth density compatible with our layer density distribution + earth_density = 0 + for i in range(len(self.numlayers)): + earth_density = earth_density + (self.radius[i + 1] ** 3 - self.radius[i] ** 3) * self.density[i] + earth_density = earth_density / self.radius[self.numlayers + 1] ** 3 + self.earth_density = earth_density elif nat == 'hydro': WriteData(fid, prefix, 'object', self, 'class', 'materials', 'fieldname', 'rho_ice', 'format', 'Double') WriteData(fid, prefix, 'object', self, 'class', 'materials', 'fieldname', 'rho_water', 'format', 'Double') - WriteData(fid, prefix, 'object', self, 'class', 'materials', 'fieldname', 'earth_density', 'format', 'Double') WriteData(fid, prefix, 'object', self, 'class', 'materials', 'fieldname', 'rho_freshwater', 'format', 'Double') else: raise RuntimeError("materials constructor error message: nature of the material not supported yet! ('ice' or 'litho' or 'hydro')") + WriteData(fid, prefix, 'data', self.earth_density, 'name', 'md.materials.earth_density', 'format', 'Double') #}}} Index: /issm/trunk-jpl/src/m/classes/matestar.m =================================================================== --- /issm/trunk-jpl/src/m/classes/matestar.m (revision 26058) +++ /issm/trunk-jpl/src/m/classes/matestar.m (revision 26059) @@ -6,24 +6,24 @@ classdef matestar properties (SetAccess=public) - rho_ice = 0.; - rho_water = 0.; - rho_freshwater = 0.; - mu_water = 0.; - heatcapacity = 0.; - latentheat = 0.; - thermalconductivity = 0.; - temperateiceconductivity = 0.; + rho_ice = 0.; + rho_water = 0.; + rho_freshwater = 0.; + mu_water = 0.; + heatcapacity = 0.; + latentheat = 0.; + thermalconductivity = 0.; + temperateiceconductivity = 0.; effectiveconductivity_averaging = 0; - meltingpoint = 0.; - beta = 0.; - mixed_layer_capacity = 0.; - thermal_exchange_velocity = 0.; - rheology_B = NaN; - rheology_Ec = NaN; - rheology_Es = NaN; - rheology_law = ''; + meltingpoint = 0.; + beta = 0.; + mixed_layer_capacity = 0.; + thermal_exchange_velocity = 0.; + rheology_B = NaN; + rheology_Ec = NaN; + rheology_Es = NaN; + rheology_law = ''; %slc - earth_density = 0; + earth_density = 0; end Index: /issm/trunk-jpl/src/m/classes/matestar.py =================================================================== --- /issm/trunk-jpl/src/m/classes/matestar.py (revision 26058) +++ /issm/trunk-jpl/src/m/classes/matestar.py (revision 26059) @@ -34,11 +34,5 @@ self.rheology_law = '' - #GIA - self.lithosphere_shear_modulus = 0. - self.lithosphere_density = 0. - self.mantle_shear_modulus = 0. - self.mantle_density = 0. - - #SLR + #slc self.earth_density = 0 @@ -66,8 +60,4 @@ s = "%s\n%s" % (s, fielddisplay(self, 'rheology_Es', 'shear enhancement factor')) s = "%s\n%s" % (s, fielddisplay(self, 'rheology_law', ['law for the temperature dependance of the rheology: ''None'', ''BuddJacka'', ''Cuffey'', ''CuffeyTemperate'', ''Paterson'', ''Arrhenius'' or ''LliboutryDuval'''])) - s = "%s\n%s" % (s, fielddisplay(self, 'lithosphere_shear_modulus', 'Lithosphere shear modulus [Pa]')) - s = "%s\n%s" % (s, fielddisplay(self, 'lithosphere_density', 'Lithosphere density [g/cm^-3]')) - s = "%s\n%s" % (s, fielddisplay(self, 'mantle_shear_modulus', 'Mantle shear modulus [Pa]')) - s = "%s\n%s" % (s, fielddisplay(self, 'mantle_density', 'Mantle density [g/cm^-3]')) s = "%s\n%s" % (s, fielddisplay(self, 'earth_density', 'Mantle density [kg/m^-3]')) @@ -113,10 +103,5 @@ #available: none, paterson and arrhenius self.rheology_law = 'Paterson' - # GIA - self.lithosphere_shear_modulus = 6.7e10 # (Pa) - self.lithosphere_density = 3.32 # (g / cm^ - 3) - self.mantle_shear_modulus = 1.45e11 # (Pa) - self.mantle_density = 3.34 # (g / cm^ - 3) - #SLR + #slc self.earth_density = 5512 # average density of the Earth, (kg / m^3) @@ -134,10 +119,4 @@ md = checkfield(md, 'fieldname', 'materials.rheology_law', 'values', ['None', 'BuddJacka', 'Cuffey', 'CuffeyTemperate', 'Paterson', 'Arrhenius', 'LliboutryDuval']) md = checkfield(md, 'fieldname', 'materials.effectiveconductivity_averaging', 'numel', [1], 'values', [0, 1, 2]) - - if 'GiaAnalysis' in analyses: - md = checkfield(md, 'fieldname', 'materials.lithosphere_shear_modulus', '>', 0, 'numel', 1) - md = checkfield(md, 'fieldname', 'materials.lithosphere_density', '>', 0, 'numel', 1) - md = checkfield(md, 'fieldname', 'materials.mantle_shear_modulus', '>', 0, 'numel', 1) - md = checkfield(md, 'fieldname', 'materials.mantle_density', '>', 0, 'numel', 1) if 'SealevelriseAnalysis' in analyses: @@ -166,8 +145,4 @@ WriteData(fid, prefix, 'object', self, 'class', 'materials', 'fieldname', 'rheology_Es', 'format', 'DoubleMat', 'mattype', 1) WriteData(fid, prefix, 'data', self.rheology_law, 'name', 'md.materials.rheology_law', 'format', 'String') - WriteData(fid, prefix, 'object', self, 'class', 'materials', 'fieldname', 'lithosphere_shear_modulus', 'format', 'Double') - WriteData(fid, prefix, 'object', self, 'class', 'materials', 'fieldname', 'lithosphere_density', 'format', 'Double', 'scale', 1.0e3) - WriteData(fid, prefix, 'object', self, 'class', 'materials', 'fieldname', 'mantle_shear_modulus', 'format', 'Double') - WriteData(fid, prefix, 'object', self, 'class', 'materials', 'fieldname', 'mantle_density', 'format', 'Double', 'scale', 1.0e3) WriteData(fid, prefix, 'object', self, 'class', 'materials', 'fieldname', 'earth_density', 'format', 'Double') # }}} Index: /issm/trunk-jpl/src/m/classes/matice.m =================================================================== --- /issm/trunk-jpl/src/m/classes/matice.m (revision 26058) +++ /issm/trunk-jpl/src/m/classes/matice.m (revision 26059) @@ -104,18 +104,16 @@ end % }}} function md = checkconsistency(self,md,solution,analyses) % {{{ - + if strcmpi(solution,'TransientSolution') & md.transient.isslc, md = checkfield(md,'fieldname','materials.earth_density','>',0,'numel',1); - return; - end - - md = checkfield(md,'fieldname','materials.rho_ice','>',0); - md = checkfield(md,'fieldname','materials.rho_water','>',0); - md = checkfield(md,'fieldname','materials.rho_freshwater','>',0); - md = checkfield(md,'fieldname','materials.mu_water','>',0); - md = checkfield(md,'fieldname','materials.rheology_B','>',0,'universal',1,'NaN',1,'Inf',1); - md = checkfield(md,'fieldname','materials.rheology_n','>',0,'universal',1,'NaN',1,'Inf',1); - md = checkfield(md,'fieldname','materials.rheology_law','values',{'None' 'BuddJacka' 'Cuffey' 'CuffeyTemperate' 'Paterson' 'Arrhenius' 'LliboutryDuval' 'NyeCO2' 'NyeH2O'}); - md = checkfield(md,'fieldname','materials.effectiveconductivity_averaging','numel',[1],'values',[0 1 2]); + else + md = checkfield(md,'fieldname','materials.rho_ice','>',0); + md = checkfield(md,'fieldname','materials.rho_water','>',0); + md = checkfield(md,'fieldname','materials.rho_freshwater','>',0); + md = checkfield(md,'fieldname','materials.mu_water','>',0); + md = checkfield(md,'fieldname','materials.rheology_B','>',0,'universal',1,'NaN',1,'Inf',1); + md = checkfield(md,'fieldname','materials.rheology_n','>',0,'universal',1,'NaN',1,'Inf',1); + md = checkfield(md,'fieldname','materials.rheology_law','values',{'None' 'BuddJacka' 'Cuffey' 'CuffeyTemperate' 'Paterson' 'Arrhenius' 'LliboutryDuval' 'NyeCO2' 'NyeH2O'}); + md = checkfield(md,'fieldname','materials.effectiveconductivity_averaging','numel',[1],'values',[0 1 2]); end % }}} Index: /issm/trunk-jpl/src/m/classes/matice.py =================================================================== --- /issm/trunk-jpl/src/m/classes/matice.py (revision 26058) +++ /issm/trunk-jpl/src/m/classes/matice.py (revision 26059) @@ -8,10 +8,9 @@ class matice(object): - ''' - MATICE class definition + """MATICE class definition Usage: matice = matice() - ''' + """ def __init__(self): #{{{ @@ -33,11 +32,5 @@ self.rheology_law = '' - #giaivins - self.lithosphere_shear_modulus = 0. - self.lithosphere_density = 0. - self.mantle_shear_modulus = 0. - self.mantle_density = 0. - - #slr + #slc self.earth_density = 0 self.setdefaultparameters() @@ -62,8 +55,4 @@ s = "%s\n%s" % (s, fielddisplay(self, "rheology_n", "Glen's flow law exponent")) s = "%s\n%s" % (s, fielddisplay(self, "rheology_law", "law for the temperature dependance of the rheology: 'None', 'BuddJacka', 'Cuffey', 'CuffeyTemperate', 'Paterson', 'Arrhenius', 'LliboutryDuval', 'NyeCO2', or 'NyeH2O'")) - s = "%s\n%s" % (s, fielddisplay(self, "lithosphere_shear_modulus", "Lithosphere shear modulus [Pa]")) - s = "%s\n%s" % (s, fielddisplay(self, "lithosphere_density", "Lithosphere density [g/cm^-3]")) - s = "%s\n%s" % (s, fielddisplay(self, "mantle_shear_modulus", "Mantle shear modulus [Pa]")) - s = "%s\n%s" % (s, fielddisplay(self, "mantle_density", "Mantle density [g/cm^-3]")) s = "%s\n%s" % (s, fielddisplay(self, "earth_density", "Mantle density [kg/m^-3]")) @@ -108,9 +97,7 @@ self.rheology_law = 'Paterson' - # GIA: - self.lithosphere_shear_modulus = 6.7e10 # (Pa) - self.lithosphere_density = 3.32 # (g/cm^-3) - self.mantle_shear_modulus = 1.45e11 # (Pa) - self.mantle_density = 3.34 # (g/cm^-3) + # Rheology for ice + self.rheology_B = 2.1 * 1e8 + self.rheology_n = 3 # SLR @@ -119,5 +106,7 @@ def checkconsistency(self, md, solution, analyses): #{{{ - if solution != 'SealevelriseSolution': + if solution == 'TransientSolution' and md.transient.isslc: + md = checkfield(md, 'fieldname', 'materials.earth_density', '>', 0, 'numel', [1]) + else: md = checkfield(md, 'fieldname', 'materials.rho_ice', '>', 0) md = checkfield(md, 'fieldname', 'materials.rho_water', '>', 0) @@ -125,16 +114,7 @@ md = checkfield(md, 'fieldname', 'materials.mu_water', '>', 0) md = checkfield(md, 'fieldname', 'materials.rheology_B', '>', 0, 'universal', 1, 'NaN', 1, 'Inf', 1) - md = checkfield(md, 'fieldname', 'materials.rheology_n', '>', 0, 'size', [md.mesh.numberofelements]) + md = checkfield(md, 'fieldname', 'materials.rheology_n', '>', 0, 'universal',1, 'NaN', 1, 'Inf', 1) md = checkfield(md, 'fieldname', 'materials.rheology_law', 'values', ['None', 'BuddJacka', 'Cuffey', 'CuffeyTemperate', 'Paterson', 'Arrhenius', 'LliboutryDuval', 'NyeCO2', 'NyeH2O']) - md = checkfield(md, 'fieldname', 'materials.effectiveconductivity_averaging', 'numel', [1], 'values', [0, 1, 2]) - - if 'GiaAnalysis' in analyses: - md = checkfield(md, 'fieldname', 'materials.lithosphere_shear_modulus', '>', 0, 'numel', [1]) - md = checkfield(md, 'fieldname', 'materials.lithosphere_density', '>', 0, 'numel', [1]) - md = checkfield(md, 'fieldname', 'materials.mantle_shear_modulus', '>', 0, 'numel', [1]) - md = checkfield(md, 'fieldname', 'materials.mantle_density', '>', 0, 'numel', [1]) - - if 'SealevelriseAnalysis' in analyses: - md = checkfield(md, 'fieldname', 'materials.earth_density', '>', 0, 'numel', [1]) + md = checkfield(md,'fieldname','materials.effectiveconductivity_averaging', 'numel', [1], 'values', [0, 1, 2]) return md @@ -159,8 +139,4 @@ WriteData(fid, prefix, 'object', self, 'class', 'materials', 'fieldname', 'rheology_n', 'format', 'DoubleMat', 'mattype', 2) WriteData(fid, prefix, 'data', self.rheology_law, 'name', 'md.materials.rheology_law', 'format', 'String') - WriteData(fid, prefix, 'object', self, 'class', 'materials', 'fieldname', 'lithosphere_shear_modulus', 'format', 'Double') - WriteData(fid, prefix, 'object', self, 'class', 'materials', 'fieldname', 'lithosphere_density', 'format', 'Double', 'scale', 10.**3.) - WriteData(fid, prefix, 'object', self, 'class', 'materials', 'fieldname', 'mantle_shear_modulus', 'format', 'Double') - WriteData(fid, prefix, 'object', self, 'class', 'materials', 'fieldname', 'mantle_density', 'format', 'Double', 'scale', 10.**3.) WriteData(fid, prefix, 'object', self, 'class', 'materials', 'fieldname', 'earth_density', 'format', 'Double') #}}} Index: /issm/trunk-jpl/src/m/classes/model.m =================================================================== --- /issm/trunk-jpl/src/m/classes/model.m (revision 26058) +++ /issm/trunk-jpl/src/m/classes/model.m (revision 26059) @@ -44,5 +44,5 @@ love = 0; esa = 0; - sampling = 0; + sampling = 0; autodiff = 0; @@ -156,5 +156,5 @@ if isa(md.frontalforcings,'double'); if(isprop('meltingrate',md.calving) & ~isnan(md.calving.meltingrate)) - disp('Warning: md.calving.meltingrate is now in md.frontalforcings'); + gia disp('Warning: md.calving.meltingrate is now in md.frontalforcings'); end md.frontalforcings=frontalforcings(md.calving); @@ -187,6 +187,6 @@ end end - end - %2021 February 17 + end + %2021 February 17 if isa(md.sampling,'double'); md.sampling=sampling(); end end% }}} @@ -242,5 +242,5 @@ md.love = fourierlove(); md.esa = esa(); - md.sampling = sampling(); + md.sampling = sampling(); md.autodiff = autodiff(); md.inversion = inversion(); @@ -615,6 +615,6 @@ elseif fieldsize(1)==numberofelements1 md2.(model_fields{i}).(object_fields{j})=field(pos_elem,:); - elseif (fieldsize(1)==numberofelements1+1) - md2.(model_fields{i}).(object_fields{j})=[field(pos_elem,:); field(end,:)]; + elseif (fieldsize(1)==numberofelements1+1) + md2.(model_fields{i}).(object_fields{j})=[field(pos_elem,:); field(end,:)]; end end @@ -628,5 +628,5 @@ elseif fieldsize(1)==numberofelements1 md2.(model_fields{i})=field(pos_elem,:); - elseif (fieldsize(1)==numberofelements1+1) + elseif (fieldsize(1)==numberofelements1+1) md2.(model_fields{i})=[field(pos_elem,:); field(end,:)]; end @@ -775,16 +775,16 @@ % loop over time steps for p=1:length(md1.results.(solutionfields{i})) - current = md1.results.(solutionfields{i})(p); - solutionsubfields=fields(current); - for j=1:length(solutionsubfields), + current = md1.results.(solutionfields{i})(p); + solutionsubfields=fields(current); + for j=1:length(solutionsubfields), field=md1.results.(solutionfields{i})(p).(solutionsubfields{j}); if length(field)==numberofvertices1, - md2.results.(solutionfields{i})(p).(solutionsubfields{j})=field(pos_node); + md2.results.(solutionfields{i})(p).(solutionsubfields{j})=field(pos_node); elseif length(field)==numberofelements1, - md2.results.(solutionfields{i})(p).(solutionsubfields{j})=field(pos_elem); + md2.results.(solutionfields{i})(p).(solutionsubfields{j})=field(pos_elem); else - md2.results.(solutionfields{i})(p).(solutionsubfields{j})=field; + md2.results.(solutionfields{i})(p).(solutionsubfields{j})=field; end - end + end end else @@ -1501,5 +1501,5 @@ 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','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')); Index: /issm/trunk-jpl/src/m/classes/model.py =================================================================== --- /issm/trunk-jpl/src/m/classes/model.py (revision 26058) +++ /issm/trunk-jpl/src/m/classes/model.py (revision 26059) @@ -55,6 +55,4 @@ from steadystate import steadystate from transient import transient -from giaivins import giaivins -from giamme import giamme from esa import esa from autodiff import autodiff @@ -114,5 +112,4 @@ self.calving = None self.frontalforcings = None - self.gia = None self.love = None self.esa = None @@ -172,5 +169,4 @@ 'calving', 'frontalforcings', - 'gia', 'love', 'esa', @@ -225,5 +221,4 @@ s = "%s\n%s" % (s, "%19s: %-22s -- %s" % ("calving", "[%s %s]" % ("1x1", obj.calving.__class__.__name__), "parameters for calving")) s = "%s\n%s" % (s, "%19s: %-22s -- %s" % ("frontalforcings", "[%s %s]" % ("1x1", obj.frontalforcings.__class__.__name__), "parameters for frontalforcings")) - s = "%s\n%s" % (s, "%19s: %-22s -- %s" % ("gia", "[%s %s]" % ("1x1", obj.gia.__class__.__name__), "parameters for gia solution")) s = "%s\n%s" % (s, '%19s: %-22s -- %s' % ("esa", "[%s %s]" % ("1x1", obj.esa.__class__.__name__), "parameters for elastic adjustment solution")) s = "%s\n%s" % (s, '%19s: %-22s -- %s' % ("sampling", "[%s %s]" % ("1x1", obj.sampling.__class__.__name__), "parameters for stochastic sampler")) @@ -272,5 +267,4 @@ self.calving = calving() self.frontalforcings = frontalforcings() - self.gia = giamme() self.love = fourierlove() self.esa = esa() @@ -852,11 +846,4 @@ md.initialization.watercolumn = project2d(md, md.initialization.watercolumn, 1) - # giaivins - if md.gia.__class__.__name__ == 'giaivins': - if not np.isnan(md.gia.mantle_viscosity).all(): - md.gia.mantle_viscosity = project2d(md, md.gia.mantle_viscosity, 1) - if not np.isnan(md.gia.lithosphere_thickness).all(): - md.gia.lithosphere_thickness = project2d(md, md.gia.lithosphere_thickness, 1) - # elementstype if not np.isnan(md.flowequation.element_equation).all(): Index: /issm/trunk-jpl/src/m/classes/sealevelmodel.m =================================================================== --- /issm/trunk-jpl/src/m/classes/sealevelmodel.m (revision 26058) +++ /issm/trunk-jpl/src/m/classes/sealevelmodel.m (revision 26059) @@ -91,13 +91,13 @@ md= slm.icecaps{i}; if ~isempty(md.solidearth.external), - error(sprintf('cannot run external forcings on an ice sheet when runing a coupling earth/ice sheet model'); - end - - end - %make sure that we have the rigth grd model for computing ouf sealevel patterns: + error('cannot run external forcings on an ice sheet when running a coupling earth/ice sheet model'); + end + + end + %make sure that we have the right grd model for computing out sealevel patterns: for i=1:length(slm.icecaps), md= slm.icecaps{i}; if md.solidearth.settings.grdmodel~=0 - error(sprintf('sealevelmodel checkconsistenty error message: ice sheets do not run GRD module, specify solidearth.settings.grdmodel=0 on ice cap %i',i)); + error(sprintf('sealevelmodel checkconsistency error message: ice sheets do not run GRD module, specify solidearth.settings.grdmodel=0 on ice cap %i',i)); end end Index: /issm/trunk-jpl/src/m/classes/sealevelmodel.py =================================================================== --- /issm/trunk-jpl/src/m/classes/sealevelmodel.py (revision 26058) +++ /issm/trunk-jpl/src/m/classes/sealevelmodel.py (revision 26059) @@ -111,4 +111,20 @@ if np.nonzero(md.slr.steric_rate - slm.earth.slr.steric_rate[slm.earth.slr.transitions[i]]) != []: raise Exception('steric rate on ice cap {} is not the same as for the earth'.format(md.miscellaneous.name)) + + # Make sure grd is the same everywhere + for i in range(len(slm.icecaps)): + md = slm.icecaps[i] + if md.solidearthsettings.isgrd != slm.earth.solidearthsettings.isgrd: + raise RuntimeError('isgrd on ice cap {} is not the same as for the earth\n'.format(md.miscellaneous.name)) + # Make sure that there is no solid earth external forcing on the basins + for i in range(len(slm.icecaps)): + md = slm.icecaps[i] + if md.solidearth.external: + raise RuntimeError('cannot run external forcings on an ice sheet when running a coupling earth/ice sheet model') + # Make sure that we have the right grd model for computing our sealevel patterns + for i in range(len(slm.icecaps)): + md = slm.icecaps[i] + if md.solidearth.settings.grdmodel != 0: + raise RuntimeError('sealevelmodel checkconsistency error message: ice sheets do not run GRD module, specify solidearth.settings.grdmodel=0 on ice cap {}'.format(i)) #}}} Index: /issm/trunk-jpl/src/m/classes/surfaceload.m =================================================================== --- /issm/trunk-jpl/src/m/classes/surfaceload.m (revision 26058) +++ /issm/trunk-jpl/src/m/classes/surfaceload.m (revision 26059) @@ -20,11 +20,11 @@ end % }}} function self = setdefaultparameters(self) % {{{ - + icethicknesschange=[]; waterheightchange=[]; otherchange=[]; - + end % }}} - function md = checkconsistency(self,md,solution,analyses) % {{{ + function md = checkconsistency(self,md,solution,analyses) % {{{ if ~ismember('SealevelchangeAnalysis',analyses) | (strcmp(solution,'TransientSolution') & md.transient.isslc==0), Index: /issm/trunk-jpl/src/m/classes/transient.m =================================================================== --- /issm/trunk-jpl/src/m/classes/transient.m (revision 26058) +++ /issm/trunk-jpl/src/m/classes/transient.m (revision 26059) @@ -8,5 +8,5 @@ issmb = 0; ismasstransport = 0; - isoceantransport = 0; + isoceantransport = 0; isstressbalance = 0; isthermal = 0; @@ -16,5 +16,5 @@ ismovingfront = 0; ishydrology = 0; - issampling = 0; + issampling = 0; isslc = 0; amr_frequency = 0; @@ -44,5 +44,5 @@ self.ismovingfront =0; self.ishydrology = 0; - self.issampling = 0; + self.issampling = 0; self.isslc = 0; self.isoceancoupling = 0; @@ -65,5 +65,5 @@ self.ismovingfront = 0; self.ishydrology = 0; - self.issampling = 0; + self.issampling = 0; self.isslc = 0; self.isoceancoupling = 0; @@ -98,5 +98,5 @@ md = checkfield(md,'fieldname','transient.isslc','numel',[1],'values',[0 1]); md = checkfield(md,'fieldname','transient.isoceancoupling','numel',[1],'values',[0 1]); - md = checkfield(md,'fieldname','transient.issampling','numel',[1],'values',[0 1]); + md = checkfield(md,'fieldname','transient.issampling','numel',[1],'values',[0 1]); md = checkfield(md,'fieldname','transient.amr_frequency','numel',[1],'>=',0,'NaN',1,'Inf',1); @@ -121,5 +121,5 @@ fielddisplay(self,'ismovingfront','indicates whether a moving front capability is used in the transient'); fielddisplay(self,'ishydrology','indicates whether an hydrology model is used'); - fielddisplay(self,'issampling','indicates whether sampling is used in the transient') + fielddisplay(self,'issampling','indicates whether sampling is used in the transient') fielddisplay(self,'isslc','indicates whether a sea-level change solution is used in the transient'); fielddisplay(self,'isoceancoupling','indicates whether a coupling with an ocean model is used in the transient'); @@ -139,5 +139,5 @@ WriteData(fid,prefix,'object',self,'fieldname','ishydrology','format','Boolean'); WriteData(fid,prefix,'object',self,'fieldname','ismovingfront','format','Boolean'); - WriteData(fid,prefix,'object',self,'fieldname','issampling','format','Boolean'); + WriteData(fid,prefix,'object',self,'fieldname','issampling','format','Boolean'); WriteData(fid,prefix,'object',self,'fieldname','isslc','format','Boolean'); WriteData(fid,prefix,'object',self,'fieldname','isoceancoupling','format','Boolean'); @@ -165,5 +165,5 @@ writejsdouble(fid,[modelname '.trans.ismovingfront'],self.ismovingfront); writejsdouble(fid,[modelname '.trans.ishydrology'],self.ishydrology); - writejsdouble(fid,[modelname '.trans.issampling'],self.issampling); + writejsdouble(fid,[modelname '.trans.issampling'],self.issampling); writejsdouble(fid,[modelname '.trans.isslc'],self.isslc); writejsdouble(fid,[modelname '.trans.isoceancoupling'],self.isoceancoupling); Index: /issm/trunk-jpl/src/m/classes/transient.py =================================================================== --- /issm/trunk-jpl/src/m/classes/transient.py (revision 26058) +++ /issm/trunk-jpl/src/m/classes/transient.py (revision 26059) @@ -14,8 +14,8 @@ self.issmb = 0 self.ismasstransport = 0 + self.isoceantransport = 0 self.isstressbalance = 0 self.isthermal = 0 self.isgroundingline = 0 - self.isgia = 0 self.isesa = 0 self.isdamageevolution = 0 @@ -24,5 +24,4 @@ self.issampling = 0 self.isslc = 0 - self.iscoupler = 0 self.amr_frequency = 0 self.isoceancoupling = 0 @@ -36,8 +35,8 @@ s += '{}\n'.format(fielddisplay(self, 'issmb', 'indicates if a surface mass balance solution is used in the transient')) s += '{}\n'.format(fielddisplay(self, 'ismasstransport', 'indicates if a masstransport solution is used in the transient')) + s += '{}\n'.format(fielddisplay(self, 'isoceantransport', 'indicates whether an ocean masstransport solution is used in the transient')) s += '{}\n'.format(fielddisplay(self, 'isstressbalance', 'indicates if a stressbalance solution is used in the transient')) s += '{}\n'.format(fielddisplay(self, 'isthermal', 'indicates if a thermal solution is used in the transient')) s += '{}\n'.format(fielddisplay(self, 'isgroundingline', 'indicates if a groundingline migration is used in the transient')) - s += '{}\n'.format(fielddisplay(self, 'isgia', 'indicates if a postglacial rebound is used in the transient')) s += '{}\n'.format(fielddisplay(self, 'isesa', 'indicates whether an elastic adjustment model is used in the transient')) s += '{}\n'.format(fielddisplay(self, 'isdamageevolution', 'indicates whether damage evolution is used in the transient')) @@ -47,5 +46,4 @@ s += '{}\n'.format(fielddisplay(self, 'isslc', 'indicates if a sea level change solution is used in the transient')) s += '{}\n'.format(fielddisplay(self, 'isoceancoupling', 'indicates whether coupling with an ocean model is used in the transient')) - s += '{}\n'.format(fielddisplay(self, 'iscoupler', 'indicates whether different models are being run with need for coupling')) s += '{}\n'.format(fielddisplay(self, 'amr_frequency', 'frequency at which mesh is refined in simulations with multiple time_steps')) s += '{}\n'.format(fielddisplay(self, 'requested_outputs', 'list of additional outputs requested')) @@ -63,8 +61,8 @@ self.issmb = 0 self.ismasstransport = 0 + self.isoceantransport = 0 self.isstressbalance = 0 self.isthermal = 0 self.isgroundingline = 0 - self.isgia = 0 self.isesa = 0 self.isdamageevolution = 0 @@ -74,5 +72,4 @@ self.isslc = 0 self.isoceancoupling = 0 - self.iscoupler = 0 self.amr_frequency = 0 @@ -86,8 +83,8 @@ self.issmb = 1 self.ismasstransport = 1 + self.isoceantransport = 0 self.isstressbalance = 1 self.isthermal = 1 self.isgroundingline = 0 - self.isgia = 0 self.isesa = 0 self.isdamageevolution = 0 @@ -97,5 +94,4 @@ self.isslc = 0 self.isoceancoupling = 0 - self.iscoupler = 0 self.amr_frequency = 0 @@ -112,8 +108,8 @@ md = checkfield(md, 'fieldname', 'transient.issmb', 'numel', [1], 'values', [0, 1]) md = checkfield(md, 'fieldname', 'transient.ismasstransport', 'numel', [1], 'values', [0, 1]) + md = checkfield(md, 'fieldname', 'transient.isocealtransport', 'numel', [1], 'values', [0, 1]) md = checkfield(md, 'fieldname', 'transient.isstressbalance', 'numel', [1], 'values', [0, 1]) md = checkfield(md, 'fieldname', 'transient.isthermal', 'numel', [1], 'values', [0, 1]) md = checkfield(md, 'fieldname', 'transient.isgroundingline', 'numel', [1], 'values', [0, 1]) - md = checkfield(md, 'fieldname', 'transient.isgia', 'numel', [1], 'values', [0, 1]) md = checkfield(md, 'fieldname', 'transient.isesa', 'numel', [1], 'values', [0, 1]) md = checkfield(md, 'fieldname', 'transient.isdamageevolution', 'numel', [1], 'values', [0, 1]) @@ -123,5 +119,4 @@ md = checkfield(md, 'fieldname', 'transient.isslc', 'numel', [1], 'values', [0, 1]) md = checkfield(md, 'fieldname', 'transient.isoceancoupling', 'numel', [1], 'values', [0, 1]) - md = checkfield(md, 'fieldname', 'transient.iscoupler', 'numel', [1], 'values', [0, 1]) md = checkfield(md, 'fieldname', 'transient.amr_frequency', 'numel', [1], '>=', 0, 'NaN', 1, 'Inf', 1) @@ -136,8 +131,8 @@ WriteData(fid, prefix, 'object', self, 'fieldname', 'issmb', 'format', 'Boolean') WriteData(fid, prefix, 'object', self, 'fieldname', 'ismasstransport', 'format', 'Boolean') + WriteData(fid, prefix, 'object', self, 'fieldname', 'isoceantransport', 'format', 'Boolean') WriteData(fid, prefix, 'object', self, 'fieldname', 'isstressbalance', 'format', 'Boolean') WriteData(fid, prefix, 'object', self, 'fieldname', 'isthermal', 'format', 'Boolean') WriteData(fid, prefix, 'object', self, 'fieldname', 'isgroundingline', 'format', 'Boolean') - WriteData(fid, prefix, 'object', self, 'fieldname', 'isgia', 'format', 'Boolean') WriteData(fid, prefix, 'object', self, 'fieldname', 'isesa', 'format', 'Boolean') WriteData(fid, prefix, 'object', self, 'fieldname', 'isdamageevolution', 'format', 'Boolean') @@ -147,5 +142,4 @@ WriteData(fid, prefix, 'object', self, 'fieldname', 'isslc', 'format', 'Boolean') WriteData(fid, prefix, 'object', self, 'fieldname', 'isoceancoupling', 'format', 'Boolean') - WriteData(fid, prefix, 'object', self, 'fieldname', 'iscoupler', 'format', 'Boolean') WriteData(fid, prefix, 'object', self, 'fieldname', 'amr_frequency', 'format', 'Integer') Index: /issm/trunk-jpl/src/m/consistency/ismodelselfconsistent.m =================================================================== --- /issm/trunk-jpl/src/m/consistency/ismodelselfconsistent.m (revision 26058) +++ /issm/trunk-jpl/src/m/consistency/ismodelselfconsistent.m (revision 26059) @@ -3,5 +3,5 @@ % % Usage: -% ismodelselfconsistent(md), +% ismodelselfconsistent(md); %initialize consistency as true Index: /issm/trunk-jpl/src/m/consistency/ismodelselfconsistent.py =================================================================== --- /issm/trunk-jpl/src/m/consistency/ismodelselfconsistent.py (revision 26058) +++ /issm/trunk-jpl/src/m/consistency/ismodelselfconsistent.py (revision 26059) @@ -1,9 +1,8 @@ def ismodelselfconsistent(md): #{{{ - ''' - ISMODELSELFCONSISTENT - check that model forms a closed form solvable problem. + """ISMODELSELFCONSISTENT - check that model forms a closed form solvable problem. - Usage: - ismodelselfconsistent(md), - ''' + Usage: + ismodelselfconsistent(md) + """ #initialize consistency as true @@ -53,4 +52,6 @@ elif solutiontype == 'MasstransportSolution': analyses = ['MasstransportAnalysis'] + elif solutiontype == 'OceantransportSolution': + analyses = ['OceantransportAnalysis'] elif solutiontype == 'BalancethicknessSolution': analyses = ['BalancethicknessAnalysis'] @@ -72,9 +73,9 @@ analyses = ['EsaAnalysis'] elif solutiontype == 'TransientSolution': - analyses = ['StressbalanceAnalysis', 'StressbalanceVerticalAnalysis', 'StressbalanceSIAAnalysis', 'L2ProjectionBaseAnalysis', 'ThermalAnalysis', 'MeltingAnalysis', 'EnthalpyAnalysis', 'MasstransportAnalysis', 'HydrologyShaktiAnalysis', 'HydrologyGladsAnalysis', 'HydrologyDCInefficientAnalysis', 'HydrologyDCEfficientAnalysis', 'SealevelriseAnalysis'] + analyses = ['StressbalanceAnalysis', 'StressbalanceVerticalAnalysis', 'StressbalanceSIAAnalysis', 'L2ProjectionBaseAnalysis', 'ThermalAnalysis', 'MeltingAnalysis', 'EnthalpyAnalysis', 'MasstransportAnalysis', 'OceantransportAnalysis', 'HydrologyShaktiAnalysis', 'HydrologyGladsAnalysis', 'HydrologyShreveAnalysis', 'HydrologyTwsAnalysis', 'HydrologyDCInefficientAnalysis', 'HydrologyDCEfficientAnalysis', 'SealevelriseAnalysis'] elif solutiontype == 'SealevelriseSolution': analyses = ['SealevelriseAnalysis'] elif solutiontype == 'HydrologySolution': - analyses = ['L2ProjectionBaseAnalysis', 'HydrologyShreveAnalysis', 'HydrologyDCInefficientAnalysis', 'HydrologyDCEfficientAnalysis'] + analyses = ['L2ProjectionBaseAnalysis', 'HydrologyShreveAnalysis', 'HydrologyDCInefficientAnalysis', 'HydrologyDCEfficientAnalysis', 'HydrologyGladsAnalysis', 'HydrologyShaktiAnalysis', 'HydrologyTwsAnalysis'] elif 'DamageEvolutionSolution': analyses = ['DamageEvolutionAnalysis'] Index: /issm/trunk-jpl/src/m/plot/radarpower.py =================================================================== --- /issm/trunk-jpl/src/m/plot/radarpower.py (revision 26059) +++ /issm/trunk-jpl/src/m/plot/radarpower.py (revision 26059) @@ -0,0 +1,16 @@ +def solveslm(slm, solutionstringi, *args): + """RADARPOWER - overlay a power radar image on an existing mesh + + This routine will overlay a power radar image on an existing mesh. + The power amplitude will be output to vel for now. + In the future, think about a field to hold this value. + + Usage: + md=radarpower(md,options); + md=radarpower(md) + + TODO: + - Translate from MATLAB API as we bring Python plotting capabilities online + """ + + return md Index: /issm/trunk-jpl/src/m/solve/solve.m =================================================================== --- /issm/trunk-jpl/src/m/solve/solve.m (revision 26058) +++ /issm/trunk-jpl/src/m/solve/solve.m (revision 26059) @@ -10,5 +10,5 @@ % - 'Stressbalance' or 'sb' % - 'Masstransport' or 'mt' -% - 'Oceantransport' or 'oceant' +% - 'Oceantransport' or 'oceant' % - 'Thermal' or 'th' % - 'Steadystate' or 'ss' Index: /issm/trunk-jpl/src/m/solve/solve.py =================================================================== --- /issm/trunk-jpl/src/m/solve/solve.py (revision 26058) +++ /issm/trunk-jpl/src/m/solve/solve.py (revision 26059) @@ -1,5 +1,4 @@ from datetime import datetime import os -import shutil from ismodelselfconsistent import ismodelselfconsistent @@ -22,4 +21,5 @@ - 'Stressbalance' or 'sb' - 'Masstransport' or 'mt' + - 'Oceantransport' or 'oceant' - 'Thermal' or 'th' - 'Steadystate' or 'ss' @@ -32,7 +32,6 @@ - 'DamageEvolution' or 'da' - 'Gia' or 'gia' + - 'Love' or 'lv' - 'Esa' or 'esa' - - 'Sealevelchange' or 'slc' - - 'Love' or 'lv' - 'Sampling' or 'smp' @@ -55,4 +54,6 @@ elif solutionstring.lower() == 'mt' or solutionstring.lower() == 'masstransport': solutionstring = 'MasstransportSolution' + elif solutionstring.lower() == 'oceant' or solutionstring.lower() == 'oceantransport': + solutionstring = 'OceantransportSolution' elif solutionstring.lower() == 'th' or solutionstring.lower() == 'thermal': solutionstring = 'ThermalSolution' @@ -79,6 +80,4 @@ elif solutionstring.lower() == 'esa': solutionstring = 'EsaSolution' - elif solutionstring.lower() == 'slc' or solutionstring.lower() == 'sealevelchange': - solutionstring = 'SealevelchangeSolution' elif solutionstring.lower() == 'smp' or solutionstring.lower() == 'sampling': solutionstring = 'SamplingSolution' @@ -151,7 +150,7 @@ if md.settings.waitonlock > 0: islock = waitonlock(md) - if islock == 0: # no results to be loaded + if islock == 0: # no results to be loaded print('The results must be loaded manually with md = loadresultsfromcluster(md).') - else: # load results + else: # load results if md.verbose.solution: print('loading results from cluster') Index: /issm/trunk-jpl/src/m/solve/solveslm.m =================================================================== --- /issm/trunk-jpl/src/m/solve/solveslm.m (revision 26058) +++ /issm/trunk-jpl/src/m/solve/solveslm.m (revision 26059) @@ -1,4 +1,4 @@ function slm=solveslm(slm,solutionstringi,varargin) -%SOLVESLR - apply solution sequence for this sealevel model +%SOLVESLM - apply solution sequence for this sealevel model % % Usage: @@ -43,5 +43,5 @@ cluster=slm.cluster; batch=0; -%now, go through icecaps, glacies and earth, and upload all the data independently: +%now, go through icecaps, glaciers and earth, and upload all the data independently: disp('solving ice caps first'); for i=1:length(slm.icecaps), @@ -51,5 +51,5 @@ slm.earth=solve(slm.earth,solutionstringi,'batch','yes'); -%Firs, build a runtime name that is unique +%First, build a runtime name that is unique c=clock; slm.private.runtimename=sprintf('%s-%02i-%02i-%04i-%02i-%02i-%02i-%i',slm.miscellaneous.name,c(2),c(3),c(1),c(4),c(5),floor(c(6)),feature('GetPid')); Index: /issm/trunk-jpl/src/m/solve/solveslm.py =================================================================== --- /issm/trunk-jpl/src/m/solve/solveslm.py (revision 26059) +++ /issm/trunk-jpl/src/m/solve/solveslm.py (revision 26059) @@ -0,0 +1,95 @@ +from datetime import datetime +import os + +import numpy as np + +from loadresultsfromcluster import loadresultsfromcluster +from pairoptions import pairoptions +from waitonlock import waitonlock + + +def solveslm(slm, solutionstringi, *args): + """SOLVESLM - apply solution sequence for this sealevel model + + Usage: + slm=solve(slm,solutionstring,varargin) + where varargin is a lit of paired arguments of string OR enums + + solution types available comprise: + - 'Transient' + + extra options: + + Examples: + slm=solve(slm,'Transient'); + """ + + # Recover and process solve options + if solutionstringi.lower() == 'tr' or solutionstringi.lower() == 'transient': + solutionstring = 'TransientSolution' + else: + raise RuntimeError('solutionstring {} not supported!'.format(solutionstringi)) + + # Default settings for debugging + valgrind = 0 + #slm.cluster.interactive = 0 + #valgrind = 1 + + # Check consistency + slm.checkconsistency(solutionstring) + + # Process options + options = pairoptions('solutionstring', solutionstring, *args) + + # Make sure we request sum of cluster processors + totalnp = 0 + for i in range(len(slm.icecaps)): + totalnp = totalnp + slm.icecaps[i].cluster.np + totalnp = totalnp + slm.earth.cluster.np + if totalnp != slm.cluster.np: + raise RuntimeError('sum of all icecaps and earch cluster processors requestes should be equal to slm.cluster.np') + + # Recover some fields + slm.private.solution = solutionstring + cluster = slm.cluster + batch = 0 + # Now, go through icecaps, glaciers and earth, and upload all the data independently + print('solving ice caps first') + for i in range(len(slm.icecaps)): + slm.icecaps[i] = solve(slm.icecaps[i], solutionastringi,'batch','yes') + print('solving earth now') + slm.earth = solve(slm.earth, solutionstringi, 'batch', 'yes') + + # First, build a runtime name that is unique + c = datetime.now() + md.private.runtimename = "%s-%02i-%02i-%04i-%02i-%02i-%02i-%i" % (md.miscellaneous.name, c.month, c.day, c.year, c.hour, c.minute, c.second, os.getpid()) + + # Write all input files + privateruntimenames = [] + miscellaneousnames = [] + nps = [] + for i in range(len(slm.icecaps)): + privateruntimenames.append(slm.icecaps[i],private.runtimename) + miscellaneousnames.append(slm.earth.miscellaneous.name) + nps.append(slm.earth.cluster.np) + + BuildQueueScriptMultipleModels(cluster, slm.private.runtimename, slm.miscellaneous.name, slm.private.solution, valgrind, privateruntimenames, miscellaneousnames, nps) + + # Upload all required files, given that each individual solution for icecaps and earth model already did + filelist = [slm.miscellaneous.name + '.queue'] + UploadQueueJob(cluster, slm.miscellaneous.name, slm.private.runtimename, filelist) + + # Launch queue job + LaunchQueueJob(cluster, slm.miscellaneous.name, slm.private.runtimename, filelist, '', batch) + + # Wait on lock + if slm.settings.waitonlock > 0: + islock = waitonlock(slm) + if islock == 0: # no results to be loaded + print('The results must be loaded manually with md = loadresultsfromcluster(md).') + else: # load results + if slm.verbose.solution: + print('loading results from cluster') + slm = loadresultsfromcluster(slm) + + return slm Index: /issm/trunk-jpl/test/NightlyRun/test2001.m =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test2001.m (revision 26058) +++ /issm/trunk-jpl/test/NightlyRun/test2001.m (revision 26059) @@ -22,6 +22,6 @@ md.timestepping.start_time=-2400000; %4,800 kyr :: EVALUATION TIME md.timestepping.time_step= 2400000; %2,400 kyr :: EVALUATION TIME -% to get rid of default final_time: make sure final_time>start_time -md.timestepping.final_time=2400000; %2,500 kyr +% to get rid of default final_time, make sure final_time > start_time +md.timestepping.final_time=2400000; %2,400 kyr md.masstransport.spcthickness=[... [md.geometry.thickness; 0],... @@ -46,5 +46,5 @@ md.verbose=verbose('11111111111'); md.verbose.solver=0; -md=solve(md,'tr'); +md=solve(md,'Transient'); %Fields and tolerances to track changes Index: /issm/trunk-jpl/test/NightlyRun/test2001.py =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test2001.py (revision 26058) +++ /issm/trunk-jpl/test/NightlyRun/test2001.py (revision 26059) @@ -4,5 +4,5 @@ import numpy as np -from giaivins import giaivins +from materials import * from model import * from parameterize import * @@ -17,34 +17,49 @@ md = parameterize(md, '../Par/SquareSheetConstrained.py') -#GIA -md.gia = giaivins() -md.gia.lithosphere_thickness = 100. * np.ones(md.mesh.numberofvertices) # in km -md.gia.mantle_viscosity = 1.0e21 * np.ones(md.mesh.numberofvertices) # in Pa.s -md.materials.lithosphere_shear_modulus = 6.7e10 # in Pa -md.materials.lithosphere_density = 3.32 # in g/cm^3 -md.materials.mantle_shear_modulus = 1.45e11 # in Pa -md.materials.mantle_density = 3.34 # in g/cm^3 +# GIA Ivins, 2 layer model +md.solidearth.settings.grdmodel = 2 -#Indicate what you want to compute -md.gia.cross_section_shape = 1 # for square-edged x-section +md.materials = materials('litho','ice') +md.materials.numlayers = 2; +md.materials.radius = [10, 6271e3, 6371e3] +md.materials.density = [3.34e3, 3.32e3] +md.materials.lame_mu = [1.45e11, 6.7e10] +md.materials.viscosity = [1e21, 0] +md.initialization.sealevel = np.zeros(md.mesh.numberofvertices) +md.solidearth.settings.cross_section_shape = 1 # for square-edged x-section +md.solidearth.settings.computesealevelchange = 0 # do not compute sea level, only deformation +md.solidearth.requested_outputs = ['Sealevel', 'SealevelUGrd'] -#Define loading history (see test2001.m for the description) -md.timestepping.start_time = 2400000 # 2, 400 kyr -md.timestepping.final_time = 2500000 # 2, 500 kyr -md.geometry.thickness = np.array([ - np.append(md.geometry.thickness * 0.0, 0.0), - np.append(md.geometry.thickness / 2.0, 0.1), - np.append(md.geometry.thickness, 0.2), - np.append(md.geometry.thickness, 1.0), - np.append(md.geometry.thickness, md.timestepping.start_time) +# Loading history +md.timestepping.start_time = -2400000 # 4,800 kyr :: EVALUATION TIME +md.timestepping.time_step = 2400000 # 2,400 kyr :: EVALUATION TIME +# To get rid of default final_time, make sure final_time > start_time +md.timestepping.final_time = 2400000 # 2,400 kyr +md.masstransport.spcthickness np.array([ + np.append(md.geometry.thickness, 0), + np.append(md.geometry.thickness, 2400000) ]).T +# Geometry at 0 initially +md.geometry.thickness = np.zeros(md.mesh.numberofvertices) +md.geometry.surface = np.zeros(md.mesh.numberofvertices) +md.geometry.base = np.zeros(md.mesh.numberofvertices) + +# Physics +md.transient.issmb = 0 +md.transient.isstressbalance = 0 +md.transient.isthermal = 0 +md.transient.ismasstransport = 0 +md.transient.isslc = 0 + #Solve for GIA deflection +md.cluster = generic('name', gethostname(), 'np', 1) md.cluster = generic('name', gethostname(), 'np', 3) -md.verbose = verbose('1111111') -md = solve(md, 'Gia') +md.verbose = verbose('11111111111') +md.verbose.solver = 0 +md = solve(md, 'Transient') #Fields and tolerances to track changes -field_names = ['UGia', 'UGiaRate'] -field_tolerances = [1e-13, 1e-13] -field_values = [md.results.GiaSolution.UGia, md.results.GiaSolution.UGiaRate] +field_names = ['UGrd'] +field_tolerances = [1e-13] +field_values = [md.results.TransientSolution[0].SealevelUGrd]