Index: /issm/trunk-jpl/src/m/archive/arch.py =================================================================== --- /issm/trunk-jpl/src/m/archive/arch.py (revision 27416) +++ /issm/trunk-jpl/src/m/archive/arch.py (revision 27417) @@ -12,6 +12,5 @@ """ - nargs = len(args) - if nargs % 2 != 0: + if len(args) % 2 != 0: raise ValueError('Incorrect number of arguments.') # open file Index: /issm/trunk-jpl/src/m/classes/SMBarma.m =================================================================== --- /issm/trunk-jpl/src/m/classes/SMBarma.m (revision 27416) +++ /issm/trunk-jpl/src/m/classes/SMBarma.m (revision 27417) @@ -100,6 +100,6 @@ end if (any(isnan(md.smb.refelevation)==0) || numel(md.smb.refelevation)>1) - md = checkfield(md,'fieldname','smb.refelevation','NaN',1,'Inf',1,'>=',0,'size',[1,md.smb.num_basins],'numel',md.smb.num_basins); - end + md = checkfield(md,'fieldname','smb.refelevation','NaN',1,'Inf',1,'>=',0,'size',[1,md.smb.num_basins],'numel',md.smb.num_basins); + end [nbas,nbins] = size(md.smb.lapserates); if (any(isnan(reshape(md.smb.lapserates,[1,nbas*nbins]))==0) || numel(md.smb.lapserates)>1) @@ -126,7 +126,7 @@ fielddisplay(self,'basin_id','basin number assigned to each element [unitless]'); fielddisplay(self,'num_breaks','number of different breakpoints in the piecewise-polynomial (separating num_breaks+1 periods)'); - fielddisplay(self,'num_params','number of different parameters in the piecewise-polynomial (1:intercept only, 2:with linear trend, 3:with quadratic trend, etc.)'); - fielddisplay(self,'polynomialparams','coefficients for the polynomial (const,trend,quadratic,etc.),dim1 for basins,dim2 for periods,dim3 for orders'); - disp(sprintf('%51s ex: polyparams=cat(3,intercepts,trendlinearcoefs,trendquadraticcoefs)',' ')); + fielddisplay(self,'num_params','number of different parameters in the piecewise-polynomial (1:intercept only, 2:with linear trend, 3:with quadratic trend, etc.)'); + fielddisplay(self,'polynomialparams','coefficients for the polynomial (const,trend,quadratic,etc.),dim1 for basins,dim2 for periods,dim3 for orders'); + disp(sprintf('%51s ex: polyparams=cat(3,intercepts,trendlinearcoefs,trendquadraticcoefs)',' ')); fielddisplay(self,'datebreaks','dates at which the breakpoints in the piecewise polynomial occur (1 row per basin) [yr]'); fielddisplay(self,'ar_order','order of the autoregressive model [unitless]'); @@ -150,6 +150,6 @@ yts=md.constants.yts; nbas = md.smb.num_basins; - nprm = md.smb.num_params; - nper = md.smb.num_breaks+1; + nprm = md.smb.num_params; + nper = md.smb.num_breaks+1; templapserates = md.smb.lapserates; @@ -180,6 +180,6 @@ % Scale the parameters % - polyparamsScaled = md.smb.polynomialparams; - polyparams2dScaled = zeros(nbas,nper*nprm); + polyparamsScaled = md.smb.polynomialparams; + polyparams2dScaled = zeros(nbas,nper*nprm); if(nprm>1) % Case 3D % Index: /issm/trunk-jpl/src/m/classes/SMBarma.py =================================================================== --- /issm/trunk-jpl/src/m/classes/SMBarma.py (revision 27416) +++ /issm/trunk-jpl/src/m/classes/SMBarma.py (revision 27417) @@ -30,6 +30,5 @@ self.requested_outputs = [] - nargs = len(args) - if nargs == 0: + if len(args) == 0: self.setdefaultparameters() else: Index: /issm/trunk-jpl/src/m/classes/SMBcomponents.py =================================================================== --- /issm/trunk-jpl/src/m/classes/SMBcomponents.py (revision 27416) +++ /issm/trunk-jpl/src/m/classes/SMBcomponents.py (revision 27417) @@ -22,6 +22,5 @@ self.requested_outputs = [] - nargs = len(args) - if nargs == 0: + if len(args) == 0: self.setdefaultparameters() else: Index: /issm/trunk-jpl/src/m/classes/SMBforcing.py =================================================================== --- /issm/trunk-jpl/src/m/classes/SMBforcing.py (revision 27416) +++ /issm/trunk-jpl/src/m/classes/SMBforcing.py (revision 27417) @@ -20,6 +20,5 @@ self.averaging = 0 - nargs = len(args) - if nargs == 0: + if len(args) == 0: self.setdefaultparameters() else: Index: /issm/trunk-jpl/src/m/classes/SMBpddSicopolis.py =================================================================== --- /issm/trunk-jpl/src/m/classes/SMBpddSicopolis.py (revision 27416) +++ /issm/trunk-jpl/src/m/classes/SMBpddSicopolis.py (revision 27417) @@ -31,6 +31,5 @@ self.requested_outputs = [] - nargs = len(args) - if nargs == 0: + if len(args) == 0: self.setdefaultparameters() else: Index: /issm/trunk-jpl/src/m/classes/calvingparameterization.py =================================================================== --- /issm/trunk-jpl/src/m/classes/calvingparameterization.py (revision 27417) +++ /issm/trunk-jpl/src/m/classes/calvingparameterization.py (revision 27417) @@ -0,0 +1,69 @@ +from checkfield import checkfield +from fielddisplay import fielddisplay +from WriteData import WriteData + + +class calvingparameterization(object): + """CALVINGPARAMETERIZATION class definition + For test calving laws and coefficients + + Usage: + calvingparameterization = calvingparameterization() + """ + + def __init__(self): # {{{ + self.stress_threshold_groundedice = 0 + self.stress_threshold_floatingice = 0 + self.min_thickness = 0 + self.use_param = 0 + self.scale_theta = 0 + self.amp_alpha = 0 + self.midp = 0 + self.nonlinearlaw = 0 + + #set defaults + self.setdefaultparameters() + + #}}} + + def __repr__(self): # {{{ + s = ' Calving VonMises parameters:\n' + s += '{}\n'.format(fielddisplay(self, 'stress_threshold_groundedice', 'sigma_max applied to grounded ice only [Pa]')) + s += '{}\n'.format(fielddisplay(self, 'stress_threshold_floatingice', 'sigma_max applied to floating ice only [Pa]')) + s += '{}\n'.format(fielddisplay(self, 'min_thickness', 'minimum thickness below which no ice is allowed [m]')) + return s + #}}} + + def extrude(self, md): # {{{ + return self + #}}} + + def setdefaultparameters(self): # {{{ + # Default sigma max + self.stress_threshold_groundedice = 1e6 + self.stress_threshold_floatingice = 150e3 + + # Turn off min_thickness by default + self.min_thickness = 0. + return self + #}}} + + def checkconsistency(self, md, solution, analyses): # {{{ + # Early return + if solution == 'TransientSolution' or not md.transient.ismovingfront: + return + + md = checkfield(md, 'fieldname', 'calving.stress_threshold_groundedice', '>', 0, 'nan', 1, 'Inf', 1, 'size', 'universal') + md = checkfield(md, 'fieldname', 'calving.stress_threshold_floatingice', '>', 0, 'nan', 1, 'Inf', 1, 'size', 'universal') + md = checkfield(md, 'fieldname', 'calving.min_thickness', '>=', 0, 'NaN', 1, 'Inf', 1, 'numel', [1]) + + return md + # }}} + + def marshall(self, prefix, md, fid): # {{{ + yts = md.constants.yts + WriteData(fid, prefix, 'name', 'md.calving.law', 'data', 2, 'format', 'Integer') + WriteData(fid, prefix, 'object', self, 'fieldname', 'stress_threshold_groundedice', 'format', 'DoubleMat', 'mattype', 1, 'timeserieslength', md.mesh.numberofvertices + 1, 'yts', yts ) + WriteData(fid, prefix, 'object', self, 'fieldname', 'stress_threshold_floatingice', 'format', 'DoubleMat', 'mattype', 1, 'timeserieslength', md.mesh.numberofvertices + 1, 'yts', yts) + WriteData(fid, prefix, 'object', self, 'fieldname', 'min_thickness', 'format', 'Double') + # }}} Index: /issm/trunk-jpl/src/m/classes/debris.m =================================================================== --- /issm/trunk-jpl/src/m/classes/debris.m (revision 27416) +++ /issm/trunk-jpl/src/m/classes/debris.m (revision 27417) @@ -6,14 +6,14 @@ classdef debris properties (SetAccess=public) - spcthickness = NaN; - min_thickness = 0; - stabilization = 0; - packingfraction = 0; - removalmodel = 0; - displacementmodel = 0; - removal_slope_threshold= 0; - removal_stress_threshold=0; - vertex_pairing = NaN; - requested_outputs = {}; + spcthickness = NaN; + min_thickness = 0; + stabilization = 0; + packingfraction = 0; + removalmodel = 0; + displacementmodel = 0; + removal_slope_threshold = 0; + removal_stress_threshold = 0; + vertex_pairing = NaN; + requested_outputs = {}; end methods (Static) @@ -23,6 +23,4 @@ % old fields must be recovered (make sure they are in the deprecated % model properties) - - if verLessThan('matlab','7.9'), @@ -89,6 +87,6 @@ %Slope threshold for removalmodel (1) - self.removal_slope_threshold=0; - + self.removal_slope_threshold=0; + %Stress threshold for removalmodel (2) self.removal_stress_threshold=0; @@ -100,5 +98,5 @@ %Early return - if ~ismember('MasstransportAnalysis',analyses) | (strcmp(solution,'TransientSolution') & md.transient.isdebris==0), return; end + if ~ismember('MasstransportAnalysis',analyses) | (strcmp(solution,'TransientSolution') & md.transient.isdebris==0), return; end md = checkfield(md,'fieldname','debris.spcthickness'); @@ -116,5 +114,5 @@ end % }}} function disp(self) % {{{ - disp(sprintf(' Debris solution parameters:')); + disp(sprintf(' debris solution parameters:')); fielddisplay(self,'spcthickness','debris thickness constraints (NaN means no constraint) [m]'); fielddisplay(self,'min_thickness','minimum debris thickness allowed [m]'); @@ -132,7 +130,4 @@ end % }}} function marshall(self,prefix,md,fid) % {{{ - - yts=md.constants.yts; - WriteData(fid,prefix,'object',self,'fieldname','spcthickness','format','DoubleMat','mattype',1,'timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts); WriteData(fid,prefix,'object',self,'fieldname','min_thickness','format','Double'); Index: /issm/trunk-jpl/src/m/classes/debris.py =================================================================== --- /issm/trunk-jpl/src/m/classes/debris.py (revision 27417) +++ /issm/trunk-jpl/src/m/classes/debris.py (revision 27417) @@ -0,0 +1,129 @@ +import numpy as np + +from checkfield import checkfield +from fielddisplay import fielddisplay +from project3d import project3d +from WriteData import WriteData + + +class debris(object): + """DEBRIS class definition + + Usage: + debris = debris() + """ + + def __init__(self, *args): # {{{ + self.spcthickness = np.nan + self.min_thickness = 0 + self.stabilization = 0 + self.packingfraction = 0 + self.removalmodel = 0 + self.displacementmodel = 0 + self.removal_slope_threshold = 0 + self.removal_stress_threshold = 0 + self.vertex_pairing = np.nan + self.requested_outputs = [] + + nargs = len(args) + if nargs == 0: + self.setdefaultparameters() + elif nargs == 1: + # TODO: Replace the following with constructor + self.setdefaultparameters() + else: + raise Exception('constructor not supported') + # }}} + + def __repr__(self): # {{{ + s = ' debris solution parameters:\n' + s += '{}\n'.format(fielddisplay(self,'spcthickness','debris thickness constraints (NaN means no constraint) [m]')) + s += '{}\n'.format(fielddisplay(self,'min_thickness','minimum debris thickness allowed [m]')) + s += '{}\n'.format(fielddisplay(self,'packingfraction','fraction of debris covered in the ice')) + s += '{}\n'.format(fielddisplay(self,'stabilization','0: no stabilization, 1: artificial diffusion, 2: streamline upwinding, 3: streamline upwind Petrov-Galerkin (SUPG)')) + s += '{}\n'.format(fielddisplay(self,'removalmodel','frontal removal of debris. 0: no removal, 1: Slope-triggered debris removal, 2: driving-stress triggered debris removal')) + s += '{}\n'.format(fielddisplay(self,'displacementmodel','debris displacement. 0: no displacement, 1: ...')) + s += '{}\n'.format(fielddisplay(self,'removal_slope_threshold','critical slope (degrees) for removalmodel (1)')) + s += '{}\n'.format(fielddisplay(self,'removal_stress_threshold','critical stress (Pa) for removalmodel (2)')) + + s += '\n {}\n'.format('Penalty options:') + s += '{}\n'.format(fielddisplay(self,'vertex_pairing','pairs of vertices that are penalized')) + s += '{}\n'.format(fielddisplay(self,'requested_outputs','additional outputs requested')) + return s + # }}} + + def defaultoutputs(self, md): # {{{ + return ['DebrisThickness'] + # }}} + + def setdefaultparameters(self): # {{{ + # Type of stabilization to use 0:nothing 1:artificial_diffusivity 3:Discontinuous Galerkin + self.stabilization = 2 + + # Minimum debris thickness that can be used + self.min_thickness = 0 + + # Fraction of debris covered in the ice + self.packingfraction = 0.01 + + # Type of frontal debris removal + self.removalmodel = 0 + + # Type of debris displacement + self.displacementmodel = 0 + + # Slope threshold for removalmodel (1) + self.removal_slope_threshold = 0 + + # Stress threshold for removalmodel (2) + self.removal_stress_threshold = 0 + + # Default output + self.requested_outputs = ['default'] + return self + # }}} + + def checkconsistency(self, md, solution, analyses): # {{{ + # Early return + if not 'MasstransportAnalysis' in analyses or solution == 'TransientSolution' and not md.transient.isdebris: + return md + + md = checkfield(md, 'fieldname', 'debris.spcthickness') + md = checkfield(md, 'fieldname', 'debris.stabilization', 'values', [0, 1, 2, 3]) + md = checkfield(md, 'fieldname', 'debris.min_thickness', '>=', 0) + md = checkfield(md, 'fieldname', 'debris.packingfraction', '>=', 0) + md = checkfield(md, 'fieldname', 'debris.removalmodel', 'values', [0, 1, 2]) + md = checkfield(md, 'fieldname', 'debris.displacementmodel', 'values', [0, 1, 2]) + md = checkfield(md, 'fieldname', 'debris.removal_slope_threshold', '>=', 0) + md = checkfield(md, 'fieldname', 'debris.removal_stress_threshold', '>=', 0) + md = checkfield(md, 'fieldname', 'debris.requested_outputs', 'stringrow', 1) + + if not np.any(np.isnan(md.stressbalance.vertex_pairing)): + md = checkfield(md, 'fieldname', 'stressbalance.vertex_pairing', '>', 0) + return md + # }}} + + def marshall(self, prefix, md, fid): # {{{ + WriteData(fid, prefix, 'object', self, 'fieldname', 'spcthickness', 'format', 'DoubleMat', 'mattype', 1, 'timeserieslength', md.mesh.numberofvertices + 1, 'yts', md.constants.yts) + WriteData(fid, prefix, 'object', self, 'fieldname', 'min_thickness', 'format', 'Double') + WriteData(fid, prefix, 'object', self, 'fieldname', 'stabilization', 'format', 'Integer') + WriteData(fid, prefix, 'object', self, 'fieldname', 'removalmodel', 'format', 'Integer') + WriteData(fid, prefix, 'object', self, 'fieldname', 'displacementmodel', 'format', 'Integer') + WriteData(fid, prefix, 'object', self, 'fieldname', 'removal_slope_threshold', 'format', 'Double') + WriteData(fid, prefix, 'object', self, 'fieldname', 'removal_stress_threshold', 'format', 'Double') + WriteData(fid, prefix, 'object', self, 'fieldname', 'packingfraction', 'format', 'Double') + WriteData(fid, prefix, 'object', self, 'fieldname', 'vertex_pairing', 'format', 'DoubleMat', 'mattype', 3) + + # Process requested outputs + outputs = self.requested_outputs + indices = [i for i, x in enumerate(outputs) if x == 'default'] + 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.debris.requested_outputs', 'format', 'StringArray') + # }}} + + def extrude(self, md): #{{{ + self.spcthickness = project3d(md, 'vector', self.spcthickness, 'type', 'node') + return + #}}} Index: /issm/trunk-jpl/src/m/classes/dslmme.py =================================================================== --- /issm/trunk-jpl/src/m/classes/dslmme.py (revision 27416) +++ /issm/trunk-jpl/src/m/classes/dslmme.py (revision 27417) @@ -18,7 +18,5 @@ 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. - nargs = len(args) - - if nargs == 0: + if len(args) == 0: self.setdefaultparameters() else: Index: /issm/trunk-jpl/src/m/classes/frontalforcings.py =================================================================== --- /issm/trunk-jpl/src/m/classes/frontalforcings.py (revision 27416) +++ /issm/trunk-jpl/src/m/classes/frontalforcings.py (revision 27417) @@ -15,22 +15,26 @@ def __init__(self): # {{{ self.meltingrate = float('NaN') + self.ablationrate = float('NaN') - #set defaults self.setdefaultparameters() #}}} def __repr__(self): # {{{ - string = ' Frontalforcings parameters:' - string = "%s\n%s" % (string, fielddisplay(self, 'meltingrate', 'melting rate at given location [m / a]')) + s = ' Frontalforcings parameters:' + s += '{}\n'.format(fielddisplay(self, 'meltingrate', 'melting rate at given location [m/a]')) + s += '{}\n'.format(fielddisplay(self, 'ablationrate', 'frontal ablation rate at given location [m/a], it contains both calving and melting')) - return string + return s #}}} def extrude(self, md): # {{{ self.meltingrate = project3d(md, 'vector', self.meltingrate, 'type', 'node') + self.ablationrate = project3d(md, 'vector', self.ablationrate, 'type', 'node') return self #}}} def setdefaultparameters(self): # {{{ + self.meltingrate = float('NaN') + self.ablationrate = float('NaN') return self #}}} @@ -42,4 +46,6 @@ md = checkfield(md, 'fieldname', 'frontalforcings.meltingrate', 'NaN', 1, 'Inf', 1, 'timeseries', 1, '>=', 0) + if not np.isnan(md.frontalforcings.ablationrate): + md = checkfield(md, 'fieldname', 'frontalforcings.ablationrate', 'NaN', 1, 'Inf', 1, 'timeseries', 1) return md # }}} @@ -49,3 +55,5 @@ WriteData(fid, prefix, 'name', 'md.frontalforcings.parameterization', 'data', 1, 'format', 'Integer') WriteData(fid, prefix, 'object', self, 'fieldname', 'meltingrate', 'format', 'DoubleMat', 'mattype', 1, 'timeserieslength', md.mesh.numberofvertices + 1, 'yts', yts, 'scale', 1. / yts) + if not np.isnan(md.frontalforcings.ablationrate): + WriteData(fid, prefix, 'object', self, 'fieldname', 'ablationrate', 'format', 'DoubleMat', 'mattype', 1, 'timeserieslength', md.mesh.numberofvertices + 1, 'yts', yts, 'scale', 1. / yts) # }}} Index: /issm/trunk-jpl/src/m/classes/levelset.py =================================================================== --- /issm/trunk-jpl/src/m/classes/levelset.py (revision 27416) +++ /issm/trunk-jpl/src/m/classes/levelset.py (revision 27417) @@ -56,5 +56,5 @@ md = checkfield(md, 'fieldname', 'levelset.spclevelset', 'Inf', 1, 'timeseries', 1) - md = checkfield(md, 'fieldname', 'levelset.stabilization', 'numel', [1], 'values', [0, 1, 2, 5]) + md = checkfield(md, 'fieldname', 'levelset.stabilization', 'numel', [1], 'values', [0, 1, 2, 5, 6]) md = checkfield(md, 'fieldname', 'levelset.kill_icebergs', 'numel', [1], 'values', [0, 1]) md = checkfield(md, 'fieldname', 'levelset.migration_max', 'numel', [1], 'NaN', 1, 'Inf', 1, '>', 0) Index: /issm/trunk-jpl/src/m/classes/linearbasalforcingsarma.py =================================================================== --- /issm/trunk-jpl/src/m/classes/linearbasalforcingsarma.py (revision 27416) +++ /issm/trunk-jpl/src/m/classes/linearbasalforcingsarma.py (revision 27417) @@ -31,6 +31,5 @@ self.geothermalflux = np.nan - nargs = len(args) - if nargs == 0: + if len(args) == 0: self.setdefaultparameters() else: Index: /issm/trunk-jpl/src/m/classes/lovenumbers.m =================================================================== --- /issm/trunk-jpl/src/m/classes/lovenumbers.m (revision 27416) +++ /issm/trunk-jpl/src/m/classes/lovenumbers.m (revision 27417) @@ -11,7 +11,7 @@ %loading love numbers: - h = []; %provided by PREM model - k = []; %idem - l = []; %idem + h = []; %provided by PREM model + k = []; %idem + l = []; %idem %tidal love numbers for computing rotational feedback: @@ -24,6 +24,6 @@ %time/frequency for visco-elastic love numbers - timefreq = []; - istime = 1; + timefreq = []; + istime = 1; end @@ -65,11 +65,11 @@ %secular fluid love number: - self.tk2secular=0.942; + self.tk2secular=0.942; self.pmtf_colinear=0.0; self.pmtf_ortho=0.0; if maxdeg>=2 - self.pmtf_colinear= (1.0+self.k(3,:))/(1.0-self.tk(3,:)/self.tk2secular); %valid only for elastic regime, not viscous. Also neglects chandler wobble - self.pmtf_ortho= 0.0; + self.pmtf_colinear=(1.0+self.k(3,:))/(1.0-self.tk(3,:)/self.tk2secular); %valid only for elastic regime, not viscous. Also neglects chandler wobble + self.pmtf_ortho=0.0; end %time: Index: /issm/trunk-jpl/src/m/classes/lovenumbers.py =================================================================== --- /issm/trunk-jpl/src/m/classes/lovenumbers.py (revision 27416) +++ /issm/trunk-jpl/src/m/classes/lovenumbers.py (revision 27417) @@ -31,6 +31,4 @@ self.pmtf_colinear = [] self.pmtf_ortho = [] - pmtf_colinear = [] - pmtf_ortho = [] # Time/frequency for visco-elastic love numbers @@ -73,9 +71,4 @@ # Secular fluid love number self.tk2secular = 0.942 - self.pmtf_colinear=0.0 - self.pmtf_ortho=0.0 - if maxdeg>=2: - self.pmtf_colinear= (1.0+self.k[3-1,:])/(1.0-self.tk[3-1,:]/self.tk2secular) #valid only for elastic regime, not viscous. Also neglects chandler wobble - self.pmtf_ortho= 0.0 self.pmtf_colinear = np.array([0.0]).reshape(-1, 1) @@ -131,9 +124,7 @@ WriteData(fid, prefix, 'object', self, 'fieldname', 'tk', 'name', 'md.solidearth.lovenumbers.tk', 'format', 'DoubleMat', 'mattype', 1) WriteData(fid, prefix, 'object', self, 'fieldname', 'tl', 'name', 'md.solidearth.lovenumbers.tl', 'format', 'DoubleMat', 'mattype', 1) - WriteData(fid, prefix, 'object', self, 'fieldname', 'pmtf_colinear', 'name', 'md.solidearth.lovenumbers.pmtf_colinear', 'format', 'DoubleMat', 'mattype', 1) - WriteData(fid, prefix, 'object', self, 'fieldname', 'pmtf_ortho', 'name', 'md.solidearth.lovenumbers.pmtf_ortho', 'format', 'DoubleMat', 'mattype', 1) + WriteData(fid, prefix, 'object', self, 'fieldname', 'pmtf_colinear','name','md.solidearth.lovenumbers.pmtf_colinear','format','DoubleMat','mattype',1) + WriteData(fid, prefix, 'object', self, 'fieldname', 'pmtf_ortho','name','md.solidearth.lovenumbers.pmtf_ortho','format','DoubleMat','mattype',1) WriteData(fid, prefix, 'object', self, 'data', self.tk2secular, 'fieldname', 'lovenumbers.tk2secular', 'format', 'Double') - WriteData(fid, prefix, 'object', self, 'fieldname', 'pmtf_colinear','name','md.solidearth.lovenumbers.pmtf_colinear','format','DoubleMat','mattype',1); - WriteData(fid, prefix, 'object', self, 'fieldname', 'pmtf_ortho','name','md.solidearth.lovenumbers.pmtf_ortho','format','DoubleMat','mattype',1); if (self.istime): @@ -144,6 +135,6 @@ WriteData(fid, prefix, 'object', self, 'fieldname', 'timefreq', 'name', 'md.solidearth.lovenumbers.timefreq', 'format', 'DoubleMat', 'mattype', 1, 'scale', scale); #}}} - + def extrude(self, md): #{{{ return -#}}} + #}}} Index: /issm/trunk-jpl/src/m/classes/mesh3dsurface.py =================================================================== --- /issm/trunk-jpl/src/m/classes/mesh3dsurface.py (revision 27416) +++ /issm/trunk-jpl/src/m/classes/mesh3dsurface.py (revision 27417) @@ -40,5 +40,5 @@ nargs = len(args) - if not nargs: + if nargs == 0: self.setdefaultparameters() elif nargs == 1: Index: /issm/trunk-jpl/src/m/classes/model.m =================================================================== --- /issm/trunk-jpl/src/m/classes/model.m (revision 27416) +++ /issm/trunk-jpl/src/m/classes/model.m (revision 27417) @@ -35,5 +35,5 @@ groundingline = 0; hydrology = 0; - debris = 0; + debris = 0; masstransport = 0; thermal = 0; @@ -214,48 +214,48 @@ %}}} function disp(self) % {{{ - disp(sprintf('%19s: %-22s -- %s','mesh' ,['[1x1 ' class(self.mesh) ']'],'mesh properties')); - disp(sprintf('%19s: %-22s -- %s','mask' ,['[1x1 ' class(self.mask) ']'],'defines grounded and floating elements')); - disp(sprintf('%19s: %-22s -- %s','geometry' ,['[1x1 ' class(self.geometry) ']'],'surface elevation, bedrock topography, ice thickness,...')); - disp(sprintf('%19s: %-22s -- %s','constants' ,['[1x1 ' class(self.constants) ']'],'physical constants')); - disp(sprintf('%19s: %-22s -- %s','smb' ,['[1x1 ' class(self.smb) ']'],'surface mass balance')); - disp(sprintf('%19s: %-22s -- %s','basalforcings' ,['[1x1 ' class(self.basalforcings) ']'],'bed forcings')); - disp(sprintf('%19s: %-22s -- %s','materials' ,['[1x1 ' class(self.materials) ']'],'material properties')); - disp(sprintf('%19s: %-22s -- %s','damage' ,['[1x1 ' class(self.damage) ']'],'parameters for damage evolution solution')); - disp(sprintf('%19s: %-22s -- %s','friction' ,['[1x1 ' class(self.friction) ']'],'basal friction/drag properties')); - disp(sprintf('%19s: %-22s -- %s','flowequation' ,['[1x1 ' class(self.flowequation) ']'],'flow equations')); - disp(sprintf('%19s: %-22s -- %s','timestepping' ,['[1x1 ' class(self.timestepping) ']'],'time stepping for transient models')); - disp(sprintf('%19s: %-22s -- %s','initialization' ,['[1x1 ' class(self.initialization) ']'],'initial guess/state')); - disp(sprintf('%19s: %-22s -- %s','rifts' ,['[1x1 ' class(self.rifts) ']'],'rifts properties')); - disp(sprintf('%19s: %-22s -- %s','solidearth' ,['[1x1 ' class(self.solidearth) ']'],'solidearth inputs and settings')); - disp(sprintf('%19s: %-22s -- %s','dsl' ,['[1x1 ' class(self.dsl) ']'],'dynamic sea-level ')); - disp(sprintf('%19s: %-22s -- %s','debug' ,['[1x1 ' class(self.debug) ']'],'debugging tools (valgrind, gprof)')); - disp(sprintf('%19s: %-22s -- %s','verbose' ,['[1x1 ' class(self.verbose) ']'],'verbosity level in solve')); - disp(sprintf('%19s: %-22s -- %s','settings' ,['[1x1 ' class(self.settings) ']'],'settings properties')); - disp(sprintf('%19s: %-22s -- %s','toolkits' ,['[1x1 ' class(self.toolkits) ']'],'PETSc options for each solution')); - disp(sprintf('%19s: %-22s -- %s','cluster' ,['[1x1 ' class(self.cluster) ']'],'cluster parameters (number of CPUs...)')); - disp(sprintf('%19s: %-22s -- %s','balancethickness',['[1x1 ' class(self.balancethickness) ']'],'parameters for balancethickness solution')); - disp(sprintf('%19s: %-22s -- %s','stressbalance' ,['[1x1 ' class(self.stressbalance) ']'],'parameters for stressbalance solution')); - disp(sprintf('%19s: %-22s -- %s','groundingline' ,['[1x1 ' class(self.groundingline) ']'],'parameters for groundingline solution')); - disp(sprintf('%19s: %-22s -- %s','hydrology' ,['[1x1 ' class(self.hydrology) ']'],'parameters for hydrology solution')); - disp(sprintf('%19s: %-22s -- %s','debris' ,['[1x1 ' class(self.debris) ']'],'parameters for debris solution')); - disp(sprintf('%19s: %-22s -- %s','masstransport' ,['[1x1 ' class(self.masstransport) ']'],'parameters for masstransport solution')); - disp(sprintf('%19s: %-22s -- %s','thermal' ,['[1x1 ' class(self.thermal) ']'],'parameters for thermal solution')); - disp(sprintf('%19s: %-22s -- %s','steadystate' ,['[1x1 ' class(self.steadystate) ']'],'parameters for steadystate solution')); - disp(sprintf('%19s: %-22s -- %s','transient' ,['[1x1 ' class(self.transient) ']'],'parameters for transient solution')); - disp(sprintf('%19s: %-22s -- %s','levelset' ,['[1x1 ' class(self.levelset) ']'],'parameters for moving boundaries (level-set method)')); - disp(sprintf('%19s: %-22s -- %s','calving' ,['[1x1 ' class(self.calving) ']'],'parameters for calving')); - disp(sprintf('%19s: %-22s -- %s','frontalforcings' ,['[1x1 ' class(self.frontalforcings) ']'],'parameters for frontalforcings')); - disp(sprintf('%19s: %-22s -- %s','esa' ,['[1x1 ' class(self.esa) ']'],'parameters for elastic adjustment solution')); - disp(sprintf('%19s: %-22s -- %s','love' ,['[1x1 ' class(self.love) ']'],'parameters for love solution')); - disp(sprintf('%19s: %-22s -- %s','sampling' ,['[1x1 ' class(self.sampling) ']'],'parameters for stochastic sampler')); - disp(sprintf('%19s: %-22s -- %s','autodiff' ,['[1x1 ' class(self.autodiff) ']'],'automatic differentiation parameters')); - disp(sprintf('%19s: %-22s -- %s','inversion' ,['[1x1 ' class(self.inversion) ']'],'parameters for inverse methods')); - disp(sprintf('%19s: %-22s -- %s','qmu' ,['[1x1 ' class(self.qmu) ']'],'Dakota properties')); - disp(sprintf('%19s: %-22s -- %s','amr' ,['[1x1 ' class(self.amr) ']'],'adaptive mesh refinement properties')); - disp(sprintf('%19s: %-22s -- %s','outputdefinition',['[1x1 ' class(self.outputdefinition) ']'],'output definition')); - disp(sprintf('%19s: %-22s -- %s','results' ,['[1x1 ' class(self.results) ']'],'model results')); - disp(sprintf('%19s: %-22s -- %s','radaroverlay' ,['[1x1 ' class(self.radaroverlay) ']'],'radar image for plot overlay')); - disp(sprintf('%19s: %-22s -- %s','miscellaneous' ,['[1x1 ' class(self.miscellaneous) ']'],'miscellaneous fields')); - disp(sprintf('%19s: %-22s -- %s','stochasticforcing',['[1x1 ' class(self.stochasticforcing) ']'],'stochasticity applied to model forcings')); + disp(sprintf('%19s: %-23s -- %s','mesh' ,['[1x1 ' class(self.mesh) ']'],'mesh properties')); + disp(sprintf('%19s: %-23s -- %s','mask' ,['[1x1 ' class(self.mask) ']'],'defines grounded and floating elements')); + disp(sprintf('%19s: %-23s -- %s','geometry' ,['[1x1 ' class(self.geometry) ']'],'surface elevation, bedrock topography, ice thickness,...')); + disp(sprintf('%19s: %-23s -- %s','constants' ,['[1x1 ' class(self.constants) ']'],'physical constants')); + disp(sprintf('%19s: %-23s -- %s','smb' ,['[1x1 ' class(self.smb) ']'],'surface mass balance')); + disp(sprintf('%19s: %-23s -- %s','basalforcings' ,['[1x1 ' class(self.basalforcings) ']'],'bed forcings')); + disp(sprintf('%19s: %-23s -- %s','materials' ,['[1x1 ' class(self.materials) ']'],'material properties')); + disp(sprintf('%19s: %-23s -- %s','damage' ,['[1x1 ' class(self.damage) ']'],'parameters for damage evolution solution')); + disp(sprintf('%19s: %-23s -- %s','friction' ,['[1x1 ' class(self.friction) ']'],'basal friction/drag properties')); + disp(sprintf('%19s: %-23s -- %s','flowequation' ,['[1x1 ' class(self.flowequation) ']'],'flow equations')); + disp(sprintf('%19s: %-23s -- %s','timestepping' ,['[1x1 ' class(self.timestepping) ']'],'time stepping for transient models')); + disp(sprintf('%19s: %-23s -- %s','initialization' ,['[1x1 ' class(self.initialization) ']'],'initial guess/state')); + disp(sprintf('%19s: %-23s -- %s','rifts' ,['[1x1 ' class(self.rifts) ']'],'rifts properties')); + disp(sprintf('%19s: %-23s -- %s','solidearth' ,['[1x1 ' class(self.solidearth) ']'],'solidearth inputs and settings')); + disp(sprintf('%19s: %-23s -- %s','dsl' ,['[1x1 ' class(self.dsl) ']'],'dynamic sea-level ')); + disp(sprintf('%19s: %-23s -- %s','debug' ,['[1x1 ' class(self.debug) ']'],'debugging tools (valgrind, gprof)')); + disp(sprintf('%19s: %-23s -- %s','verbose' ,['[1x1 ' class(self.verbose) ']'],'verbosity level in solve')); + disp(sprintf('%19s: %-23s -- %s','settings' ,['[1x1 ' class(self.settings) ']'],'settings properties')); + disp(sprintf('%19s: %-23s -- %s','toolkits' ,['[1x1 ' class(self.toolkits) ']'],'PETSc options for each solution')); + disp(sprintf('%19s: %-23s -- %s','cluster' ,['[1x1 ' class(self.cluster) ']'],'cluster parameters (number of CPUs...)')); + disp(sprintf('%19s: %-23s -- %s','balancethickness',['[1x1 ' class(self.balancethickness) ']'],'parameters for balancethickness solution')); + disp(sprintf('%19s: %-23s -- %s','stressbalance' ,['[1x1 ' class(self.stressbalance) ']'],'parameters for stressbalance solution')); + disp(sprintf('%19s: %-23s -- %s','groundingline' ,['[1x1 ' class(self.groundingline) ']'],'parameters for groundingline solution')); + disp(sprintf('%19s: %-23s -- %s','hydrology' ,['[1x1 ' class(self.hydrology) ']'],'parameters for hydrology solution')); + disp(sprintf('%19s: %-23s -- %s','debris' ,['[1x1 ' class(self.debris) ']'],'parameters for debris solution')); + disp(sprintf('%19s: %-23s -- %s','masstransport' ,['[1x1 ' class(self.masstransport) ']'],'parameters for masstransport solution')); + disp(sprintf('%19s: %-23s -- %s','thermal' ,['[1x1 ' class(self.thermal) ']'],'parameters for thermal solution')); + disp(sprintf('%19s: %-23s -- %s','steadystate' ,['[1x1 ' class(self.steadystate) ']'],'parameters for steadystate solution')); + disp(sprintf('%19s: %-23s -- %s','transient' ,['[1x1 ' class(self.transient) ']'],'parameters for transient solution')); + disp(sprintf('%19s: %-23s -- %s','levelset' ,['[1x1 ' class(self.levelset) ']'],'parameters for moving boundaries (level-set method)')); + disp(sprintf('%19s: %-23s -- %s','calving' ,['[1x1 ' class(self.calving) ']'],'parameters for calving')); + disp(sprintf('%19s: %-23s -- %s','frontalforcings' ,['[1x1 ' class(self.frontalforcings) ']'],'parameters for frontalforcings')); + disp(sprintf('%19s: %-23s -- %s','esa' ,['[1x1 ' class(self.esa) ']'],'parameters for elastic adjustment solution')); + disp(sprintf('%19s: %-23s -- %s','love' ,['[1x1 ' class(self.love) ']'],'parameters for love solution')); + disp(sprintf('%19s: %-23s -- %s','sampling' ,['[1x1 ' class(self.sampling) ']'],'parameters for stochastic sampler')); + disp(sprintf('%19s: %-23s -- %s','autodiff' ,['[1x1 ' class(self.autodiff) ']'],'automatic differentiation parameters')); + disp(sprintf('%19s: %-23s -- %s','inversion' ,['[1x1 ' class(self.inversion) ']'],'parameters for inverse methods')); + disp(sprintf('%19s: %-23s -- %s','qmu' ,['[1x1 ' class(self.qmu) ']'],'Dakota properties')); + disp(sprintf('%19s: %-23s -- %s','amr' ,['[1x1 ' class(self.amr) ']'],'adaptive mesh refinement properties')); + disp(sprintf('%19s: %-23s -- %s','outputdefinition',['[1x1 ' class(self.outputdefinition) ']'],'output definition')); + disp(sprintf('%19s: %-23s -- %s','results' ,['[1x1 ' class(self.results) ']'],'model results')); + disp(sprintf('%19s: %-23s -- %s','radaroverlay' ,['[1x1 ' class(self.radaroverlay) ']'],'radar image for plot overlay')); + disp(sprintf('%19s: %-23s -- %s','miscellaneous' ,['[1x1 ' class(self.miscellaneous) ']'],'miscellaneous fields')); + disp(sprintf('%19s: %-23s -- %s','stochasticforcing',['[1x1 ' class(self.stochasticforcing) ']'],'stochasticity applied to model forcings')); end % }}} function md = setdefaultparameters(md,planet) % {{{ @@ -286,5 +286,5 @@ md.stressbalance = stressbalance(); md.hydrology = hydrologyshreve(); - md.debris = debris(); + md.debris = debris(); md.masstransport = masstransport(); md.thermal = thermal(); @@ -959,5 +959,5 @@ md2.mesh.vertexonboundary = zeros(md2.mesh.numberofvertices,1); md2.mesh.vertexonboundary(md2.mesh.segments(:,1:2)) = 1; - %Deal with boudary + %Deal with boundary md2.mesh.vertexonboundary = [md.mesh.vertexonboundary;sum(md.mesh.vertexonboundary(edges),2)==2]; md2.mesh.elementconnectivity=bamgmesh_out.ElementConnectivity; Index: /issm/trunk-jpl/src/m/classes/model.py =================================================================== --- /issm/trunk-jpl/src/m/classes/model.py (revision 27416) +++ /issm/trunk-jpl/src/m/classes/model.py (revision 27417) @@ -45,4 +45,5 @@ from hydrologypism import hydrologypism from hydrologyshakti import hydrologyshakti +from debris import debris from masstransport import masstransport from thermal import thermal @@ -97,5 +98,4 @@ self.dsl = None self.solidearth = None - self.debug = None self.verbose = None @@ -103,9 +103,9 @@ self.toolkits = None self.cluster = None - self.balancethickness = None self.stressbalance = None self.groundingline = None self.hydrology = None + self.debris = None self.masstransport = None self.thermal = None @@ -118,5 +118,4 @@ self.esa = None self.sampling = None - self.autodiff = None self.inversion = None @@ -143,47 +142,47 @@ # already converted .__repr__ method for examples) # - s = "%19s: %-22s -- %s" % ("mesh", "[%s %s]" % ("1x1", obj.mesh.__class__.__name__), "mesh properties") - s = "%s\n%s" % (s, "%19s: %-22s -- %s" % ("mask", "[%s %s]" % ("1x1", obj.mask.__class__.__name__), "defines grounded and floating elements")) - s = "%s\n%s" % (s, "%19s: %-22s -- %s" % ("geometry", "[%s %s]" % ("1x1", obj.geometry.__class__.__name__), "surface elevation, bedrock topography, ice thickness, ...")) - s = "%s\n%s" % (s, "%19s: %-22s -- %s" % ("constants", "[%s %s]" % ("1x1", obj.constants.__class__.__name__), "physical constants")) - s = "%s\n%s" % (s, "%19s: %-22s -- %s" % ("smb", "[%s %s]" % ("1x1", obj.smb.__class__.__name__), "surface mass balance")) - s = "%s\n%s" % (s, "%19s: %-22s -- %s" % ("basalforcings", "[%s %s]" % ("1x1", obj.basalforcings.__class__.__name__), "bed forcings")) - s = "%s\n%s" % (s, "%19s: %-22s -- %s" % ("materials", "[%s %s]" % ("1x1", obj.materials.__class__.__name__), "material properties")) - s = "%s\n%s" % (s, "%19s: %-22s -- %s" % ("damage", "[%s %s]" % ("1x1", obj.damage.__class__.__name__), "damage propagation laws")) - s = "%s\n%s" % (s, "%19s: %-22s -- %s" % ("friction", "[%s %s]" % ("1x1", obj.friction.__class__.__name__), "basal friction / drag properties")) - s = "%s\n%s" % (s, "%19s: %-22s -- %s" % ("flowequation", "[%s %s]" % ("1x1", obj.flowequation.__class__.__name__), "flow equations")) - s = "%s\n%s" % (s, "%19s: %-22s -- %s" % ("timestepping", "[%s %s]" % ("1x1", obj.timestepping.__class__.__name__), "time stepping for transient models")) - s = "%s\n%s" % (s, "%19s: %-22s -- %s" % ("initialization", "[%s %s]" % ("1x1", obj.initialization.__class__.__name__), "initial guess / state")) - s = "%s\n%s" % (s, "%19s: %-22s -- %s" % ("rifts", "[%s %s]" % ("1x1", obj.rifts.__class__.__name__), "rifts properties")) - s = "%s\n%s" % (s, "%19s: %-22s -- %s" % ("solidearth", "[%s %s]" % ("1x1", obj.solidearth.__class__.__name__), "solidearth inputs and settings")) - s = "%s\n%s" % (s, "%19s: %-22s -- %s" % ("dsl", "[%s %s]" % ("1x1", obj.dsl.__class__.__name__), "dynamic sea level")) - s = "%s\n%s" % (s, "%19s: %-22s -- %s" % ("debug", "[%s %s]" % ("1x1", obj.debug.__class__.__name__), "debugging tools (valgrind, gprof)")) - s = "%s\n%s" % (s, "%19s: %-22s -- %s" % ("verbose", "[%s %s]" % ("1x1", obj.verbose.__class__.__name__), "verbosity level in solve")) - s = "%s\n%s" % (s, "%19s: %-22s -- %s" % ("settings", "[%s %s]" % ("1x1", obj.settings.__class__.__name__), "settings properties")) - s = "%s\n%s" % (s, "%19s: %-22s -- %s" % ("toolkits", "[%s %s]" % ("1x1", obj.toolkits.__class__.__name__), "PETSc options for each solution")) - s = "%s\n%s" % (s, "%19s: %-22s -- %s" % ("cluster", "[%s %s]" % ("1x1", obj.cluster.__class__.__name__), "cluster parameters (number of CPUs...)")) - s = "%s\n%s" % (s, "%19s: %-22s -- %s" % ("balancethickness", "[%s %s]" % ("1x1", obj.balancethickness.__class__.__name__), "parameters for balancethickness solution")) - s = "%s\n%s" % (s, "%19s: %-22s -- %s" % ("stressbalance", "[%s %s]" % ("1x1", obj.stressbalance.__class__.__name__), "parameters for stressbalance solution")) - s = "%s\n%s" % (s, "%19s: %-22s -- %s" % ("groundingline", "[%s %s]" % ("1x1", obj.groundingline.__class__.__name__), "parameters for groundingline solution")) - s = "%s\n%s" % (s, "%19s: %-22s -- %s" % ("hydrology", "[%s %s]" % ("1x1", obj.hydrology.__class__.__name__), "parameters for hydrology solution")) - s = "%s\n%s" % (s, "%19s: %-22s -- %s" % ("masstransport", "[%s %s]" % ("1x1", obj.masstransport.__class__.__name__), "parameters for masstransport solution")) - s = "%s\n%s" % (s, "%19s: %-22s -- %s" % ("thermal", "[%s %s]" % ("1x1", obj.thermal.__class__.__name__), "parameters for thermal solution")) - s = "%s\n%s" % (s, "%19s: %-22s -- %s" % ("steadystate", "[%s %s]" % ("1x1", obj.steadystate.__class__.__name__), "parameters for steadystate solution")) - s = "%s\n%s" % (s, "%19s: %-22s -- %s" % ("transient", "[%s %s]" % ("1x1", obj.transient.__class__.__name__), "parameters for transient solution")) - s = "%s\n%s" % (s, "%19s: %-22s -- %s" % ("levelset", "[%s %s]" % ("1x1", obj.levelset.__class__.__name__), "parameters for moving boundaries (level-set method)")) - 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' % ("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")) - s = "%s\n%s" % (s, '%19s: %-22s -- %s' % ("love", "[%s %s]" % ("1x1", obj.love.__class__.__name__), "parameters for love solution")) - s = "%s\n%s" % (s, "%19s: %-22s -- %s" % ("autodiff", "[%s %s]" % ("1x1", obj.autodiff.__class__.__name__), "automatic differentiation parameters")) - s = "%s\n%s" % (s, "%19s: %-22s -- %s" % ("inversion", "[%s %s]" % ("1x1", obj.inversion.__class__.__name__), "parameters for inverse methods")) - s = "%s\n%s" % (s, "%19s: %-22s -- %s" % ("qmu", "[%s %s]" % ("1x1", obj.qmu.__class__.__name__), "Dakota properties")) - s = "%s\n%s" % (s, "%19s: %-22s -- %s" % ("amr", "[%s %s]" % ("1x1", obj.amr.__class__.__name__), "adaptive mesh refinement properties")) - s = "%s\n%s" % (s, "%19s: %-22s -- %s" % ("outputdefinition", "[%s %s]" % ("1x1", obj.outputdefinition.__class__.__name__), "output definition")) - s = "%s\n%s" % (s, "%19s: %-22s -- %s" % ("results", "[%s %s]" % ("1x1", obj.results.__class__.__name__), "model results")) - s = "%s\n%s" % (s, "%19s: %-22s -- %s" % ("radaroverlay", "[%s %s]" % ("1x1", obj.radaroverlay.__class__.__name__), "radar image for plot overlay")) - s = "%s\n%s" % (s, "%19s: %-22s -- %s" % ("miscellaneous", "[%s %s]" % ("1x1", obj.miscellaneous.__class__.__name__), "miscellaneous fields")) - s = "%s\n%s" % (s, "%19s: %-22s -- %s" % ("stochasticforcing", "[%s %s]" % ("1x1", obj.stochasticforcing.__class__.__name__), "stochasticity applied to model forcings")) + s = '%19s: %-23s -- %s' % ('mesh', '[%s %s]' % ('1x1', obj.mesh.__class__.__name__), 'mesh properties') + s = '%s\n%s' % (s, '%19s: %-23s -- %s' % ('mask', '[%s %s]' % ('1x1', obj.mask.__class__.__name__), 'defines grounded and floating elements')) + s = '%s\n%s' % (s, '%19s: %-23s -- %s' % ('geometry', '[%s %s]' % ('1x1', obj.geometry.__class__.__name__), 'surface elevation, bedrock topography, ice thickness, ...')) + s = '%s\n%s' % (s, '%19s: %-23s -- %s' % ('constants', '[%s %s]' % ('1x1', obj.constants.__class__.__name__), 'physical constants')) + s = '%s\n%s' % (s, '%19s: %-23s -- %s' % ('smb', '[%s %s]' % ('1x1', obj.smb.__class__.__name__), 'surface mass balance')) + s = '%s\n%s' % (s, '%19s: %-23s -- %s' % ('basalforcings', '[%s %s]' % ('1x1', obj.basalforcings.__class__.__name__), 'bed forcings')) + s = '%s\n%s' % (s, '%19s: %-23s -- %s' % ('materials', '[%s %s]' % ('1x1', obj.materials.__class__.__name__), 'material properties')) + s = '%s\n%s' % (s, '%19s: %-23s -- %s' % ('damage', '[%s %s]' % ('1x1', obj.damage.__class__.__name__), 'damage propagation laws')) + s = '%s\n%s' % (s, '%19s: %-23s -- %s' % ('friction', '[%s %s]' % ('1x1', obj.friction.__class__.__name__), 'basal friction / drag properties')) + s = '%s\n%s' % (s, '%19s: %-23s -- %s' % ('flowequation', '[%s %s]' % ('1x1', obj.flowequation.__class__.__name__), 'flow equations')) + s = '%s\n%s' % (s, '%19s: %-23s -- %s' % ('timestepping', '[%s %s]' % ('1x1', obj.timestepping.__class__.__name__), 'time stepping for transient models')) + s = '%s\n%s' % (s, '%19s: %-23s -- %s' % ('initialization', '[%s %s]' % ('1x1', obj.initialization.__class__.__name__), 'initial guess / state')) + s = '%s\n%s' % (s, '%19s: %-23s -- %s' % ('rifts', '[%s %s]' % ('1x1', obj.rifts.__class__.__name__), 'rifts properties')) + s = '%s\n%s' % (s, '%19s: %-23s -- %s' % ('solidearth', '[%s %s]' % ('1x1', obj.solidearth.__class__.__name__), 'solidearth inputs and settings')) + s = '%s\n%s' % (s, '%19s: %-23s -- %s' % ('dsl', '[%s %s]' % ('1x1', obj.dsl.__class__.__name__), 'dynamic sea level')) + s = '%s\n%s' % (s, '%19s: %-23s -- %s' % ('debug', '[%s %s]' % ('1x1', obj.debug.__class__.__name__), 'debugging tools (valgrind, gprof)')) + s = '%s\n%s' % (s, '%19s: %-23s -- %s' % ('verbose', '[%s %s]' % ('1x1', obj.verbose.__class__.__name__), 'verbosity level in solve')) + s = '%s\n%s' % (s, '%19s: %-23s -- %s' % ('settings', '[%s %s]' % ('1x1', obj.settings.__class__.__name__), 'settings properties')) + s = '%s\n%s' % (s, '%19s: %-23s -- %s' % ('toolkits', '[%s %s]' % ('1x1', obj.toolkits.__class__.__name__), 'PETSc options for each solution')) + s = '%s\n%s' % (s, '%19s: %-23s -- %s' % ('cluster', '[%s %s]' % ('1x1', obj.cluster.__class__.__name__), 'cluster parameters (number of CPUs...)')) + s = '%s\n%s' % (s, '%19s: %-23s -- %s' % ('balancethickness', '[%s %s]' % ('1x1', obj.balancethickness.__class__.__name__), 'parameters for balancethickness solution')) + s = '%s\n%s' % (s, '%19s: %-23s -- %s' % ('stressbalance', '[%s %s]' % ('1x1', obj.stressbalance.__class__.__name__), 'parameters for stressbalance solution')) + s = '%s\n%s' % (s, '%19s: %-23s -- %s' % ('groundingline', '[%s %s]' % ('1x1', obj.groundingline.__class__.__name__), 'parameters for groundingline solution')) + s = '%s\n%s' % (s, '%19s: %-23s -- %s' % ('hydrology', '[%s %s]' % ('1x1', obj.hydrology.__class__.__name__), 'parameters for hydrology solution')) + s = '%s\n%s' % (s, '%19s: %-23s -- %s' % ('masstransport', '[%s %s]' % ('1x1', obj.masstransport.__class__.__name__), 'parameters for masstransport solution')) + s = '%s\n%s' % (s, '%19s: %-23s -- %s' % ('thermal', '[%s %s]' % ('1x1', obj.thermal.__class__.__name__), 'parameters for thermal solution')) + s = '%s\n%s' % (s, '%19s: %-23s -- %s' % ('steadystate', '[%s %s]' % ('1x1', obj.steadystate.__class__.__name__), 'parameters for steadystate solution')) + s = '%s\n%s' % (s, '%19s: %-23s -- %s' % ('transient', '[%s %s]' % ('1x1', obj.transient.__class__.__name__), 'parameters for transient solution')) + s = '%s\n%s' % (s, '%19s: %-23s -- %s' % ('levelset', '[%s %s]' % ('1x1', obj.levelset.__class__.__name__), 'parameters for moving boundaries (level-set method)')) + s = '%s\n%s' % (s, '%19s: %-23s -- %s' % ('calving', '[%s %s]' % ('1x1', obj.calving.__class__.__name__), 'parameters for calving')) + s = '%s\n%s' % (s, '%19s: %-23s -- %s' % ('frontalforcings', '[%s %s]' % ('1x1', obj.frontalforcings.__class__.__name__), 'parameters for frontalforcings')) + s = '%s\n%s' % (s, '%19s: %-23s -- %s' % ('esa', '[%s %s]' % ('1x1', obj.esa.__class__.__name__), 'parameters for elastic adjustment solution')) + s = '%s\n%s' % (s, '%19s: %-23s -- %s' % ('sampling', '[%s %s]' % ('1x1', obj.sampling.__class__.__name__), 'parameters for stochastic sampler')) + s = '%s\n%s' % (s, '%19s: %-23s -- %s' % ('love', '[%s %s]' % ('1x1', obj.love.__class__.__name__), 'parameters for love solution')) + s = '%s\n%s' % (s, '%19s: %-23s -- %s' % ('autodiff', '[%s %s]' % ('1x1', obj.autodiff.__class__.__name__), 'automatic differentiation parameters')) + s = '%s\n%s' % (s, '%19s: %-23s -- %s' % ('inversion', '[%s %s]' % ('1x1', obj.inversion.__class__.__name__), 'parameters for inverse methods')) + s = '%s\n%s' % (s, '%19s: %-23s -- %s' % ('qmu', '[%s %s]' % ('1x1', obj.qmu.__class__.__name__), 'Dakota properties')) + s = '%s\n%s' % (s, '%19s: %-23s -- %s' % ('amr', '[%s %s]' % ('1x1', obj.amr.__class__.__name__), 'adaptive mesh refinement properties')) + s = '%s\n%s' % (s, '%19s: %-23s -- %s' % ('outputdefinition', '[%s %s]' % ('1x1', obj.outputdefinition.__class__.__name__), 'output definition')) + s = '%s\n%s' % (s, '%19s: %-23s -- %s' % ('results', '[%s %s]' % ('1x1', obj.results.__class__.__name__), 'model results')) + s = '%s\n%s' % (s, '%19s: %-23s -- %s' % ('radaroverlay', '[%s %s]' % ('1x1', obj.radaroverlay.__class__.__name__), 'radar image for plot overlay')) + s = '%s\n%s' % (s, '%19s: %-23s -- %s' % ('miscellaneous', '[%s %s]' % ('1x1', obj.miscellaneous.__class__.__name__), 'miscellaneous fields')) + s = '%s\n%s' % (s, '%19s: %-23s -- %s' % ('stochasticforcing', '[%s %s]' % ('1x1', obj.stochasticforcing.__class__.__name__), 'stochasticity applied to model forcings')) return s #}}} @@ -216,4 +215,5 @@ 'groundingline', 'hydrology', + 'debris', 'masstransport', 'thermal', @@ -264,4 +264,5 @@ self.stressbalance = stressbalance() self.hydrology = hydrologyshreve() + self.debris = debris() self.masstransport = masstransport() self.thermal = thermal() @@ -597,5 +598,5 @@ #some checks on list of arguments if len(args) > 3 or len(args) < 1: - raise RuntimeError("extrude error message") + raise RuntimeError('extrude error message') #Extrude the mesh @@ -603,5 +604,5 @@ clist = args[0] if any(clist < 0) or any(clist > 1): - raise TypeError("extrusioncoefficients must be between 0 and 1") + raise TypeError('extrusioncoefficients must be between 0 and 1') clist.extend([0., 1.]) clist.sort() @@ -611,5 +612,5 @@ elif len(args) == 2: #one polynomial law if args[1] <= 0: - raise TypeError("extrusionexponent must be >= 0") + raise TypeError('extrusionexponent must be >= 0') numlayers = args[0] extrusionlist = (np.arange(0., float(numlayers - 1) + 1., 1.) / float(numlayers - 1))**args[1] @@ -621,5 +622,5 @@ if args[1] <= 0 or args[2] <= 0: - raise TypeError("lower and upper extrusionexponents must be >= 0") + raise TypeError('lower and upper extrusionexponents must be >= 0') lowerextrusionlist = (np.arange(0., 1. + 2. / float(numlayers - 1), 2. / float(numlayers - 1)))**lowerexp / 2. @@ -628,7 +629,7 @@ if numlayers < 2: - raise TypeError("number of layers should be at least 2") + raise TypeError('number of layers should be at least 2') if md.mesh.__class__.__name__ == 'mesh3dprisms': - raise TypeError("Cannot extrude a 3d mesh (extrude cannot be called more than once)") + raise TypeError('Cannot extrude a 3d mesh (extrude cannot be called more than once)') #Initialize with the 2d mesh @@ -732,4 +733,5 @@ md.frontalforcings.extrude(md) md.hydrology.extrude(md) + md.debris.extrude(md) md.solidearth.extrude(md) md.dsl.extrude(md) @@ -774,5 +776,5 @@ # Check that the model is really a 3d model if md.mesh.elementtype() != 'Penta': - raise Exception("collapse error message: only a 3d mesh can be collapsed") + raise Exception('collapse error message: only a 3d mesh can be collapsed') # Start with changing all the fields from the 3d mesh Index: /issm/trunk-jpl/src/m/classes/modellist.m =================================================================== --- /issm/trunk-jpl/src/m/classes/modellist.m (revision 27416) +++ /issm/trunk-jpl/src/m/classes/modellist.m (revision 27417) @@ -298,5 +298,5 @@ % BuildMultipleQueueingScript(executionpath,codepath) -disp('building queueing script'); +disp('building queuing script'); %First try and figure out if there is a special script for this particular cluster @@ -363,12 +363,12 @@ end% }}} function md=LaunchMultipleQueueJobgemini(cluster,name,executionpath)% {{{ -%LAUNCHMULTIPLEQUEUEJOBGEMINI - Launch multiple queueing script on Gemini cluster +%LAUNCHMULTIPLEQUEUEJOBGEMINI - Launch multiple queuing script on Gemini cluster % % Usage: % LaunchMultipleQueueJobgemini(cluster,name,executionpath) -%first, check we have the binary file and the queueing script +%first, check we have the binary file and the queuing script if ~exist([ name '.queue'],'file'), - error('LaunchMultipleQueueJobgemini error message: queueing script issing, cannot go forward'); + error('LaunchMultipleQueueJobgemini error message: queuing script issing, cannot go forward'); end @@ -378,5 +378,5 @@ %upload both files to cluster -disp('uploading input file, queueing script and variables script'); +disp('uploading input file, queuing script and variables script'); eval(['!scp ModelList.tar.gz ' name '.queue ' cluster ':' executionpath]); Index: /issm/trunk-jpl/src/m/classes/sampling.m =================================================================== --- /issm/trunk-jpl/src/m/classes/sampling.m (revision 27416) +++ /issm/trunk-jpl/src/m/classes/sampling.m (revision 27417) @@ -75,5 +75,6 @@ if(md.sampling.robin) md = checkfield(md,'fieldname','sampling.beta','NaN',1,'Inf',1,'size',[md.mesh.numberofvertices 1],'>',0); - end + end + md = checkfield(md,'fieldname','sampling.phi','NaN',1,'Inf',1,'numel',1,'>=',0); md = checkfield(md,'fieldname','sampling.alpha','NaN',1,'Inf',1,'numel',1,'>',0); md = checkfield(md,'fieldname','sampling.seed','NaN',1,'Inf',1,'numel',1); Index: /issm/trunk-jpl/src/m/classes/sampling.py =================================================================== --- /issm/trunk-jpl/src/m/classes/sampling.py (revision 27416) +++ /issm/trunk-jpl/src/m/classes/sampling.py (revision 27417) @@ -20,5 +20,5 @@ self.tau = 0 self.beta = np.nan - self.phi = 0 + self.phi = np.nan self.alpha = 0 self.robin = 0 @@ -54,15 +54,10 @@ def setdefaultparameters(self): #{{{ - # Scaling coefficient - self.tau = 1 # Apply Robin boundary conditions self.robin = 0 - # Temporal correlation factor - self.phi = 0 - # Exponent in fraction SPDE (default: 2, biLaplacian covariance operator) - self.alpha = 2 + self.alpha = 2 # Default # Seed for pseudorandom number generator (default: -1, for random seed) @@ -99,7 +94,7 @@ def marshall(self, prefix, md, fid): #{{{ WriteData(fid, prefix, 'object', self, 'fieldname', 'kappa', 'format', 'DoubleMat', 'mattype', 1) - WriteData(fid, prefix, 'object', self, 'fieldname', 'tau', 'format', 'Double') + WriteData(fid, prefix, 'object', self, 'fieldname', 'tau', 'format', 'DoubleMat', 'mattype', 1) WriteData(fid, prefix, 'object', self, 'fieldname', 'beta', 'format', 'DoubleMat', 'mattype', 1) - WriteData(fid, prefix, 'object', self, 'fieldname', 'phi', 'format', 'Double') + WriteData(fid, prefix, 'object', self, 'fieldname', 'phi', 'format', 'DoubleMat', 'mattype', 1) WriteData(fid, prefix, 'object', self, 'fieldname', 'alpha', 'format', 'Integer') WriteData(fid, prefix, 'object', self, 'fieldname', 'robin', 'format', 'Boolean') Index: /issm/trunk-jpl/src/m/classes/solidearthsettings.m =================================================================== --- /issm/trunk-jpl/src/m/classes/solidearthsettings.m (revision 27416) +++ /issm/trunk-jpl/src/m/classes/solidearthsettings.m (revision 27417) @@ -88,5 +88,5 @@ self.cross_section_shape=1; %square as default (see iedge in GiaDeflectionCorex) - %no grd model by default: + %grd model by default self.grdmodel=1; @@ -171,5 +171,4 @@ WriteData(fid,prefix,'object',self,'fieldname','viscous','name','md.solidearth.settings.viscous','format','Boolean'); WriteData(fid,prefix,'object',self,'fieldname','rotation','name','md.solidearth.settings.rotation','format','Boolean'); - WriteData(fid,prefix,'object',self,'fieldname','rotation','name','md.solidearth.settings.satellitegravity','format','Boolean'); WriteData(fid,prefix,'object',self,'fieldname','grdocean','name','md.solidearth.settings.grdocean','format','Boolean'); WriteData(fid,prefix,'object',self,'fieldname','ocean_area_scaling','name','md.solidearth.settings.ocean_area_scaling','format','Boolean'); Index: /issm/trunk-jpl/src/m/classes/solidearthsettings.py =================================================================== --- /issm/trunk-jpl/src/m/classes/solidearthsettings.py (revision 27416) +++ /issm/trunk-jpl/src/m/classes/solidearthsettings.py (revision 27417) @@ -93,6 +93,6 @@ self.cross_section_shape = 1 # Square as default (see iedge in GiaDeflectionCorex) - # No GRD model by default - self.grdmodel = 0 + # GRD model by default + self.grdmodel = 1 # }}} Index: /issm/trunk-jpl/src/m/classes/transient.m =================================================================== --- /issm/trunk-jpl/src/m/classes/transient.m (revision 27416) +++ /issm/trunk-jpl/src/m/classes/transient.m (revision 27417) @@ -144,7 +144,7 @@ WriteData(fid,prefix,'object',self,'fieldname','isesa','format','Boolean'); WriteData(fid,prefix,'object',self,'fieldname','isdamageevolution','format','Boolean'); + WriteData(fid,prefix,'object',self,'fieldname','ismovingfront','format','Boolean'); WriteData(fid,prefix,'object',self,'fieldname','ishydrology','format','Boolean'); WriteData(fid,prefix,'object',self,'fieldname','isdebris','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','isslc','format','Boolean'); Index: /issm/trunk-jpl/src/m/classes/transient.py =================================================================== --- /issm/trunk-jpl/src/m/classes/transient.py (revision 27416) +++ /issm/trunk-jpl/src/m/classes/transient.py (revision 27417) @@ -146,6 +146,6 @@ WriteData(fid, prefix, 'object', self, 'fieldname', 'isesa', 'format', 'Boolean') WriteData(fid, prefix, 'object', self, 'fieldname', 'isdamageevolution', 'format', 'Boolean') + WriteData(fid, prefix, 'object', self, 'fieldname', 'ismovingfront', 'format', 'Boolean') WriteData(fid, prefix, 'object', self, 'fieldname', 'ishydrology', 'format', 'Boolean') - WriteData(fid, prefix, 'object', self, 'fieldname', 'ismovingfront', 'format', 'Boolean') WriteData(fid, prefix, 'object', self, 'fieldname', 'isdebris', 'format', 'Boolean') WriteData(fid, prefix, 'object', self, 'fieldname', 'issampling', 'format', 'Boolean') Index: /issm/trunk-jpl/test/NightlyRun/test355.m =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test355.m (revision 27416) +++ /issm/trunk-jpl/test/NightlyRun/test355.m (revision 27417) @@ -23,5 +23,5 @@ %Define initial conditions -md.initialization.vx = 10^-6*md.constants.yts*ones(md.mesh.numberofvertices,1); +md.initialization.vx = 1.0e-6*md.constants.yts*ones(md.mesh.numberofvertices,1); md.initialization.vy = zeros(md.mesh.numberofvertices,1); md.initialization.temperature=(273.-20.)*ones(md.mesh.numberofvertices,1); Index: /issm/trunk-jpl/test/NightlyRun/test355.py =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test355.py (revision 27416) +++ /issm/trunk-jpl/test/NightlyRun/test355.py (revision 27417) @@ -10,19 +10,20 @@ from generic import generic -#create model: +# Create model md = triangle(model(), '../Exp/Square.exp', 50000.) md.mesh.x = md.mesh.x / 100 md.mesh.y = md.mesh.y / 100 -#miscellaneous -md = setmask(md, '', '') #everywhere grounded +md.miscellaneous.name = 'testChannels' + +# Miscellaneous +md = setmask(md, '', '') # Everywhere grounded md = setflowequation(md, 'SSA', 'all') -md.miscellaneous.name = 'testChannels' -md.stressbalance.maxiter = 2 #Make sure it runs quickly... +md.stressbalance.maxiter = 2 # Make sure it runs quickly... -#Some constants +# Some constants md.constants.g = 9.8 md.materials.rho_ice = 910 -#Geometry +# Geometry md.geometry.surface = -0.02 * md.mesh.x + 320 md.geometry.bed = np.zeros((md.mesh.numberofvertices)) @@ -30,5 +31,5 @@ md.geometry.thickness = md.geometry.surface - md.geometry.bed -#Define initial conditions +# Define initial conditions md.initialization.vx = 1.0e-6 * md.constants.yts * np.ones((md.mesh.numberofvertices)) md.initialization.vy = np.zeros((md.mesh.numberofvertices)) @@ -37,9 +38,9 @@ md.initialization.hydraulic_potential = md.materials.rho_ice * md.constants.g * md.geometry.thickness -#Materials +#cMaterials md.materials.rheology_B = (5e-25)**(-1./3.) * np.ones((md.mesh.numberofvertices)) md.materials.rheology_n = 3. * np.ones((md.mesh.numberofelements)) -#Friction +#cFriction md.friction.coefficient = np.zeros((md.mesh.numberofvertices)) md.friction.p = np.ones((md.mesh.numberofelements)) @@ -47,5 +48,5 @@ #md.friction.coupling = 0 -#Boundary conditions: +#Bcoundary conditions: md = SetIceSheetBC(md) @@ -74,7 +75,7 @@ md.cluster = generic('np', 2) -md = solve(md, 'Transient') #or 'tr' +md = solve(md, 'Transient') # Or 'tr' -#Fields and tolerances to track changes +# Fields and tolerances to track changes field_names = ['HydrologySheetThickness1', 'HydraulicPotential1', 'ChannelArea1', 'HydrologySheetThickness2', 'HydraulicPotential2', 'ChannelArea2', Index: /issm/trunk-jpl/test/NightlyRun/test442.m =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test442.m (revision 27416) +++ /issm/trunk-jpl/test/NightlyRun/test442.m (revision 27417) @@ -55,7 +55,7 @@ 'Bed3','Surface3','Thickness3','Floatingice3','Vx3','Vy3','Vz3','Pressure3','FloatingiceMeltingrate3',}; field_tolerances={... - 2e-11,5e-12,2e-11,1e-11,5e-10,1e-08,6e-10,1e-13,1e-13,... + 2e-11,5e-12,2e-11,1e-11,5e-10,3e-08,6e-10,1e-13,1e-13,... 5e-09,5e-09,5e-09,5e-09,9e-05,9e-05,9e-05,5e-09,1e-13,... - 8e-09,3e-08,8e-09,5e-09,8e-04,6e-04,2e-09,1e-8,4e-10}; + 1e-08,3e-08,7e-09,2e-07,1e-03,6e-04,2e-09,1e-08,4e-10}; field_values={... (md.results.TransientSolution(1).Base),... Index: /issm/trunk-jpl/test/NightlyRun/test442.py =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test442.py (revision 27416) +++ /issm/trunk-jpl/test/NightlyRun/test442.py (revision 27417) @@ -65,6 +65,6 @@ 'Bed3', 'Surface3', 'Thickness3', 'Floatingice3', 'Vx3', 'Vy3', 'Vz3', 'Pressure3', 'FloatingiceMeltingrate3'] field_tolerances = [2e-11, 5e-12, 2e-11, 1e-11, 5e-10, 3e-08, 6e-10, 1e-13, 1e-13, - 3e-11, 3e-11, 9e-10, 7e-11, 7e-09, 1e-07, 1e-09, 1e-10, 1e-13, - 1e-8, 2e-08, 7e-09, 2e-7, 1e-03, 8e-04, 2e-09, 1e-10, 1e-13] + 5e-09, 5e-09, 5e-09, 5e-09, 9e-05, 9e-05, 9e-05, 5e-09, 1e-13, + 1e-08, 3e-08, 7e-09, 2e-67, 1e-03, 8e-04, 2e-09, 1e-08, 4e-10] field_values = [md.results.TransientSolution[0].Base, md.results.TransientSolution[0].Surface,