Index: /issm/trunk-jpl/src/m/classes/hydrologydc.py =================================================================== --- /issm/trunk-jpl/src/m/classes/hydrologydc.py (revision 18585) +++ /issm/trunk-jpl/src/m/classes/hydrologydc.py (revision 18586) @@ -1,2 +1,3 @@ +import numpy from fielddisplay import fielddisplay from EnumDefinitions import * @@ -5,12 +6,12 @@ class hydrologydc(object): - """ - Hydrologydc class definition - - Usage: - hydrologydc=hydrologydc(); - """ + """ + Hydrologydc class definition - def __init__(self): # {{{ + Usage: + hydrologydc=hydrologydc(); + """ + + def __init__(self): # {{{ self.water_compressibility = 0 self.isefficientlayer = 0 @@ -23,5 +24,5 @@ self.transfer_flag = 0 self.leakage_factor = 0 - self.basal_moulin_input = float('NaN') + self.basal_moulin_input = float('NaN') self.spcsediment_head = float('NaN') @@ -39,152 +40,146 @@ self.epl_max_thickness = 0 self.epl_conductivity = 0 - + #set defaults self.setdefaultparameters() + #}}} - #}}} - def __repr__(self): # {{{ - string=' hydrology Dual Porous Continuum Equivalent parameters:' - string=' - general parameters' - string="%s\n%s"%(string,fielddisplay(self,'water_compressibility','compressibility of water [Pa^-1]')) - string="%s\n%s"%(string,fielddisplay(self,'isefficientlayer','do we use an efficient drainage system [1: true 0: false]')) - string="%s\n%s"%(string,fielddisplay(self,'penalty_factor','exponent of the value used in the penalisation method [dimensionless]')) - string="%s\n%s"%(string,fielddisplay(self,'penalty_lock','stabilize unstable constraints that keep zigzagging after n iteration (default is 0, no stabilization)')) - 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,'sedimentlimit_flag','what kind of upper limit is applied for the inefficient layer')) - string='%55s 0: no limit',' ' - string='%55s 1: user defined: %s',' ','sedimentlimit' - string='%55s 2: hydrostatic pressure',' ' - string='%55s 3: normal stress',' ' + string=' hydrology Dual Porous Continuum Equivalent parameters:' + string=' - general parameters' + string="%s\n%s"%(string,fielddisplay(self,'water_compressibility','compressibility of water [Pa^-1]')) + string="%s\n%s"%(string,fielddisplay(self,'isefficientlayer','do we use an efficient drainage system [1: true 0: false]')) + string="%s\n%s"%(string,fielddisplay(self,'penalty_factor','exponent of the value used in the penalisation method [dimensionless]')) + string="%s\n%s"%(string,fielddisplay(self,'penalty_lock','stabilize unstable constraints that keep zigzagging after n iteration (default is 0, no stabilization)')) + 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,'sedimentlimit_flag','what kind of upper limit is applied for the inefficient layer')) + string='%55s 0: no limit',' ' + string='%55s 1: user defined: %s',' ','sedimentlimit' + string='%55s 2: hydrostatic pressure',' ' + string='%55s 3: normal stress',' ' - if self.sedimentlimit_flag==1: - string="%s\n%s"%(string,fielddisplay(self,'sedimentlimit','user defined upper limit for the inefficient layer [m]')) - 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,'transfer_flag','what kind of transfer method is applied between the layers')) - string='%55s 0: no transfer',' ' - string='%55s 1: constant leakage factor: %s',' ','leakage_factor' - - if self.transfer_flag is 1: - string="%s\n%s"%(string,fielddisplay(self,'leakage_factor','user defined leakage factor [m]')) - string=' - for the sediment layer' - string="%s\n%s"%(string,fielddisplay(self,'spcsediment_head','sediment water head constraints (NaN means no constraint) [m above MSL]')) - string="%s\n%s"%(string,fielddisplay(self,'sediment_compressibility','sediment compressibility [Pa^-1]')) - string="%s\n%s"%(string,fielddisplay(self,'sediment_porosity','sediment [dimensionless]')) - string="%s\n%s"%(string,fielddisplay(self,'sediment_thickness','sediment thickness [m]')) - string="%s\n%s"%(string,fielddisplay(self,'sediment_transmitivity','sediment transmitivity [m^2/s]')) - - if self.isefficientlayer==1: - string=' - for the epl layer' - string="%s\n%s"%(string,fielddisplay(self,'spcepl_head','epl water head constraints (NaN means no constraint) [m above MSL]')) - string="%s\n%s"%(string,fielddisplay(self,'mask_eplactive_node','active (1) or not (0) EPL')) - string="%s\n%s"%(string,fielddisplay(self,'epl_compressibility','epl compressibility [Pa^-1]')) - string="%s\n%s"%(string,fielddisplay(self,'epl_porosity','epl [dimensionless]')) - string="%s\n%s"%(string,fielddisplay(self,'epl_initial_thickness','epl initial thickness [m]')) - string="%s\n%s"%(string,fielddisplay(self,'epl_conductivity','epl conductivity [m^2/s]')) - - #}}} - def setdefaultparameters(self): #{{{ - - #Parameters from de Fleurian 2014 - self.water_compressibility = 5.04e-10 - self.isefficientlayer = 1 - self.penalty_factor = 3 - self.rel_tol = 1.0e-06 - self.max_iter = 100 - self.sedimentlimit_flag = 0 - self.sedimentlimit = 0 - self.transfer_flag = 0 - self.leakage_factor = 10.0 - - self.sediment_compressibility = 1.0e-08 - self.sediment_porosity = 0.4 - self.sediment_thickness = 20.0 - self.sediment_transmitivity = 8.0e-04 - - self.epl_compressibility = 1.0e-08 - self.epl_porosity = 0.4 - self.epl_initial_thickness = 1.0 - self.epl_max_thickness = 5.0 - self.epl_conductivity = 8.0e-02 - - return self - # }}} - - def initialize(self): # {{{ - if numpy.all(numpy.isnan(self.basal_moulin_input): - self.basal_moulin_input=numpy.zeros(md.mesh.numberofvertices,1) - print" no hydrology.basal_moulin_input specified: values set as zero" + if self.sedimentlimit_flag==1: + string="%s\n%s"%(string,fielddisplay(self,'sedimentlimit','user defined upper limit for the inefficient layer [m]')) + 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,'transfer_flag','what kind of transfer method is applied between the layers')) + string='%55s 0: no transfer',' ' + string='%55s 1: constant leakage factor: %s',' ','leakage_factor' + + if self.transfer_flag is 1: + string="%s\n%s"%(string,fielddisplay(self,'leakage_factor','user defined leakage factor [m]')) + string=' - for the sediment layer' + string="%s\n%s"%(string,fielddisplay(self,'spcsediment_head','sediment water head constraints (NaN means no constraint) [m above MSL]')) + string="%s\n%s"%(string,fielddisplay(self,'sediment_compressibility','sediment compressibility [Pa^-1]')) + string="%s\n%s"%(string,fielddisplay(self,'sediment_porosity','sediment [dimensionless]')) + string="%s\n%s"%(string,fielddisplay(self,'sediment_thickness','sediment thickness [m]')) + string="%s\n%s"%(string,fielddisplay(self,'sediment_transmitivity','sediment transmitivity [m^2/s]')) - return self + if self.isefficientlayer==1: + string=' - for the epl layer' + string="%s\n%s"%(string,fielddisplay(self,'spcepl_head','epl water head constraints (NaN means no constraint) [m above MSL]')) + string="%s\n%s"%(string,fielddisplay(self,'mask_eplactive_node','active (1) or not (0) EPL')) + string="%s\n%s"%(string,fielddisplay(self,'epl_compressibility','epl compressibility [Pa^-1]')) + string="%s\n%s"%(string,fielddisplay(self,'epl_porosity','epl [dimensionless]')) + string="%s\n%s"%(string,fielddisplay(self,'epl_initial_thickness','epl initial thickness [m]')) + string="%s\n%s"%(string,fielddisplay(self,'epl_conductivity','epl conductivity [m^2/s]')) + #}}} + def setdefaultparameters(self): #{{{ - # }}} - def checkconsistency(self,md,solution,analyses): #{{{ - - #Early return - if HydrologyDCInefficientAnalysisEnum() not in analyses: - return md + #Parameters from de Fleurian 2014 + self.water_compressibility = 5.04e-10 + self.isefficientlayer = 1 + self.penalty_factor = 3 + self.rel_tol = 1.0e-06 + self.max_iter = 100 + self.sedimentlimit_flag = 0 + self.sedimentlimit = 0 + self.transfer_flag = 0 + self.leakage_factor = 10.0 - md = checkfield(md,'fieldname','hydrology.water_compressibility','>',0,'numel',1) - md = checkfield(md,'fieldname','hydrology.isefficientlayer','numel',[1],'values',[0 1]) - md = checkfield(md,'fieldname','hydrology.penalty_factor','>',0,'numel',1) - 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.sedimentlimit_flag','numel',[1],'values',[0 1 2 3]) - md = checkfield(md,'fieldname','hydrology.transfer_flag','numel',[1],'values',[0 1]) - - if self.sedimentlimit_flag==1: - md = checkfield(md,'fieldname','hydrology.sedimentlimit','>',0,'numel',1) - - if self.transfer_flag==1: - md = checkfield(md,'fieldname','hydrology.leakage_factor','>',0,'numel',1) + self.sediment_compressibility = 1.0e-08 + self.sediment_porosity = 0.4 + self.sediment_thickness = 20.0 + self.sediment_transmitivity = 8.0e-04 - md = checkfield(md,'fieldname','hydrology.basal_moulin_input','NaN',1,'forcing',1) - md = checkfield(md,'fieldname','hydrology.spcsediment_head','forcing',1) - md = checkfield(md,'fieldname','hydrology.sediment_compressibility','>',0,'numel',1) - md = checkfield(md,'fieldname','hydrology.sediment_porosity','>',0,'numel',1) - md = checkfield(md,'fieldname','hydrology.sediment_thickness','>',0,'numel',1) - md = checkfield(md,'fieldname','hydrology.sediment_transmitivity','>=',0,'size',[md.mesh.numberofvertices 1]) - if self.isefficientlayer==1: - md = checkfield(md,'fieldname','hydrology.spcepl_head','forcing',1) - md = checkfield(md,'fieldname','hydrology.mask_eplactive_node','size',[md.mesh.numberofvertices 1],'values',[0 1]) - md = checkfield(md,'fieldname','hydrology.epl_compressibility','>',0,'numel',1) - md = checkfield(md,'fieldname','hydrology.epl_porosity','>',0,'numel',1) - md = checkfield(md,'fieldname','hydrology.epl_initial_thickness','>',0,'numel',1) - md = checkfield(md,'fieldname','hydrology.epl_conductivity','>',0,'numel',1) + self.epl_compressibility = 1.0e-08 + self.epl_porosity = 0.4 + self.epl_initial_thickness = 1.0 + self.epl_max_thickness = 5.0 + self.epl_conductivity = 8.0e-02 - # }}} - def marshall(self,md,fid): #{{{ - WriteData(fid,'enum',HydrologyModelEnum(),'data',HydrologydcEnum(),'format','Integer') - WriteData(fid,'object',self,'fieldname','water_compressibility','format','Double') - WriteData(fid,'object',self,'fieldname','isefficientlayer','format','Boolean') - WriteData(fid,'object',self,'fieldname','penalty_factor','format','Double') - WriteData(fid,'object',self,'fieldname','penalty_lock','format','Integer') - WriteData(fid,'object',self,'fieldname','rel_tol','format','Double') - WriteData(fid,'object',self,'fieldname','max_iter','format','Integer') - WriteData(fid,'object',self,'fieldname','sedimentlimit_flag','format','Integer') - WriteData(fid,'object',self,'fieldname','transfer_flag','format','Integer') + return self + # }}} - if self.sedimentlimit_flag==1: - WriteData(fid,'object',self,'fieldname','sedimentlimit','format','Double') + def initialize(self): # {{{ + if numpy.all(numpy.isnan(self.basal_moulin_input)): + self.basal_moulin_input=numpy.zeros(md.mesh.numberofvertices,1) + print" no hydrology.basal_moulin_input specified: values set as zero" - if self.transfer_flag==1: - WriteData(fid,'object',self,'fieldname','leakage_factor','format','Double') + return self + # }}} + def checkconsistency(self,md,solution,analyses): #{{{ - WriteData(fid,'object',self,'fieldname','basal_moulin_input','format','DoubleMat','mattype',1,'forcinglength',md.mesh.numberofvertices+1) - WriteData(fid,'object',self,'fieldname','spcsediment_head','format','DoubleMat','mattype',1,'forcinglength',md.mesh.numberofvertices+1) - WriteData(fid,'object',self,'fieldname','sediment_compressibility','format','Double') - WriteData(fid,'object',self,'fieldname','sediment_porosity','format','Double') - WriteData(fid,'object',self,'fieldname','sediment_thickness','format','Double') - WriteData(fid,'object',self,'fieldname','sediment_transmitivity','format','DoubleMat','mattype',1) + #Early return + if HydrologyDCInefficientAnalysisEnum() not in analyses: + return md - if self.isefficientlayer==1: - WriteData(fid,'object',self,'fieldname','spcepl_head','format','DoubleMat','mattype',1,'forcinglength',md.mesh.numberofvertices+1) - WriteData(fid,'object',self,'fieldname','mask_eplactive_node','format','DoubleMat','mattype',1) - WriteData(fid,'object',self,'fieldname','epl_compressibility','format','Double') - WriteData(fid,'object',self,'fieldname','epl_porosity','format','Double') - WriteData(fid,'object',self,'fieldname','epl_initial_thickness','format','Double') - WriteData(fid,'object',self,'fieldname','epl_conductivity','format','Double') -# }}} + md = checkfield(md,'fieldname','hydrology.water_compressibility','>',0,'numel',1) + md = checkfield(md,'fieldname','hydrology.isefficientlayer','numel',[1],'values',[0,1]) + md = checkfield(md,'fieldname','hydrology.penalty_factor','>',0,'numel',1) + 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.sedimentlimit_flag','numel',[1],'values',[0,1,2,3]) + md = checkfield(md,'fieldname','hydrology.transfer_flag','numel',[1],'values',[0,1]) + if self.sedimentlimit_flag==1: + md = checkfield(md,'fieldname','hydrology.sedimentlimit','>',0,'numel',1) + + if self.transfer_flag==1: + md = checkfield(md,'fieldname','hydrology.leakage_factor','>',0,'numel',1) + + md = checkfield(md,'fieldname','hydrology.basal_moulin_input','NaN',1,'forcing',1) + md = checkfield(md,'fieldname','hydrology.spcsediment_head','forcing',1) + md = checkfield(md,'fieldname','hydrology.sediment_compressibility','>',0,'numel',1) + md = checkfield(md,'fieldname','hydrology.sediment_porosity','>',0,'numel',1) + md = checkfield(md,'fieldname','hydrology.sediment_thickness','>',0,'numel',1) + md = checkfield(md,'fieldname','hydrology.sediment_transmitivity','>=',0,'size',[md.mesh.numberofvertices,1]) + if self.isefficientlayer==1: + md = checkfield(md,'fieldname','hydrology.spcepl_head','forcing',1) + md = checkfield(md,'fieldname','hydrology.mask_eplactive_node','size',[md.mesh.numberofvertices,1],'values',[0,1]) + md = checkfield(md,'fieldname','hydrology.epl_compressibility','>',0,'numel',1) + md = checkfield(md,'fieldname','hydrology.epl_porosity','>',0,'numel',1) + md = checkfield(md,'fieldname','hydrology.epl_initial_thickness','>',0,'numel',1) + md = checkfield(md,'fieldname','hydrology.epl_conductivity','>',0,'numel',1) + # }}} + def marshall(self,md,fid): #{{{ + WriteData(fid,'enum',HydrologyModelEnum(),'data',HydrologydcEnum(),'format','Integer') + WriteData(fid,'object',self,'fieldname','water_compressibility','format','Double') + WriteData(fid,'object',self,'fieldname','isefficientlayer','format','Boolean') + WriteData(fid,'object',self,'fieldname','penalty_factor','format','Double') + WriteData(fid,'object',self,'fieldname','penalty_lock','format','Integer') + WriteData(fid,'object',self,'fieldname','rel_tol','format','Double') + WriteData(fid,'object',self,'fieldname','max_iter','format','Integer') + WriteData(fid,'object',self,'fieldname','sedimentlimit_flag','format','Integer') + WriteData(fid,'object',self,'fieldname','transfer_flag','format','Integer') + if self.sedimentlimit_flag==1: + WriteData(fid,'object',self,'fieldname','sedimentlimit','format','Double') + + if self.transfer_flag==1: + WriteData(fid,'object',self,'fieldname','leakage_factor','format','Double') + + WriteData(fid,'object',self,'fieldname','basal_moulin_input','format','DoubleMat','mattype',1,'forcinglength',md.mesh.numberofvertices+1) + WriteData(fid,'object',self,'fieldname','spcsediment_head','format','DoubleMat','mattype',1,'forcinglength',md.mesh.numberofvertices+1) + WriteData(fid,'object',self,'fieldname','sediment_compressibility','format','Double') + WriteData(fid,'object',self,'fieldname','sediment_porosity','format','Double') + WriteData(fid,'object',self,'fieldname','sediment_thickness','format','Double') + WriteData(fid,'object',self,'fieldname','sediment_transmitivity','format','DoubleMat','mattype',1) + + if self.isefficientlayer==1: + WriteData(fid,'object',self,'fieldname','spcepl_head','format','DoubleMat','mattype',1,'forcinglength',md.mesh.numberofvertices+1) + WriteData(fid,'object',self,'fieldname','mask_eplactive_node','format','DoubleMat','mattype',1) + WriteData(fid,'object',self,'fieldname','epl_compressibility','format','Double') + WriteData(fid,'object',self,'fieldname','epl_porosity','format','Double') + WriteData(fid,'object',self,'fieldname','epl_initial_thickness','format','Double') + WriteData(fid,'object',self,'fieldname','epl_conductivity','format','Double') + # }}}