Index: ../trunk-jpl/src/m/classes/matdamageice.py =================================================================== --- ../trunk-jpl/src/m/classes/matdamageice.py (revision 25168) +++ ../trunk-jpl/src/m/classes/matdamageice.py (revision 25169) @@ -129,6 +129,16 @@ md = checkfield(md, 'fieldname', 'materials.mantle_shear_modulus', '>', 0, 'numel', [1]) 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: + md = checkfield(md, 'fieldname', 'materials.earth_density', '>', 0, 'numel', 1) + return md # }}} Index: ../trunk-jpl/src/m/classes/toolkits.m =================================================================== --- ../trunk-jpl/src/m/classes/toolkits.m (revision 25168) +++ ../trunk-jpl/src/m/classes/toolkits.m (revision 25169) @@ -4,170 +4,176 @@ % self=toolkits(); classdef toolkits < dynamicprops - properties (SetAccess=public) - DefaultAnalysis = struct(); - RecoveryAnalysis = struct(); - %The other properties are dynamic - end - methods (Static) - function self = loadobj(self) % {{{ - % This function is directly called by matlab when a model object is - % loaded. Update old properties here + properties (SetAccess=public) + DefaultAnalysis = struct(); + RecoveryAnalysis = struct(); + %The other properties are dynamic + end + methods (Static) + function self = loadobj(self) % {{{ + % This function is directly called by matlab when a model object is + % loaded. Update old properties here - if isempty(fieldnames(self.RecoveryAnalysis)); - disp('WARNING: updating toolkits (RecoveryAnalysis not set)'); - self.RecoveryAnalysis = self.DefaultAnalysis; - end - end% }}} - end - methods - function self = toolkits(varargin) % {{{ - switch nargin - case 0 - self=setdefaultparameters(self); - case 1 - self=structtoobj(self,varargin{1}); - otherwise - error('constructor not supported'); - end - end % }}} - function self = addoptions(self,analysis,varargin) % {{{ - % Usage example: - % md.toolkits=addoptions(md.toolkits,'StressbalanceAnalysis',FSoptions()); - % md.toolkits=addoptions(md.toolkits,'StressbalanceAnalysis'); + if isempty(fieldnames(self.RecoveryAnalysis)); + disp('WARNING: updating toolkits (RecoveryAnalysis not set)'); + self.RecoveryAnalysis = self.DefaultAnalysis; + end + end% }}} + end + methods + function self = toolkits(varargin) % {{{ + switch nargin + case 0 + self=setdefaultparameters(self); + case 1 + self=structtoobj(self,varargin{1}); + otherwise + error('constructor not supported'); + end + end % }}} + function self = addoptions(self,analysis,varargin) % {{{ + %ADDOPTIONS - add analysis to md.toolkits.analysis + % + % Optional third parameter adds toolkits options to analysis. + % + % Usage: + % md.toolkits=addoptions(md.toolkits,'StressbalanceAnalysis',FSoptions()); + % md.toolkits=addoptions(md.toolkits,'StressbalanceAnalysis'); - %Create dynamic property if property does not exist yet - if ~ismember(analysis,properties(self)), - self.addprop(analysis); - end + %Create dynamic property if property does not exist yet + if ~ismember(analysis,properties(self)), + self.addprop(analysis); + end - %Add toolkits options to analysis - if nargin==3, self.(analysis) = varargin{1}; end - end - %}}} - function self = setdefaultparameters(self) % {{{ + %Add toolkits options to analysis + if nargin==3, + self.(analysis) = varargin{1}; + end + end + %}}} + function self = setdefaultparameters(self) % {{{ - %default toolkits: - if IssmConfig('_HAVE_PETSC_'), - %MUMPS is the default toolkits - if IssmConfig('_HAVE_MUMPS_'), - self.DefaultAnalysis = mumpsoptions(); - else - self.DefaultAnalysis = iluasmoptions(); - end - else - if IssmConfig('_HAVE_MUMPS_'), - self.DefaultAnalysis = issmmumpssolver(); - elseif IssmConfig('_HAVE_GSL_'), - self.DefaultAnalysis = issmgslsolver(); - else - disp('WARNING: Need at least Mumps or Gsl to define an issm solver type, no default solver assigned'); - end - end + %default toolkits: + if IssmConfig('_HAVE_PETSC_'), + %MUMPS is the default toolkits + if IssmConfig('_HAVE_MUMPS_'), + self.DefaultAnalysis = mumpsoptions(); + else + self.DefaultAnalysis = iluasmoptions(); + end + else + if IssmConfig('_HAVE_MUMPS_'), + self.DefaultAnalysis = issmmumpssolver(); + elseif IssmConfig('_HAVE_GSL_'), + self.DefaultAnalysis = issmgslsolver(); + else + disp('WARNING: Need at least Mumps or Gsl to define an issm solver type, no default solver assigned'); + end + end - %Use same solver for Recovery mode - self.RecoveryAnalysis = self.DefaultAnalysis; + %Use same solver for Recovery mode + self.RecoveryAnalysis = self.DefaultAnalysis; - end % }}} - function disp(self) % {{{ - analyses=properties(self); - disp(sprintf('List of toolkits options per analysis:\n')); - for i=1:numel(analyses), - analysis=analyses{i}; - disp([analysis ':']); - disp(self.(analysis)); - end - end % }}} - function md = checkconsistency(self,md,solution,analyses) % {{{ - analyses=properties(self); - for i=1:numel(analyses), - switch analyses{i} - case 'DefaultAnalysis' - case 'RecoveryAnalysis' - case 'StressbalanceAnalysis' - case 'GLheightadvectionAnalysis' - case 'MasstransportAnalysis' - case 'ThermalAnalysis' - case 'EnthalpyAnalysis' - case 'AdjointBalancethicknessAnalysis' - case 'BalancethicknessAnalysis' - case 'Balancethickness2Analysis' - case 'BalancethicknessSoftAnalysis' - case 'BalancevelocityAnalysis' - case 'DamageEvolutionAnalysis' - case 'LoveAnalysis' - case 'EsaAnalysis' - case 'SealevelriseAnalysis' - otherwise + end % }}} + function disp(self) % {{{ + analyses=properties(self); + disp(sprintf('List of toolkits options per analysis:\n')); + for i=1:numel(analyses), + analysis=analyses{i}; + disp([analysis ':']); + disp(self.(analysis)); + end + end % }}} + function md = checkconsistency(self,md,solution,analyses) % {{{ + analyses=properties(self); + for i=1:numel(analyses), + switch analyses{i} + case 'DefaultAnalysis' + case 'RecoveryAnalysis' + case 'StressbalanceAnalysis' + case 'GLheightadvectionAnalysis' + case 'MasstransportAnalysis' + case 'ThermalAnalysis' + case 'EnthalpyAnalysis' + case 'AdjointBalancethicknessAnalysis' + case 'BalancethicknessAnalysis' + case 'Balancethickness2Analysis' + case 'BalancethicknessSoftAnalysis' + case 'BalancevelocityAnalysis' + case 'DamageEvolutionAnalysis' + case 'LoveAnalysis' + case 'EsaAnalysis' + case 'SealevelriseAnalysis' + otherwise md = checkmessage(md,['md.toolkits.' analyses{i} ' not supported yet']); - end - if isempty(fieldnames(self.(analyses{i}))) - md = checkmessage(md,['md.toolkits.' analyses{i} ' is empty']); - end - end - end % }}} - function ToolkitsFile(toolkits,filename) % {{{ - %TOOLKITSFILE - build toolkits file - % - % Build a Petsc compatible options file, from the toolkits model field + return options string. - % This file will also be used when the toolkit used is 'issm' instead of 'petsc' - % - % Usage: ToolkitsFile(toolkits,filename); + end + if isempty(fieldnames(self.(analyses{i}))) + md = checkmessage(md,['md.toolkits.' analyses{i} ' is empty']); + end + end + end % }}} + function ToolkitsFile(toolkits,filename) % {{{ + %TOOLKITSFILE - build toolkits file + % + % Build a Petsc compatible options file, from the toolkits model field + return options string. + % This file will also be used when the toolkit used is 'issm' instead of 'petsc' + % + % Usage: ToolkitsFile(toolkits,filename); - %open file for writing - fid=fopen(filename,'w'); - if fid==-1, - error(['ToolkitsFile error: could not open ' filename ' for writing']); - end + %open file for writing + fid=fopen(filename,'w'); + if fid==-1, + error(['ToolkitsFile error: could not open ' filename ' for writing']); + end - %write header - fprintf(fid,'%s%s%s\n','%Toolkits options file: ',filename,' written from Matlab toolkits array'); + %write header + fprintf(fid,'%s%s%s\n','%Toolkits options file: ',filename,' written from Matlab toolkits array'); - %start writing options - analyses=properties(toolkits); - for i=1:numel(analyses), - analysis=analyses{i}; - options=toolkits.(analysis); + %start writing options + analyses=properties(toolkits); + for i=1:numel(analyses), + analysis=analyses{i}; + options=toolkits.(analysis); - %first write analysis: - fprintf(fid,'\n+%s\n',analysis); %append a + to recognize it's an analysis enum + %first write analysis: + fprintf(fid,'\n+%s\n',analysis); %append a + to recognize it's an analysis enum - %now, write options - optionslist=fieldnames(options); - for j=1:numel(optionslist), - optionname=optionslist{j}; - optionvalue=options.(optionname); + %now, write options + optionslist=fieldnames(options); + for j=1:numel(optionslist), + optionname=optionslist{j}; + optionvalue=options.(optionname); - if isempty(optionvalue), - %this option has only one argument - fprintf(fid,'-%s\n',optionname); - else - %option with value. value can be string or scalar - if isnumeric(optionvalue), - fprintf(fid,'-%s %g\n',optionname,optionvalue); - elseif ischar(optionvalue), - fprintf(fid,'-%s %s\n',optionname,optionvalue); - else - error(['ToolkitsFile error: option ' optionname ' is not well formatted']); - end - end - end - end + if isempty(optionvalue), + %this option has only one argument + fprintf(fid,'-%s\n',optionname); + else + %option with value. value can be string or scalar + if isnumeric(optionvalue), + fprintf(fid,'-%s %g\n',optionname,optionvalue); + elseif ischar(optionvalue), + fprintf(fid,'-%s %s\n',optionname,optionvalue); + else + error(['ToolkitsFile error: option ' optionname ' is not well formatted']); + end + end + end + end - fclose(fid); - end %}}} - function savemodeljs(self,fid,modelname) % {{{ + fclose(fid); + end %}}} + function savemodeljs(self,fid,modelname) % {{{ - analyses=properties(self); - for i=1:numel(analyses), - if isempty(fieldnames(self.(analyses{i}))) - error(['md.toolkits.' analyses{i} ' is empty']); - else - writejsstruct(fid,[modelname '.toolkits.' analyses{i}],self.(analyses{i})); - end - end + analyses=properties(self); + for i=1:numel(analyses), + if isempty(fieldnames(self.(analyses{i}))) + error(['md.toolkits.' analyses{i} ' is empty']); + else + writejsstruct(fid,[modelname '.toolkits.' analyses{i}],self.(analyses{i})); + end + end - end % }}} - end + end % }}} + end end Index: ../trunk-jpl/src/m/classes/giamme.py =================================================================== --- ../trunk-jpl/src/m/classes/giamme.py (revision 25168) +++ ../trunk-jpl/src/m/classes/giamme.py (revision 25169) @@ -55,8 +55,8 @@ WriteData(fid, prefix, 'name', 'md.gia.model', 'data', 3, 'format', 'Integer') if isinstance(self.Ngia, list) or self.Ngia.ndim == 1: #list or 1D numpy.ndarray - WriteData(fid, prefix, 'data', np.zeros(md.mesh.numberofvertices), 'format', 'DoubleMat', 'mattype', 1, 'name', 'md.gia.Ngia') - WriteData(fid, prefix, 'data', np.zeros(md.mesh.numberofvertices), 'format', 'DoubleMat', 'mattype', 1, 'name', 'md.gia.Ugia') + WriteData(fid, prefix, 'data', np.zeros((md.mesh.numberofvertices, 1)), 'format', 'DoubleMat', 'mattype', 1, 'name', 'md.gia.Ngia') + WriteData(fid, prefix, 'data', np.zeros((md.mesh.numberofvertices, 1)), 'format', 'DoubleMat', 'mattype', 1, 'name', 'md.gia.Ugia') WriteData(fid, prefix, 'data', 1, 'name', 'md.gia.modelid', 'format', 'Double') WriteData(fid, prefix, 'name', 'md.gia.nummodels', 'data', 1, 'format', 'Integer') else: #2D numpy.ndarray Index: ../trunk-jpl/src/m/classes/dslmme.py =================================================================== --- ../trunk-jpl/src/m/classes/dslmme.py (nonexistent) +++ ../trunk-jpl/src/m/classes/dslmme.py (revision 25169) @@ -0,0 +1,80 @@ +import numpy as np + +from checkfield import checkfield +from fielddisplay import fielddisplay +from project3d import project3d +from WriteData import WriteData + + +class dslmme(object): + ''' + DSLMME class definition + + Usage: + dsl = dslmme() #dynamic sea level class based on a multi-model ensemble of CMIP5 outputs + ''' + + 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) + + if nargin == 0: + self.setdefaultparameters() + else: + raise Exception('constructor not supported') + #}}} + + def __repr__(self): # {{{ + s = ' gia 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, 'Ngia', 'geoid (mm/yr).')) + s += '{}\n'.format(fielddisplay(self, 'Ugia', 'uplift (mm/yr).')) + + return s + #}}} + + def setdefaultparameters(self): # {{{ + return + #}}} + + def checkconsistency(self, md, solution, analyses): # {{{ + if ('SealevelriseAnalysis' not in analyses) or (solution == 'TransientSolution' and md.transient.isslr == 0): + 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 md.slr.geodetic==0, + raise Exception('DSL checkconsistency error message: if bottom pressure fingerprints computations are requested, slr class should have geodetic flag on') + + return md + #}}} + + 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', 'Integer') + 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 / md.constants.yts) + 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) + #}}} + + 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) + + return self + #}}} Index: ../trunk-jpl/src/m/classes/geometry.py =================================================================== --- ../trunk-jpl/src/m/classes/geometry.py (revision 25168) +++ ../trunk-jpl/src/m/classes/geometry.py (revision 25169) @@ -1,41 +1,41 @@ import numpy as np + +from checkfield import checkfield +from fielddisplay import fielddisplay from project3d import project3d -from fielddisplay import fielddisplay -from checkfield import checkfield from WriteData import WriteData class geometry(object): - """ + ''' GEOMETRY class definition - Usage: - geometry = geometry() - """ + Usage: + geometry = geometry() + ''' - def __init__(self): # {{{ - self.surface = float('NaN') - self.thickness = float('NaN') - self.base = float('NaN') - self.bed = float('NaN') - self.hydrostatic_ratio = float('NaN') + def __init__(self): #{{{ + self.surface = np.nan + self.thickness = np.nan + self.base = np.nan + self.bed = np.nan + self.hydrostatic_ratio = np.nan - #set defaults + #set defaults self.setdefaultparameters() - #}}} - def __repr__(self): # {{{ - - string = " geometry parameters:" - string = "%s\n%s" % (string, fielddisplay(self, 'surface', 'ice upper surface elevation [m]')) - string = "%s\n%s" % (string, fielddisplay(self, 'thickness', 'ice thickness [m]')) - string = "%s\n%s" % (string, fielddisplay(self, 'base', 'ice base elevation [m]')) - string = "%s\n%s" % (string, fielddisplay(self, 'bed', 'bed elevation [m]')) - return string + def __repr__(self): #{{{ + s = " geometry parameters:" + s = "%s\n%s" % (s, fielddisplay(self, 'surface', 'ice upper surface elevation [m]')) + s = "%s\n%s" % (s, fielddisplay(self, 'thickness', 'ice thickness [m]')) + s = "%s\n%s" % (s, fielddisplay(self, 'base', 'ice base elevation [m]')) + s = "%s\n%s" % (s, fielddisplay(self, 'bed', 'bed elevation [m]')) + + return s #}}} - def extrude(self, md): # {{{ + def extrude(self, md): #{{{ self.surface = project3d(md, 'vector', self.surface, 'type', 'node') self.thickness = project3d(md, 'vector', self.thickness, 'type', 'node') self.hydrostatic_ratio = project3d(md, 'vector', self.hydrostatic_ratio, 'type', 'node') @@ -44,13 +44,15 @@ return self #}}} - def setdefaultparameters(self): # {{{ + def setdefaultparameters(self): #{{{ return self #}}} - def checkconsistency(self, md, solution, analyses): # {{{ + def checkconsistency(self, md, solution, analyses): #{{{ if (solution == 'TransientSolution' and md.transient.isgia) or (solution == 'GiaSolution'): - md = checkfield(md, 'fieldname', 'geometry.thickness', 'NaN', 1, 'Inf', 1, 'timeseries', 1) + 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 else: @@ -71,7 +73,7 @@ return md # }}} - def marshall(self, prefix, md, fid): # {{{ + def marshall(self, prefix, md, fid): #{{{ 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) Index: ../trunk-jpl/src/m/classes/geometry.m =================================================================== --- ../trunk-jpl/src/m/classes/geometry.m (revision 25168) +++ ../trunk-jpl/src/m/classes/geometry.m (revision 25169) @@ -60,7 +60,7 @@ elseif strcmpi(solution,'LoveSolution'), return; else - md = checkfield(md,'fieldname','geometry.surface' ,'NaN',1,'Inf',1,'size',[md.mesh.numberofvertices 1]); + md = checkfield(md,'fieldname','geometry.surface' ,'NaN',1,'Inf',1,'size',[md.mesh.numberofvertices 1]); md = checkfield(md,'fieldname','geometry.base' ,'NaN',1,'Inf',1,'size',[md.mesh.numberofvertices 1]); md = checkfield(md,'fieldname','geometry.thickness','NaN',1,'Inf',1,'size',[md.mesh.numberofvertices 1],'>',0); if any(abs(self.thickness-self.surface+self.base)>10^-9), Index: ../trunk-jpl/src/m/classes/matice.m =================================================================== --- ../trunk-jpl/src/m/classes/matice.m (revision 25168) +++ ../trunk-jpl/src/m/classes/matice.m (revision 25169) @@ -67,7 +67,7 @@ self.rho_freshwater=1000.; %water viscosity (N.s/m^2) - self.mu_water=0.001787; + self.mu_water=0.001787; %ice heat capacity cp (J/kg/K) self.heatcapacity=2093.; @@ -177,7 +177,6 @@ WriteData(fid,prefix,'object',self,'class','materials','fieldname','rheology_B','format','DoubleMat','mattype',1,'timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts); 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'); Index: ../trunk-jpl/src/m/classes/dsl.py =================================================================== --- ../trunk-jpl/src/m/classes/dsl.py (revision 25168) +++ ../trunk-jpl/src/m/classes/dsl.py (revision 25169) @@ -1,19 +1,20 @@ import numpy as np + +from checkfield import checkfield from fielddisplay import fielddisplay -from checkfield import checkfield +from project3d import project3d from WriteData import WriteData -from project3d import project3d class dsl(object): - """ - dsl Class definition + ''' + DSL - slass definition - Usage: - dsl = dsl() - """ + Usage: + dsl = dsl() + ''' - def __init__(self): # {{{ + 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) @@ -20,28 +21,29 @@ self.compute_fingerprints = 0; #do we use the sea water pressure change to compute fingerprints and correct sea_surface_height_change_above_geoid #}}} - def __repr__(self): # {{{ - string = " dsl parameters:" - string = "%s\n%s" % (string, fielddisplay(self, 'global_average_thermosteric_sea_level_change', 'corresponds to zostoga field in CMIP5 archives. Specified as a temporally variable global rate (mm/yr)')) - string = "%s\n%s" % (string, 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)')) - string = "%s\n%s" % (string, 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)')) - string = "%s\n%s" % (string, fielddisplay(self, 'compute_fingerprints', 'do we use the sea water pressure change to compute fingerprints and correct sea_surface_height_change_above_geoid')) - return string + def __repr__(self): #{{{ + s = " dsl parameters:" + s = "%s\n%s" % (s, 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 = "%s\n%s" % (s, 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 = "%s\n%s" % (s, 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 = "%s\n%s" % (s, fielddisplay(self, 'compute_fingerprints', 'do we use the sea water pressure change to compute fingerprints and correct sea_surface_height_change_above_geoid')) + + return s #}}} - def extrude(self, md): # {{{ + 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') return self #}}} - def defaultoutputs(self, md): # {{{ + def defaultoutputs(self, md): #{{{ return [] #}}} - def checkconsistency(self, md, solution, analyses): # {{{ + def checkconsistency(self, md, solution, analyses): #{{{ # Early retun - if not 'SealevelriseAnalysis' in analyses: + if 'SealevelriseAnalysis' not in analyses: return md if solution == 'TransientSolution' and md.transient.isslr == 0: @@ -59,7 +61,7 @@ return md # }}} - def marshall(self, prefix, md, fid): # {{{ + def marshall(self, prefix, md, fid): #{{{ yts = md.constants.yts WriteData(fid, prefix, 'name', 'md.dsl.model', 'data', 1, 'format', 'Integer') WriteData(fid, prefix, 'object', self, 'class', 'dsl', 'fieldname', 'compute_fingerprints', 'format', 'Integer') Index: ../trunk-jpl/src/m/classes/nodalvalue.py =================================================================== --- ../trunk-jpl/src/m/classes/nodalvalue.py (nonexistent) +++ ../trunk-jpl/src/m/classes/nodalvalue.py (revision 25169) @@ -0,0 +1,69 @@ +import numpy as np + +from checkfield import checkfield +from fielddisplay import fielddisplay +from WriteData import WriteData + + +class dslmme(object): + ''' + NODALVALUE class definition + + Usage: + nodalvalue=nodalvalue() + nodalvalue=nodalvalue( + 'name', 'SealevelriseSNodalValue', + 'definitionstring', 'Outputdefinition1', + 'model_string', 'SealevelriseS', + 'node', 1 + ) + ''' + + def __init__(self, *args): #{{{ + self.name = '' + self.definitionstring = '' #string that identifies this output definition uniquely, from 'Outputdefinition[1-10]' + self.model_string = '' #string for field that is being retrieved + self.node = np.nan #for which node are we retrieving the value? + + #use provided options to change fields + options = pairoptions(*args) + + #get name + self.name = options.getfieldvalue(options, 'name', '') + self.definitionstring = options.getfieldvalue(options, 'definitionstring', '') + self.model_string = options.getfieldvalue(options, 'model_string', '') + self.node = options.getfieldvalue(options, 'node', '') + #}}} + + def __repr__(self): # {{{ + s = ' Nodalvalue:\n' + s += '{}\n'.format(fielddisplay(self, 'name', 'identifier for this nodalvalue response')) + s += '{}\n'.format(fielddisplay(self, 'definitionstring', 'string that identifies this output definition uniquely, from \'Outputdefinition[1-10]\'')) + s += '{}\n'.format(fielddisplay(self, 'model_string', 'string for field that is being retrieved')) + s += '{}\n'.format(fielddisplay(self, 'node', 'vertex index at which we retrieve the value')) + + return s + #}}} + + def setdefaultparameters(self): # {{{ + return + #}}} + + def checkconsistency(self, md, solution, analyses): # {{{ + if not isinstance(self.name, basestring): + raise Exception('nodalvalue error message: \'name\' field should be a string!') + OutputdefinitionStringArray = [] + for i in range(100): + OutputdefinitionStringArray[i] = 'Outputdefinition{}'.format(i) + md = checkfield(md, 'fieldname', 'self.definitionstring', 'field', self.definitionstring, 'values', OutputdefinitionStringArray) + md = checkfield(md, 'fieldname', 'self.node', 'field', self.node, 'values', range(md.mesh.numberofvertices)) + + return md + #}}} + + def marshall(self, prefix, md, fid): #{{{ + WriteData(fid, prefix, 'data', self.name, 'name', 'md.nodalvalue.name', 'format', 'String') + WriteData(fid, prefix, 'data', self.definitionstring, 'name', 'md.nodalvalue.definitionenum', 'format', 'String') + WriteData(fid, prefix, 'data', self.model_string, 'name', 'md.nodalvalue.model_enum', 'format', 'String') + WriteData(fid, prefix, 'data', self.node, 'name', 'md.nodalvalue.node', 'format', 'Integer') + #}}} Index: ../trunk-jpl/src/m/classes/solidearth.py =================================================================== --- ../trunk-jpl/src/m/classes/solidearth.py (revision 25168) +++ ../trunk-jpl/src/m/classes/solidearth.py (revision 25169) @@ -23,7 +23,7 @@ self.sealevel = np.nan self.settings = solidearthsettings() self.surfaceload = surfaceload() - self.love = lovenumbers() + self.lovenumbers = lovenumbers() self.rotational = rotational() self.planetradius = planetradius('earth') self.requested_outputs = ['default'] @@ -64,7 +64,7 @@ self.settings.checkconsistency(md, solution, analyses) self.surfaceload.checkconsistency(md, solution, analyses) - self.love.checkconsistency(md, solution, analyses) + self.lovenumbers.checkconsistency(md, solution, analyses) self.rotational.checkconsistency(md, solution, analyses) return md @@ -81,7 +81,7 @@ self.settings.marshall(prefix, md, fid) self.surfaceload.marshall(prefix, md, fid) - self.love.marshall(prefix, md, fid) + self.lovenumbers.marshall(prefix, md, fid) self.rotational.marshall(prefix, md, fid) #process requested outputs Index: ../trunk-jpl/src/m/classes/toolkits.py =================================================================== --- ../trunk-jpl/src/m/classes/toolkits.py (revision 25168) +++ ../trunk-jpl/src/m/classes/toolkits.py (revision 25169) @@ -7,14 +7,14 @@ class toolkits(object): - """ + ''' TOOLKITS class definition - Usage: - self = toolkits() - """ + Usage: + self = toolkits() + ''' - def __init__(self): # {{{ + def __init__(self): #{{{ #default toolkits if IssmConfig('_HAVE_PETSC_')[0]: #MUMPS is the default toolkits @@ -33,18 +33,18 @@ #Use same solver for Recovery mode self.RecoveryAnalysis = self.DefaultAnalysis - #The other properties are dynamic - # }}} + #The other properties are dynamic + #}}} - def __repr__(self): # {{{ + def __repr__(self): #{{{ s = "List of toolkits options per analysis:\n\n" for analysis in list(vars(self).keys()): - s += "%s\n" % fielddisplay(self, analysis, '') + s += "{}\n".format(fielddisplay(self, analysis, '')) - return s - # }}} + return s + #}}} - def addoptions(self, analysis, *args): # {{{ + def addoptions(self, analysis, *args): #{{{ # Usage example: # md.toolkits = addoptions(md.toolkits, 'StressbalanceAnalysis', FSoptions()) # md.toolkits = addoptions(md.toolkits, 'StressbalanceAnalysis') @@ -58,43 +58,49 @@ setattr(self, analysis, args[0]) return self - # }}} + #}}} - def checkconsistency(self, md, solution, analyses): # {{{ + def checkconsistency(self, md, solution, analyses): #{{{ + # TODO + # - Implement something closer to a switch as in + # src/m/classes/toolkits.m? + # for analysis in list(vars(self).keys()): if not getattr(self, analysis): - md.checkmessage("md.toolkits.%s is empty" % analysis) + md.checkmessage("md.toolkits.{} is empty".format(analysis)) return md - # }}} + #}}} - def ToolkitsFile(self, filename): # {{{ - """ + def ToolkitsFile(self, filename): #{{{ + ''' TOOLKITSFILE - build toolkits file - Build a Petsc compatible options file, from the toolkits model field + return options string - This file will also be used when the toolkit used is 'issm' instead of 'petsc' + Build a Petsc compatible options file, from the toolkits model + field + return options string. + This file will also be used when the toolkit used is 'issm' instead + of 'petsc'.s + Usage: + ToolkitsFile(toolkits, filename) + ''' - Usage: ToolkitsFile(toolkits, filename) - """ - - #open file for writing + #open file for writing try: fid = open(filename, 'w') except IOError as e: raise IOError("ToolkitsFile error: could not open {}' for writing due to".format(filename), e) - #write header + #write header fid.write("%s%s%s\n" % ('%Toolkits options file: ', filename, ' written from Python toolkits array')) - #start writing options + #start writing options for analysis in list(vars(self).keys()): options = getattr(self, analysis) - #first write analysis: + #first write analysis: fid.write("\n+{}\n".format(analysis)) #append a + to recognize it's an analysis enum - #now, write options + #now, write options for optionname, optionvalue in list(options.items()): if not optionvalue: @@ -110,4 +116,4 @@ raise TypeError("ToolkitsFile error: option '{}' is not well formatted.".format(optionname)) fid.close() - # }}} + #}}} Index: ../trunk-jpl/src/m/classes/matdamageice.m =================================================================== --- ../trunk-jpl/src/m/classes/matdamageice.m (revision 25168) +++ ../trunk-jpl/src/m/classes/matdamageice.m (revision 25169) @@ -119,7 +119,7 @@ md = checkfield(md,'fieldname','materials.rheology_n','>',0,'size',[md.mesh.numberofelements 1]); 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 ismember('GiaAnalysis',analyses), md = checkfield(md,'fieldname','materials.lithosphere_shear_modulus','>',0,'numel',1); md = checkfield(md,'fieldname','materials.lithosphere_density','>',0,'numel',1); Index: ../trunk-jpl/src/m/classes/dslmme.m =================================================================== --- ../trunk-jpl/src/m/classes/dslmme.m (revision 25168) +++ ../trunk-jpl/src/m/classes/dslmme.m (revision 25169) @@ -14,12 +14,6 @@ end methods - function self = extrude(self,md) % {{{ - for i=1:length(self.global_average_thermosteric_sea_level_change), - self.sea_surface_height_change_above_geoid{i}=project3d(md,'vector',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); - end - end % }}} function self = dslmme(varargin) % {{{ switch nargin case 0 @@ -79,6 +73,12 @@ 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); end % }}} + function self = extrude(self,md) % {{{ + for i=1:length(self.global_average_thermosteric_sea_level_change), + self.sea_surface_height_change_above_geoid{i}=project3d(md,'vector',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); + end + end % }}} function savemodeljs(self,fid,modelname) % {{{ writejsdouble(fid,[modelname '.dsl.modelid'],self.modelid); Index: ../trunk-jpl/src/m/solve/solve.py =================================================================== --- ../trunk-jpl/src/m/solve/solve.py (revision 25168) +++ ../trunk-jpl/src/m/solve/solve.py (revision 25169) @@ -1,13 +1,13 @@ import datetime import os import shutil -from pairoptions import pairoptions + from ismodelselfconsistent import ismodelselfconsistent +from loadresultsfromcluster import loadresultsfromcluster from marshall import marshall +from pairoptions import pairoptions +from preqmu import * from waitonlock import waitonlock -from loadresultsfromcluster import loadresultsfromcluster -from preqmu import * -#from MatlabFuncs import * def solve(md, solutionstring, *args): @@ -14,35 +14,37 @@ ''' SOLVE - apply solution sequence for this model - Usage: - md = solve(md, solutionstring, varargin) - where varargin is a list of paired arguments of string OR enums + Usage: + md = solve(md, solutionstring, varargin) + where varargin is a list of paired arguments of string OR enums solution types available comprise: - - 'Stressbalance' or 'sb' - - 'Masstransport' or 'mt' - - 'Thermal' or 'th' - - 'Steadystate' or 'ss' - - 'Transient' or 'tr' - - 'Balancethickness' or 'mc' - - 'Balancevelocity' or 'bv' - - 'BedSlope' or 'bsl' - - 'SurfaceSlope' or 'ssl' - - 'Hydrology' or 'hy' - - 'DamageEvolution' or 'da' - - 'Gia' or 'gia' - - 'Esa' or 'esa' - - 'Sealevelrise' or 'slr' - - 'Love' or 'lv' + - 'Stressbalance' or 'sb' + - 'Masstransport' or 'mt' + - 'Thermal' or 'th' + - 'Steadystate' or 'ss' + - 'Transient' or 'tr' + - 'Balancethickness' or 'mc' + - 'Balancevelocity' or 'bv' + - 'BedSlope' or 'bsl' + - 'SurfaceSlope' or 'ssl' + - 'Hydrology' or 'hy' + - 'DamageEvolution' or 'da' + - 'Gia' or 'gia' + - 'Esa' or 'esa' + - 'Sealevelrise' or 'slr' + - 'Love' or 'lv' - extra options: - - loadonly : does not solve. only load results - - checkconsistency : 'yes' or 'no' (default is 'yes'), ensures checks on consistency of model - - restart: 'directory name (relative to the execution directory) where the restart file is located. + extra options: + - loadonly : does not solve. only load results + - checkconsistency : 'yes' or 'no' (default is 'yes'), ensures + checks on consistency of model + - restart : 'directory name (relative to the execution + directory) where the restart file is located - Examples: - md = solve(md, 'Stressbalance') - md = solve(md, 'sb') + Examples: + md = solve(md, 'Stressbalance') + md = solve(md, 'sb') ''' #recover and process solve options Index: ../trunk-jpl/src/m/consistency/checkfield.py =================================================================== --- ../trunk-jpl/src/m/consistency/checkfield.py (revision 25168) +++ ../trunk-jpl/src/m/consistency/checkfield.py (revision 25169) @@ -7,32 +7,35 @@ def checkfield(md, *args): - """ + ''' CHECKFIELD - check field consistency - Used to check model consistency., - Requires: - 'field' or 'fieldname' option. If 'fieldname' is provided, it will retrieve it from the model md. (md.(fieldname)) - If 'field' is provided, it will assume the argument following 'field' is a numeric array. + Used to check model consistency + + Requires: + 'field' or 'fieldname' option. If 'fieldname' is provided, it will retrieve it from the model md. (md.(fieldname)) + If 'field' is provided, it will assume the argument following + 'field' is a numeric array. - Available options: - - NaN: 1 if check that there is no NaN - - size: [lines cols], NaN for non checked dimensions, or 'universal' for any input type (nodal, element, time series, etc) - -> : greater than provided value - ->= : greater or equal to provided value - - < : smallerthan provided value - - <=: smaller or equal to provided value - - < vec: smallerthan provided values on each vertex - - timeseries: 1 if check time series consistency (size and time) - - values: cell of strings or vector of acceptable values - - numel: list of acceptable number of elements - - cell: 1 if check that is cell - - empty: 1 if check that non empty - - message: overloaded error message + Available options: + - NaN: 1 if check that there is no NaN + - size: [lines cols], NaN for non checked dimensions, or 'universal' + for any input type (nodal, element, time series, etc) + - > : greater than provided value + - >= : greater or equal to provided value + - < : smallerthan provided value + - <=: smaller or equal to provided value + - < vec: smallerthan provided values on each vertex + - timeseries: 1 if check time series consistency (size and time) + - values: list of acceptable values + - numel: list of acceptable number of elements + - cell: 1 if check that is cell + - empty: 1 if check that non empty + - message: overloaded error message - Usage: - md = checkfield(md, fieldname, options) - """ + Usage: + md = checkfield(md, fieldname, options) + ''' #get options options = pairoptions(*args) @@ -142,7 +145,7 @@ #check cell if options.getfieldvalue('cell', 0): if not isinstance(field, (tuple, list, dict)): - md = md.checkmessage(options.getfieldvalue('message', "field '{}' should be a cell".format(fieldname))) + md = md.checkmessage(options.getfieldvalue('message', "field '{}' should be a tuple, list, or dict".format(fieldname))) #check values if options.exist('values'): Index: ../trunk-jpl/test/NightlyRun/test2002.py =================================================================== --- ../trunk-jpl/test/NightlyRun/test2002.py (revision 25168) +++ ../trunk-jpl/test/NightlyRun/test2002.py (revision 25169) @@ -38,7 +38,7 @@ #mask: {{{ mask = gmtmask(md.mesh.lat, md.mesh.long) -icemask = np.ones(md.mesh.numberofvertices) +icemask = np.ones((md.mesh.numberofvertices, 1)) pos = np.where(mask == 0)[0] icemask[pos] = -1 pos = np.where(np.sum(mask[md.mesh.elements - 1], axis=1) < 3) @@ -47,8 +47,8 @@ md.mask.ocean_levelset = -icemask #make sure that the ice level set is all inclusive: -md.mask.land_levelset = np.zeros((md.mesh.numberofvertices)) -md.mask.ocean_levelset = -np.ones((md.mesh.numberofvertices)) +md.mask.land_levelset = np.zeros((md.mesh.numberofvertices, 1)) +md.mask.ocean_levelset = -np.ones((md.mesh.numberofvertices, 1)) #make sure that the elements that have loads are fully grounded pos = np.nonzero(md.solidearth.surfaceload.icethicknesschange)[0] Index: ../trunk-jpl/test/NightlyRun/test2002.m =================================================================== --- ../trunk-jpl/test/NightlyRun/test2002.m (revision 25168) +++ ../trunk-jpl/test/NightlyRun/test2002.m (revision 25169) @@ -64,7 +64,9 @@ md.solidearth.settings.maxiter=10; %eustatic run: -md.solidearth.settings.rigid=0; md.solidearth.settings.elastic=0;md.solidearth.settings.rotation=0; +md.solidearth.settings.rigid=0; +md.solidearth.settings.elastic=0; +md.solidearth.settings.rotation=0; md=solve(md,'Sealevelrise'); Seustatic=md.results.SealevelriseSolution.Sealevel; Index: ../trunk-jpl/src/m/classes/matice.py =================================================================== --- ../trunk-jpl/src/m/classes/matice.py (revision 25168) +++ ../trunk-jpl/src/m/classes/matice.py (revision 25169) @@ -1,3 +1,5 @@ +import numpy as np + from fielddisplay import fielddisplay from project3d import project3d from checkfield import checkfield @@ -5,14 +7,14 @@ class matice(object): - """ + ''' MATICE class definition - Usage: - matice = matice() - """ + Usage: + matice = matice() + ''' - def __init__(self): # {{{ + def __init__(self): #{{{ self.rho_ice = 0. self.rho_water = 0. self.rho_freshwater = 0. @@ -26,80 +28,80 @@ self.beta = 0. self.mixed_layer_capacity = 0. self.thermal_exchange_velocity = 0. - self.rheology_B = float('NaN') - self.rheology_n = float('NaN') + self.rheology_B = np.nan + self.rheology_n = np.nan self.rheology_law = '' - #giaivins: + #giaivins self.lithosphere_shear_modulus = 0. self.lithosphere_density = 0. self.mantle_shear_modulus = 0. self.mantle_density = 0. - #SLR - self.earth_density = 5512 + #slr + self.earth_density = 0 self.setdefaultparameters() #}}} - def __repr__(self): # {{{ - string = " Materials:" + def __repr__(self): #{{{ + s = " Materials:" + s = "%s\n%s" % (s, fielddisplay(self, "rho_ice", "ice density [kg/m^3]")) + s = "%s\n%s" % (s, fielddisplay(self, "rho_water", "water density [kg/m^3]")) + s = "%s\n%s" % (s, fielddisplay(self, "rho_freshwater", "fresh water density [kg/m^3]")) + s = "%s\n%s" % (s, fielddisplay(self, "mu_water", "water viscosity [Ns/m^2]")) + s = "%s\n%s" % (s, fielddisplay(self, "heatcapacity", "heat capacity [J/kg/K]")) + s = "%s\n%s" % (s, fielddisplay(self, "thermalconductivity", "ice thermal conductivity [W/m/K]")) + s = "%s\n%s" % (s, fielddisplay(self, "temperateiceconductivity", "temperate ice thermal conductivity [W/m/K]")) + s = "%s\n%s" % (s, fielddisplay(self, "effectiveconductivity_averaging", "computation of effectiveconductivity: (0) arithmetic mean, (1) harmonic mean, (2) geometric mean (default)")) + s = "%s\n%s" % (s, fielddisplay(self, "meltingpoint", "melting point of ice at 1atm in K")) + s = "%s\n%s" % (s, fielddisplay(self, "latentheat", "latent heat of fusion [J/m^3]")) + s = "%s\n%s" % (s, fielddisplay(self, "beta", "rate of change of melting point with pressure [K/Pa]")) + s = "%s\n%s" % (s, fielddisplay(self, "mixed_layer_capacity", "mixed layer capacity [W/kg/K]")) + s = "%s\n%s" % (s, fielddisplay(self, "thermal_exchange_velocity", "thermal exchange velocity [m/s]")) + s = "%s\n%s" % (s, fielddisplay(self, "rheology_B", "flow law parameter [Pa s^(1/n)]")) + 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]")) - string = "%s\n%s" % (string, fielddisplay(self, "rho_ice", "ice density [kg / m^3]")) - string = "%s\n%s" % (string, fielddisplay(self, "rho_water", "water density [kg / m^3]")) - string = "%s\n%s" % (string, fielddisplay(self, "rho_freshwater", "fresh water density [kg / m^3]")) - string = "%s\n%s" % (string, fielddisplay(self, "mu_water", "water viscosity [N s / m^2]")) - string = "%s\n%s" % (string, fielddisplay(self, "heatcapacity", "heat capacity [J / kg / K]")) - string = "%s\n%s" % (string, fielddisplay(self, "thermalconductivity", "ice thermal conductivity [W / m / K]")) - string = "%s\n%s" % (string, fielddisplay(self, "temperateiceconductivity", "temperate ice thermal conductivity [W / m / K]")) - string = "%s\n%s" % (string, fielddisplay(self, "effectiveconductivity_averaging", "computation of effectiveconductivity: (0) arithmetic mean, (1) harmonic mean, (2) geometric mean (default)")) - string = "%s\n%s" % (string, fielddisplay(self, "meltingpoint", "melting point of ice at 1atm in K")) - string = "%s\n%s" % (string, fielddisplay(self, "latentheat", "latent heat of fusion [J / m^3]")) - string = "%s\n%s" % (string, fielddisplay(self, "beta", "rate of change of melting point with pressure [K / Pa]")) - string = "%s\n%s" % (string, fielddisplay(self, "mixed_layer_capacity", "mixed layer capacity [W / kg / K]")) - string = "%s\n%s" % (string, fielddisplay(self, "thermal_exchange_velocity", "thermal exchange velocity [m / s]")) - string = "%s\n%s" % (string, fielddisplay(self, "rheology_B", "flow law parameter [Pa s^(1 / n)]")) - 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', 'LliboutryDuval', 'NyeCO2', or 'NyeH2O'")) - 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 + return s #}}} - def extrude(self, md): # {{{ + def extrude(self, md): #{{{ self.rheology_B = project3d(md, 'vector', self.rheology_B, 'type', 'node') self.rheology_n = project3d(md, 'vector', self.rheology_n, 'type', 'element') return self #}}} - def setdefaultparameters(self): # {{{ - #ice density (kg / m^3) + def setdefaultparameters(self): #{{{ + #ice density (kg/m^3) self.rho_ice = 917. - #ocean water density (kg / m^3) + #ocean water density (kg/m^3) self.rho_water = 1023. - #fresh water density (kg / m^3) + #fresh water density (kg/m^3) self.rho_freshwater = 1000. - #water viscosity (N.s / m^2) + #water viscosity (N.s/m^2) self.mu_water = 0.001787 - #ice heat capacity cp (J / kg / K) + #ice heat capacity cp (J/kg/K) self.heatcapacity = 2093. - #ice latent heat of fusion L (J / kg) + #ice latent heat of fusion L (J/kg) self.latentheat = 3.34e5 - #ice thermal conductivity (W / m / K) + #ice thermal conductivity (W/m/K) self.thermalconductivity = 2.4 + #temperate ice thermal conductivity (W/m/K) + self.temperateiceconductivity = 0.24 #computation of effective conductivity self.effectiveconductivity_averaging = 1 - #temperate ice thermal conductivity (W / m / K) - self.temperateiceconductivity = 0.24 #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) + #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) + #mixed layer (ice-water interface) heat capacity (J/kg/K) self.mixed_layer_capacity = 3974. - #thermal exchange velocity (ice-water interface) (m / s) + #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 @@ -106,34 +108,39 @@ self.rheology_law = 'Paterson' # GIA: - self.lithosphere_shear_modulus = 6.7e10 # (Pa) - self.lithosphere_density = 3.32 # (g / cm^ - 3) + 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) + self.mantle_density = 3.34 # (g/cm^-3) - #SLR - self.earth_density = 5512 # average density of the Earth, (kg / m^3) - return self + # SLR + self.earth_density = 5512 # average density of the Earth, (kg/m^3) #}}} - def checkconsistency(self, md, solution, analyses): # {{{ - 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, 'timeseries', 1, 'NaN', 1, 'Inf', 1) - md = checkfield(md, 'fieldname', 'materials.rheology_n', '>', 0, 'size', [md.mesh.numberofelements]) - 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]) - 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]) - md = checkfield(md, 'fieldname', 'materials.earth_density', '>', 0, 'numel', [1]) + def checkconsistency(self, md, solution, analyses): #{{{ + if solution != 'SealevelriseSolution': + 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, 'timeseries', 1, 'NaN', 1, 'Inf', 1) + md = checkfield(md, 'fieldname', 'materials.rheology_n', '>', 0, 'size', [md.mesh.numberofelements]) + 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]) + return md - # }}} + #}}} - def marshall(self, prefix, md, fid): # {{{ + def marshall(self, prefix, md, fid): #{{{ WriteData(fid, prefix, 'name', 'md.materials.type', 'data', 3, 'format', 'Integer') WriteData(fid, prefix, 'object', self, 'class', 'materials', 'fieldname', 'rho_ice', 'format', 'Double') WriteData(fid, prefix, 'object', self, 'class', 'materials', 'fieldname', 'rho_water', 'format', 'Double') @@ -150,12 +157,10 @@ WriteData(fid, prefix, 'object', self, 'class', 'materials', 'fieldname', 'thermal_exchange_velocity', 'format', 'Double') WriteData(fid, prefix, 'object', self, 'class', 'materials', 'fieldname', 'rheology_B', 'format', 'DoubleMat', 'mattype', 1, 'timeserieslength', md.mesh.numberofvertices + 1, 'yts', md.constants.yts) 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, 'data', self.rheology_law, 'name', 'md.materials.rheology_law', 'format', 's') 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: ../trunk-jpl/src/m/classes/slr.py =================================================================== --- ../trunk-jpl/src/m/classes/slr.py (revision 25168) +++ ../trunk-jpl/src/m/classes/slr.py (revision 25169) @@ -49,34 +49,34 @@ #}}} def __repr__(self): # {{{ - string = ' slr parameters:' - string = "%s\n%s" % (string, fielddisplay(self, 'deltathickness', 'thickness change: ice height equivalent [m]')) - string = "%s\n%s" % (string, fielddisplay(self, 'sealevel', 'current sea level (prior to computation) [m]')) - string = "%s\n%s" % (string, fielddisplay(self, 'spcthickness', 'thickness constraints (NaN means no constraint) [m]')) - string = "%s\n%s" % (string, fielddisplay(self, 'reltol', 'sea level rise relative convergence criterion, (NaN: not applied)')) - string = "%s\n%s" % (string, fielddisplay(self, 'abstol', 'sea level rise absolute convergence criterion, (default, NaN: not applied)')) - string = "%s\n%s" % (string, fielddisplay(self, 'maxiter', 'maximum number of nonlinear iterations')) - string = "%s\n%s" % (string, fielddisplay(self, 'love_h', 'load Love number for radial displacement')) - string = "%s\n%s" % (string, fielddisplay(self, 'love_k', 'load Love number for gravitational potential perturbation')) - string = "%s\n%s" % (string, fielddisplay(self, 'love_l', 'load Love number for horizontal displaements')) - string = "%s\n%s" % (string, fielddisplay(self, 'tide_love_k', 'tidal load Love number (degree 2)')) - string = "%s\n%s" % (string, fielddisplay(self, 'tide_love_h', 'tidal load Love number (degree 2)')) - string = "%s\n%s" % (string, fielddisplay(self, 'fluid_love', 'secular fluid Love number')) - string = "%s\n%s" % (string, fielddisplay(self, 'equatorial_moi', 'mean equatorial moment of inertia [kg m^2]')) - string = "%s\n%s" % (string, fielddisplay(self, 'polar_moi', 'polar moment of inertia [kg m^2]')) - string = "%s\n%s" % (string, fielddisplay(self, 'angular_velocity', 'mean rotational velocity of earth [per second]')) - string = "%s\n%s" % (string, fielddisplay(self, 'ocean_area_scaling', 'correction for model representation of ocean area [default: No correction]')) - string = "%s\n%s" % (string, fielddisplay(self, 'hydro_rate', 'rate of hydrological expansion [mm / yr]')) - string = "%s\n%s" % (string, fielddisplay(self, 'geodetic', 'compute geodetic SLR? (in addition to steric?) default 0')) - string = "%s\n%s" % (string, fielddisplay(self, 'geodetic_run_frequency', 'how many time steps we skip before we run SLR solver during transient (default: 1)')) - string = "%s\n%s" % (string, fielddisplay(self, 'rigid', 'rigid earth graviational potential perturbation')) - string = "%s\n%s" % (string, fielddisplay(self, 'elastic', 'elastic earth graviational potential perturbation')) - string = "%s\n%s" % (string, fielddisplay(self, 'rotation', 'earth rotational potential perturbation')) - string = "%s\n%s" % (string, fielddisplay(self, 'degacc', 'accuracy (default .01 deg) for numerical discretization of the Green''s functions')) - string = "%s\n%s" % (string, fielddisplay(self, 'transitions', 'indices into parts of the mesh that will be icecaps')) - string = "%s\n%s" % (string, fielddisplay(self, 'requested_outputs', 'additional outputs requested')) + s = ' slr parameters:' + s = "%s\n%s" % (s, fielddisplay(self, 'deltathickness', 'thickness change: ice height equivalent [m]')) + s = "%s\n%s" % (s, fielddisplay(self, 'sealevel', 'current sea level (prior to computation) [m]')) + s = "%s\n%s" % (s, fielddisplay(self, 'spcthickness', 'thickness constraints (NaN means no constraint) [m]')) + s = "%s\n%s" % (s, fielddisplay(self, 'reltol', 'sea level rise relative convergence criterion, (NaN: not applied)')) + s = "%s\n%s" % (s, fielddisplay(self, 'abstol', 'sea level rise absolute convergence criterion, (default, NaN: not applied)')) + s = "%s\n%s" % (s, fielddisplay(self, 'maxiter', 'maximum number of nonlinear iterations')) + s = "%s\n%s" % (s, fielddisplay(self, 'love_h', 'load Love number for radial displacement')) + s = "%s\n%s" % (s, fielddisplay(self, 'love_k', 'load Love number for gravitational potential perturbation')) + s = "%s\n%s" % (s, fielddisplay(self, 'love_l', 'load Love number for horizontal displaements')) + s = "%s\n%s" % (s, fielddisplay(self, 'tide_love_k', 'tidal load Love number (degree 2)')) + s = "%s\n%s" % (s, fielddisplay(self, 'tide_love_h', 'tidal load Love number (degree 2)')) + s = "%s\n%s" % (s, fielddisplay(self, 'fluid_love', 'secular fluid Love number')) + s = "%s\n%s" % (s, fielddisplay(self, 'equatorial_moi', 'mean equatorial moment of inertia [kg m^2]')) + s = "%s\n%s" % (s, fielddisplay(self, 'polar_moi', 'polar moment of inertia [kg m^2]')) + s = "%s\n%s" % (s, fielddisplay(self, 'angular_velocity', 'mean rotational velocity of earth [per second]')) + s = "%s\n%s" % (s, fielddisplay(self, 'ocean_area_scaling', 'correction for model representation of ocean area [default: No correction]')) + s = "%s\n%s" % (s, fielddisplay(self, 'hydro_rate', 'rate of hydrological expansion [mm / yr]')) + s = "%s\n%s" % (s, fielddisplay(self, 'geodetic', 'compute geodetic SLR? (in addition to steric?) default 0')) + s = "%s\n%s" % (s, fielddisplay(self, 'geodetic_run_frequency', 'how many time steps we skip before we run SLR solver during transient (default: 1)')) + s = "%s\n%s" % (s, fielddisplay(self, 'rigid', 'rigid earth graviational potential perturbation')) + s = "%s\n%s" % (s, fielddisplay(self, 'elastic', 'elastic earth graviational potential perturbation')) + s = "%s\n%s" % (s, fielddisplay(self, 'rotation', 'earth rotational potential perturbation')) + s = "%s\n%s" % (s, fielddisplay(self, 'degacc', 'accuracy (default .01 deg) for numerical discretization of the Green''s functions')) + s = "%s\n%s" % (s, fielddisplay(self, 'transitions', 'indices into parts of the mesh that will be icecaps')) + s = "%s\n%s" % (s, fielddisplay(self, 'requested_outputs', 'additional outputs requested')) - return string + return s # }}} def setdefaultparameters(self): # {{{ @@ -142,7 +142,7 @@ md = checkfield(md, 'fieldname', 'slr.geodetic_run_frequency', 'size', [1, 1], '>=', 1) md = checkfield(md, 'fieldname', 'slr.hydro_rate', 'NaN', 1, 'Inf', 1, 'size', [md.mesh.numberofvertices]) md = checkfield(md, 'fieldname', 'slr.degacc', 'size', [1, 1], '>=', 1e-10) - md = checkfield(md, 'fieldname', 'slr.requested_outputs', 'stringrow', 1) + md = checkfield(md, 'fieldname', 'slr.requested_outputs', 'srow', 1) md = checkfield(md, 'fieldname', 'slr.horiz', 'NaN', 1, 'Inf', 1, 'values', [0, 1]) #check that love numbers are provided at the same level of accuracy: @@ -201,5 +201,5 @@ if len(indices) > 0: outputscopy = outputs[0:max(0, indices[0] - 1)] + self.defaultoutputs(md) + outputs[indices[0] + 1:] outputs = outputscopy - WriteData(fid, prefix, 'data', outputs, 'name', 'md.slr.requested_outputs', 'format', 'StringArray') + WriteData(fid, prefix, 'data', outputs, 'name', 'md.slr.requested_outputs', 'format', 'sArray') # }}}