Changeset 18586


Ignore:
Timestamp:
10/07/14 09:24:40 (10 years ago)
Author:
Mathieu Morlighem
Message:

BUG: fixed python version of hydrologydc

File:
1 edited

Legend:

Unmodified
Added
Removed

  • issm/trunk-jpl/src/m/classes/hydrologydc.py

    r18585 r18586  
     1import numpy
    12from fielddisplay import fielddisplay
    23from EnumDefinitions import *
    … …  
    56
    67class hydrologydc(object):
    7     """
    8     Hydrologydc class definition
    9    
    10     Usage:
    11     hydrologydc=hydrologydc();
    12     """
     8        """
     9        Hydrologydc class definition
    1310
    14     def __init__(self): # {{{
     11        Usage:
     12                hydrologydc=hydrologydc();
     13        """
     14
     15        def __init__(self): # {{{
    1516                self.water_compressibility    = 0
    1617                self.isefficientlayer         = 0
    … …  
    2324                self.transfer_flag            = 0
    2425                self.leakage_factor           = 0
    25                 self.basal_moulin_input       = float('NaN')
     26                self.basal_moulin_input       = float('NaN')
    2627
    2728                self.spcsediment_head         = float('NaN')
    … …  
    3940                self.epl_max_thickness        = 0
    4041                self.epl_conductivity         = 0
    41                
     42                               
    4243                #set defaults
    4344                self.setdefaultparameters()
     45        #}}}
    4446
    45                 #}}}
    46                
    4747        def __repr__(self): # {{{
    48             string='   hydrology Dual Porous Continuum Equivalent parameters:'
    49             string='   - general parameters'
    50             string="%s\n%s"%(string,fielddisplay(self,'water_compressibility','compressibility of water [Pa^-1]'))
    51             string="%s\n%s"%(string,fielddisplay(self,'isefficientlayer','do we use an efficient drainage system [1: true 0: false]'))
    52             string="%s\n%s"%(string,fielddisplay(self,'penalty_factor','exponent of the value used in the penalisation method [dimensionless]'))
    53             string="%s\n%s"%(string,fielddisplay(self,'penalty_lock','stabilize unstable constraints that keep zigzagging after n iteration (default is 0, no stabilization)'))
    54             string="%s\n%s"%(string,fielddisplay(self,'rel_tol','tolerance of the nonlinear iteration for the transfer between layers [dimensionless]'))
    55             string="%s\n%s"%(string,fielddisplay(self,'max_iter','maximum number of nonlinear iteration'))
    56             string="%s\n%s"%(string,fielddisplay(self,'sedimentlimit_flag','what kind of upper limit is applied for the inefficient layer'))
    57             string='%55s  0: no limit',' '
    58             string='%55s  1: user defined: %s',' ','sedimentlimit'
    59             string='%55s  2: hydrostatic pressure',' '
    60             string='%55s  3: normal stress',' '
     48                string='   hydrology Dual Porous Continuum Equivalent parameters:'
     49                string='   - general parameters'
     50                string="%s\n%s"%(string,fielddisplay(self,'water_compressibility','compressibility of water [Pa^-1]'))
     51                string="%s\n%s"%(string,fielddisplay(self,'isefficientlayer','do we use an efficient drainage system [1: true 0: false]'))
     52                string="%s\n%s"%(string,fielddisplay(self,'penalty_factor','exponent of the value used in the penalisation method [dimensionless]'))
     53                string="%s\n%s"%(string,fielddisplay(self,'penalty_lock','stabilize unstable constraints that keep zigzagging after n iteration (default is 0, no stabilization)'))
     54                string="%s\n%s"%(string,fielddisplay(self,'rel_tol','tolerance of the nonlinear iteration for the transfer between layers [dimensionless]'))
     55                string="%s\n%s"%(string,fielddisplay(self,'max_iter','maximum number of nonlinear iteration'))
     56                string="%s\n%s"%(string,fielddisplay(self,'sedimentlimit_flag','what kind of upper limit is applied for the inefficient layer'))
     57                string='%55s  0: no limit',' '
     58                string='%55s  1: user defined: %s',' ','sedimentlimit'
     59                string='%55s  2: hydrostatic pressure',' '
     60                string='%55s  3: normal stress',' '
    6161
    62             if self.sedimentlimit_flag==1:
    63                 string="%s\n%s"%(string,fielddisplay(self,'sedimentlimit','user defined upper limit for the inefficient layer [m]'))
    64                 string="%s\n%s"%(string,fielddisplay(self,'basal_moulin_input','water flux at a given point [m3 s-1]'))
    65                 string="%s\n%s"%(string,fielddisplay(self,'transfer_flag','what kind of transfer method is applied between the layers'))
    66                 string='%55s  0: no transfer',' '
    67                 string='%55s  1: constant leakage factor: %s',' ','leakage_factor'
    68                
    69             if self.transfer_flag is 1:
    70                 string="%s\n%s"%(string,fielddisplay(self,'leakage_factor','user defined leakage factor [m]'))
    71                 string='   - for the sediment layer'
    72                 string="%s\n%s"%(string,fielddisplay(self,'spcsediment_head','sediment water head constraints (NaN means no constraint) [m above MSL]'))
    73                 string="%s\n%s"%(string,fielddisplay(self,'sediment_compressibility','sediment compressibility [Pa^-1]'))
    74                 string="%s\n%s"%(string,fielddisplay(self,'sediment_porosity','sediment [dimensionless]'))
    75                 string="%s\n%s"%(string,fielddisplay(self,'sediment_thickness','sediment thickness [m]'))
    76                 string="%s\n%s"%(string,fielddisplay(self,'sediment_transmitivity','sediment transmitivity [m^2/s]'))
    77            
    78             if self.isefficientlayer==1:
    79                 string='   - for the epl layer'
    80                 string="%s\n%s"%(string,fielddisplay(self,'spcepl_head','epl water head constraints (NaN means no constraint) [m above MSL]'))
    81                 string="%s\n%s"%(string,fielddisplay(self,'mask_eplactive_node','active (1) or not (0) EPL'))
    82                 string="%s\n%s"%(string,fielddisplay(self,'epl_compressibility','epl compressibility [Pa^-1]'))
    83                 string="%s\n%s"%(string,fielddisplay(self,'epl_porosity','epl [dimensionless]'))
    84                 string="%s\n%s"%(string,fielddisplay(self,'epl_initial_thickness','epl initial thickness [m]'))
    85                 string="%s\n%s"%(string,fielddisplay(self,'epl_conductivity','epl conductivity [m^2/s]'))
    86            
    87             #}}}
    88     def setdefaultparameters(self): #{{{
    89        
    90         #Parameters from de Fleurian 2014
    91         self.water_compressibility    = 5.04e-10
    92         self.isefficientlayer         = 1
    93         self.penalty_factor           = 3
    94         self.rel_tol                  = 1.0e-06
    95         self.max_iter                 = 100
    96         self.sedimentlimit_flag       = 0
    97         self.sedimentlimit            = 0
    98         self.transfer_flag            = 0
    99         self.leakage_factor           = 10.0
    100        
    101         self.sediment_compressibility = 1.0e-08
    102         self.sediment_porosity        = 0.4
    103         self.sediment_thickness       = 20.0
    104         self.sediment_transmitivity   = 8.0e-04
    105        
    106         self.epl_compressibility      = 1.0e-08
    107         self.epl_porosity             = 0.4
    108         self.epl_initial_thickness    = 1.0
    109         self.epl_max_thickness        = 5.0
    110         self.epl_conductivity         = 8.0e-02
    111        
    112         return self
    113     # }}}
    114    
    115     def initialize(self): # {{{
    116         if numpy.all(numpy.isnan(self.basal_moulin_input):
    117             self.basal_moulin_input=numpy.zeros(md.mesh.numberofvertices,1)
    118             print"      no hydrology.basal_moulin_input specified: values set as zero"
     62                if self.sedimentlimit_flag==1:
     63                        string="%s\n%s"%(string,fielddisplay(self,'sedimentlimit','user defined upper limit for the inefficient layer [m]'))
     64                        string="%s\n%s"%(string,fielddisplay(self,'basal_moulin_input','water flux at a given point [m3 s-1]'))
     65                        string="%s\n%s"%(string,fielddisplay(self,'transfer_flag','what kind of transfer method is applied between the layers'))
     66                        string='%55s  0: no transfer',' '
     67                        string='%55s  1: constant leakage factor: %s',' ','leakage_factor'
     68                         
     69                if self.transfer_flag is 1:
     70                        string="%s\n%s"%(string,fielddisplay(self,'leakage_factor','user defined leakage factor [m]'))
     71                        string='   - for the sediment layer'
     72                        string="%s\n%s"%(string,fielddisplay(self,'spcsediment_head','sediment water head constraints (NaN means no constraint) [m above MSL]'))
     73                        string="%s\n%s"%(string,fielddisplay(self,'sediment_compressibility','sediment compressibility [Pa^-1]'))
     74                        string="%s\n%s"%(string,fielddisplay(self,'sediment_porosity','sediment [dimensionless]'))
     75                        string="%s\n%s"%(string,fielddisplay(self,'sediment_thickness','sediment thickness [m]'))
     76                        string="%s\n%s"%(string,fielddisplay(self,'sediment_transmitivity','sediment transmitivity [m^2/s]'))
    11977
    120         return self
     78                if self.isefficientlayer==1:
     79                        string='   - for the epl layer'
     80                        string="%s\n%s"%(string,fielddisplay(self,'spcepl_head','epl water head constraints (NaN means no constraint) [m above MSL]'))
     81                        string="%s\n%s"%(string,fielddisplay(self,'mask_eplactive_node','active (1) or not (0) EPL'))
     82                        string="%s\n%s"%(string,fielddisplay(self,'epl_compressibility','epl compressibility [Pa^-1]'))
     83                        string="%s\n%s"%(string,fielddisplay(self,'epl_porosity','epl [dimensionless]'))
     84                        string="%s\n%s"%(string,fielddisplay(self,'epl_initial_thickness','epl initial thickness [m]'))
     85                        string="%s\n%s"%(string,fielddisplay(self,'epl_conductivity','epl conductivity [m^2/s]'))
     86        #}}}
     87        def setdefaultparameters(self): #{{{
    12188
    122     # }}}
    123     def checkconsistency(self,md,solution,analyses): #{{{
    124        
    125         #Early return
    126         if HydrologyDCInefficientAnalysisEnum() not in analyses:
    127             return md
     89                #Parameters from de Fleurian 2014
     90                self.water_compressibility    = 5.04e-10
     91                self.isefficientlayer         = 1
     92                self.penalty_factor           = 3
     93                self.rel_tol                  = 1.0e-06
     94                self.max_iter                 = 100
     95                self.sedimentlimit_flag       = 0
     96                self.sedimentlimit            = 0
     97                self.transfer_flag            = 0
     98                self.leakage_factor           = 10.0
    12899
    129         md = checkfield(md,'fieldname','hydrology.water_compressibility','>',0,'numel',1)
    130         md = checkfield(md,'fieldname','hydrology.isefficientlayer','numel',[1],'values',[0 1])
    131         md = checkfield(md,'fieldname','hydrology.penalty_factor','>',0,'numel',1)
    132         md = checkfield(md,'fieldname','hydrology.rel_tol','>',0,'numel',1)
    133         md = checkfield(md,'fieldname','hydrology.max_iter','>',0,'numel',1)
    134         md = checkfield(md,'fieldname','hydrology.sedimentlimit_flag','numel',[1],'values',[0 1 2 3])
    135         md = checkfield(md,'fieldname','hydrology.transfer_flag','numel',[1],'values',[0 1])
    136        
    137         if self.sedimentlimit_flag==1:
    138             md = checkfield(md,'fieldname','hydrology.sedimentlimit','>',0,'numel',1)
    139    
    140         if self.transfer_flag==1:
    141             md = checkfield(md,'fieldname','hydrology.leakage_factor','>',0,'numel',1)
     100                self.sediment_compressibility = 1.0e-08
     101                self.sediment_porosity        = 0.4
     102                self.sediment_thickness       = 20.0
     103                self.sediment_transmitivity   = 8.0e-04
    142104
    143         md = checkfield(md,'fieldname','hydrology.basal_moulin_input','NaN',1,'forcing',1)
    144         md = checkfield(md,'fieldname','hydrology.spcsediment_head','forcing',1)
    145         md = checkfield(md,'fieldname','hydrology.sediment_compressibility','>',0,'numel',1)
    146         md = checkfield(md,'fieldname','hydrology.sediment_porosity','>',0,'numel',1)
    147         md = checkfield(md,'fieldname','hydrology.sediment_thickness','>',0,'numel',1)
    148         md = checkfield(md,'fieldname','hydrology.sediment_transmitivity','>=',0,'size',[md.mesh.numberofvertices 1])
    149         if self.isefficientlayer==1:
    150             md = checkfield(md,'fieldname','hydrology.spcepl_head','forcing',1)
    151             md = checkfield(md,'fieldname','hydrology.mask_eplactive_node','size',[md.mesh.numberofvertices 1],'values',[0 1])
    152             md = checkfield(md,'fieldname','hydrology.epl_compressibility','>',0,'numel',1)
    153             md = checkfield(md,'fieldname','hydrology.epl_porosity','>',0,'numel',1)
    154             md = checkfield(md,'fieldname','hydrology.epl_initial_thickness','>',0,'numel',1)
    155             md = checkfield(md,'fieldname','hydrology.epl_conductivity','>',0,'numel',1)
     105                self.epl_compressibility      = 1.0e-08
     106                self.epl_porosity             = 0.4
     107                self.epl_initial_thickness    = 1.0
     108                self.epl_max_thickness        = 5.0
     109                self.epl_conductivity         = 8.0e-02
    156110
    157         # }}}
    158     def marshall(self,md,fid): #{{{
    159         WriteData(fid,'enum',HydrologyModelEnum(),'data',HydrologydcEnum(),'format','Integer')
    160         WriteData(fid,'object',self,'fieldname','water_compressibility','format','Double')
    161         WriteData(fid,'object',self,'fieldname','isefficientlayer','format','Boolean')
    162         WriteData(fid,'object',self,'fieldname','penalty_factor','format','Double')
    163         WriteData(fid,'object',self,'fieldname','penalty_lock','format','Integer')
    164         WriteData(fid,'object',self,'fieldname','rel_tol','format','Double')
    165         WriteData(fid,'object',self,'fieldname','max_iter','format','Integer')
    166         WriteData(fid,'object',self,'fieldname','sedimentlimit_flag','format','Integer')
    167         WriteData(fid,'object',self,'fieldname','transfer_flag','format','Integer')
     111                return self
     112        # }}}
    168113
    169         if self.sedimentlimit_flag==1:
    170             WriteData(fid,'object',self,'fieldname','sedimentlimit','format','Double')
     114        def initialize(self): # {{{
     115                if numpy.all(numpy.isnan(self.basal_moulin_input)):
     116                        self.basal_moulin_input=numpy.zeros(md.mesh.numberofvertices,1)
     117                        print"      no hydrology.basal_moulin_input specified: values set as zero"
    171118
    172         if self.transfer_flag==1:
    173             WriteData(fid,'object',self,'fieldname','leakage_factor','format','Double')
     119                return self
     120        # }}}
     121        def checkconsistency(self,md,solution,analyses): #{{{
    174122
    175         WriteData(fid,'object',self,'fieldname','basal_moulin_input','format','DoubleMat','mattype',1,'forcinglength',md.mesh.numberofvertices+1)
    176         WriteData(fid,'object',self,'fieldname','spcsediment_head','format','DoubleMat','mattype',1,'forcinglength',md.mesh.numberofvertices+1)
    177         WriteData(fid,'object',self,'fieldname','sediment_compressibility','format','Double')
    178         WriteData(fid,'object',self,'fieldname','sediment_porosity','format','Double')                 
    179         WriteData(fid,'object',self,'fieldname','sediment_thickness','format','Double')
    180         WriteData(fid,'object',self,'fieldname','sediment_transmitivity','format','DoubleMat','mattype',1)             
     123                #Early return
     124                if HydrologyDCInefficientAnalysisEnum() not in analyses:
     125                        return md
    181126
    182         if self.isefficientlayer==1:   
    183             WriteData(fid,'object',self,'fieldname','spcepl_head','format','DoubleMat','mattype',1,'forcinglength',md.mesh.numberofvertices+1) 
    184             WriteData(fid,'object',self,'fieldname','mask_eplactive_node','format','DoubleMat','mattype',1)
    185             WriteData(fid,'object',self,'fieldname','epl_compressibility','format','Double')                   
    186             WriteData(fid,'object',self,'fieldname','epl_porosity','format','Double')                   
    187             WriteData(fid,'object',self,'fieldname','epl_initial_thickness','format','Double')
    188             WriteData(fid,'object',self,'fieldname','epl_conductivity','format','Double')
    189 # }}}
     127                md = checkfield(md,'fieldname','hydrology.water_compressibility','>',0,'numel',1)
     128                md = checkfield(md,'fieldname','hydrology.isefficientlayer','numel',[1],'values',[0,1])
     129                md = checkfield(md,'fieldname','hydrology.penalty_factor','>',0,'numel',1)
     130                md = checkfield(md,'fieldname','hydrology.rel_tol','>',0,'numel',1)
     131                md = checkfield(md,'fieldname','hydrology.max_iter','>',0,'numel',1)
     132                md = checkfield(md,'fieldname','hydrology.sedimentlimit_flag','numel',[1],'values',[0,1,2,3])
     133                md = checkfield(md,'fieldname','hydrology.transfer_flag','numel',[1],'values',[0,1])
    190134
     135                if self.sedimentlimit_flag==1:
     136                        md = checkfield(md,'fieldname','hydrology.sedimentlimit','>',0,'numel',1)
     137
     138                if self.transfer_flag==1:
     139                        md = checkfield(md,'fieldname','hydrology.leakage_factor','>',0,'numel',1)
     140
     141                md = checkfield(md,'fieldname','hydrology.basal_moulin_input','NaN',1,'forcing',1)
     142                md = checkfield(md,'fieldname','hydrology.spcsediment_head','forcing',1)
     143                md = checkfield(md,'fieldname','hydrology.sediment_compressibility','>',0,'numel',1)
     144                md = checkfield(md,'fieldname','hydrology.sediment_porosity','>',0,'numel',1)
     145                md = checkfield(md,'fieldname','hydrology.sediment_thickness','>',0,'numel',1)
     146                md = checkfield(md,'fieldname','hydrology.sediment_transmitivity','>=',0,'size',[md.mesh.numberofvertices,1])
     147                if self.isefficientlayer==1:
     148                        md = checkfield(md,'fieldname','hydrology.spcepl_head','forcing',1)
     149                        md = checkfield(md,'fieldname','hydrology.mask_eplactive_node','size',[md.mesh.numberofvertices,1],'values',[0,1])
     150                        md = checkfield(md,'fieldname','hydrology.epl_compressibility','>',0,'numel',1)
     151                        md = checkfield(md,'fieldname','hydrology.epl_porosity','>',0,'numel',1)
     152                        md = checkfield(md,'fieldname','hydrology.epl_initial_thickness','>',0,'numel',1)
     153                        md = checkfield(md,'fieldname','hydrology.epl_conductivity','>',0,'numel',1)
     154        # }}}
     155        def marshall(self,md,fid): #{{{
     156                WriteData(fid,'enum',HydrologyModelEnum(),'data',HydrologydcEnum(),'format','Integer')
     157                WriteData(fid,'object',self,'fieldname','water_compressibility','format','Double')
     158                WriteData(fid,'object',self,'fieldname','isefficientlayer','format','Boolean')
     159                WriteData(fid,'object',self,'fieldname','penalty_factor','format','Double')
     160                WriteData(fid,'object',self,'fieldname','penalty_lock','format','Integer')
     161                WriteData(fid,'object',self,'fieldname','rel_tol','format','Double')
     162                WriteData(fid,'object',self,'fieldname','max_iter','format','Integer')
     163                WriteData(fid,'object',self,'fieldname','sedimentlimit_flag','format','Integer')
     164                WriteData(fid,'object',self,'fieldname','transfer_flag','format','Integer')
     165                if self.sedimentlimit_flag==1:
     166                        WriteData(fid,'object',self,'fieldname','sedimentlimit','format','Double')
     167
     168                if self.transfer_flag==1:
     169                        WriteData(fid,'object',self,'fieldname','leakage_factor','format','Double')
     170
     171                WriteData(fid,'object',self,'fieldname','basal_moulin_input','format','DoubleMat','mattype',1,'forcinglength',md.mesh.numberofvertices+1)
     172                WriteData(fid,'object',self,'fieldname','spcsediment_head','format','DoubleMat','mattype',1,'forcinglength',md.mesh.numberofvertices+1)
     173                WriteData(fid,'object',self,'fieldname','sediment_compressibility','format','Double')
     174                WriteData(fid,'object',self,'fieldname','sediment_porosity','format','Double')                 
     175                WriteData(fid,'object',self,'fieldname','sediment_thickness','format','Double')
     176                WriteData(fid,'object',self,'fieldname','sediment_transmitivity','format','DoubleMat','mattype',1)             
     177
     178                if self.isefficientlayer==1:   
     179                        WriteData(fid,'object',self,'fieldname','spcepl_head','format','DoubleMat','mattype',1,'forcinglength',md.mesh.numberofvertices+1)     
     180                        WriteData(fid,'object',self,'fieldname','mask_eplactive_node','format','DoubleMat','mattype',1)
     181                        WriteData(fid,'object',self,'fieldname','epl_compressibility','format','Double')                       
     182                        WriteData(fid,'object',self,'fieldname','epl_porosity','format','Double')                       
     183                        WriteData(fid,'object',self,'fieldname','epl_initial_thickness','format','Double')
     184                        WriteData(fid,'object',self,'fieldname','epl_conductivity','format','Double')
     185        # }}}
Note: See TracChangeset for help on using the changeset viewer.