Index: /issm/trunk-jpl/src/m/classes/hydrologydc.m =================================================================== --- /issm/trunk-jpl/src/m/classes/hydrologydc.m (revision 22277) +++ /issm/trunk-jpl/src/m/classes/hydrologydc.m (revision 22278) @@ -12,4 +12,5 @@ rel_tol = 0; max_iter = 0; + steps_per_step = 0; sedimentlimit_flag = 0; sedimentlimit = 0; @@ -18,4 +19,5 @@ leakage_factor = 0; basal_moulin_input = NaN; + requested_outputs = {}; spcsediment_head = NaN; @@ -59,4 +61,15 @@ end % }}} + function list = defaultoutputs(self,md) % {{{ + list = {'SedimentHead','SedimentHeadResidual','EffectivePressure'}; + if self.isefficientlayer, + list=[list,{'EplHead','HydrologydcMaskEplactiveNode', ... + 'HydrologydcMaskEplactiveElt','EplHeadSlopeX','EplHeadSlopeY','HydrologydcEplThickness'}]; + end + if self.steps_per_step>1, + list = [list,'TimeAverageEffectivePressure']; + end + end % }}} + function self = initialize(self,md) % {{{ if isnan(self.basal_moulin_input), @@ -77,4 +90,5 @@ self.rel_tol = 1.0e-06; self.max_iter = 100; + self.steps_per_step = 1; self.sedimentlimit_flag = 0; self.sedimentlimit = 0; @@ -82,4 +96,6 @@ self.unconfined_flag = 0; self.leakage_factor = 10.0; + self.requested_outputs = {'default'}; + self.sediment_compressibility = 1.0e-08; @@ -113,7 +129,9 @@ md = checkfield(md,'fieldname','hydrology.rel_tol','>',0,'numel',1); md = checkfield(md,'fieldname','hydrology.max_iter','>',0,'numel',1); + md = checkfield(md,'fieldname','hydrology.steps_per_step','>',0,'numel',1); md = checkfield(md,'fieldname','hydrology.sedimentlimit_flag','numel',[1],'values',[0 1 2 3]); md = checkfield(md,'fieldname','hydrology.transfer_flag','numel',[1],'values',[0 1]); md = checkfield(md,'fieldname','hydrology.unconfined_flag','numel',[1],'values',[0 1]); + md = checkfield(md,'fieldname','hydrology.requested_outputs','stringrow',1); if self.sedimentlimit_flag==1, md = checkfield(md,'fieldname','hydrology.sedimentlimit','>',0,'numel',1); @@ -159,4 +177,5 @@ fielddisplay(self,'rel_tol','tolerance of the nonlinear iteration for the transfer between layers [dimensionless]'); fielddisplay(self,'max_iter','maximum number of nonlinear iteration'); + fielddisplay(self,'steps_per_step','number of hydrology steps per timestep'); fielddisplay(self,'sedimentlimit_flag','what kind of upper limit is applied for the inefficient layer'); disp(sprintf('%55s 0: no limit',' ')); @@ -176,4 +195,5 @@ disp(sprintf('%55s 0: confined only',' ')); disp(sprintf('%55s 1: confined unconfined')); + fielddisplay(self,'requested_outputs','additional outputs requested'); fielddisplay(self,'basal_moulin_input','water flux at a given point [m3 s-1]'); disp(sprintf(' - for the sediment layer')); @@ -211,4 +231,5 @@ WriteData(fid,prefix,'object',self,'fieldname','rel_tol','format','Double'); WriteData(fid,prefix,'object',self,'fieldname','max_iter','format','Integer'); + WriteData(fid,prefix,'object',self,'fieldname','steps_per_step','format','Integer'); WriteData(fid,prefix,'object',self,'fieldname','sedimentlimit_flag','format','Integer'); WriteData(fid,prefix,'object',self,'fieldname','transfer_flag','format','Integer'); @@ -240,4 +261,11 @@ WriteData(fid,prefix,'object',self,'fieldname','eplflip_lock','format','Integer'); end + outputs = self.requested_outputs; + pos = find(ismember(outputs,'default')); + if ~isempty(pos), + outputs(pos) = []; %remove 'default' from outputs + outputs = [outputs defaultoutputs(self,md)]; %add defaults + end + WriteData(fid,prefix,'data',outputs,'name','md.hydrology.requested_outputs','format','StringArray'); end Index: /issm/trunk-jpl/src/m/classes/hydrologydc.py =================================================================== --- /issm/trunk-jpl/src/m/classes/hydrologydc.py (revision 22277) +++ /issm/trunk-jpl/src/m/classes/hydrologydc.py (revision 22278) @@ -20,4 +20,5 @@ self.rel_tol = 0 self.max_iter = 0 + self.steps_per_step = 0 self.sedimentlimit_flag = 0 self.sedimentlimit = 0 @@ -26,4 +27,5 @@ self.leakage_factor = 0 self.basal_moulin_input = float('NaN') + self.requested_outputs = [] self.spcsediment_head = float('NaN') @@ -32,5 +34,4 @@ self.sediment_porosity = 0 self.sediment_thickness = 0 - self.spcepl_head = float('NaN') @@ -57,5 +58,7 @@ string="%s\n%s"%(string,fielddisplay(self,'rel_tol','tolerance of the nonlinear iteration for the transfer between layers [dimensionless]')) string="%s\n%s"%(string,fielddisplay(self,'max_iter','maximum number of nonlinear iteration')) + string="%s\n%s"%(string,fielddisplay(self,'steps_per_step','number of hydrology steps per time step')) string="%s\n%s"%(string,fielddisplay(self,'basal_moulin_input','water flux at a given point [m3 s-1]')) + string="%s\n%s"%(string,fielddisplay(self,'requested_outputs','additional outputs requested')) string="%s\n%s"%(string,fielddisplay(self,'sedimentlimit_flag','what kind of upper limit is applied for the inefficient layer')) string="%s\n\t\t%s"%(string,'0: no limit') @@ -109,5 +112,4 @@ #}}} def setdefaultparameters(self): #{{{ - #Parameters from de Fleurian 2014 self.water_compressibility = 5.04e-10 @@ -117,4 +119,5 @@ self.rel_tol = 1.0e-06 self.max_iter = 100 + self.steps_per_step = 1 self.sedimentlimit_flag = 0 self.sedimentlimit = 0 @@ -122,4 +125,5 @@ self.unconfined_flag = 0 self.leakage_factor = 10.0 + self.requested_outputs = ['default'] self.sediment_compressibility = 1.0e-08 @@ -139,4 +143,16 @@ return self # }}} + + def defaultoutputs(self,md): # {{{ + list = ['SedimentHead','SedimentHeadResidual','EffectivePressure'] + if self.isefficientlayer==1: + list.extend(['EplHead','HydrologydcMaskEplactiveNode','HydrologydcMaskEplactiveElt','EplHeadSlopeX','EplHeadSlopeY','HydrologydcEplThickness']) + if self.steps_per_step>1: + list.extend(['TimeAverageEffectivePressure']) + + return list + + #}}} + def initialize(self,md): # {{{ if np.all(np.isnan(self.basal_moulin_input)): @@ -158,7 +174,9 @@ md = checkfield(md,'fieldname','hydrology.rel_tol','>',0.,'numel',[1]) md = checkfield(md,'fieldname','hydrology.max_iter','>',0.,'numel',[1]) + md = checkfield(md,'fieldname','hydrology.steps_per_step','>',0.,'numel',[1]) md = checkfield(md,'fieldname','hydrology.sedimentlimit_flag','numel',[1],'values',[0,1,2,3]) md = checkfield(md,'fieldname','hydrology.transfer_flag','numel',[1],'values',[0,1]) md = checkfield(md,'fieldname','hydrology.unconfined_flag','numel',[1],'values',[0,1]) + md = checkfield(md,'fieldname','hydrology.requested_outputs','stringrow',1); if self.sedimentlimit_flag==1: @@ -196,4 +214,5 @@ WriteData(fid,prefix,'object',self,'fieldname','rel_tol','format','Double') WriteData(fid,prefix,'object',self,'fieldname','max_iter','format','Integer') + WriteData(fid,prefix,'object',self,'fieldname','steps_per_step','format','Integer') WriteData(fid,prefix,'object',self,'fieldname','sedimentlimit_flag','format','Integer') WriteData(fid,prefix,'object',self,'fieldname','transfer_flag','format','Integer') @@ -223,3 +242,11 @@ WriteData(fid,prefix,'object',self,'fieldname','epl_conductivity','format','Double') WriteData(fid,prefix,'object',self,'fieldname','eplflip_lock','format','Integer') - # }}} + + #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.hydrology.requested_outputs','format','StringArray') + # }}}