source: issm/oecreview/Archive/25834-26739/ISSM-26300-26301.diff

Last change on this file was 26740, checked in by Mathieu Morlighem, 3 years ago

CHG: added 25834-26739
File size: 167.9 KB

  • ../trunk-jpl/src/m/boundaryconditions/getlovenumbers.py

     
    44
    55
    66def getlovenumbers(*args): #{{{
    7     '''
    8     GETLOVENUMBERS - provide love numbers retrieved from:
     7    """GETLOVENUMBERS - provide love numbers retrieved from:
    98    http://www.srosat.com/iag-jsg/loveNb.php in a chosen reference frame
    109
    11         Usage:
    12             series = love_numbers('type', 'loadingverticaldisplacement', 'referenceframe', 'CM', 'maxdeg', 1001)
     10    Usage:
     11        series = love_numbers('type', 'loadingverticaldisplacement', 'referenceframe', 'CM', 'maxdeg', 1001)
    1312
    14         - type = one of 'loadingverticaldisplacement',
    15         'loadinggravitationalpotential', 'loadinghorizontaldisplacement',
    16         'tidalverticaldisplacement', 'tidalgravitationalpotential',
    17         'tidalhorizontaldisplacement'
    18         - reference_frame = one of 'CM' (default) and 'CF'
    19         - maxdeg = default 1000
     13    - type = one of 'loadingverticaldisplacement',
     14    'loadinggravitationalpotential', 'loadinghorizontaldisplacement',
     15    'tidalverticaldisplacement', 'tidalgravitationalpotential',
     16    'tidalhorizontaldisplacement'
     17    - reference_frame = one of 'CM' (default) and 'CF'
     18    - maxdeg = default 1000
    2019
    21         Example:
    22             h = love_number
    23     '''
     20    Example:
     21        h = getlovenumbers('type', 'loadingverticaldisplacement', 'referenceframe', 'CM', 'maxdeg', maxdeg)
     22        k = getlovenumbers('type', 'loadinggravitationalpotential', 'referenceframe', 'CM', 'maxdeg', maxdeg)
     23        l = getlovenumbers('type', 'loadinghorizontaldisplacement', 'referenceframe', 'CM', 'maxdeg', maxdeg)
     24        th = getlovenumbers('type', 'tidalverticaldisplacement', 'referenceframe', 'CM', 'maxdeg', maxdeg)
     25        tk = getlovenumbers('type', 'tidalgravitationalpotential', 'referenceframe', 'CM', 'maxdeg', maxdeg)
     26        tl = getlovenumbers('type', 'tidalhorizontaldisplacement', 'referenceframe', 'CM', 'maxdeg', maxdeg)
     27    """
    2428
    2529    TYPES = [
    2630        'loadingverticaldisplacement',
    … …  
    3135        'tidalhorizontaldisplacement'
    3236    ]
    3337
    34     #recover options
     38    # Recover options
    3539    options = pairoptions(*args)
    3640    type = options.getfieldvalue('type')
    3741    frame = options.getfieldvalue('referenceframe', 'CM')
    3842    maxdeg = options.getfieldvalue('maxdeg', 1000)
    3943
     44    if (maxdeg > 10000):
     45        raise Exception('PREM love numbers computed only for deg < 10,000. Request lower maxdeg')
     46
    4047    if type not in TYPES:
    4148        raise Exception('type should be one of {}'.format(TYPES))
    4249
    … …  
    1004410051        [-6.27342778, -0.00030945, 0.00018956, 0.00031100, -0.00000116, 0.00000000]
    1004510052    ])
    1004610053
    10047     #cut off series at maxdeg
     10054    # Cut off series at maxdeg
    1004810055    love_numbers = np.delete(love_numbers, range(maxdeg + 1, len(love_numbers)), axis=0)
    1004910056
    10050     #retrieve right type
     10057    # Retrieve right type
    1005110058    if type == 'loadingverticaldisplacement':
    1005210059        series = love_numbers[:, 0]
    1005310060    elif type == 'loadinggravitationalpotential':
    … …  
    1006110068    elif type == 'tidalhorizontaldisplacement':
    1006210069        series = love_numbers[:, 5]
    1006310070    else:
    10064         raise Exception("love_numbers error message: unknown type: {}".format(type))
     10071        raise Exception('love_numbers error message: unknown type: {}'.format(type))
    1006510072
    1006610073    #choose degree 1 term for CF reference system
    1006710074    if frame == 'CM':
    … …  
    1007410081        elif type == 'loadinghorizontaldisplacement':
    1007510082            series[1] = 0.134
    1007610083    else:
    10077         raise Exception("love_numbers error message: unknown reference frame: {}".format(frame))
     10084        raise Exception('love_numbers error message: unknown reference frame: {}'.format(frame))
    1007810085
    1007910086    return series
    1008010087#}}}

  • ../trunk-jpl/src/m/classes/amr.m

     
    8383                        %control of element lengths
    8484                        self.gradation=1.5;
    8585
    86                         %other criterias
     86                        %other criteria
    8787                        self.groundingline_resolution=500.;
    8888                        self.groundingline_distance=0.;
    8989                        self.icefront_resolution=500.;

  • ../trunk-jpl/src/m/classes/amr.py

     
    44
    55
    66class amr(object):
    7     """
    8     AMR Class definition
     7    """AMR Class definition
    98
    109    Usage:
    1110        amr = amr()
    1211    """
    1312
    14     def __init__(self):  # {{{
    15         self.hmin = 0.
    16         self.hmax = 0.
     13    def __init__(self): #{{{
     14        self.hmin = 0
     15        self.hmax = 0
    1716        self.fieldname = ''
    18         self.err = 0.
    19         self.keepmetric = 0.
    20         self.gradation = 0.
    21         self.groundingline_resolution = 0.
    22         self.groundingline_distance = 0.
    23         self.icefront_resolution = 0.
    24         self.icefront_distance = 0.
    25         self.thicknesserror_resolution = 0.
    26         self.thicknesserror_threshold = 0.
    27         self.thicknesserror_groupthreshold = 0.
    28         self.thicknesserror_maximum = 0.
    29         self.deviatoricerror_resolution = 0.
    30         self.deviatoricerror_threshold = 0.
    31         self.deviatoricerror_groupthreshold = 0.
    32         self.deviatoricerror_maximum = 0.
    33         self.restart = 0.
    34     #set defaults
     17        self.err = 0
     18        self.keepmetric = 0
     19        self.gradation = 0
     20        self.groundingline_resolution = 0
     21        self.groundingline_distance = 0
     22        self.icefront_resolution = 0
     23        self.icefront_distance = 0
     24        self.thicknesserror_resolution = 0
     25        self.thicknesserror_threshold = 0
     26        self.thicknesserror_groupthreshold = 0
     27        self.thicknesserror_maximum = 0
     28        self.deviatoricerror_resolution = 0
     29        self.deviatoricerror_threshold = 0
     30        self.deviatoricerror_groupthreshold = 0
     31        self.deviatoricerror_maximum = 0
     32        self.restart = 0
     33
    3534        self.setdefaultparameters()
    3635    #}}}
    3736
    38     def __repr__(self):  # {{{
    39         string = "   amr parameters:"
    40         string = "%s\n%s" % (string, fielddisplay(self, "hmin", "minimum element length"))
    41         string = "%s\n%s" % (string, fielddisplay(self, "hmax", "maximum element length"))
    42         string = "%s\n%s" % (string, fielddisplay(self, "fieldname", "name of input that will be used to compute the metric (should be an input of FemModel)"))
    43         string = "%s\n%s" % (string, fielddisplay(self, "keepmetric", "indicates whether the metric should be kept every remeshing time"))
    44         string = "%s\n%s" % (string, fielddisplay(self, "gradation", "maximum ratio between two adjacent edges"))
    45         string = "%s\n%s" % (string, fielddisplay(self, "groundingline_resolution", "element length near the grounding line"))
    46         string = "%s\n%s" % (string, fielddisplay(self, "groundingline_distance", "distance around the grounding line which elements will be refined"))
    47         string = "%s\n%s" % (string, fielddisplay(self, "icefront_resolution", "element length near the ice front"))
    48         string = "%s\n%s" % (string, fielddisplay(self, "icefront_distance", "distance around the ice front which elements will be refined"))
    49         string = "%s\n%s" % (string, fielddisplay(self, "thicknesserror_resolution", "element length when thickness error estimator is used"))
    50         string = "%s\n%s" % (string, fielddisplay(self, "thicknesserror_threshold", "maximum threshold thickness error permitted"))
    51         string = "%s\n%s" % (string, fielddisplay(self, "thicknesserror_groupthreshold", "maximum group threshold thickness error permitted"))
    52         string = "%s\n%s" % (string, fielddisplay(self, "thicknesserror_maximum", "maximum thickness error permitted"))
    53         string = "%s\n%s" % (string, fielddisplay(self, "deviatoricerror_resolution", "element length when deviatoric stress error estimator is used"))
    54         string = "%s\n%s" % (string, fielddisplay(self, "deviatoricerror_threshold", "maximum threshold deviatoricstress error permitted"))
    55         string = "%s\n%s" % (string, fielddisplay(self, "deviatoricerror_groupthreshold", "maximum group threshold deviatoric stress error permitted"))
    56         string = "%s\n%s" % (string, fielddisplay(self, "deviatoricerror_maximum", "maximum deviatoricstress error permitted"))
    57         string = "%s\n%s" % (string, fielddisplay(self, "restart", "indicates if ReMesh() will call before first time step"))
    58         return string
     37    def __repr__(self): #{{{
     38        s = '   amr parameters:\n'
     39        s += '{}\n'.format(fielddisplay(self, 'hmin', 'minimum element length'))
     40        s += '{}\n'.format(fielddisplay(self, 'hmax', 'maximum element length'))
     41        s += '{}\n'.format(fielddisplay(self, 'fieldname', 'name of input that will be used to compute the metric (should be an input of FemModel)'))
     42        s += '{}\n'.format(fielddisplay(self, 'keepmetric', 'indicates whether the metric should be kept every remeshing time'))
     43        s += '{}\n'.format(fielddisplay(self, 'gradation', 'maximum ratio between two adjacent edges'))
     44        s += '{}\n'.format(fielddisplay(self, 'groundingline_resolution', 'element length near the grounding line'))
     45        s += '{}\n'.format(fielddisplay(self, 'groundingline_distance', 'distance around the grounding line which elements will be refined'))
     46        s += '{}\n'.format(fielddisplay(self, 'icefront_resolution', 'element length near the ice front'))
     47        s += '{}\n'.format(fielddisplay(self, 'icefront_distance', 'distance around the ice front which elements will be refined'))
     48        s += '{}\n'.format(fielddisplay(self, 'thicknesserror_resolution', 'element length when thickness error estimator is used'))
     49        s += '{}\n'.format(fielddisplay(self, 'thicknesserror_threshold', 'maximum threshold thickness error permitted'))
     50        s += '{}\n'.format(fielddisplay(self, 'thicknesserror_groupthreshold', 'maximum group threshold thickness error permitted'))
     51        s += '{}\n'.format(fielddisplay(self, 'thicknesserror_maximum', 'maximum thickness error permitted'))
     52        s += '{}\n'.format(fielddisplay(self, 'deviatoricerror_resolution', 'element length when deviatoric stress error estimator is used'))
     53        s += '{}\n'.format(fielddisplay(self, 'deviatoricerror_threshold', 'maximum threshold deviatoricstress error permitted'))
     54        s += '{}\n'.format(fielddisplay(self, 'deviatoricerror_groupthreshold', 'maximum group threshold deviatoric stress error permitted'))
     55        s += '{}\n'.format(fielddisplay(self, 'deviatoricerror_maximum', 'maximum deviatoricstress error permitted'))
     56        s += '{}\n'.format(fielddisplay(self, 'restart', 'indicates if ReMesh() will call before first time step'))
     57        return s
    5958    #}}}
    6059
    61     def setdefaultparameters(self):  # {{{
    62         self.hmin = 100.
    63         self.hmax = 100.e3
     60    def setdefaultparameters(self): #{{{
     61        self.hmin = 100
     62        self.hmax = 100e3
     63
     64        # Fields
    6465        self.fieldname = 'Vel'
    65         self.err = 3.
     66        self.err = 3
     67
     68        # Keep metric?
    6669        self.keepmetric = 1
     70
     71        # Control of element lengths
    6772        self.gradation = 1.5
    68         self.groundingline_resolution = 500.
     73
     74        # Other criteria
     75        self.groundingline_resolution = 500
    6976        self.groundingline_distance = 0
    70         self.icefront_resolution = 500.
     77        self.icefront_resolution = 500
    7178        self.icefront_distance = 0
    72         self.thicknesserror_resolution = 500.
     79        self.thicknesserror_resolution = 500
    7380        self.thicknesserror_threshold = 0
    7481        self.thicknesserror_groupthreshold = 0
    7582        self.thicknesserror_maximum = 0
    76         self.deviatoricerror_resolution = 500.
     83        self.deviatoricerror_resolution = 500
    7784        self.deviatoricerror_threshold = 0
    7885        self.deviatoricerror_groupthreshold = 0
    7986        self.deviatoricerror_maximum = 0
    80         self.restart = 0.
     87
     88        # Is restart? This calls femmodel->ReMesh() before first time step.
     89        self.restart = 0
    8190        return self
    8291    #}}}
    8392
    84     def checkconsistency(self, md, solution, analyses):  # {{{
     93    def checkconsistency(self, md, solution, analyses): #{{{
    8594        md = checkfield(md, 'fieldname', 'amr.hmax', 'numel', [1], '>', 0, 'NaN', 1)
    8695        md = checkfield(md, 'fieldname', 'amr.hmin', 'numel', [1], '>', 0, '<', self.hmax, 'NaN', 1)
    8796        md = checkfield(md, 'fieldname', 'amr.keepmetric', 'numel', [1], '>=', 0, '<=', 1, 'NaN', 1)
    … …  
    100109        md = checkfield(md, 'fieldname', 'amr.deviatoricerror_maximum', 'numel', [1], '>=', 0, 'NaN', 1, 'Inf', 1)
    101110        md = checkfield(md, 'fieldname', 'amr.restart', 'numel', [1], '>=', 0, '<=', 1, 'NaN', 1)
    102111        return md
    103     # }}}
     112   # }}}
    104113
    105     def marshall(self, prefix, md, fid):  # {{{
     114    def marshall(self, prefix, md, fid): #{{{
    106115        WriteData(fid, prefix, 'name', 'md.amr.type', 'data', 1, 'format', 'Integer')
    107116        WriteData(fid, prefix, 'object', self, 'fieldname', 'hmin', 'format', 'Double')
    108117        WriteData(fid, prefix, 'object', self, 'fieldname', 'hmax', 'format', 'Double')
    … …  
    123132        WriteData(fid, prefix, 'object', self, 'fieldname', 'deviatoricerror_groupthreshold', 'format', 'Double')
    124133        WriteData(fid, prefix, 'object', self, 'fieldname', 'deviatoricerror_maximum', 'format', 'Double')
    125134        WriteData(fid, prefix, 'object', self, 'class', 'amr', 'fieldname', 'restart', 'format', 'Integer')
    126     # }}}
     135    #}}}

  • ../trunk-jpl/src/m/classes/dsl.m

     
    1212
    1313        end
    1414        methods
    15                 function self = extrude(self,md) % {{{
    16                         self.sea_surface_height_above_geoid=project3d(md,'vector',self.sea_surface_height_above_geoid,'type','node','layer',1);
    17                         self.sea_water_pressure_at_sea_floor=project3d(md,'vector',self.sea_water_pressure_at_sea_floor,'type','node','layer',1);
    18                 end % }}}
    1915                function self = dsl(varargin) % {{{
    2016                        switch nargin
    2117                                case 0
    … …  
    3329                        self.sea_water_pressure_at_sea_floor=NaN;
    3430
    3531                end % }}}
     32                function disp(self) % {{{
     33
     34                        disp(sprintf('   dsl parameters:'));
     35                        fielddisplay(self,'global_average_thermosteric_sea_level','Corresponds to zostoga field in CMIP5 archives. Specified as a temporally variable quantity (in m).');
     36                        fielddisplay(self,'sea_surface_height_above_geoid','Corresponds to zos field in CMIP5 archives. Spatial average is 0. Specified as a spatio-temporally variable quantity (in m).');
     37                        fielddisplay(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!).');
     38
     39                end % }}}
    3640                function md = checkconsistency(self,md,solution,analyses) % {{{
    3741
    3842                        %Early return
    … …  
    4852                        end
    4953                       
    5054                end % }}}
    51                 function disp(self) % {{{
    52 
    53                         disp(sprintf('   dsl parameters:'));
    54                         fielddisplay(self,'global_average_thermosteric_sea_level','Corresponds to zostoga field in CMIP5 archives. Specified as a temporally variable quantity (in m).');
    55                         fielddisplay(self,'sea_surface_height_above_geoid','Corresponds to zos field in CMIP5 archives. Spatial average is 0. Specified as a spatio-temporally variable quantity (in m).');
    56                         fielddisplay(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!).');
    57 
    58                 end % }}}
    5955                function marshall(self,prefix,md,fid) % {{{
    60 
     56                        yts=md.constants.yts;
    6157                        WriteData(fid,prefix,'name','md.dsl.model','data',1,'format','Integer');
    62                         WriteData(fid,prefix,'object',self,'fieldname','global_average_thermosteric_sea_level','format','DoubleMat','mattype',2,'timeseries',1,'timeserieslength',2,'yts',md.constants.yts); %mattype 2, because we are sending a GMSL value identical everywhere on each element.
    63                         WriteData(fid,prefix,'object',self,'fieldname','sea_surface_height_above_geoid','format','DoubleMat','mattype',1,'timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts); %mattype 1 because we specify DSL at vertex locations.
    64                         WriteData(fid,prefix,'object',self,'fieldname','sea_water_pressure_at_sea_floor','format','DoubleMat','mattype',1,'timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts); %mattype 1 because we specify bottom pressure at vertex locations.
     58                        WriteData(fid,prefix,'object',self,'fieldname','global_average_thermosteric_sea_level','format','DoubleMat','mattype',2,'timeseries',1,'timeserieslength',2,'yts',yts); %mattype 2, because we are sending a GMSL value identical everywhere on each element.
     59                        WriteData(fid,prefix,'object',self,'fieldname','sea_surface_height_above_geoid','format','DoubleMat','mattype',1,'timeserieslength',md.mesh.numberofvertices+1,'yts',yts); %mattype 1 because we specify DSL at vertex locations.
     60                        WriteData(fid,prefix,'object',self,'fieldname','sea_water_pressure_at_sea_floor','format','DoubleMat','mattype',1,'timeserieslength',md.mesh.numberofvertices+1,'yts',yts); %mattype 1 because we specify bottom pressure at vertex locations.
    6561
    6662                end % }}}
    67                 function savemodeljs(self,fid,modelname) % {{{
    68 
    69                         writejs1Darray(fid,[modelname '.dsl.global_average_thermosteric_sea_level'],self.global_average_thermosteric_sea_level);
    70                         writejs1Darray(fid,[modelname '.dsl.sea_surface_height_above_geoid'],self.sea_surface_height_above_geoid);
    71                         writejs1Darray(fid,[modelname '.dsl.sea_water_pressure_at_sea_floor'],self.sea_water_pressure_at_sea_floor);
    72 
     63                function self = extrude(self,md) % {{{
     64                        self.sea_surface_height_above_geoid=project3d(md,'vector',self.sea_surface_height_above_geoid,'type','node','layer',1);
     65                        self.sea_water_pressure_at_sea_floor=project3d(md,'vector',self.sea_water_pressure_at_sea_floor,'type','node','layer',1);
    7366                end % }}}
    7467                function self = initialize(self,md) % {{{
    7568
    … …  
    8679                                disp('      no dsl.sea_water_pressure_at_sea_floor specified: transient values set at zero');
    8780                        end
    8881                end % }}}
    89        
     82                function savemodeljs(self,fid,modelname) % {{{
     83
     84                        writejs1Darray(fid,[modelname '.dsl.global_average_thermosteric_sea_level'],self.global_average_thermosteric_sea_level);
     85                        writejs1Darray(fid,[modelname '.dsl.sea_surface_height_above_geoid'],self.sea_surface_height_above_geoid);
     86                        writejs1Darray(fid,[modelname '.dsl.sea_water_pressure_at_sea_floor'],self.sea_water_pressure_at_sea_floor);
     87
     88                end % }}}
    9089        end
    9190end

  • ../trunk-jpl/src/m/classes/fourierlove.m

     
    55
    66classdef fourierlove
    77        properties (SetAccess=public)
    8                 nfreq                      = 0;
    9                 frequencies                = 0;
    10                 sh_nmax                    = 0;
    11                 sh_nmin                    = 0;
    12                 g0                         = 0;
    13                 r0                         = 0;
    14                 mu0                        = 0;
    15                 Gravitational_Constant     = 0;
    16                 allow_layer_deletion       = 0;
    17                 underflow_tol              = 0;
    18                 integration_steps_per_layer= 0;
    19                 istemporal                 = 0;
    20                 n_temporal_iterations      = 0;
    21                 time                               = 0;
    22                 love_kernels               = 0;
    23                 forcing_type               = 0;
    24                 inner_core_boundary            = 0;
    25                 core_mantle_boundary       = 0;
    26                 complex_computation        = 0;
     8                nfreq                       = 0;
     9                frequencies                 = 0;
     10                sh_nmax                     = 0;
     11                sh_nmin                     = 0;
     12                g0                          = 0;
     13                r0                          = 0;
     14                mu0                         = 0;
     15                Gravitational_Constant      = 0;
     16                allow_layer_deletion        = 0;
     17                underflow_tol               = 0;
     18                integration_steps_per_layer = 0;
     19                istemporal                  = 0;
     20                n_temporal_iterations       = 0;
     21                time                        = 0;
     22                love_kernels                = 0;
     23                forcing_type                = 0;
     24                inner_core_boundary         = 0;
     25                core_mantle_boundary        = 0;
     26                complex_computation         = 0;
    2727
    2828        end
    2929        methods (Static)
    … …  
    3333                end% }}}
    3434        end
    3535        methods
    36                 function self = extrude(self,md) % {{{
    37                 end % }}}
    3836                function self = fourierlove(varargin) % {{{
    3937                        switch nargin
    4038                                case 0
    … …  
    5250                        % work on matlab script for computing g0 for given Earth's structure.
    5351                        self.g0=9.81; % m/s^2;
    5452                        self.r0=6371*1e3; %m;
    55                         self.mu0=10^11; % Pa
     53                        self.mu0=1e11; % Pa
    5654                        self.Gravitational_Constant=6.67259e-11; % m^3 kg^-1 s^-2
    5755                        self.allow_layer_deletion=1;
    5856                        self.underflow_tol=1e-16; %threshold of deep to surface love number ratio to trigger the deletion of layer
    … …  
    6765                        self.complex_computation=0;
    6866                end % }}}
    6967                function disp(self) % {{{
    70                         fielddisplay(self,'nfreq','number of frequencies sampled (default 1, elastic) [Hz]');
     68                        fielddisplay(self,'nfreq','number of frequencies sampled (default: 1, elastic) [Hz]');
    7169                        fielddisplay(self,'frequencies','frequencies sampled (convention defaults to 0 for the elastic case) [Hz]');
    72                         fielddisplay(self,'sh_nmax','maximum spherical harmonic degree (default 256, .35 deg, or 40 km at equator)');
    73                         fielddisplay(self,'sh_nmin','minimum spherical harmonic degree (default 1)');
    74                         fielddisplay(self,'g0','adimensioning constant for gravity (default 10) [m/s^2]');
    75                         fielddisplay(self,'r0','adimensioning constant for radius (default 6371*10^3) [m]');
    76                         fielddisplay(self,'mu0','adimensioning constant for stress (default 10^11) [Pa]');
    77                         fielddisplay(self,'Gravitational_Constant','Newtonian constant of gravitation (default 6.67259e-11 [m^3 kg^-1 s^-2])');
    78                         fielddisplay(self,'allow_layer_deletion','allow for migration of the integration boundary with increasing spherical harmonics degree (default 1)');
    79                         fielddisplay(self,'underflow_tol','threshold of deep to surface love number ratio to trigger the deletion of layers (default 2.2204460492503131E-016)');
    80                         fielddisplay(self,'integration_steps_per_layer','number of radial steps to propagate the yi system from the bottom to the top of each layer (default 100)');
    81                         fielddisplay(self,'istemporal',{'1 for time-dependent love numbers, 0 for frequency-dependent or elastic love numbers (default 0)', 'If 1: use fourierlove function build_frequencies_from_time to meet consistency'});
    82                         fielddisplay(self,'n_temporal_iterations','max number of iterations in the inverse Laplace transform. Also the number of spectral samples per time step requested (default 8)');
    83                         fielddisplay(self,'time','time vector for deformation if istemporal (default 0) [s]');
    84                         fielddisplay(self,'love_kernels','compute love numbers at depth? (default 0)');
    85                         fielddisplay(self,'forcing_type',{'integer indicating the nature and depth of the forcing for the Love number calculation (default 11) :','1:  Inner core boundary -- Volumic Potential','2:  Inner core boundary -- Pressure','3:  Inner core boundary -- Loading','4:  Inner core boundary -- Tangential traction','5:  Core mantle boundary -- Volumic Potential','6:  Core mantle boundary -- Pressure','7:  Core mantle boundary -- Loading','8:  Core mantle boundary -- Tangential traction','9:  Surface -- Volumic Potential','10: Surface -- Pressure','11: Surface -- Loading','12: Surface -- Tangential traction '});
    86                         fielddisplay(self,'inner_core_boundary','interface index in materials.radius locating forcing. Only used for forcing_type 1--4 (default 1)');
    87                         fielddisplay(self,'core_mantle_boundary','interface index in materials.radius locating forcing. Only used for forcing_type 5--8 (default 2)');
     70                        fielddisplay(self,'sh_nmax','maximum spherical harmonic degree (default: 256, .35 deg, or 40 km at equator)');
     71                        fielddisplay(self,'sh_nmin','minimum spherical harmonic degree (default: 1)');
     72                        fielddisplay(self,'g0','adimensioning constant for gravity (default: 10) [m/s^2]');
     73                        fielddisplay(self,'r0','adimensioning constant for radius (default: 6371*10^3) [m]');
     74                        fielddisplay(self,'mu0','adimensioning constant for stress (default: 10^11) [Pa]');
     75                        fielddisplay(self,'Gravitational_Constant','Newtonian constant of gravitation (default: 6.67259e-11 [m^3 kg^-1 s^-2])');
     76                        fielddisplay(self,'allow_layer_deletion','allow for migration of the integration boundary with increasing spherical harmonics degree (default: 1)');
     77                        fielddisplay(self,'underflow_tol','threshold of deep to surface love number ratio to trigger the deletion of layers (default: 1e-16)');
     78                        fielddisplay(self,'integration_steps_per_layer','number of radial steps to propagate the yi system from the bottom to the top of each layer (default: 100)');
     79                        fielddisplay(self,'istemporal',{'1 for time-dependent love numbers, 0 for frequency-dependent or elastic love numbers (default: 0)', 'If 1: use fourierlove function build_frequencies_from_time to meet consistency'});
     80                        fielddisplay(self,'n_temporal_iterations','max number of iterations in the inverse Laplace transform. Also the number of spectral samples per time step requested (default: 8)');
     81                        fielddisplay(self,'time','time vector for deformation if istemporal (default: 0) [s]');
     82                        fielddisplay(self,'love_kernels','compute love numbers at depth? (default: 0)');
     83                        fielddisplay(self,'forcing_type',{'integer indicating the nature and depth of the forcing for the Love number calculation (default: 11):','1:  Inner core boundary -- Volumic Potential','2:  Inner core boundary -- Pressure','3:  Inner core boundary -- Loading','4:  Inner core boundary -- Tangential traction','5:  Core mantle boundary -- Volumic Potential','6:  Core mantle boundary -- Pressure','7:  Core mantle boundary -- Loading','8:  Core mantle boundary -- Tangential traction','9:  Surface -- Volumic Potential','10: Surface -- Pressure','11: Surface -- Loading','12: Surface -- Tangential traction '});
     84                        fielddisplay(self,'inner_core_boundary','interface index in materials.radius locating forcing. Only used for forcing_type 1--4 (default: 1)');
     85                        fielddisplay(self,'core_mantle_boundary','interface index in materials.radius locating forcing. Only used for forcing_type 5--8 (default: 2)');
    8886
    8987                end % }}}
    9088                function md = checkconsistency(self,md,solution,analyses) % {{{
    … …  
    111109                                md = checkfield(md,'fieldname','love.n_temporal_iterations','NaN',1,'Inf',1,'numel',1,'>',0);
    112110                                md = checkfield(md,'fieldname','love.time','NaN',1,'Inf',1,'numel',md.love.nfreq/2/md.love.n_temporal_iterations);
    113111                        end
    114                         if md.love.sh_nmin<=1 & (md.love.forcing_type==9 || md.love.forcing_type==5 || md.love.forcing_type==1)
    115                                 error('Degree 1 not supported for Volumetric Potential forcing. Use sh_min>=2 for this kind of calculation.')
     112                        if md.love.sh_nmin<=1 & (md.love.forcing_type==1 || md.love.forcing_type==5 || md.love.forcing_type==9)
     113                                error(['Degree 1 not supported for forcing type ' num2str(md.love.forcing_type) '. Use sh_min>=2 for this kind of calculation.'])
    116114                        end
    117115
    118116                        %need 'litho' material:
    … …  
    129127
    130128                end % }}}
    131129                function marshall(self,prefix,md,fid) % {{{
    132                
     130
    133131                        WriteData(fid,prefix,'object',self,'fieldname','nfreq','format','Integer');
    134132                        WriteData(fid,prefix,'object',self,'fieldname','frequencies','format','DoubleMat','mattype',3);
    135133                        WriteData(fid,prefix,'object',self,'fieldname','sh_nmax','format','Integer');
    … …  
    151149                        WriteData(fid,prefix,'object',self,'fieldname','core_mantle_boundary','format','Integer');
    152150
    153151                end % }}}
     152                function self = extrude(self,md) % {{{
     153                end % }}}
    154154                function savemodeljs(self,fid,modelname) % {{{
    155155                        error('not implemented yet!');
    156156                end % }}}

  • ../trunk-jpl/src/m/classes/fourierlove.py

     
     1import numpy as np
     2
    13from fielddisplay import fielddisplay
    24from checkfield import checkfield
    35from WriteData import WriteData
    … …  
    46
    57
    68class fourierlove(object):
    7     """FOURIERLOVE - Fourier Love Number class definition
     9    """FOURIERLOVE - class definition
    810
    911    Usage:
    10         fourierlove = fourierlove()
     12        md.love = fourierlove()
    1113    """
    12     def __init__(self):  # {{{
     14
     15    def __init__(self): #{{{
    1316        self.nfreq = 0
    1417        self.frequencies = 0
    1518        self.sh_nmax = 0
    … …  
    1720        self.g0 = 0
    1821        self.r0 = 0
    1922        self.mu0 = 0
     23        self.Gravitational_Constant = 0
    2024        self.allow_layer_deletion = 0
     25        self.underflow_tol = 0
     26        self.integration_steps_per_layer = 0
     27        self.istemporal = 0
     28        self.n_temporal_iterations = 0
     29        self.time = 0
    2130        self.love_kernels = 0
    2231        self.forcing_type = 0
     32        self.inner_core_boundary = 0
     33        self.core_mantle_boundary = 0
     34        self.complex_computation = 0
    2335
    2436        self.setdefaultparameters()
    2537    #}}}
    2638
    27     def __repr__(self):  # {{{
    28         # TODO:
    29         # - Convert all formatting to calls to <string>.format (see any
    30         #   already converted <class>.__repr__ method for examples)
    31         #
    32         string = '   Fourier Love class:'
    33         string = "%s\n%s" % (string, fielddisplay(self, 'nfreq', 'number of frequencies sampled (default 1, elastic) [Hz]'))
    34         string = "%s\n%s" % (string, fielddisplay(self, 'frequencies', 'frequencies sampled (convention defaults to 0 for the elastic case) [Hz]'))
    35         string = "%s\n%s" % (string, fielddisplay(self, 'sh_nmax', 'maximum spherical harmonic degree (default 256, .35 deg, or 40 km at equator)'))
    36         string = "%s\n%s" % (string, fielddisplay(self, 'sh_nmin', 'minimum spherical harmonic degree (default 1)'))
    37         string = "%s\n%s" % (string, fielddisplay(self, 'g0', 'adimensioning constant for gravity (default 10) [m / s^2]'))
    38         string = "%s\n%s" % (string, fielddisplay(self, 'r0', 'adimensioning constant for radius (default 6378e3) [m]'))
    39         string = "%s\n%s" % (string, fielddisplay(self, 'mu0', 'adimensioning constant for stress (default 1.0e11) [Pa]'))
    40         string = "%s\n%s" % (string, fielddisplay(self, 'allow_layer_deletion', 'allow for migration of the integration boundary with increasing spherical harmonics degree (default 1)'))
    41         string = "%s\n%s" % (string, fielddisplay(self, 'love_kernels', 'compute love numbers at depth? (default 0)'))
    42         string = "%s\n%s" % (string, fielddisplay(self, 'forcing_type', 'integer indicating the nature and depth of the forcing for the Love number calculation (default 11) :'))
    43         string = "%s\n%s" % (string, '                                                     1:  Inner core boundary -- Volumic Potential')
    44         string = "%s\n%s" % (string, '                                                     2:  Inner core boundary --  Pressure')
    45         string = "%s\n%s" % (string, '                                                     3:  Inner core boundary --  Loading')
    46         string = "%s\n%s" % (string, '                                                     4:  Inner core boundary --  Tangential traction')
    47         string = "%s\n%s" % (string, '                                                     5:  Core mantle boundary --  Volumic Potential')
    48         string = "%s\n%s" % (string, '                                                     6:  Core mantle boundary --  Pressure')
    49         string = "%s\n%s" % (string, '                                                     7:  Core mantle boundary --  Loading')
    50         string = "%s\n%s" % (string, '                                                     8:  Core mantle boundary --  Tangential traction')
    51         string = "%s\n%s" % (string, '                                                     9:  Surface--  Volumic Potential')
    52         string = "%s\n%s" % (string, '                                                     10: Surface--  Pressure')
    53         string = "%s\n%s" % (string, '                                                     11: Surface--  Loading')
    54         string = "%s\n%s" % (string, '                                                     12: Surface--  Tangential traction ')
     39    def __repr__(self): #{{{
     40        s = '   Fourier Love class:\n'
     41        s += '{}\n'.format(fielddisplay(self, 'nfreq', 'number of frequencies sampled (default: 1, elastic) [Hz]'))
     42        s += '{}\n'.format(fielddisplay(self, 'frequencies', 'frequencies sampled (convention defaults to 0 for the elastic case) [Hz]'))
     43        s += '{}\n'.format(fielddisplay(self, 'sh_nmax', 'maximum spherical harmonic degree (default: 256, .35 deg, or 40 km at equator)'))
     44        s += '{}\n'.format(fielddisplay(self, 'sh_nmin', 'minimum spherical harmonic degree (default: 1)'))
     45        s += '{}\n'.format(fielddisplay(self, 'g0', 'adimensioning constant for gravity (default: 10) [m / s^2]'))
     46        s += '{}\n'.format(fielddisplay(self, 'r0', 'adimensioning constant for radius (default: 6378e3) [m]'))
     47        s += '{}\n'.format(fielddisplay(self, 'mu0', 'adimensioning constant for stress (default: 1.0e11) [Pa]'))
     48        s += '{}\n'.format(fielddisplay(self, 'allow_layer_deletion', 'allow for migration of the integration boundary with increasing spherical harmonics degree (default: 1)'))
     49        s += '{}\n'.format(fielddisplay(self, 'Gravitational_Constant', 'Newtonian constant of gravitation (default: 6.67259e-11 [m^3 kg^-1 s^-2])'))
     50        s += '{}\n'.format(fielddisplay(self, 'allow_layer_deletion', 'allow for migration of the integration boundary with increasing spherical harmonics degree (default: 1)'))
     51        s += '{}\n'.format(fielddisplay(self, 'underflow_tol', 'threshold of deep to surface love number ratio to trigger the deletion of layers (default: 1e-16)'))
     52        s += '{}\n'.format(fielddisplay(self, 'integration_steps_per_layer', 'number of radial steps to propagate the yi system from the bottom to the top of each layer (default: 100)'))
     53        s += '{}\n'.format(fielddisplay(self, 'istemporal', {'1 for time-dependent love numbers, 0 for frequency-dependent or elastic love numbers (default: 0)', 'If 1: use fourierlove function build_frequencies_from_time to meet consistency'}))
     54        s += '{}\n'.format(fielddisplay(self, 'n_temporal_iterations', 'max number of iterations in the inverse Laplace transform. Also the number of spectral samples per time step requested (default: 8)'))
     55        s += '{}\n'.format(fielddisplay(self, 'time', 'time vector for deformation if istemporal (default: 0) [s]'))
     56        s += '{}\n'.format(fielddisplay(self, 'love_kernels', 'compute love numbers at depth? (default: 0)'))
     57        s += '{}\n'.format(fielddisplay(self, 'forcing_type', 'integer indicating the nature and depth of the forcing for the Love number calculation (default: 11):'))
     58        s += '{}\n'.format('                                                     1:  Inner core boundary -- Volumic Potential')
     59        s += '{}\n'.format('                                                     2:  Inner core boundary --  Pressure')
     60        s += '{}\n'.format('                                                     3:  Inner core boundary --  Loading')
     61        s += '{}\n'.format('                                                     4:  Inner core boundary --  Tangential traction')
     62        s += '{}\n'.format('                                                     5:  Core mantle boundary --  Volumic Potential')
     63        s += '{}\n'.format('                                                     6:  Core mantle boundary --  Pressure')
     64        s += '{}\n'.format('                                                     7:  Core mantle boundary --  Loading')
     65        s += '{}\n'.format('                                                     8:  Core mantle boundary --  Tangential traction')
     66        s += '{}\n'.format('                                                     9:  Surface--  Volumic Potential')
     67        s += '{}\n'.format('                                                     10: Surface--  Pressure')
     68        s += '{}\n'.format('                                                     11: Surface--  Loading')
     69        s += '{}\n'.format('                                                     12: Surface--  Tangential traction ')
     70        s += '{}\n'.format(fielddisplay(self, 'inner_core_boundary', 'interface index in materials.radius locating forcing. Only used for forcing_type 1--4 (default: 1)'))
     71        s += '{}\n'.format(fielddisplay(self, 'core_mantle_boundary', 'interface index in materials.radius locating forcing. Only used for forcing_type 5--8 (default: 2)'))
    5572
    56         return string
     73        return s
    5774    #}}}
    5875
    59     def extrude(self, md):  # {{{
    60         return self
    61     #}}}
    62 
    63     def setdefaultparameters(self):  # {{{
    64         #we setup an elastic love number computation by default.
     76    def setdefaultparameters(self): #{{{
     77        # We setup an elastic love number computation by default
    6578        self.nfreq = 1
    66         self.frequencies = [0]  #Hz
    67         self.sh_nmax = 256  # .35 degree, 40 km at the equator.
     79        self.frequencies = [0] # Hz
     80        self.sh_nmax = 256 # .35 degree, 40 km at the equator
    6881        self.sh_nmin = 1
     82        # Work on matlab script for computing g0 for given Earth's structure
    6983        self.g0 = 9.81 # m/s^2
    70         self.r0 = 6371e3 #m
     84        self.r0 = 6371 * 1e3 # m
    7185        self.mu0 = 1e11 # Pa
     86        self.Gravitational_Constant = 6.67259e-11 # m^3 kg^-1 s^-2
    7287        self.allow_layer_deletion = 1
     88        self.underflow_tol = 1e-16 # Threshold of deep to surface love number ratio to trigger the deletion of layer
     89        self.integration_steps_per_layer = 100
     90        self.istemporal = 0
     91        self.n_temporal_iterations = 8
     92        self.time = [0] # s
    7393        self.love_kernels = 0
    74         self.forcing_type = 11
    75 
    76         return self
     94        self.forcing_type = 11 # Surface loading
     95        self.inner_core_boundary = 1
     96        self.core_mantle_boundary = 2
     97        self.complex_computation = 0
    7798    #}}}
    7899
    79     def checkconsistency(self, md, solution, analyses):  # {{{
     100    def checkconsistency(self, md, solution, analyses): #{{{
    80101        if 'LoveAnalysis' not in analyses:
    81102            return md
    82         md = checkfield(md, 'fieldname', 'love.nfreq', 'NaN', 1, 'Inf', 1, 'numel', [1], '>', 0)
    83         md = checkfield(md, 'fieldname', 'love.frequencies', 'NaN', 1, 'Inf', 1, 'numel', [md.love.nfreq])
    84         md = checkfield(md, 'fieldname', 'love.sh_nmax', 'NaN', 1, 'Inf', 1, 'numel', [1], '>', 0)
    85         md = checkfield(md, 'fieldname', 'love.sh_nmin', 'NaN', 1, 'Inf', 1, 'numel', [1], '>', 0)
    86         md = checkfield(md, 'fieldname', 'love.g0', 'NaN', 1, 'Inf', 1, 'numel', [1], '>', 0)
    87         md = checkfield(md, 'fieldname', 'love.r0', 'NaN', 1, 'Inf', 1, 'numel', [1], '>', 0)
    88         md = checkfield(md, 'fieldname', 'love.mu0', 'NaN', 1, 'Inf', 1, 'numel', [1], '>', 0)
     103
     104        md = checkfield(md, 'fieldname', 'love.nfreq', 'NaN', 1, 'Inf', 1, 'numel', 1, '>', 0)
     105        md = checkfield(md, 'fieldname', 'love.frequencies', 'NaN', 1, 'Inf', 1, 'numel', md.love.nfreq)
     106        md = checkfield(md, 'fieldname', 'love.sh_nmax', 'NaN', 1, 'Inf', 1, 'numel', 1, '>', 0)
     107        md = checkfield(md, 'fieldname', 'love.sh_nmin', 'NaN', 1, 'Inf', 1, 'numel', 1, '>', 0)
     108        md = checkfield(md, 'fieldname', 'love.g0', 'NaN', 1, 'Inf', 1, 'numel', 1, '>', 0)
     109        md = checkfield(md, 'fieldname', 'love.r0', 'NaN', 1, 'Inf', 1, 'numel', 1, '>', 0)
     110        md = checkfield(md, 'fieldname', 'love.mu0', 'NaN', 1, 'Inf', 1, 'numel', 1, '>', 0)
     111        md = checkfield(md, 'fieldname', 'love.Gravitational_Constant', 'NaN', 1, 'Inf', 1, 'numel', 1, '>', 0)
    89112        md = checkfield(md, 'fieldname', 'love.allow_layer_deletion', 'values', [0, 1])
     113        md = checkfield(md, 'fieldname', 'love.underflow_tol', 'NaN', 1, 'Inf', 1, 'numel', 1, '>', 0)
     114        md = checkfield(md, 'fieldname', 'love.integration_steps_per_layer', 'NaN', 1, 'Inf', 1, 'numel', 1, '>', 0)
    90115        md = checkfield(md, 'fieldname', 'love.love_kernels', 'values', [0, 1])
    91         md = checkfield(md, 'fieldname', 'love.forcing_type', 'NaN', 1, 'Inf', 1, 'numel', [1], '>', 0, '<=', 12)
    92         if md.love.sh_nmin <= 1 and md.love.forcing_type == 9:
    93             raise RuntimeError("Degree 1 not supported for Volumetric Potential forcing. Use sh_min >= 2 for this kind of calculation.")
     116        md = checkfield(md, 'fieldname', 'love.forcing_type', 'NaN', 1, 'Inf', 1, 'numel', 1, '>', 0, '<=', 12)
     117        md = checkfield(md, 'fieldname', 'love.complex_computation', 'NaN', 1, 'Inf', 1, 'numel', 1, 'values', [0, 1])
    94118
    95         # need 'litho' material
     119        md = checkfield(md, 'fieldname', 'love.istemporal', 'values', [0, 1])
     120        if md.love.istemporal:
     121            md = checkfield(md, 'fieldname', 'love.n_temporal_iterations', 'NaN', 1, 'Inf', 1, 'numel', 1, '>', 0)
     122            md = checkfield(md, 'fieldname', 'love.time', 'NaN', 1, 'Inf', 1, 'numel', md.love.nfreq / 2 / md.love.n_temporal_iterations)
     123        if md.love.sh_nmin <= 1 and (md.love.forcing_type == 1 or md.love.forcing_type == 5 or md.love.forcing_type == 9):
     124            raise RuntimeError('Degree 1 not supported for forcing type {}. Use sh_min >= 2 for this kind of calculation.'.format(md.love.forcing_type))
     125
     126        # Need 'litho' material
    96127        if not hasattr(md.materials, 'materials') or 'litho' not in md.materials.nature:
    97128            raise RuntimeError('Need a \'litho\' material to run a Fourier Love number analysis')
     129
     130        mat = np.where(np.array(nature) == 'litho')[0]
     131        if md.love.forcing_type <= 4:
     132            md = checkfield(md, 'fieldname', 'love.inner_core_boundary', 'NaN', 1, 'Inf', 1, 'numel', 1, '>', 0, '<=', md.materials[mat].numlayers)
     133        elif md.love.forcing_type <= 8:
     134            md = checkfield(md, 'fieldname', 'love.core_mantle_boundary', 'NaN', 1, 'Inf', 1, 'numel', 1, '>', 0, '<=', md.materials[mat].numlayers)
     135
    98136        return md
    99     # }}}
     137    #}}}
    100138
    101     def marshall(self, prefix, md, fid):  # {{{
     139    def marshall(self, prefix, md, fid): #{{{
    102140        WriteData(fid, prefix, 'object', self, 'fieldname', 'nfreq', 'format', 'Integer')
    103         WriteData(fid, prefix, 'object', self, 'fieldname', 'frequencies', 'format', 'DoubleMat', 'mattype', 3)
     141        WriteData(fid, prefix, 'object', self, 'fieldname', 'frequencies', 'format', 'DoubleMat', 'mattype',3)
    104142        WriteData(fid, prefix, 'object', self, 'fieldname', 'sh_nmax', 'format', 'Integer')
    105143        WriteData(fid, prefix, 'object', self, 'fieldname', 'sh_nmin', 'format', 'Integer')
    106144        WriteData(fid, prefix, 'object', self, 'fieldname', 'g0', 'format', 'Double')
    107145        WriteData(fid, prefix, 'object', self, 'fieldname', 'r0', 'format', 'Double')
    108146        WriteData(fid, prefix, 'object', self, 'fieldname', 'mu0', 'format', 'Double')
     147        WriteData(fid, prefix, 'object', self, 'fieldname', 'Gravitational_Constant', 'format', 'Double')
    109148        WriteData(fid, prefix, 'object', self, 'fieldname', 'allow_layer_deletion', 'format', 'Boolean')
     149        WriteData(fid, prefix, 'object', self, 'fieldname', 'underflow_tol', 'format', 'Double')
     150        WriteData(fid, prefix, 'object', self, 'fieldname', 'integration_steps_per_layer', 'format', 'Integer')
     151        WriteData(fid, prefix, 'object', self, 'fieldname', 'istemporal', 'format', 'Boolean')
     152        WriteData(fid, prefix, 'object', self, 'fieldname', 'n_temporal_iterations', 'format', 'Integer')
     153        WriteData(fid, prefix, 'object', self, 'fieldname', 'complex_computation', 'format', 'Boolean')
     154        # Note: no need to marshall the time vector, we have frequencies
    110155        WriteData(fid, prefix, 'object', self, 'fieldname', 'love_kernels', 'format', 'Boolean')
    111156        WriteData(fid, prefix, 'object', self, 'fieldname', 'forcing_type', 'format', 'Integer')
    112     # }}}
     157        WriteData(fid, prefix, 'object', self, 'fieldname', 'inner_core_boundary', 'format', 'Integer')
     158        WriteData(fid, prefix, 'object', self, 'fieldname', 'core_mantle_boundary', 'format', 'Integer')
     159    #}}}
     160
     161    def extrude(self, md): #{{{
     162        return self
     163    #}}}
     164
     165    def build_frequencies_from_time(self): #{{{
     166        if not self.istemporal:
     167            raise RuntimeError('cannot build frequencies for temporal love numbers if love.istemporal==0')
     168        print('Temporal love numbers: Overriding md.love.nfreq and md.love.frequencies')
     169        self.nfreq = len(self.time) * 2 * self.n_temporal_iterations
     170        self.frequencies = np.zeros((self.nfreq, 1))
     171        for i in range(len(self.time)):
     172            for j in range(2 * self.n_temporal_iterations):
     173                self.frequencies[(i - 1) * 2 * self.n_temporal_iterations + j] = j * np.log(2) / self.time[i] / 2 / np.pi
     174    #}}}

  • ../trunk-jpl/src/m/classes/geometry.py

     
    3636    #}}}
    3737
    3838    def setdefaultparameters(self): #{{{
    39         return self
     39        return
    4040    #}}}
    4141
    4242    def checkconsistency(self, md, solution, analyses): #{{{

  • ../trunk-jpl/src/m/classes/hydrologyshreve.py

     
     1import numpy as np
     2
    13from fielddisplay import fielddisplay
    24from checkfield import checkfield
    35from WriteData import WriteData
    … …  
    1012        hydrologyshreve = hydrologyshreve()
    1113    """
    1214
    13     def __init__(self):  # {{{
    14         self.spcwatercolumn = float('NaN')
     15    def __init__(self, *args): #{{{
     16        self.spcwatercolumn = np.nan
    1517        self.stabilization = 0
    1618        self.requested_outputs = []
    1719
    18         self.setdefaultparameters()
     20        nargs = len(args)
     21        if nargs == 0:
     22            self.setdefaultparameters()
     23        elif nargs == 1:
     24            self.setdefaultparameters(args)
     25        else:
     26            raise RuntimeError('constructor not supported')
    1927    #}}}
    2028
    21     def __repr__(self):  # {{{
    22         # TODO:
    23         # - Convert all formatting to calls to <string>.format (see any
    24         #   already converted <class>.__repr__ method for examples)
    25         #
    26         string = '   hydrologyshreve solution parameters:'
    27         string = "%s\n%s" % (string, fielddisplay(self, 'spcwatercolumn', 'water thickness constraints (NaN means no constraint) [m]'))
    28         string = "%s\n%s" % (string, fielddisplay(self, 'stabilization', 'artificial diffusivity (default is 1). can be more than 1 to increase diffusivity.'))
    29         string = "%s\n%s" % (string, fielddisplay(self, 'requested_outputs', 'additional outputs requested'))
    30         return string
     29    def __repr__(self): #{{{
     30        s = '   hydrologyshreve solution parameters:\n'
     31        s += '{}\n'.format(fielddisplay(self, 'spcwatercolumn', 'water thickness constraints (NaN means no constraint) [m]'))
     32        s += '{}\n'.format(fielddisplay(self, 'stabilization', 'artificial diffusivity (default: 1). can be more than 1 to increase diffusivity.'))
     33        s += '{}\n'.format(fielddisplay(self, 'requested_outputs', 'additional outputs requested'))
     34        return s
    3135    #}}}
    3236
    33     def extrude(self, md):  # {{{
    34         return self
    35     #}}}
    36 
    37     def setdefaultparameters(self):  # {{{
    38         #Type of stabilization to use 0:nothing 1:artificial_diffusivity
     37    def setdefaultparameters(self): #{{{
     38        # Type of stabilization to use 0:nothing 1:artificial_diffusivity
    3939        self.stabilization = 1
    4040        self.requested_outputs = ['default']
    41         return self
    4241    #}}}
    4342
    44     def defaultoutputs(self, md):  # {{{
    45         list = ['Watercolumn', 'HydrologyWaterVx', 'HydrologyWaterVy']
    46         return list
    47     #}}}
    48 
    49     def checkconsistency(self, md, solution, analyses):  # {{{
     43    def checkconsistency(self, md, solution, analyses): #{{{
    5044        #Early return
    5145        if 'HydrologyShreveAnalysis' not in analyses or (solution == 'TransientSolution' and not md.transient.ishydrology):
    5246            return md
    … …  
    5347
    5448        md = checkfield(md, 'fieldname', 'hydrology.spcwatercolumn', 'Inf', 1, 'timeseries', 1)
    5549        md = checkfield(md, 'fieldname', 'hydrology.stabilization', '>=', 0)
    56         md = checkfield(md, 'fieldname', 'hydrology.requested_outputs', 'stringrow', 1)
    5750        return md
    5851    # }}}
    5952
    60     def marshall(self, prefix, md, fid):  # {{{
     53    def defaultoutputs(self, md): #{{{
     54        return ['Watercolumn', 'HydrologyWaterVx', 'HydrologyWaterVy']
     55    #}}}
     56
     57    def marshall(self, prefix, md, fid): #{{{
    6158        WriteData(fid, prefix, 'name', 'md.hydrology.model', 'data', 2, 'format', 'Integer')
    6259        WriteData(fid, prefix, 'object', self, 'fieldname', 'spcwatercolumn', 'format', 'DoubleMat', 'mattype', 1, 'timeserieslength', md.mesh.numberofvertices + 1, 'yts', md.constants.yts)
    6360        WriteData(fid, prefix, 'object', self, 'fieldname', 'stabilization', 'format', 'Double')
    64         #process requested outputs
    6561        outputs = self.requested_outputs
    6662        indices = [i for i, x in enumerate(outputs) if x == 'default']
    67         if len(indices) > 0:
     63        if len(indices):
    6864            outputscopy = outputs[0:max(0, indices[0] - 1)] + self.defaultoutputs(md) + outputs[indices[0] + 1:]
    6965            outputs = outputscopy
    7066        WriteData(fid, prefix, 'data', outputs, 'name', 'md.hydrology.requested_outputs', 'format', 'StringArray')
     67    # }}}
    7168
    72     # }}}
     69    def extrude(self, md): #{{{
     70        return self
     71    #}}}

  • ../trunk-jpl/src/m/classes/initialization.m

     
    2626                sample              = NaN;
    2727        end
    2828        methods
    29                 function self = extrude(self,md) % {{{
    30                         self.vx=project3d(md,'vector',self.vx,'type','node');
    31                         self.vy=project3d(md,'vector',self.vy,'type','node');
    32                         self.vz=project3d(md,'vector',self.vz,'type','node');
    33                         self.vel=project3d(md,'vector',self.vel,'type','node');
    34                         self.temperature=project3d(md,'vector',self.temperature,'type','node');
    35                         self.enthalpy=project3d(md,'vector',self.enthalpy,'type','node');
    36                         self.waterfraction=project3d(md,'vector',self.waterfraction,'type','node');
    37                         self.watercolumn=project3d(md,'vector',self.watercolumn,'type','node','layer',1);
    38                         self.sediment_head=project3d(md,'vector',self.sediment_head,'type','node','layer',1);
    39                         self.epl_head=project3d(md,'vector',self.epl_head,'type','node','layer',1);
    40                         self.epl_thickness=project3d(md,'vector',self.epl_thickness,'type','node','layer',1);
    41                         self.sealevel=project3d(md,'vector',self.sealevel,'type','node','layer',1);
    42                         self.bottompressure=project3d(md,'vector',self.bottompressure,'type','node','layer',1);
    43                         self.dsl=project3d(md,'vector',self.dsl,'type','node','layer',1);
    44                         self.str=project3d(md,'vector',self.str,'type','node','layer',1);
    45 
    46                         %Lithostatic pressure by default
    47                         self.pressure=md.constants.g*md.materials.rho_ice*(md.geometry.surface-md.mesh.z);
    48                 end % }}}
    4929                function self = initialization(varargin) % {{{
    5030                        switch nargin
    5131                                case 0
    … …  
    200180                                WriteData(fid,prefix,'data',enthalpy,'format','DoubleMat','mattype',1,'name','md.initialization.enthalpy');
    201181                        end
    202182                end % }}}
     183                function self = extrude(self,md) % {{{
     184                        self.vx=project3d(md,'vector',self.vx,'type','node');
     185                        self.vy=project3d(md,'vector',self.vy,'type','node');
     186                        self.vz=project3d(md,'vector',self.vz,'type','node');
     187                        self.vel=project3d(md,'vector',self.vel,'type','node');
     188                        self.temperature=project3d(md,'vector',self.temperature,'type','node');
     189                        self.enthalpy=project3d(md,'vector',self.enthalpy,'type','node');
     190                        self.waterfraction=project3d(md,'vector',self.waterfraction,'type','node');
     191                        self.watercolumn=project3d(md,'vector',self.watercolumn,'type','node','layer',1);
     192                        self.sediment_head=project3d(md,'vector',self.sediment_head,'type','node','layer',1);
     193                        self.epl_head=project3d(md,'vector',self.epl_head,'type','node','layer',1);
     194                        self.epl_thickness=project3d(md,'vector',self.epl_thickness,'type','node','layer',1);
     195                        self.sealevel=project3d(md,'vector',self.sealevel,'type','node','layer',1);
     196                        self.bottompressure=project3d(md,'vector',self.bottompressure,'type','node','layer',1);
     197                        self.dsl=project3d(md,'vector',self.dsl,'type','node','layer',1);
     198                        self.str=project3d(md,'vector',self.str,'type','node','layer',1);
     199
     200                        %Lithostatic pressure by default
     201                        self.pressure=md.constants.g*md.materials.rho_ice*(md.geometry.surface-md.mesh.z);
     202                end % }}}
    203203                function savemodeljs(self,fid,modelname) % {{{
    204204
    205205                        writejs1Darray(fid,[modelname '.initialization.vx'],self.vx);

  • ../trunk-jpl/src/m/classes/lovenumbers.m

     
    1515                l           = []; %idem
    1616               
    1717                %tidal love numbers for computing rotational feedback:
    18                 th          = []; 
    19                 tk          = []; 
    20                 tl          = []; 
    21                 tk2secular  = 0;  %deg 2 secular number.
     18                th          = [];
     19                tk          = [];
     20                tl          = [];
     21                tk2secular  = 0; %deg 2 secular number.
    2222
    2323                %time/frequency for visco-elastic love numbers
    2424                timefreq    = [];
    25                 istime      = 1; 
     25                istime      = 1;
    2626
    2727        end
    2828        methods
    … …  
    3232                        referenceframe=getfieldvalue(options,'referenceframe','CM');
    3333                        self=setdefaultparameters(self,maxdeg,referenceframe);
    3434                end % }}}
     35                function disp(self) % {{{
     36                        disp(sprintf('   lovenumbers parameters:'));
     37
     38                        fielddisplay(self,'h','load Love number for radial displacement');
     39                        fielddisplay(self,'k','load Love number for gravitational potential perturbation');
     40                        fielddisplay(self,'l','load Love number for horizontal displacements');
     41
     42                        fielddisplay(self,'th','tidal load Love number (deg 2)');
     43                        fielddisplay(self,'tk','tidal load Love number (deg 2)');
     44                        fielddisplay(self,'tl','tidal load Love number (deg 2)');
     45                        fielddisplay(self,'tk2secular','secular fluid Love number');
     46
     47                        fielddisplay(self,'istime','time (default: 1) or frequency love numbers (0)');
     48                        fielddisplay(self,'timefreq','time/frequency vector (yr or 1/yr)');
     49
     50                end % }}}
    3551                function self = setdefaultparameters(self,maxdeg,referenceframe) % {{{
    36                
     52
    3753                        %initialize love numbers:
    3854                        self.h=getlovenumbers('type','loadingverticaldisplacement','referenceframe',referenceframe,'maxdeg',maxdeg);
    3955                        self.k=getlovenumbers('type','loadinggravitationalpotential','referenceframe',referenceframe,'maxdeg',maxdeg);
    … …  
    4460
    4561                        %secular fluid love number:
    4662                        self.tk2secular=0.942;
    47                        
     63
    4864                        %time:
    4965                        self.istime=1; %temporal love numbers by default
    5066                        self.timefreq=0; %elastic case by default.
    … …  
    7288                        if (size(self.h,1)~=size(self.k,1) | size(self.h,1)~=size(self.l,1)),
    7389                                error('lovenumbers error message: love numbers should be provided at the same level of accuracy');
    7490                        end
    75                        
     91
    7692                        ntf=length(self.timefreq);
    7793                        if( size(self.h,2) ~= ntf | size(self.k,2) ~= ntf | size(self.l,2) ~= ntf | size(self.th,2) ~= ntf | size(self.tk,2) ~= ntf | size(self.tl,2) ~= ntf ),
    7894                                error('lovenumbers error message: love numbers should have as many time/frequency steps as the time/frequency vector');
    … …  
    8298                function list=defaultoutputs(self,md) % {{{
    8399                        list = {};
    84100                end % }}}
    85                 function disp(self) % {{{
    86                         disp(sprintf('   lovenumbers parameters:'));
    87 
    88                         fielddisplay(self,'h','load Love number for radial displacement');
    89                         fielddisplay(self,'k','load Love number for gravitational potential perturbation');
    90                         fielddisplay(self,'l','load Love number for horizontal displacements');
    91 
    92                         fielddisplay(self,'th','tidal load Love number (deg 2)');
    93                         fielddisplay(self,'tk','tidal load Love number (deg 2)');
    94                         fielddisplay(self,'tl','tidal load Love number (deg 2)');
    95                         fielddisplay(self,'tk2secular','secular fluid Love number');
    96                        
    97                         fielddisplay(self,'istime','time (default=1) or frequency love numbers (0)');
    98                         fielddisplay(self,'timefreq','time/frequency vector (yr or 1/yr)');
    99 
    100                 end % }}}
    101101                function marshall(self,prefix,md,fid) % {{{
    102102
    103103                        WriteData(fid,prefix,'object',self,'fieldname','h','name','md.solidearth.lovenumbers.h','format','DoubleMat','mattype',1);

  • ../trunk-jpl/src/m/classes/materials.m

     
    203203                                        fielddisplay(self,'ebm_alpha','array describing each layer''s exponent parameter controlling the shape of shear modulus curve between taul and tauh, only for EBM rheology (numlayers)');
    204204                                        fielddisplay(self,'ebm_delta','array describing each layer''s amplitude of the transient relaxation (ratio between elastic rigity to pre-maxwell relaxation rigity), only for EBM rheology (numlayers)');
    205205                                        fielddisplay(self,'ebm_taul','array describing each layer''s starting period for transient relaxation, only for EBM rheology  (numlayers) [s]');
    206                                         fielddisplay(self,'ebm_tauh','array describing each layer''s array describing each layer''s end period for transient relaxation, only for Burgers rheology  (numlayers) [s]');
     206                                        fielddisplay(self,'ebm_tauh','array describing each layer''s array describing each layer''s end period for transient relaxation, only for Burgers rheology (numlayers) [s]');
    207207
    208208
    209209                                        fielddisplay(self,'rheologymodel','array describing whether we adopt a Maxwell (0), Burgers (1) or EBM (2) rheology (default: 0)');
    … …  
    255255                                                if md.materials.rheologymodel(i)==1 & (isnan(md.materials.burgers_viscosity(i) | isnan(md.materials.burgers_mu(i)))),
    256256                                                        error('materials checkconsistency error message: Litho burgers_viscosity or burgers_mu has NaN values, inconsistent with rheologymodel choice');
    257257                                                end
    258                                                 if md.materials.rheologymodel(i)==2 & (isnan(md.materials.ebm_alpha(i)) | isnan(md.materials.ebm_delta(i)) | isnan(md.materials.ebm_taul(i)) | isnan(md.materials.ebm_tauh(i)) ),
     258                                                if md.materials.rheologymodel(i)==2 & (isnan(md.materials.ebm_alpha(i)) | isnan(md.materials.ebm_delta(i)) | isnan(md.materials.ebm_taul(i)) | isnan(md.materials.ebm_tauh(i))),
    259259                                                        error('materials checkconsistency error message: Litho ebm_alpha, ebm_delta, ebm_taul or ebm_tauh has NaN values, inconsistent with rheologymodel choice');
    260260                                                end
    261261                                        end
    … …  
    283283
    284284                        %1: MatdamageiceEnum 2: MatestarEnum 3: MaticeEnum 4: MatenhancediceEnum 5: MaterialsEnum
    285285                        WriteData(fid,prefix,'name','md.materials.nature','data',naturetointeger(self.nature),'format','IntMat','mattype',3);
    286                         WriteData(fid,prefix,'name','md.materials.type','data',5,'format','Integer'); %DANGER, this can evolve if you have classes.
     286                        WriteData(fid,prefix,'name','md.materials.type','data',5,'format','Integer'); %DANGER: this can evolve if you have classes.
    287287                        for i=1:length(self.nature),
    288288                                nat=self.nature{i};
    289289                                switch nat
    … …  
    325325                                                earth_density=earth_density + (self.radius(i+1)^3-self.radius(i)^3)*self.density(i);
    326326                                        end
    327327                                        earth_density=earth_density/self.radius(self.numlayers+1)^3;
     328                                        disp(earth_density)
    328329                                        self.earth_density=earth_density;
    329330                                case 'hydro'
    330331                                        WriteData(fid,prefix,'object',self,'class','materials','fieldname','rho_ice','format','Double');
    … …  
    540541                                t4 = s(i,4)/((ra^3)*6.);
    541542                                vs =  vs + t1*(r2^3) + t2*(r2^4) + t3*(r2^5) + t4*(r2^6) - ( t1*(r1^3) + t2*(r1^4) + t3*(r1^5) + t4*(r1^6) );
    542543
    543                         end 
     544                        end
    544545                        ro = ro*3 / (rad(j+1)^3-rad(j)^3);
    545546                        vp = vp*3 /(rad(j+1)^3-rad(j)^3);
    546547                        vs = vs*3 / (rad(j+1)^3-rad(j)^3);
    … …  
    556557                  end
    557558
    558559                end % }}}
    559    
    560560        end
    561561end
    562562
    … …  
    583583                end
    584584        end
    585585end % }}}
    586 
    587    

  • ../trunk-jpl/src/m/classes/model.m

     
    3838                thermal          = 0;
    3939                steadystate      = 0;
    4040                transient        = 0;
    41                 levelset             = 0;
     41                levelset         = 0;
    4242                calving          = 0;
    4343                frontalforcings  = 0;
    44                 love                    = 0;
     44                love            = 0;
    4545                esa              = 0;
    4646                sampling         = 0;
    4747
    … …  
    4848                autodiff         = 0;
    4949                inversion        = 0;
    5050                qmu              = 0;
    51                 amr                              = 0;
     51                amr              = 0;
    5252                results          = 0;
    5353                outputdefinition = 0;
    5454                radaroverlay     = 0;
    … …  
    206206
    207207                end
    208208                %}}}
     209                function disp(self) % {{{
     210                        disp(sprintf('%19s: %-22s -- %s','mesh'            ,['[1x1 ' class(self.mesh) ']'],'mesh properties'));
     211                        disp(sprintf('%19s: %-22s -- %s','mask'            ,['[1x1 ' class(self.mask) ']'],'defines grounded and floating elements'));
     212                        disp(sprintf('%19s: %-22s -- %s','geometry'        ,['[1x1 ' class(self.geometry) ']'],'surface elevation, bedrock topography, ice thickness,...'));
     213                        disp(sprintf('%19s: %-22s -- %s','constants'       ,['[1x1 ' class(self.constants) ']'],'physical constants'));
     214                        disp(sprintf('%19s: %-22s -- %s','smb'             ,['[1x1 ' class(self.smb) ']'],'surface mass balance'));
     215                        disp(sprintf('%19s: %-22s -- %s','basalforcings'   ,['[1x1 ' class(self.basalforcings) ']'],'bed forcings'));
     216                        disp(sprintf('%19s: %-22s -- %s','materials'       ,['[1x1 ' class(self.materials) ']'],'material properties'));
     217                        disp(sprintf('%19s: %-22s -- %s','damage'          ,['[1x1 ' class(self.damage) ']'],'parameters for damage evolution solution'));
     218                        disp(sprintf('%19s: %-22s -- %s','friction'        ,['[1x1 ' class(self.friction) ']'],'basal friction/drag properties'));
     219                        disp(sprintf('%19s: %-22s -- %s','flowequation'    ,['[1x1 ' class(self.flowequation) ']'],'flow equations'));
     220                        disp(sprintf('%19s: %-22s -- %s','timestepping'    ,['[1x1 ' class(self.timestepping) ']'],'time stepping for transient models'));
     221                        disp(sprintf('%19s: %-22s -- %s','initialization'  ,['[1x1 ' class(self.initialization) ']'],'initial guess/state'));
     222                        disp(sprintf('%19s: %-22s -- %s','rifts'           ,['[1x1 ' class(self.rifts) ']'],'rifts properties'));
     223                        disp(sprintf('%19s: %-22s -- %s','solidearth'      ,['[1x1 ' class(self.solidearth) ']'],'solidearth inputs and settings'));
     224                        disp(sprintf('%19s: %-22s -- %s','dsl'             ,['[1x1 ' class(self.dsl) ']'],'dynamic sea-level '));
     225                        disp(sprintf('%19s: %-22s -- %s','debug'           ,['[1x1 ' class(self.debug) ']'],'debugging tools (valgrind, gprof)'));
     226                        disp(sprintf('%19s: %-22s -- %s','verbose'         ,['[1x1 ' class(self.verbose) ']'],'verbosity level in solve'));
     227                        disp(sprintf('%19s: %-22s -- %s','settings'        ,['[1x1 ' class(self.settings) ']'],'settings properties'));
     228                        disp(sprintf('%19s: %-22s -- %s','toolkits'        ,['[1x1 ' class(self.toolkits) ']'],'PETSc options for each solution'));
     229                        disp(sprintf('%19s: %-22s -- %s','cluster'         ,['[1x1 ' class(self.cluster) ']'],'cluster parameters (number of CPUs...)'));
     230                        disp(sprintf('%19s: %-22s -- %s','balancethickness',['[1x1 ' class(self.balancethickness) ']'],'parameters for balancethickness solution'));
     231                        disp(sprintf('%19s: %-22s -- %s','stressbalance'   ,['[1x1 ' class(self.stressbalance) ']'],'parameters for stressbalance solution'));
     232                        disp(sprintf('%19s: %-22s -- %s','groundingline'   ,['[1x1 ' class(self.groundingline) ']'],'parameters for groundingline solution'));
     233                        disp(sprintf('%19s: %-22s -- %s','hydrology'       ,['[1x1 ' class(self.hydrology) ']'],'parameters for hydrology solution'));
     234                        disp(sprintf('%19s: %-22s -- %s','masstransport'   ,['[1x1 ' class(self.masstransport) ']'],'parameters for masstransport solution'));
     235                        disp(sprintf('%19s: %-22s -- %s','thermal'         ,['[1x1 ' class(self.thermal) ']'],'parameters for thermal solution'));
     236                        disp(sprintf('%19s: %-22s -- %s','steadystate'     ,['[1x1 ' class(self.steadystate) ']'],'parameters for steadystate solution'));
     237                        disp(sprintf('%19s: %-22s -- %s','transient'       ,['[1x1 ' class(self.transient) ']'],'parHwoameters for transient solution'));
     238                        disp(sprintf('%19s: %-22s -- %s','levelset'        ,['[1x1 ' class(self.levelset) ']'],'parameters for moving boundaries (level-set method)'));
     239                        disp(sprintf('%19s: %-22s -- %s','calving'         ,['[1x1 ' class(self.calving) ']'],'parameters for calving'));
     240                        disp(sprintf('%19s: %-22s -- %s','frontalforcings' ,['[1x1 ' class(self.frontalforcings) ']'],'parameters for frontalforcings'));
     241                        disp(sprintf('%19s: %-22s -- %s','esa'             ,['[1x1 ' class(self.esa) ']'],'parameters for elastic adjustment solution'));
     242                        disp(sprintf('%19s: %-22s -- %s','love'            ,['[1x1 ' class(self.love) ']'],'parameters for love solution'));
     243                        disp(sprintf('%19s: %-22s -- %s','sampling'        ,['[1x1 ' class(self.sampling) ']'],'parameters for stochastic sampler'));
     244                        disp(sprintf('%19s: %-22s -- %s','autodiff'        ,['[1x1 ' class(self.autodiff) ']'],'automatic differentiation parameters'));
     245                        disp(sprintf('%19s: %-22s -- %s','inversion'       ,['[1x1 ' class(self.inversion) ']'],'parameters for inverse methods'));
     246                        disp(sprintf('%19s: %-22s -- %s','qmu'             ,['[1x1 ' class(self.qmu) ']'],'Dakota properties'));
     247                        disp(sprintf('%19s: %-22s -- %s','amr'             ,['[1x1 ' class(self.amr) ']'],'adaptive mesh refinement properties'));
     248                        disp(sprintf('%19s: %-22s -- %s','outputdefinition',['[1x1 ' class(self.outputdefinition) ']'],'output definition'));
     249                        disp(sprintf('%19s: %-22s -- %s','results'         ,['[1x1 ' class(self.results) ']'],'model results'));
     250                        disp(sprintf('%19s: %-22s -- %s','radaroverlay'    ,['[1x1 ' class(self.radaroverlay) ']'],'radar image for plot overlay'));
     251                        disp(sprintf('%19s: %-22s -- %s','miscellaneous'   ,['[1x1 ' class(self.miscellaneous) ']'],'miscellaneous fields'));
     252                end % }}}
    209253                function md = setdefaultparameters(md,planet) % {{{
    210254
    211255                        %initialize subclasses
    … …  
    11321176                        if md.mesh.average_vertex_connectivity<=25,
    11331177                                md.mesh.average_vertex_connectivity=100;
    11341178                        end
    1135                         end % }}}
     1179                end % }}}
    11361180                function md = structtomodel(md,structmd) % {{{
    11371181
    11381182                        if ~isstruct(structmd) error('input model is not a structure'); end
    … …  
    15831627                                md.mesh.average_vertex_connectivity=100;
    15841628                        end
    15851629                end % }}}
    1586                 function disp(self) % {{{
    1587                         disp(sprintf('%19s: %-22s -- %s','mesh'            ,['[1x1 ' class(self.mesh) ']'],'mesh properties'));
    1588                         disp(sprintf('%19s: %-22s -- %s','mask'            ,['[1x1 ' class(self.mask) ']'],'defines grounded and floating elements'));
    1589                         disp(sprintf('%19s: %-22s -- %s','geometry'        ,['[1x1 ' class(self.geometry) ']'],'surface elevation, bedrock topography, ice thickness,...'));
    1590                         disp(sprintf('%19s: %-22s -- %s','constants'       ,['[1x1 ' class(self.constants) ']'],'physical constants'));
    1591                         disp(sprintf('%19s: %-22s -- %s','smb'             ,['[1x1 ' class(self.smb) ']'],'surface mass balance'));
    1592                         disp(sprintf('%19s: %-22s -- %s','basalforcings'   ,['[1x1 ' class(self.basalforcings) ']'],'bed forcings'));
    1593                         disp(sprintf('%19s: %-22s -- %s','materials'       ,['[1x1 ' class(self.materials) ']'],'material properties'));
    1594                         disp(sprintf('%19s: %-22s -- %s','damage'          ,['[1x1 ' class(self.damage) ']'],'parameters for damage evolution solution'));
    1595                         disp(sprintf('%19s: %-22s -- %s','friction'        ,['[1x1 ' class(self.friction) ']'],'basal friction/drag properties'));
    1596                         disp(sprintf('%19s: %-22s -- %s','flowequation'    ,['[1x1 ' class(self.flowequation) ']'],'flow equations'));
    1597                         disp(sprintf('%19s: %-22s -- %s','timestepping'    ,['[1x1 ' class(self.timestepping) ']'],'time stepping for transient models'));
    1598                         disp(sprintf('%19s: %-22s -- %s','initialization'  ,['[1x1 ' class(self.initialization) ']'],'initial guess/state'));
    1599                         disp(sprintf('%19s: %-22s -- %s','rifts'           ,['[1x1 ' class(self.rifts) ']'],'rifts properties'));
    1600                         disp(sprintf('%19s: %-22s -- %s','solidearth'      ,['[1x1 ' class(self.solidearth) ']'],'solidearth inputs and settings'));
    1601                         disp(sprintf('%19s: %-22s -- %s','dsl'             ,['[1x1 ' class(self.dsl) ']'],'dynamic sea-level '));
    1602                         disp(sprintf('%19s: %-22s -- %s','debug'           ,['[1x1 ' class(self.debug) ']'],'debugging tools (valgrind, gprof)'));
    1603                         disp(sprintf('%19s: %-22s -- %s','verbose'         ,['[1x1 ' class(self.verbose) ']'],'verbosity level in solve'));
    1604                         disp(sprintf('%19s: %-22s -- %s','settings'        ,['[1x1 ' class(self.settings) ']'],'settings properties'));
    1605                         disp(sprintf('%19s: %-22s -- %s','toolkits'        ,['[1x1 ' class(self.toolkits) ']'],'PETSc options for each solution'));
    1606                         disp(sprintf('%19s: %-22s -- %s','cluster'         ,['[1x1 ' class(self.cluster) ']'],'cluster parameters (number of CPUs...)'));
    1607                         disp(sprintf('%19s: %-22s -- %s','balancethickness',['[1x1 ' class(self.balancethickness) ']'],'parameters for balancethickness solution'));
    1608                         disp(sprintf('%19s: %-22s -- %s','stressbalance'   ,['[1x1 ' class(self.stressbalance) ']'],'parameters for stressbalance solution'));
    1609                         disp(sprintf('%19s: %-22s -- %s','groundingline'   ,['[1x1 ' class(self.groundingline) ']'],'parameters for groundingline solution'));
    1610                         disp(sprintf('%19s: %-22s -- %s','hydrology'       ,['[1x1 ' class(self.hydrology) ']'],'parameters for hydrology solution'));
    1611                         disp(sprintf('%19s: %-22s -- %s','masstransport'   ,['[1x1 ' class(self.masstransport) ']'],'parameters for masstransport solution'));
    1612                         disp(sprintf('%19s: %-22s -- %s','thermal'         ,['[1x1 ' class(self.thermal) ']'],'parameters for thermal solution'));
    1613                         disp(sprintf('%19s: %-22s -- %s','steadystate'     ,['[1x1 ' class(self.steadystate) ']'],'parameters for steadystate solution'));
    1614                         disp(sprintf('%19s: %-22s -- %s','transient'       ,['[1x1 ' class(self.transient) ']'],'parHwoameters for transient solution'));
    1615                         disp(sprintf('%19s: %-22s -- %s','levelset'        ,['[1x1 ' class(self.levelset) ']'],'parameters for moving boundaries (level-set method)'));
    1616                         disp(sprintf('%19s: %-22s -- %s','calving'         ,['[1x1 ' class(self.calving) ']'],'parameters for calving'));
    1617                         disp(sprintf('%19s: %-22s -- %s','frontalforcings' ,['[1x1 ' class(self.frontalforcings) ']'],'parameters for frontalforcings'));
    1618                         disp(sprintf('%19s: %-22s -- %s','esa'             ,['[1x1 ' class(self.esa) ']'],'parameters for elastic adjustment solution'));
    1619                         disp(sprintf('%19s: %-22s -- %s','love'            ,['[1x1 ' class(self.love) ']'],'parameters for love solution'));
    1620                         disp(sprintf('%19s: %-22s -- %s','sampling'        ,['[1x1 ' class(self.sampling) ']'],'parameters for stochastic sampler'));
    1621                         disp(sprintf('%19s: %-22s -- %s','autodiff'        ,['[1x1 ' class(self.autodiff) ']'],'automatic differentiation parameters'));
    1622                         disp(sprintf('%19s: %-22s -- %s','inversion'       ,['[1x1 ' class(self.inversion) ']'],'parameters for inverse methods'));
    1623                         disp(sprintf('%19s: %-22s -- %s','qmu'             ,['[1x1 ' class(self.qmu) ']'],'Dakota properties'));
    1624                         disp(sprintf('%19s: %-22s -- %s','amr'             ,['[1x1 ' class(self.amr) ']'],'adaptive mesh refinement properties'));
    1625                         disp(sprintf('%19s: %-22s -- %s','outputdefinition',['[1x1 ' class(self.outputdefinition) ']'],'output definition'));
    1626                         disp(sprintf('%19s: %-22s -- %s','results'         ,['[1x1 ' class(self.results) ']'],'model results'));
    1627                         disp(sprintf('%19s: %-22s -- %s','radaroverlay'    ,['[1x1 ' class(self.radaroverlay) ']'],'radar image for plot overlay'));
    1628                         disp(sprintf('%19s: %-22s -- %s','miscellaneous'   ,['[1x1 ' class(self.miscellaneous) ']'],'miscellaneous fields'));
    1629                 end % }}}
    16301630                function memory(self) % {{{
    16311631
    1632                 disp(sprintf('\nMemory imprint:\n'));
     1632                        disp(sprintf('\nMemory imprint:\n'));
    16331633
    1634                 fields=properties('model');
    1635                 mem=0;
     1634                        fields=properties('model');
     1635                        mem=0;
    16361636
    1637                 for i=1:length(fields),
    1638                         field=self.(fields{i});
    1639                         s=whos('field');
    1640                         mem=mem+s.bytes/1e6;
    1641                         disp(sprintf('%19s: %6.2f Mb',fields{i},s.bytes/1e6));
     1637                        for i=1:length(fields),
     1638                                field=self.(fields{i});
     1639                                s=whos('field');
     1640                                mem=mem+s.bytes/1e6;
     1641                                disp(sprintf('%19s: %6.2f Mb',fields{i},s.bytes/1e6));
     1642                        end
     1643                        disp(sprintf('%19s--%10s','--------------','--------------'));
     1644                        disp(sprintf('%19s: %g Mb','Total',mem));
    16421645                end
    1643                 disp(sprintf('%19s--%10s','--------------','--------------'));
    1644                 disp(sprintf('%19s: %g Mb','Total',mem));
    1645                 end % }}}
     1646                % }}}
    16461647                function netcdf(self,filename) % {{{
    1647                 %NETCDF - save model as netcdf
    1648                 %
    1649                 %   Usage:
    1650                 %      netcdf(md,filename)
    1651                 %
    1652                 %   Example:
    1653                 %      netcdf(md,'model.nc');
     1648                        %NETCDF - save model as netcdf
     1649                        %
     1650                        %   Usage:
     1651                        %      netcdf(md,filename)
     1652                        %
     1653                        %   Example:
     1654                        %      netcdf(md,'model.nc');
    16541655
    1655                 disp('Saving model as NetCDF');
    1656                 %1. Create NetCDF file
    1657                 ncid=netcdf.create(filename,'CLOBBER');
    1658                 netcdf.putAtt(ncid,netcdf.getConstant('NC_GLOBAL'),'Conventions','CF-1.4');
    1659                 netcdf.putAtt(ncid,netcdf.getConstant('NC_GLOBAL'),'Title',['ISSM model (' self.miscellaneous.name ')']);
    1660                 netcdf.putAtt(ncid,netcdf.getConstant('NC_GLOBAL'),'Author',getenv('USER'));
    1661                 netcdf.putAtt(ncid,netcdf.getConstant('NC_GLOBAL'),'Date',datestr(now));
     1656                        disp('Saving model as NetCDF');
     1657                        %1. Create NetCDF file
     1658                        ncid=netcdf.create(filename,'CLOBBER');
     1659                        netcdf.putAtt(ncid,netcdf.getConstant('NC_GLOBAL'),'Conventions','CF-1.4');
     1660                        netcdf.putAtt(ncid,netcdf.getConstant('NC_GLOBAL'),'Title',['ISSM model (' self.miscellaneous.name ')']);
     1661                        netcdf.putAtt(ncid,netcdf.getConstant('NC_GLOBAL'),'Author',getenv('USER'));
     1662                        netcdf.putAtt(ncid,netcdf.getConstant('NC_GLOBAL'),'Date',datestr(now));
    16621663
    1663                 %Preallocate variable id, needed to write variables in netcdf file
    1664                 var_id=zeros(1000,1);%preallocate
     1664                        %Preallocate variable id, needed to write variables in netcdf file
     1665                        var_id=zeros(1000,1);%preallocate
    16651666
    1666                 for step=1:2,
    1667                         counter=0;
    1668                         [var_id,counter]=structtonc(ncid,'md',self,0,var_id,counter,step);
    1669                         if step==1, netcdf.endDef(ncid); end
    1670                 end
     1667                        for step=1:2,
     1668                                counter=0;
     1669                                [var_id,counter]=structtonc(ncid,'md',self,0,var_id,counter,step);
     1670                                if step==1, netcdf.endDef(ncid); end
     1671                        end
    16711672
    1672                 if counter>1000,
    1673                         warning(['preallocation of var_id need to be updated from ' num2str(1000) ' to ' num2str(counter)]);
    1674                 end
     1673                        if counter>1000,
     1674                                warning(['preallocation of var_id need to be updated from ' num2str(1000) ' to ' num2str(counter)]);
     1675                        end
    16751676
    1676                 netcdf.close(ncid)
     1677                        netcdf.close(ncid)
    16771678                end % }}}
    16781679                function xylim(self) % {{{
    16791680
    … …  
    16811682                        ylim([min(self.mesh.y) max(self.mesh.y)])
    16821683                end % }}}
    16831684                function md=upload(md) % {{{
    1684                 %the goal of this routine is to upload the model onto a server, and to empty it.
    1685                 %So first, save the model with a unique name and upload the file to the server:
    1686                 random_part=fix(rand(1)*10000);
    1687                 id=[md.miscellaneous.name '-' regexprep(datestr(now),'[^\w'']','') '-' num2str(random_part)  '-' getenv('USER') '-' oshostname() '.upload'];
    1688                 eval(['save ' id ' md']);
     1685                        %the goal of this routine is to upload the model onto a server, and to empty it.
     1686                        %So first, save the model with a unique name and upload the file to the server:
     1687                        random_part=fix(rand(1)*10000);
     1688                        id=[md.miscellaneous.name '-' regexprep(datestr(now),'[^\w'']','') '-' num2str(random_part)  '-' getenv('USER') '-' oshostname() '.upload'];
     1689                        eval(['save ' id ' md']);
    16891690
    1690                 %Now, upload the file:
    1691                 issmscpout(md.settings.upload_server,md.settings.upload_path,md.settings.upload_login,md.settings.upload_port,{id},1);
     1691                        %Now, upload the file:
     1692                        issmscpout(md.settings.upload_server,md.settings.upload_path,md.settings.upload_login,md.settings.upload_port,{id},1);
    16921693
    1693                 %Now, empty this model of everything except settings, and record name of file we just uploaded!
    1694                 settings_back=md.settings;
    1695                 md=model();
    1696                 md.settings=settings_back;
    1697                 md.settings.upload_filename=id;
     1694                        %Now, empty this model of everything except settings, and record name of file we just uploaded!
     1695                        settings_back=md.settings;
     1696                        md=model();
     1697                        md.settings=settings_back;
     1698                        md.settings.upload_filename=id;
    16981699
    1699                 %get locally rid of file that was uploaded
    1700                 eval(['delete ' id]);
     1700                        %get locally rid of file that was uploaded
     1701                        eval(['delete ' id]);
    17011702
    17021703                end % }}}
    17031704                function md=download(md) % {{{
    17041705
    1705                 %the goal of this routine is to download the internals of the current model from a server, because
    1706                 %this model is empty, except for the settings which tell us where to go and find this model!
     1706                        %the goal of this routine is to download the internals of the current model from a server, because
     1707                        %this model is empty, except for the settings which tell us where to go and find this model!
    17071708
    1708                 %Download the file:
    1709                 issmscpin(md.settings.upload_server, md.settings.upload_login, md.settings.upload_port, md.settings.upload_path, {md.settings.upload_filename});
     1709                        %Download the file:
     1710                        issmscpin(md.settings.upload_server, md.settings.upload_login, md.settings.upload_port, md.settings.upload_path, {md.settings.upload_filename});
    17101711
    1711                 name=md.settings.upload_filename;
     1712                        name=md.settings.upload_filename;
    17121713
    1713                 %Now, load this model:
    1714                 md=loadmodel(md.settings.upload_filename);
     1714                        %Now, load this model:
     1715                        md=loadmodel(md.settings.upload_filename);
    17151716
    1716                 %get locally rid of file that was downloaded
    1717                 eval(['delete ' name]);
     1717                        %get locally rid of file that was downloaded
     1718                        eval(['delete ' name]);
    17181719
    17191720                end % }}}
    17201721                function savemodeljs(md,modelname,websiteroot,varargin) % {{{

  • ../trunk-jpl/src/m/classes/model.py

     
    7575
    7676
    7777class model(object):
    78     #properties
     78    """MODEL - class definition
     79
     80    Usage:
     81        md = model()
     82    """
     83
    7984    def __init__(self, *args): #{{{
    8085        self.mesh = None
    8186        self.mask = None
     87
     88        self.geometry = None
    8289        self.constants = None
    83         self.geometry = None
    84         self.initialization = None
    8590        self.smb = None
    8691        self.basalforcings = None
     92        self.materials = None
     93        self.damage = None
    8794        self.friction = None
     95        self.flowequation = None
     96        self.timestepping = None
     97        self.initialization = None
    8898        self.rifts = None
     99        self.dsl = None
    89100        self.solidearth = None
    90         self.dsl = None
    91         self.timestepping = None
    92         self.groundingline = None
    93         self.materials = None
    94         self.damage = None
    95         self.flowequation = None
     101
    96102        self.debug = None
    97103        self.verbose = None
    98104        self.settings = None
    99105        self.toolkits = None
    100106        self.cluster = None
     107
    101108        self.balancethickness = None
    102109        self.stressbalance = None
     110        self.groundingline = None
    103111        self.hydrology = None
    104112        self.masstransport = None
    105113        self.thermal = None
    … …  
    111119        self.love = None
    112120        self.esa = None
    113121        self.sampling = None
     122
    114123        self.autodiff = None
    115124        self.inversion = None
    116125        self.qmu = None
    117126        self.amr = None
    118         self.radaroverlay = None
    119127        self.results = None
    120128        self.outputdefinition = None
     129        self.radaroverlay = None
    121130        self.miscellaneous = None
    122131        self.private = None
    123132
    124         nargs = len(args)
    125 
    126         if nargs == 0:
     133        if len(args) == 0:
    127134            self.setdefaultparameters('earth')
    128135        else:
    129             self.setdefaultparameters(args[0])
    130136            options = pairoptions(*args)
    131137            planet = options.getfieldvalue('planet', 'earth')
    132138            self.setdefaultparameters(planet)
    133 #}}}
     139    #}}}
    134140
    135     def properties(self):  # {{{
    136         # ordered list of properties since vars(self) is random
    137         return ['mesh',
    138                 'mask',
    139                 'constants',
    140                 'geometry',
    141                 'initialization',
    142                 'smb',
    143                 'basalforcings',
    144                 'friction',
    145                 'rifts',
    146                 'solidearth',
    147                 'dsl',
    148                 'timestepping',
    149                 'groundingline',
    150                 'materials',
    151                 'damage',
    152                 'flowequation',
    153                 'debug',
    154                 'verbose',
    155                 'settings',
    156                 'toolkits',
    157                 'cluster',
    158                 'balancethickness',
    159                 'stressbalance',
    160                 'hydrology',
    161                 'masstransport',
    162                 'thermal',
    163                 'steadystate',
    164                 'transient',
    165                 'levelset',
    166                 'calving',
    167                 'frontalforcings',
    168                 'love',
    169                 'esa',
    170                 'sampling',
    171                 'autodiff',
    172                 'inversion',
    173                 'qmu',
    174                 'amr',
    175                 'radaroverlay',
    176                 'results',
    177                 'outputdefinition',
    178                 'miscellaneous',
    179                 'private']
    180     # }}}
    181 
    182141    def __repr__(obj):  #{{{
    183142        # TODO:
    184143        # - Convert all formatting to calls to <string>.format (see any
    185144        #   already converted <class>.__repr__ method for examples)
    186145        #
    187 
    188         #print "Here %s the number: %d" % ("is", 37)
    189146        s = "%19s: %-22s -- %s" % ("mesh", "[%s %s]" % ("1x1", obj.mesh.__class__.__name__), "mesh properties")
    190147        s = "%s\n%s" % (s, "%19s: %-22s -- %s" % ("mask", "[%s %s]" % ("1x1", obj.mask.__class__.__name__), "defines grounded and floating elements"))
    191148        s = "%s\n%s" % (s, "%19s: %-22s -- %s" % ("geometry", "[%s %s]" % ("1x1", obj.geometry.__class__.__name__), "surface elevation, bedrock topography, ice thickness, ..."))
    … …  
    229186        s = "%s\n%s" % (s, "%19s: %-22s -- %s" % ("radaroverlay", "[%s %s]" % ("1x1", obj.radaroverlay.__class__.__name__), "radar image for plot overlay"))
    230187        s = "%s\n%s" % (s, "%19s: %-22s -- %s" % ("miscellaneous", "[%s %s]" % ("1x1", obj.miscellaneous.__class__.__name__), "miscellaneous fields"))
    231188        return s
    232     # }}}
     189    #}}}
    233190
     191    def properties(self): #{{{
     192        # ordered list of properties since vars(self) is random
     193        return ['mesh',
     194                'mask',
     195                'geometry',
     196                'constants',
     197                'smb',
     198                'basalforcings',
     199                'materials',
     200                'damage',
     201                'friction',
     202                'flowequation',
     203                'timestepping',
     204                'initialization',
     205                'rifts',
     206                'dsl',
     207                'solidearth',
     208                'debug',
     209                'verbose',
     210                'settings',
     211                'toolkits',
     212                'cluster',
     213                'balancethickness',
     214                'stressbalance',
     215                'groundingline',
     216                'hydrology',
     217                'masstransport',
     218                'thermal',
     219                'steadystate',
     220                'transient',
     221                'levelset',
     222                'calving',
     223                'frontalforcings',
     224                'love',
     225                'esa',
     226                'sampling',
     227                'autodiff',
     228                'inversion',
     229                'qmu',
     230                'amr',
     231                'results',
     232                'outputdefinition',
     233                'radaroverlay',
     234                'miscellaneous',
     235                'private']
     236    #}}}
     237
    234238    def setdefaultparameters(self, planet): #{{{
    235239        self.mesh = mesh2d()
    236240        self.mask = mask()
    … …  
    277281        self.private = private()
    278282    #}}}
    279283
    280     def checkmessage(self, string):  # {{{
     284    def checkmessage(self, string): #{{{
    281285        print("model not consistent: {}".format(string))
    282286        self.private.isconsistent = False
    283287        return self
    284     # }}}
     288    #}}}
    285289    #@staticmethod
    286290
    287     def extract(self, area):  # {{{
     291    def extract(self, area): #{{{
    288292        """EXTRACT - extract a model according to an Argus contour or flag list
    289293
    290294        This routine extracts a submodel from a bigger model with respect to a given contour
    … …  
    561565        md2.mesh.extractedelements = pos_elem + 1
    562566
    563567        return md2
    564     # }}}
     568    #}}}
    565569
    566     def extrude(md, *args):  # {{{
     570    def extrude(md, *args): #{{{
    567571        """EXTRUDE - vertically extrude a 2d mesh
    568572
    569573        vertically extrude a 2d mesh and create corresponding 3d mesh.
    … …  
    748752            md.mesh.average_vertex_connectivity = 100
    749753
    750754        return md
    751     # }}}
     755    #}}}
    752756
    753757    def collapse(md): #{{{
    754758        """COLLAPSE - collapses a 3d mesh into a 2d mesh

  • ../trunk-jpl/src/m/classes/sampling.m

     
    55
    66classdef sampling
    77        properties (SetAccess=public)
    8          kappa          = NaN;
    9                  tau            = 0;
    10                  beta           = NaN;
    11          phi            = 0;
    12          alpha          = 0;
    13          robin          = 0;
    14          seed           = 0;
    15                  requested_outputs      = {};
    16     end
     8                kappa             = NaN;
     9                tau               = 0;
     10                beta              = NaN;
     11                phi               = 0;
     12                alpha             = 0;
     13                robin             = 0;
     14                seed              = 0;
     15                requested_outputs = {};
     16        end
    1717        methods
    1818                function self = sampling(varargin) % {{{
    19             switch nargin
     19                        switch nargin
    2020                                case 0
    2121                                        self=setdefaultparameters(self);
    2222                                otherwise
    2323                                        error('constructor not supported');
    24             end
     24                        end
    2525                end % }}}
    26                 function list = defaultoutputs(self,md) % {{{
     26                function disp(self) % {{{
    2727
    28             list = {};
     28                        disp(sprintf('   Sampling parameters:'));
    2929
    30                 end % }}}
     30                        disp(sprintf('\n      %s','Parameters of PDE operator (kappa^2 I-Laplacian)^(alpha/2)(tau):'));
     31                        fielddisplay(self,'kappa','coefficient of the identity operator');
     32                        fielddisplay(self,'tau','scaling coefficient of the solution (default: 1.0)');
     33                        fielddisplay(self,'alpha','exponent in PDE operator, (default: 2.0, BiLaplacian covariance operator)');
     34
     35                        disp(sprintf('\n      %s','Parameters of Robin boundary conditions nabla () \cdot normvec + beta ():'));
     36                        fielddisplay(self,'robin','Apply Robin boundary conditions (1 if applied and 0 for homogenous Neumann boundary conditions) (default: 0)');
     37                        fielddisplay(self,'beta','Coefficient in Robin boundary conditions (to be defined for robin = 1)');
     38
     39                        disp(sprintf('\n      %s','Parameters for first-order autoregressive process (X_t = phi X_{t-1} + noise) (if transient):'));
     40                        fielddisplay(self,'phi','Temporal correlation factor (|phi|<1 for stationary process, phi = 1 for random walk process) (default 0)');
     41
     42                        disp(sprintf('\n      %s','Other parameters of stochastic sampler:'));
     43                        fielddisplay(self,'seed','Seed for pseudorandom number generator (given seed if >=0 and random seed if <0) (default: -1)');
     44                        fielddisplay(self,'requested_outputs','additional outputs requested (not implemented yet)');
     45
     46                end % }}}','
    3147                function self = setdefaultparameters(self) % {{{
    32            
    33             %Scaling coefficient
    34                         self.tau=1; 
    3548
    36             %Apply Robin boundary conditions
    37                         self.robin=0;   
    38            
    39             %Temporal correlation factor
    40                         self.phi=0; 
    41            
    42                         %Exponent in fraction SPDE (default=2, biLaplacian covariance
     49                        %Scaling coefficient
     50                        self.tau=1;
     51
     52                        %Apply Robin boundary conditions
     53                        self.robin=0;
     54
     55                        %Temporal correlation factor
     56                        self.phi=0;
     57
     58                        %Exponent in fraction SPDE (default: 2, biLaplacian covariance
    4359                        %operator)
    4460                        self.alpha=2; % Default
    45            
    46             %Seed for pseudorandom number generator (default -1 for random seed)
    47             self.seed=-1;
    48            
     61
     62                        %Seed for pseudorandom number generator (default: -1, for random seed)
     63                        self.seed=-1;
     64
    4965                        %default output
    5066                        self.requested_outputs={'default'};
    5167
    5268                end % }}}
    53         function md = setparameters(self,md,lc,sigma) % {{{
    54            
    55             nu = self.alpha-1;
    56             KAPPA = sqrt(8*nu)/lc;
    57             TAU = sqrt(gamma(nu)/(gamma(self.alpha)*(4*pi)*KAPPA^(2*nu)*sigma^2));
    58             md.sampling.kappa = KAPPA*ones(md.mesh.numberofvertices,1);
    59             md.sampling.tau = TAU;
    60        
     69                function list = defaultoutputs(self,md) % {{{
     70
     71                        list = {};
     72
    6173                end % }}}
    6274                function md = checkconsistency(self,md,solution,analyses) % {{{
    63            
    64             if ~ismember('SamplingAnalysis',analyses), return; end
    65                        
    66             md = checkfield(md,'fieldname','sampling.kappa','NaN',1,'Inf',1,'size',[md.mesh.numberofvertices 1],'>',0);
    67             md = checkfield(md,'fieldname','sampling.tau','NaN',1,'Inf',1,'numel',1,'>',0);
    68             md = checkfield(md,'fieldname','sampling.robin','numel',1,'values',[0 1]);
    69             if(md.sampling.robin)
    70                 md = checkfield(md,'fieldname','sampling.beta','NaN',1,'Inf',1,'size',[md.mesh.numberofvertices 1],'>',0);
    71             end   
    72             md = checkfield(md,'fieldname','sampling.phi','NaN',1,'Inf',1,'numel',1,'>=',0);
    73             md = checkfield(md,'fieldname','sampling.alpha','NaN',1,'Inf',1,'numel',1,'>',0);
    74             md = checkfield(md,'fieldname','sampling.seed','NaN',1,'Inf',1,'numel',1);
    75             md = checkfield(md,'fieldname','sampling.requested_outputs','stringrow',1);
    7675
    77                 end % }}}
    78                 function disp(self) % {{{
     76                        if ~ismember('SamplingAnalysis',analyses), return; end
    7977
    80             disp(sprintf('   Sampling parameters:'));
     78                        md = checkfield(md,'fieldname','sampling.kappa','NaN',1,'Inf',1,'size',[md.mesh.numberofvertices 1],'>',0);
     79                        md = checkfield(md,'fieldname','sampling.tau','NaN',1,'Inf',1,'numel',1,'>',0);
     80                        md = checkfield(md,'fieldname','sampling.robin','numel',1,'values',[0 1]);
     81                        if(md.sampling.robin)
     82                                md = checkfield(md,'fieldname','sampling.beta','NaN',1,'Inf',1,'size',[md.mesh.numberofvertices 1],'>',0);
     83                        end
     84                        md = checkfield(md,'fieldname','sampling.phi','NaN',1,'Inf',1,'numel',1,'>=',0);
     85                        md = checkfield(md,'fieldname','sampling.alpha','NaN',1,'Inf',1,'numel',1,'>',0);
     86                        md = checkfield(md,'fieldname','sampling.seed','NaN',1,'Inf',1,'numel',1);
     87                        md = checkfield(md,'fieldname','sampling.requested_outputs','stringrow',1);
    8188
    82                         disp(sprintf('\n      %s','Parameters of PDE operator (kappa^2 I-Laplacian)^(alpha/2)(tau):'));
    83             fielddisplay(self,'kappa','coefficient of the identity operator');
    84                         fielddisplay(self,'tau','scaling coefficient of the solution (default 1.0)');
    85                         fielddisplay(self,'alpha','exponent in PDE operator, (default 2.0, BiLaplacian covariance operator)');
    86          
    87                         disp(sprintf('\n      %s','Parameters of Robin boundary conditions nabla () \cdot normvec + beta ():'));
    88                         fielddisplay(self,'robin','Apply Robin boundary conditions (1 if applied and 0 for homogenous Neumann boundary conditions) (default 0)');
    89                         fielddisplay(self,'beta','Coefficient in Robin boundary conditions (to be defined for robin = 1)');         
    90            
    91             disp(sprintf('\n      %s','Parameters for first-order autoregressive process (X_t = phi X_{t-1} + noise) (if transient):'));
    92                         fielddisplay(self,'phi','Temporal correlation factor (|phi|<1 for stationary process, phi = 1 for random walk process) (default 0)');
    93            
    94                         disp(sprintf('\n      %s','Other parameters of stochastic sampler:'));
    95             fielddisplay(self,'seed','Seed for pseudorandom number generator (given seed if >=0 and random seed if <0) (default -1)');
    96                         fielddisplay(self,'requested_outputs','additional outputs requested (not implemented yet)');
    97  
    9889                end % }}}
    9990                function marshall(self,prefix,md,fid) % {{{
    10091
    101             WriteData(fid,prefix,'object',self,'fieldname','kappa','format','DoubleMat','mattype',1);
    102             WriteData(fid,prefix,'object',self,'fieldname','tau','format','Double');
    103             WriteData(fid,prefix,'object',self,'fieldname','beta','format','DoubleMat','mattype',1);
    104             WriteData(fid,prefix,'object',self,'fieldname','phi','format','Double');
    105             WriteData(fid,prefix,'object',self,'fieldname','alpha','format','Integer');
    106             WriteData(fid,prefix,'object',self,'fieldname','robin','format','Boolean');
    107             WriteData(fid,prefix,'object',self,'fieldname','seed','format','Integer');
    108            
     92                        WriteData(fid,prefix,'object',self,'fieldname','kappa','format','DoubleMat','mattype',1);
     93                        WriteData(fid,prefix,'object',self,'fieldname','tau','format','Double');
     94                        WriteData(fid,prefix,'object',self,'fieldname','beta','format','DoubleMat','mattype',1);
     95                        WriteData(fid,prefix,'object',self,'fieldname','phi','format','Double');
     96                        WriteData(fid,prefix,'object',self,'fieldname','alpha','format','Integer');
     97                        WriteData(fid,prefix,'object',self,'fieldname','robin','format','Boolean');
     98                        WriteData(fid,prefix,'object',self,'fieldname','seed','format','Integer');
     99
    109100                        %process requested outputs
    110101                        outputs = self.requested_outputs;
    111102                        pos  = find(ismember(outputs,'default'));
    … …  
    115106                        end
    116107                        WriteData(fid,prefix,'data',outputs,'name','md.sampling.requested_outputs','format','StringArray');
    117108                end % }}}
     109                function md = setparameters(self,md,lc,sigma) % {{{
     110
     111                        nu = self.alpha-1;
     112                        KAPPA = sqrt(8*nu)/lc;
     113                        TAU = sqrt(gamma(nu)/(gamma(self.alpha)*(4*pi)*KAPPA^(2*nu)*sigma^2));
     114                        md.sampling.kappa = KAPPA*ones(md.mesh.numberofvertices,1);
     115                        md.sampling.tau = TAU;
     116
     117                end % }}}
    118118                function savemodeljs(self,fid,modelname) % {{{
    119            
    120             writejsdouble(fid,[modelname '.sampling.kappa'],self.kappa);
    121             writejsdouble(fid,[modelname '.sampling.tau'],self.tau);
    122             writejsdouble(fid,[modelname '.sampling.beta'],self.beta);
    123             writejsdouble(fid,[modelname '.sampling.phi'],self.beta);
    124             writejsdouble(fid,[modelname '.sampling.alpha'],self.alpha);
    125             writejsdouble(fid,[modelname '.sampling.robin'],self.robin);
    126             writejsdouble(fid,[modelname '.sampling.seed'],self.seed);
     119
     120                        writejsdouble(fid,[modelname '.sampling.kappa'],self.kappa);
     121                        writejsdouble(fid,[modelname '.sampling.tau'],self.tau);
     122                        writejsdouble(fid,[modelname '.sampling.beta'],self.beta);
     123                        writejsdouble(fid,[modelname '.sampling.phi'],self.beta);
     124                        writejsdouble(fid,[modelname '.sampling.alpha'],self.alpha);
     125                        writejsdouble(fid,[modelname '.sampling.robin'],self.robin);
     126                        writejsdouble(fid,[modelname '.sampling.seed'],self.seed);
    127127                        writejscellstring(fid,[modelname '.sampling.requested_outputs'],self.requested_outputs);
    128128
    129129                end % }}}
    130130        end
    131 end
    132  No newline at end of file
     131end

  • ../trunk-jpl/src/m/classes/sampling.py

     
    11import numpy as np
    22
     3from math import *
     4
    35from checkfield import checkfield
    46from fielddisplay import fielddisplay
    57from project3d import project3d
    … …  
    1315        sampling = sampling()
    1416    """
    1517
    16     def __init__(self):  # {{{
    17         self.kappa = float('NaN')
     18    def __init__(self, *args): #{{{
     19        self.kappa = np.nan
    1820        self.tau = 0
    19         self.beta = float('NaN')
    20         self.hydrostatic_adjustment = 0
     21        self.beta = np.nan
    2122        self.phi = 0
    2223        self.alpha = 0
    2324        self.robin = 0
    … …  
    2425        self.seed = 0
    2526        self.requested_outputs = []
    2627
    27         # Set defaults
    28         self.setdefaultparameters()
    29 
     28        if len(args) == 0:
     29            self.setdefaultparameters()
     30        else:
     31            raise RuntimeError('constructor not supported')
    3032    #}}}
    3133
    32     def __repr__(self):  # {{{
     34    def __repr__(self): #{{{
    3335        s = '   Sampling parameters::\n'
    34         s += '{}\n'.format(fielddisplay(self,'kappa','coefficient of the identity operator'))
    35         s += '{}\n'.format(fielddisplay(self,'tau','scaling coefficient of the solution (default 1.0)'))
    36         s += '{}\n'.format(fielddisplay(self,'alpha','exponent in PDE operator, (default 2.0, BiLaplacian covariance operator)'))
    37         s += '{}\n'.format(disp(sprintf('\n      %s','Parameters of Robin boundary conditions nabla () \cdot normvec + beta ():')))
    38         s += '{}\n'.format(fielddisplay(self,'beta','Coefficient in Robin boundary conditions (to be defined for robin = 1)'))
    39         s += '{}\n'.format(fielddisplay(self,'phi','Temporal correlation factor (|phi|<1 for stationary process, phi = 1 for random walk process) (default 0)'))
    40         s += '{}\n'.format(fielddisplay(self,'seed','Seed for pseudorandom number generator (given seed if >=0 and random seed if <0) (default -1)'))
     36        s += '      Parameters of PDE operator (kappa^2 I-Laplacian)^(alpha/2)(tau):\n'
     37        s += '{}\n'.format(fielddisplay(self, 'kappa', 'coefficient of the identity operator'))
     38        s += '{}\n'.format(fielddisplay(self, 'tau', 'scaling coefficient of the solution (default: 1.0)'))
     39        s += '{}\n'.format(fielddisplay(self, 'alpha', 'exponent in PDE operator, (default: 2.0, BiLaplacian covariance operator)'))
     40
     41        s += '      Parameters of Robin boundary conditions nabla () \cdot normvec + beta ():\n'
     42        s += '{}\n'.format(fielddisplay(self, 'robin', 'Apply Robin boundary conditions (1 if applied and 0 for homogenous Neumann boundary conditions) (default: 0)'))
     43        s += '{}\n'.format(fielddisplay(self, 'beta', 'Coefficient in Robin boundary conditions (to be defined for robin = 1)'))
     44
     45        s += '      Parameters for first-order autoregressive process (X_t = phi X_{t-1} + noise) (if transient):\n'
     46        s += '{}\n'.format(fielddisplay(self, 'phi', 'Temporal correlation factor (|phi|<1 for stationary process, phi = 1 for random walk process) (default 0)'))
     47
     48        s += '      Other parameters of stochastic sampler:\n'
     49        s += '{}\n'.format(fielddisplay(self, 'seed', 'Seed for pseudorandom number generator (given seed if >=0 and random seed if <0) (default: -1)'))
    4150        s += '{}\n'.format(fielddisplay(self, 'requested_outputs', 'additional outputs requested (not implemented yet)'))
     51
    4252        return s
    4353    #}}}
    4454
    45     def defaultoutputs(self, md):  # {{{
    46         return []
    47 
    48     #}}}
    49     def setdefaultparameters(self):  # {{{
     55    def setdefaultparameters(self): #{{{
    5056        # Scaling coefficient
    5157        self.tau = 1
     58
    5259        # Apply Robin boundary conditions
    5360        self.robin = 0
     61
    5462        # Temporal correlation factor
    5563        self.phi = 0
    56         # Exponent in fraction SPDE (default=2, biLaplacian covariance operator)
     64
     65        # Exponent in fraction SPDE (default: 2, biLaplacian covariance operator)
    5766        self.alpha = 2
    58         # Seed for pseudorandom number generator (default -1 for random seed)
    59         self.alpha = -1
     67
     68        # Seed for pseudorandom number generator (default: -1, for random seed)
     69        self.seed = -1
     70
    6071        # Default output
    6172        self.requested_outputs = ['default']
     73
    6274        return self
    6375    #}}}
    6476
    65     def checkconsistency(self, md, solution, analyses):  # {{{
    66         # Early return
     77    def defaultoutputs(self, md): #{{{
     78        return []
     79    #}}}
     80
     81    def checkconsistency(self, md, solution, analyses): #{{{
    6782        if ('SamplingAnalysis' not in analyses):
    6883            return md
    6984
    70         md = checkfield(md,'fieldname','sampling.kappa','NaN',1,'Inf',1,'size',[md.mesh.numberofvertices],'>',0)
    71         md = checkfield(md,'fieldname','sampling.tau','NaN',1,'Inf',1,'numel',1,'>',0)
    72         md = checkfield(md,'fieldname','sampling.robin','numel',1,'values',[0, 1])
     85        md = checkfield(md, 'fieldname', 'sampling.kappa', 'NaN', 1, 'Inf', 1, 'size', [md.mesh.numberofvertices], '>', 0)
     86        md = checkfield(md, 'fieldname', 'sampling.tau', 'NaN', 1, 'Inf', 1, 'numel', 1, '>', 0)
     87        md = checkfield(md, 'fieldname', 'sampling.robin', 'numel', 1, 'values', [0, 1])
    7388        if md.sampling.robin:
    74             md = checkfield(md,'fieldname','sampling.beta','NaN',1,'Inf',1,'size',[md.mesh.numberofvertices],'>',0)
    75         md = checkfield(md,'fieldname','sampling.phi','NaN',1,'Inf',1,'numel',1,'>=',0)
    76         md = checkfield(md,'fieldname','sampling.alpha','NaN',1,'Inf',1,'numel',1,'>',0)
    77         md = checkfield(md,'fieldname','sampling.seed','NaN',1,'Inf',1,'numel',1)
    78         md = checkfield(md,'fieldname','sampling.requested_outputs','stringrow',1)
     89            md = checkfield(md, 'fieldname', 'sampling.beta', 'NaN', 1, 'Inf', 1, 'size', [md.mesh.numberofvertices], '>', 0)
     90        end
     91        md = checkfield(md, 'fieldname', 'sampling.phi', 'NaN', 1, 'Inf', 1, 'numel', 1, '>=', 0)
     92        md = checkfield(md, 'fieldname', 'sampling.alpha', 'NaN', 1, 'Inf', 1, 'numel', 1, '>', 0)
     93        md = checkfield(md, 'fieldname', 'sampling.seed', 'NaN', 1, 'Inf', 1, 'numel', 1)
     94        md = checkfield(md, 'fieldname', 'sampling.requested_outputs', 'stringrow', 1)
    7995
    8096        return md
    81     # }}}
     97    #}}}
    8298
    83     def marshall(self, prefix, md, fid):  # {{{
    84         WriteData(fid,prefix,'object',self,'fieldname','kappa','format','DoubleMat','mattype',1)
    85         WriteData(fid,prefix,'object',self,'fieldname','tau','format','Double')
    86         WriteData(fid,prefix,'object',self,'fieldname','beta','format','DoubleMat','mattype',1)
    87         WriteData(fid,prefix,'object',self,'fieldname','phi','format','Double')
    88         WriteData(fid,prefix,'object',self,'fieldname','alpha','format','Integer')
    89         WriteData(fid,prefix,'object',self,'fieldname','robin','format','Boolean')
    90         WriteData(fid,prefix,'object',self,'fieldname','seed','format','Integer')
     99    def marshall(self, prefix, md, fid): #{{{
     100        WriteData(fid, prefix, 'object', self, 'fieldname', 'kappa', 'format', 'DoubleMat', 'mattype', 1)
     101        WriteData(fid, prefix, 'object', self, 'fieldname', 'tau', 'format', 'Double')
     102        WriteData(fid, prefix, 'object', self, 'fieldname', 'beta', 'format', 'DoubleMat', 'mattype', 1)
     103        WriteData(fid, prefix, 'object', self, 'fieldname', 'phi', 'format', 'Double')
     104        WriteData(fid, prefix, 'object', self, 'fieldname', 'alpha', 'format', 'Integer')
     105        WriteData(fid, prefix, 'object', self, 'fieldname', 'robin', 'format', 'Boolean')
     106        WriteData(fid, prefix, 'object', self, 'fieldname', 'seed', 'format', 'Integer')
    91107
    92     #process requested outputs
     108        # Process requested outputs
    93109        outputs = self.requested_outputs
    94110        indices = [i for i, x in enumerate(outputs) if x == 'default']
    95111        if len(indices) > 0:
    … …  
    96112            outputscopy = outputs[0:max(0, indices[0] - 1)] + self.defaultoutputs(md) + outputs[indices[0] + 1:]
    97113            outputs = outputscopy
    98114        WriteData(fid, prefix, 'data', outputs, 'name', 'md.sampling.requested_outputs', 'format', 'StringArray')
     115    #}}}
    99116
    100     # }}}
     117    def setparameters(self, md, lc, sigma): #{{{
     118        nu = self.alpha - 1
     119        KAPPA = pow((8 * nu), 0.5) / lc
     120        TAU = pow((math.gamma(nu) / math.gamma(self.alpha) * (4 *  np.pi) * pow(KAPPA, 2 * nu) * pow(sigma, 2)), 0.5)
     121        md.sampling.kappa = KAPPA * np.ones((md.mesh.numberofvertices, 1))
     122        md.sampling.tau = TAU
     123
     124        return md
     125    #}}}

  • ../trunk-jpl/src/m/classes/solidearthsettings.m

     
    9191                self.grdmodel=0;
    9292
    9393                end % }}}
     94                function disp(self) % {{{
     95                        disp(sprintf('   solidearth settings:'));
     96
     97                        fielddisplay(self,'reltol','sea level change relative convergence criterion (default, NaN: not applied)');
     98                        fielddisplay(self,'abstol','sea level change absolute convergence criterion(default, NaN: not applied)');
     99                        fielddisplay(self,'maxiter','maximum number of nonlinear iterations');
     100                        fielddisplay(self,'grdocean','does this planet have an ocean, if set to 1: global water mass is conserved in GRD module (default: 1)');
     101                        fielddisplay(self,'ocean_area_scaling','correction for model representation of ocean area (default: No correction)');
     102                        fielddisplay(self,'sealevelloading','enables surface loading from sea-level change (default: 1)');
     103                        fielddisplay(self,'isgrd','compute GRD patterns (default: 1)');
     104                        fielddisplay(self,'compute_bp_grd','compute GRD patterns for bottom pressure loads (default: 1)');
     105                        fielddisplay(self,'runfrequency','how many time steps we skip before we run solidearthsettings solver during transient (default: 1)');
     106                        fielddisplay(self,'selfattraction','enables surface mass load to perturb the gravity field');
     107                        fielddisplay(self,'elastic','enables elastic deformation from surface loading');
     108                        fielddisplay(self,'viscous','enables viscous deformation from surface loading');
     109                        fielddisplay(self,'rotation','enables polar motion to feedback on the GRD fields');
     110                        fielddisplay(self,'degacc','accuracy (default: .01 deg) for numerical discretization of the Green''s functions');
     111                        fielddisplay(self,'timeacc','time accuracy (default: 1 yr) for numerical discretization of the Green''s functions');
     112                        fielddisplay(self,'grdmodel','type of deformation model, 0 for no GRD, 1 for spherical GRD model (SESAW model), 2 for half-space planar GRD (visco-elastic model from Ivins)');
     113                        fielddisplay(self,'cross_section_shape','1: square-edged (default). 2: elliptical. See iedge in GiaDeflectionCore');
     114                end % }}}
    94115                function md = checkconsistency(self,md,solution,analyses) % {{{
    95116
    96117                        if ~ismember('SealevelchangeAnalysis',analyses) | (strcmp(solution,'TransientSolution') & md.transient.isslc==0),
    … …  
    134155                        end
    135156
    136157                end % }}}
    137                 function disp(self) % {{{
    138                         disp(sprintf('   solidearth settings:'));
    139 
    140                         fielddisplay(self,'reltol','sea level change relative convergence criterion (default, NaN: not applied)');
    141                         fielddisplay(self,'abstol','sea level change absolute convergence criterion(default, NaN: not applied)');
    142                         fielddisplay(self,'maxiter','maximum number of nonlinear iterations');
    143                         fielddisplay(self,'grdocean','does this planet have an ocean, if set to 1: global water mass is conserved in GRD module (default: 1)');
    144                         fielddisplay(self,'ocean_area_scaling','correction for model representation of ocean area (default: No correction)');
    145                         fielddisplay(self,'sealevelloading','enables surface loading from sea-level change (default: 1)');
    146                         fielddisplay(self,'isgrd','compute GRD patterns (default: 1)');
    147                         fielddisplay(self,'compute_bp_grd','compute GRD patterns for bottom pressure loads (default: 1)');
    148                         fielddisplay(self,'runfrequency','how many time steps we skip before we run solidearthsettings solver during transient (default: 1)');
    149                         fielddisplay(self,'selfattraction','enables surface mass load to perturb the gravity field');
    150                         fielddisplay(self,'elastic','enables elastic deformation from surface loading');
    151                         fielddisplay(self,'viscous','enables viscous deformation from surface loading');
    152                         fielddisplay(self,'rotation','enables polar motion to feedback on the GRD fields');
    153                         fielddisplay(self,'degacc','accuracy (default: .01 deg) for numerical discretization of the Green''s functions');
    154                         fielddisplay(self,'timeacc','time accuracy (default: 1 yr) for numerical discretization of the Green''s functions');
    155                         fielddisplay(self,'grdmodel','type of deformation model, 0 for no GRD, 1 for spherical GRD model (SESAW model), 2 for half-space planar GRD (visco-elastic model from Ivins)');
    156                         fielddisplay(self,'cross_section_shape','1: square-edged (default). 2: elliptical. See iedge in GiaDeflectionCore');
    157                 end % }}}
    158158                function marshall(self,prefix,md,fid) % {{{
    159159                        WriteData(fid,prefix,'object',self,'fieldname','reltol','name','md.solidearth.settings.reltol','format','Double');
    160160                        WriteData(fid,prefix,'object',self,'fieldname','abstol','name','md.solidearth.settings.abstol','format','Double');

  • ../trunk-jpl/src/m/miscellaneous/MatlabFuncs.py

     
    1 def oshostname():
    2     import socket
     1def acosd(X): #{{{
     2    """ function acosd - Inverse cosine in degrees
    33
    4     return socket.gethostname()
     4    Usage:
     5        Y = acosd(X)
     6    """
     7    import numpy as np
    58
     9    return np.degrees(np.arccos(X))
     10#}}}
    611
    7 def ispc():
    8     import platform
     12def asind(X): #{{{
     13    """ function asind - Inverse sine in degrees
    914
    10     if 'Windows' in platform.system():
    11         return True
    12     else:
    13         return False
     15    Usage:
     16        Y = asind(X)
     17    """
     18    import numpy as np
    1419
     20    return np.degrees(np.arcsin(X))
     21#}}}
    1522
    16 def ismac():
    17     import platform
     23def atand(X): #{{{
     24    """ function atand - Inverse tangent in degrees
    1825
    19     if 'Darwin' in platform.system():
    20         return True
    21     else:
    22         return False
     26    Usage:
     27        Y = atand(X)
     28    """
     29    import numpy as np
    2330
     31    return np.degrees(np.arctan(X))
     32#}}}
    2433
    25 def strcmp(s1, s2):
    2634
    27     if s1 == s2:
    28         return True
    29     else:
    30         return False
     35def atan2d(Y, X): #{{{
     36    """ function atan2d - Four-quadrant inverse tangent in degrees
    3137
     38    Usage:
     39        D = atan2d(Y, X)
     40    """
     41    import numpy as np
    3242
    33 def strncmp(s1, s2, n):
     43    return np.degrees(np.arctan2(Y, X))
     44#}}}
    3445
    35     if s1[0:n] == s2[0:n]:
    36         return True
     46def det(a): #{{{
     47    if a.shape == (1, ):
     48        return a[0]
     49    elif a.shape == (1, 1):
     50        return a[0, 0]
     51    elif a.shape == (2, 2):
     52        return a[0, 0] * a[1, 1] - a[0, 1] * a[1, 0]
    3753    else:
    38         return False
     54        raise TypeError("MatlabFunc.det only implemented for shape (2, 2), not for shape %s." % str(a.shape))
     55#}}}
    3956
     57def heaviside(x): #{{{
     58    import numpy as np
    4059
    41 def strcmpi(s1, s2):
     60    y = np.zeros_like(x)
     61    y[np.nonzero(x > 0.)] = 1.
     62    y[np.nonzero(x == 0.)] = 0.5
    4263
    43     if s1.lower() == s2.lower():
    44         return True
    45     else:
    46         return False
     64    return y
     65#}}}
    4766
     67def ismac(): #{{{
     68    import platform
    4869
    49 def strncmpi(s1, s2, n):
    50 
    51     if s1.lower()[0:n] == s2.lower()[0:n]:
     70    if 'Darwin' in platform.system():
    5271        return True
    5372    else:
    5473        return False
     74#}}}
    5575
    56 
    57 def ismember(a, s):
     76def ismember(a, s): #{{{
    5877    import numpy as np
    5978
    6079    if not isinstance(s, (tuple, list, dict, np.ndarray)):
    … …  
    6584
    6685    if not isinstance(a, np.ndarray):
    6786        b = [item in s for item in a]
    68 
    6987    else:
    7088        if not isinstance(s, np.ndarray):
    7189            b = np.empty_like(a)
    7290            for i, item in enumerate(a.flat):
    7391                b.flat[i] = item in s
    74 
    7592        else:
    7693            b = np.in1d(a.flat, s.flat).reshape(a.shape)
    7794
    7895    return b
     96#}}}
    7997
     98def ispc(): #{{{
     99    import platform
    80100
    81 def det(a):
    82 
    83     if a.shape == (1, ):
    84         return a[0]
    85     elif a.shape == (1, 1):
    86         return a[0, 0]
    87     elif a.shape == (2, 2):
    88         return a[0, 0] * a[1, 1] - a[0, 1] * a[1, 0]
     101    if 'Windows' in platform.system():
     102        return True
    89103    else:
    90         raise TypeError("MatlabFunc.det only implemented for shape (2, 2), not for shape %s." % str(a.shape))
     104        return False
     105#}}}
    91106
     107def oshostname(): #{{{
     108    import socket
    92109
    93 def sparse(ivec, jvec, svec, m=0, n=0, nzmax=0):
     110    return socket.gethostname()
     111#}}}
     112
     113def sparse(ivec, jvec, svec, m=0, n=0, nzmax=0): #{{{
    94114    import numpy as np
    95115
    96116    if not m:
    … …  
    104124        a[i - 1, j - 1] += s
    105125
    106126    return a
     127#}}}
    107128
     129def strcmp(s1, s2): #{{{
     130    if s1 == s2:
     131        return True
     132    else:
     133        return False
     134#}}}
    108135
    109 def heaviside(x):
    110     import numpy as np
     136def strcmpi(s1, s2): #{{{
     137    if s1.lower() == s2.lower():
     138        return True
     139    else:
     140        return False
     141#}}}
    111142
    112     y = np.zeros_like(x)
    113     y[np.nonzero(x > 0.)] = 1.
    114     y[np.nonzero(x == 0.)] = 0.5
     143def strncmp(s1, s2, n): #{{{
     144    if s1[0:n] == s2[0:n]:
     145        return True
     146    else:
     147        return False
     148#}}}
    115149
    116     return y
     150def strncmpi(s1, s2, n): #{{{
     151    if s1.lower()[0:n] == s2.lower()[0:n]:
     152        return True
     153    else:
     154        return False
     155#}}}

  • ../trunk-jpl/src/m/miscellaneous/PythonFuncs.py

     
    11import numpy as np
    22
    33
    4 def logical_and_n(*arg):
    5 
     4def logical_and_n(*arg): #{{{
    65    if len(arg):
    76        result = arg[0]
    87        for item in arg[1:]:
    98            result = np.logical_and(result, item)
    109        return result
    11 
    1210    else:
    1311        return None
     12#}}}
    1413
    15 
    16 def logical_or_n(*arg):
    17 
     14def logical_or_n(*arg): #{{{
    1815    if len(arg):
    1916        result = arg[0]
    2017        for item in arg[1:]:
    2118            result = np.logical_or(result, item)
    2219        return result
    23 
    2420    else:
    2521        return None
     22#}}}

  • ../trunk-jpl/src/m/solve/marshall.m

     
    4444st=fclose(fid);
    4545
    4646% Uncomment the following to make a copy of the binary input file for debugging
    47 % purposes (can be fed into scripts/ReadBin.py).
     47% purposes (can be fed into scripts/BinRead.py).
    4848% copyfile([md.miscellaneous.name '.bin'], [md.miscellaneous.name '.m.bin'])
    4949
    5050if st==-1,

  • ../trunk-jpl/src/m/solve/marshall.py

     
    4545        fid.close()
    4646
    4747        # Uncomment the following to make a copy of the binary input file for
    48         # debugging purposes (can be fed into scripts/ReadBin.py).
     48        # debugging purposes (can be fed into scripts/BinRead.py).
    4949        # copy_cmd = "cp {}.bin {}.py.bin".format(md.miscellaneous.name, md.miscellaneous.name)
    5050        # subprocess.call(copy_cmd, shell=True)
    5151

  • ../trunk-jpl/src/m/classes/dsl.py

     
    77
    88
    99class dsl(object):
    10     """DSL - slass definition
     10    """DSL - class definition
    1111
    1212    Usage:
    13         dsl = dsl()
     13        dsl = dsl() # dynamic sea level class, based on CMIP5 outputs
    1414    """
    1515
    16     def __init__(self): #{{{
     16    def __init__(self, *args): #{{{
    1717        self.global_average_thermosteric_sea_level = np.nan # Corresponds to zostoga field in CMIP5 archives. Specified as a temporally variable quantity (in m).
    1818        self.sea_surface_height_above_geoid        = np.nan # Corresponds to zos field in CMIP5 archives. Spatial average is 0. Specified as a spatio-temporally variable quantity (in m).
    1919        self.sea_water_pressure_at_sea_floor       = np.nan # Corresponds to bpo field in CMIP5 archives. Specified as a spatio-temporally variable quantity (in m equivalent, not in Pa!).
     20
     21        if len(args) == 0:
     22            self.setdefaultparameters()
     23        else:
     24            raise Exception('constructor not supported')
    2025    #}}}
    2126
    2227    def __repr__(self): #{{{
    … …  
    2732        return s
    2833    #}}}
    2934
    30     def defaultoutputs(self, md): #{{{
    31         return []
     35    def setdefaultparameters(self): #{{{
     36        self.global_average_thermosteric_sea_level = np.nan
     37        self.sea_surface_height_above_geoid        = np.nan
     38        self.sea_water_pressure_at_sea_floor       = np.nan
    3239    #}}}
    3340
    34     def initialize(self, md): #{{{
    35         if np.isnan(self.global_average_thermosteric_sea_level):
    36             self.global_average_thermosteric_sea_level = np.array([0, 0]).reshape(-1, 1)
    37             print('      no dsl.global_average_thermosteric_sea_level specified: transient values set at zero')
    38 
    39         if np.isnan(self.sea_surface_height_above_geoid):
    40             self.sea_surface_height_above_geoid = np.array(np.zeros(md.mesh.numberofvertices)).reshape(-1, 1)
    41             print('      no dsl.sea_surface_height_above_geoid specified: transient values set at zero')
    42 
    43         if np.isnan(self.sea_water_pressure_at_sea_floor):
    44             self.sea_water_pressure_at_sea_floor = np.array(np.zeros(md.mesh.numberofvertices)).reshape(-1, 1)
    45             print('      no dsl.sea_water_pressure_at_sea_floor specified: transient values set at zero')
    46     #}}}
    47 
    4841    def checkconsistency(self, md, solution, analyses): #{{{
    4942        # Early return
    50         if ('SealevelriseAnalysis' not in analyses) or (solution == 'TransientSolution' and not md.transient.isslc):
     43        if ('SealevelchangeAnalysis' not in analyses) or (solution == 'TransientSolution' and not md.transient.isslc):
    5144            return md
    5245
    5346        md = checkfield(md, 'fieldname', 'dsl.global_average_thermosteric_sea_level', 'NaN', 1, 'Inf', 1)
    … …  
    6053        return md
    6154    # }}}
    6255
     56    def marshall(self, prefix, md, fid):   #{{{
     57        yts = md.constants.yts
     58        WriteData(fid, prefix, 'name', 'md.dsl.model', 'data', 1, 'format', 'Integer')
     59        WriteData(fid, prefix, 'object', self, 'fieldname', 'global_average_thermosteric_sea_level', 'format', 'DoubleMat', 'mattype', 2, 'timeseries', 1, 'yts', yts) # mattype 2, because we are sending a GMSL value identical everywhere on each element.
     60        WriteData(fid, prefix, 'object', self, 'fieldname', 'sea_surface_height_above_geoid', 'format', 'DoubleMat', 'mattype', 1, 'timeserieslength', md.mesh.numberofvertices + 1, 'yts', yts) # mattype 1 because we specify DSL at vertex locations.
     61        WriteData(fid, prefix, 'object', self, 'fieldname', 'sea_water_pressure_at_sea_floor', 'format', 'DoubleMat', 'mattype', 1, 'timeserieslength', md.mesh.numberofvertices + 1, 'yts', yts) # mattype 1 because we specify bottom pressure at vertex locations.
     62    # }}}
     63
    6364    def extrude(self, md): #{{{
    64         self.sea_surface_height_above_geoid = project3d(md, 'vector', self.sea_surface_height_above_geoid, 'type', 'node')
    65         self.sea_water_pressure_at_sea_floor = project3d(md, 'vector', self.sea_water_pressure_at_sea_floor, 'type', 'node')
     65        self.sea_surface_height_above_geoid = project3d(md, 'vector', self.sea_surface_height_above_geoid, 'type', 'node', 'layer', 1)
     66        self.sea_water_pressure_at_sea_floor = project3d(md, 'vector', self.sea_water_pressure_at_sea_floor, 'type', 'node', 'layer', 1)
    6667        return self
    6768    #}}}
    6869
    69     def marshall(self, prefix, md, fid):   #{{{
    70         yts = md.constants.yts
    71         WriteData(fid, prefix, 'name', 'md.dsl.model', 'data', 1, 'format', 'Integer')
    72         WriteData(fid, prefix, 'object', self, 'class', 'dsl', 'fieldname', 'global_average_thermosteric_sea_level', 'format', 'DoubleMat', 'mattype', 2, 'timeserieslength', 1+1, 'yts', md.constants.yts, 'scale', 1e-3/md.constants.yts) # mattype 2, because we are sending a GMSL value identical everywhere on each element.
    73         WriteData(fid, prefix, 'object', self, 'class', 'dsl', 'fieldname', 'sea_water_pressure_at_sea_floor', 'format', 'DoubleMat', 'mattype', 1, 'timeserieslength', md.mesh.numberofvertices+1, 'yts', md.constants.yts, 'scale', 1e-3/md.constants.yts) # mattype 1 because we specify DSL at vertex locations.
    74         WriteData(fid, prefix, 'object', self, 'class', 'dsl', 'fieldname', 'sea_surface_height_above_geoid', 'format', 'DoubleMat', 'mattype', 1, 'timeserieslength', md.mesh.numberofvertices+1, 'yts', md.constants.yts) # mattype 1 because we specify bottom pressure at vertex locations.
    75     # }}}
     70
     71    def initialize(self, md): #{{{
     72        if np.isnan(self.global_average_thermosteric_sea_level):
     73            self.global_average_thermosteric_sea_level = np.array([0, 0]).reshape(-1, 1)
     74            print('      no dsl.global_average_thermosteric_sea_level specified: transient values set at zero')
     75
     76        if np.isnan(self.sea_surface_height_above_geoid):
     77            self.sea_surface_height_above_geoid = np.append(np.zeros((md.mesh.numberofvertices, 1)), 0).reshape(-1, 1)
     78            print('      no dsl.sea_surface_height_above_geoid specified: transient values set at zero')
     79
     80        if np.isnan(self.sea_water_pressure_at_sea_floor):
     81            self.sea_water_pressure_at_sea_floor = np.append(np.zeros((md.mesh.numberofvertices, 1)), 0).reshape(-1, 1)
     82            print('      no dsl.sea_water_pressure_at_sea_floor specified: transient values set at zero')
     83    #}}}

  • ../trunk-jpl/src/m/classes/hydrologyshreve.m

     
    1010                requested_outputs  = {};
    1111        end
    1212        methods
    13                 function self = extrude(self,md) % {{{
    14                 end % }}}
    1513                function self = hydrologyshreve(varargin) % {{{
    1614                        switch nargin
    1715                                case 0
    … …  
    2220                                        error('constructor not supported');
    2321                        end
    2422                end % }}}
    25                 function list = defaultoutputs(self,md) % {{{
    26                         list = {'Watercolumn','HydrologyWaterVx','HydrologyWaterVy'};
    27                 end % }}}   
     23                function disp(self) % {{{
     24                        disp(sprintf('   hydrologyshreve solution parameters:'));
     25                        fielddisplay(self,'spcwatercolumn','water thickness constraints (NaN means no constraint) [m]');
     26                        fielddisplay(self,'stabilization','artificial diffusivity (default: 1). can be more than 1 to increase diffusivity.');
     27                        fielddisplay(self,'requested_outputs','additional outputs requested');
    2828
     29                end % }}}
    2930                function self = setdefaultparameters(self) % {{{
    3031
    3132                        %Type of stabilization to use 0:nothing 1:artificial_diffusivity
    … …  
    3334                        self.requested_outputs = {'default'};
    3435                end % }}}
    3536                function md = checkconsistency(self,md,solution,analyses) % {{{
    36                        
     37
    3738                        %Early return
    3839                        if ~ismember('HydrologyShreveAnalysis',analyses) |  (strcmp(solution,'TransientSolution') & md.transient.ishydrology==0), return; end
    3940
    … …  
    4041                        md = checkfield(md,'fieldname','hydrology.spcwatercolumn','Inf',1,'timeseries',1);
    4142                        md = checkfield(md,'fieldname','hydrology.stabilization','>=',0);
    4243                end % }}}
    43                 function disp(self) % {{{
    44                         disp(sprintf('   hydrologyshreve solution parameters:'));
    45                         fielddisplay(self,'spcwatercolumn','water thickness constraints (NaN means no constraint) [m]');
    46                         fielddisplay(self,'stabilization','artificial diffusivity (default is 1). can be more than 1 to increase diffusivity.');
    47       fielddisplay(self,'requested_outputs','additional outputs requested');
    48 
     44                function list = defaultoutputs(self,md) % {{{
     45                        list = {'Watercolumn','HydrologyWaterVx','HydrologyWaterVy'};
    4946                end % }}}
    5047                function marshall(self,prefix,md,fid) % {{{
    5148                        WriteData(fid,prefix,'name','md.hydrology.model','data',2,'format','Integer');
    5249                        WriteData(fid,prefix,'object',self,'fieldname','spcwatercolumn','format','DoubleMat','mattype',1,'timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts);
    5350                        WriteData(fid,prefix,'object',self,'fieldname','stabilization','format','Double');
    54       outputs = self.requested_outputs;
    55       pos = find(ismember(outputs,'default'));
    56       if ~isempty(pos),
    57         outputs(pos) = []; %remove 'default' from outputs
    58         outputs      = [outputs defaultoutputs(self,md)]; %add defaults
    59       end
    60       WriteData(fid,prefix,'data',outputs,'name','md.hydrology.requested_outputs','format','StringArray');
     51                        outputs = self.requested_outputs;
     52                        pos = find(ismember(outputs,'default'));
     53                        if ~isempty(pos),
     54                                outputs(pos) = []; %remove 'default' from outputs
     55                                outputs      = [outputs defaultoutputs(self,md)]; %add defaults
     56                        end
     57                        WriteData(fid,prefix,'data',outputs,'name','md.hydrology.requested_outputs','format','StringArray');
    6158                end % }}}
     59                function self = extrude(self,md) % {{{
     60                end % }}}
    6261                function savemodeljs(self,fid,modelname) % {{{
    6362               
    6463                        writejs1Darray(fid,[modelname '.hydrology.spcwatercolumn'],self.spcwatercolumn);

  • ../trunk-jpl/src/m/classes/initialization.py

     
    1010    """INITIALIZATION class definition
    1111
    1212    Usage:
    13     initialization = initialization()
     13        initialization = initialization()
    1414    """
    1515
    16     def __init__(self):  # {{{
     16    def __init__(self): #{{{
    1717        self.vx = np.nan
    1818        self.vy = np.nan
    1919        self.vz = np.nan
    2020        self.vel = np.nan
    21         self.enthalpy = np.nan
    2221        self.pressure = np.nan
    2322        self.temperature = np.nan
     23        self.enthalpy = np.nan
    2424        self.waterfraction = np.nan
    25         self.watercolumn = np.nan
    2625        self.sediment_head = np.nan
    2726        self.epl_head = np.nan
    2827        self.epl_thickness = np.nan
     28        self.watercolumn = np.nan
    2929        self.hydraulic_potential = np.nan
    3030        self.channelarea = np.nan
    3131        self.sealevel = np.nan
    3232        self.bottompressure = np.nan
     33        self.dsl = np.nan
     34        self.str = np.nan
    3335        self.sample = np.nan
    3436
    35         #set defaults
    3637        self.setdefaultparameters()
    3738    #}}}
    38     def __repr__(self):  # {{{
     39    def __repr__(self): #{{{
    3940        s = '   initial field values:\n'
    4041        s += '{}\n'.format(fielddisplay(self, 'vx', 'x component of velocity [m / yr]'))
    4142        s += '{}\n'.format(fielddisplay(self, 'vy', 'y component of velocity [m / yr]'))
    … …  
    5152        s += '{}\n'.format(fielddisplay(self, 'epl_thickness', 'thickness of the epl [m]'))
    5253        s += '{}\n'.format(fielddisplay(self, 'hydraulic_potential', 'Hydraulic potential (for GlaDS) [Pa]'))
    5354        s += '{}\n'.format(fielddisplay(self, 'channelarea', 'subglaciale water channel area (for GlaDS) [m2]'))
    54         s += '{}\n'.format(fielddisplay(self,'sample','Realization of a Gaussian random field'))
     55        s += '{}\n'.format(fielddisplay(self,'sample', 'Realization of a Gaussian random field'))
    5556        return s
    5657    #}}}
    57     def extrude(self, md):  # {{{
    58         self.vx = project3d(md, 'vector', self.vx, 'type', 'node')
    59         self.vy = project3d(md, 'vector', self.vy, 'type', 'node')
    60         self.vz = project3d(md, 'vector', self.vz, 'type', 'node')
    61         self.vel = project3d(md, 'vector', self.vel, 'type', 'node')
    62         self.temperature = project3d(md, 'vector', self.temperature, 'type', 'node')
    63         self.enthalpy = project3d(md, 'vector', self.enthalpy, 'type', 'node')
    64         self.waterfraction = project3d(md, 'vector', self.waterfraction, 'type', 'node')
    65         self.watercolumn = project3d(md, 'vector', self.watercolumn, 'type', 'node')
    66         self.sediment_head = project3d(md, 'vector', self.sediment_head, 'type', 'node', 'layer', 1)
    67         self.epl_head = project3d(md, 'vector', self.epl_head, 'type', 'node', 'layer', 1)
    68         self.epl_thickness = project3d(md, 'vector', self.epl_thickness, 'type', 'node', 'layer', 1)
    69         self.sealevel = project3d(md, 'vector', self.sealevel, 'type', 'node', 'layer', 1)
    70         self.bottompressure = project3d(md, 'vector', self.bottompressure, 'type', 'node', 'layer', 1)
    71 
    72         #Lithostatic pressure by default
    73         if np.ndim(md.geometry.surface) == 2:
    74             print('Reshaping md.geometry.surface for your convenience but you should fix it in your model set up')
    75             self.pressure = md.constants.g * md.materials.rho_ice * (md.geometry.surface.reshape(-1, ) - md.mesh.z)
    76         else:
    77             self.pressure = md.constants.g * md.materials.rho_ice * (md.geometry.surface - md.mesh.z)
    78 
    79         return self
     58    def setdefaultparameters(self): #{{{
     59        return
    8060    #}}}
    81     def setdefaultparameters(self):  # {{{
    82         return self
    83     #}}}
    84     def checkconsistency(self, md, solution, analyses):  # {{{
     61    def checkconsistency(self, md, solution, analyses): #{{{
    8562        if 'StressbalanceAnalysis' in analyses and not solution == 'TransientSolution' and not md.transient.isstressbalance:
    8663            if not np.any(np.logical_or(np.isnan(md.initialization.vx), np.isnan(md.initialization.vy))):
    8764                md = checkfield(md, 'fieldname', 'initialization.vx', 'NaN', 1, 'Inf', 1, 'size', [md.mesh.numberofvertices])
    … …  
    129106                md = checkfield(md, 'fieldname', 'initialization.channelarea', 'NaN', 1, 'Inf', 1, '>=', 0, 'size', [md.mesh.numberofelements])
    130107        if 'SamplingAnalysis' in analyses and not solution == 'TransientSolution' and not md.transient.issampling:
    131108            if np.any(np.isnan(md.initialization.sample)):
    132                 md = checkfield(md,'fieldname','initialization.sample','NaN',1,'Inf',1,'size',[md.mesh.numberofvertices])
     109                md = checkfield(md, 'fieldname', 'initialization.sample', 'NaN', 1,'Inf', 1, 'size', [md.mesh.numberofvertices])
    133110        return md
    134     # }}}
    135     def marshall(self, prefix, md, fid):  # {{{
     111    #}}}
     112    def marshall(self, prefix, md, fid): #{{{
    136113        yts = md.constants.yts
    137         WriteData(fid, prefix, 'object', self, 'fieldname', 'vx', 'format', 'DoubleMat', 'mattype', 1, 'scale', 1. / yts)
    138         WriteData(fid, prefix, 'object', self, 'fieldname', 'vy', 'format', 'DoubleMat', 'mattype', 1, 'scale', 1. / yts)
    139         WriteData(fid, prefix, 'object', self, 'fieldname', 'vz', 'format', 'DoubleMat', 'mattype', 1, 'scale', 1. / yts)
     114
     115        WriteData(fid, prefix, 'object', self, 'fieldname', 'vx', 'format', 'DoubleMat', 'mattype', 1, 'scale', 1 / yts)
     116        WriteData(fid, prefix, 'object', self, 'fieldname', 'vy', 'format', 'DoubleMat', 'mattype', 1, 'scale', 1 / yts)
     117        WriteData(fid, prefix, 'object', self, 'fieldname', 'vz', 'format', 'DoubleMat', 'mattype', 1, 'scale', 1 / yts)
    140118        WriteData(fid, prefix, 'object', self, 'fieldname', 'pressure', 'format', 'DoubleMat', 'mattype', 1)
    141         WriteData(fid, prefix, 'object', self, 'fieldname', 'sealevel', 'format', 'DoubleMat', 'mattype', 1, 'timeserieslength', md.mesh.numberofvertices + 1, 'yts', md.constants.yts)
     119        WriteData(fid, prefix, 'object', self, 'fieldname', 'sealevel', 'format', 'DoubleMat', 'mattype', 1,'timeserieslength', md.mesh.numberofvertices + 1, 'yts', yts)
    142120        WriteData(fid, prefix, 'object', self, 'fieldname', 'bottompressure', 'format', 'DoubleMat', 'mattype', 1)
     121        WriteData(fid, prefix, 'object', self, 'fieldname', 'str', 'format', 'DoubleMat', 'mattype', 1)
     122        WriteData(fid, prefix, 'object', self, 'fieldname', 'dsl', 'format', 'DoubleMat', 'mattype', 1)
    143123        WriteData(fid, prefix, 'object', self, 'fieldname', 'temperature', 'format', 'DoubleMat', 'mattype', 1)
    144124        WriteData(fid, prefix, 'object', self, 'fieldname', 'waterfraction', 'format', 'DoubleMat', 'mattype', 1)
    145125        WriteData(fid, prefix, 'object', self, 'fieldname', 'sediment_head', 'format', 'DoubleMat', 'mattype', 1)
    … …  
    148128        WriteData(fid, prefix, 'object', self, 'fieldname', 'watercolumn', 'format', 'DoubleMat', 'mattype', 1)
    149129        WriteData(fid, prefix, 'object', self, 'fieldname', 'channelarea', 'format', 'DoubleMat', 'mattype', 1)
    150130        WriteData(fid, prefix, 'object', self, 'fieldname', 'hydraulic_potential', 'format', 'DoubleMat', 'mattype', 1)
    151         WriteData(fid,prefix,'object',self,'fieldname','sample','format','DoubleMat','mattype',1)
     131        WriteData(fid, prefix, 'object', self, 'fieldname', 'sample', 'format', 'DoubleMat', 'mattype', 1)
     132
    152133        if md.thermal.isenthalpy:
    153134            if (np.size(self.enthalpy) <= 1):
     135                # Reconstruct enthalpy
    154136                tpmp = md.materials.meltingpoint - md.materials.beta * md.initialization.pressure
    155                 pos = np.nonzero(md.initialization.waterfraction > 0.)[0]
     137                pos = np.where(md.initialization.waterfraction > 0)[0]
    156138                self.enthalpy = md.materials.heatcapacity * (md.initialization.temperature - md.constants.referencetemperature)
    157                 self.enthalpy[pos] = md.materials.heatcapacity * (tpmp[pos].reshape(-1, ) - md.constants.referencetemperature) + md.materials.latentheat * md.initialization.waterfraction[pos].reshape(-1, )
     139                self.enthalpy[pos] = md.materials.heatcapacity * (tpmp[pos].reshape(-1, 1) - md.constants.referencetemperature) + md.materials.latentheat * md.initialization.waterfraction[pos].reshape(-1, 1)
    158140
    159141            WriteData(fid, prefix, 'data', self.enthalpy, 'format', 'DoubleMat', 'mattype', 1, 'name', 'md.initialization.enthalpy')
    160     # }}}
     142    #}}}
     143    def extrude(self, md): #{{{
     144        self.vx = project3d(md, 'vector', self.vx, 'type', 'node')
     145        self.vy = project3d(md, 'vector', self.vy, 'type', 'node')
     146        self.vz = project3d(md, 'vector', self.vz, 'type', 'node')
     147        self.vel = project3d(md, 'vector', self.vel, 'type', 'node')
     148        self.temperature = project3d(md, 'vector', self.temperature, 'type', 'node')
     149        self.enthalpy = project3d(md, 'vector', self.enthalpy, 'type', 'node')
     150        self.waterfraction = project3d(md, 'vector', self.waterfraction, 'type', 'node')
     151        self.watercolumn = project3d(md, 'vector', self.watercolumn, 'type', 'node')
     152        self.sediment_head = project3d(md, 'vector', self.sediment_head, 'type', 'node', 'layer', 1)
     153        self.epl_head = project3d(md, 'vector', self.epl_head, 'type', 'node', 'layer', 1)
     154        self.epl_thickness = project3d(md, 'vector', self.epl_thickness, 'type', 'node', 'layer', 1)
     155        self.sealevel = project3d(md, 'vector', self.sealevel, 'type', 'node', 'layer', 1)
     156        self.bottompressure = project3d(md, 'vector', self.bottompressure, 'type', 'node', 'layer', 1)
     157
     158        # Lithostatic pressure by default
     159        if np.ndim(md.geometry.surface) == 2:
     160            print('Reshaping md.geometry.surface for your convenience but you should fix it in your model set up')
     161            self.pressure = md.constants.g * md.materials.rho_ice * (md.geometry.surface.reshape(-1, 1) - md.mesh.z)
     162        else:
     163            self.pressure = md.constants.g * md.materials.rho_ice * (md.geometry.surface - md.mesh.z)
     164
     165        return self
     166    #}}}

  • ../trunk-jpl/src/m/classes/lovenumbers.py

     
    99    """LOVENUMBERS class definition
    1010
    1111    Usage:
    12         lovenumbers = lovenumbers()  #will setup love numbers deg 1001 by default
    13         lovenumbers = lovenumbers('maxdeg', 10001);   #supply numbers of degrees required (here, 10001)
     12        lovenumbers = lovenumbers()
     13        lovenumbers = lovenumbers('maxdeg', 10000, 'referenceframe', 'CF');
     14
     15    Choose numbers of degrees required (1000 by default) and reference frame
     16    (between CF and CM; CM by default)
    1417    """
    1518
    1619    def __init__(self, *args):  #{{{
    … …  
    2528        self.tl = []
    2629        self.tk2secular = 0  # deg 2 secular number
    2730
     31        # Time/frequency for visco-elastic love numbers
     32        self.timefreq = []
     33        self.istime = 1
     34
    2835        options = pairoptions(*args)
    2936        maxdeg = options.getfieldvalue('maxdeg', 1000)
    3037        referenceframe = options.getfieldvalue('referenceframe', 'CM')
    … …  
    3138        self.setdefaultparameters(maxdeg, referenceframe)
    3239    #}}}
    3340
     41    def __repr__(self):  #{{{
     42        s = '   lovenumbers parameters:\n'
     43        s += '{}\n'.format(fielddisplay(self, 'h', 'load Love number for radial displacement'))
     44        s += '{}\n'.format(fielddisplay(self, 'k', 'load Love number for gravitational potential perturbation'))
     45        s += '{}\n'.format(fielddisplay(self, 'l', 'load Love number for horizontal displacements'))
     46        s += '{}\n'.format(fielddisplay(self, 'th', 'tidal load Love number (deg 2)'))
     47        s += '{}\n'.format(fielddisplay(self, 'tk', 'tidal load Love number (deg 2)'))
     48        s += '{}\n'.format(fielddisplay(self, 'tl', 'tidal load Love number (deg 2)'))
     49        s += '{}\n'.format(fielddisplay(self, 'tk2secular', 'secular fluid Love number'))
     50        s += '{}\n'.format(fielddisplay(self, 'istime', 'time (default: 1) or frequency love numbers (0)'))
     51        s += '{}\n'.format(fielddisplay(self, 'timefreq', 'time/frequency vector (yr or 1/yr)'))
     52        return s
     53    #}}}
     54
    3455    def setdefaultparameters(self, maxdeg, referenceframe):  #{{{
    3556        # Initialize love numbers
    3657        self.h = getlovenumbers('type', 'loadingverticaldisplacement', 'referenceframe', referenceframe, 'maxdeg', maxdeg)
    … …  
    4263
    4364        # Secular fluid love number
    4465        self.tk2secular = 0.942
     66
     67        # Time
     68        self.istime = 1 # Temporal love numbers by default
     69        self.timefreq = 0 # Elastic case by default
    4570        return self
    4671    #}}}
    4772
    4873    def checkconsistency(self, md, solution, analyses):  #{{{
    49         if ('SealevelriseAnalysis' not in analyses) or (solution == 'TransientSolution' and not md.transient.isslc):
     74        if ('SealevelchangeAnalysis' not in analyses) or (solution == 'TransientSolution' and not md.transient.isslc):
    5075            return
    5176        md = checkfield(md, 'fieldname', 'solidearth.lovenumbers.h', 'NaN', 1, 'Inf', 1)
    5277        md = checkfield(md, 'fieldname', 'solidearth.lovenumbers.k', 'NaN', 1, 'Inf', 1)
    … …  
    5580        md = checkfield(md, 'fieldname', 'solidearth.lovenumbers.th', 'NaN', 1, 'Inf', 1)
    5681        md = checkfield(md, 'fieldname', 'solidearth.lovenumbers.tk', 'NaN', 1, 'Inf', 1)
    5782        md = checkfield(md, 'fieldname', 'solidearth.lovenumbers.tk2secular', 'NaN', 1, 'Inf', 1)
     83        md = checkfield(md, 'fieldname', 'solidearth.lovenumbers.timefreq', 'NaN', 1, 'Inf', 1)
     84        md = checkfield(md, 'fieldname', 'solidearth.lovenumbers.istime', 'NaN', 1, 'Inf', 1, 'values', [0, 1])
     85
    5886        # Check that love numbers are provided at the same level of accuracy
    5987        if (self.h.shape[0] != self.k.shape[0]) or (self.h.shape[0] != self.l.shape[0]):
    6088            raise ValueError('lovenumbers error message: love numbers should be provided at the same level of accuracy')
     89
     90        ntf = len(self.timefreq)
     91        if (self.h.shape[1] != ntf or self.k.shape[1] != ntf or self.l.shape[1] != ntf or self.th.shape[1] != ntf or self.tk.shape[1] != ntf or self.tl.shape[1] != ntf):
     92            raise ValueError('lovenumbers error message: love numbers should have as many time/frequency steps as the time/frequency vector')
     93
    6194        return md
    6295    #}}}
    6396
    … …  
    6598        return[]
    6699    #}}}
    67100
    68     def __repr__(self):  #{{{
    69         s = '   lovenumbers parameters:\n'
    70         s += '{}\n'.format(fielddisplay(self, 'h', 'load Love number for radial displacement'))
    71         s += '{}\n'.format(fielddisplay(self, 'k', 'load Love number for gravitational potential perturbation'))
    72         s += '{}\n'.format(fielddisplay(self, 'l', 'load Love number for horizontal displacements'))
    73         s += '{}\n'.format(fielddisplay(self, 'th', 'tidal load Love number (deg 2)'))
    74         s += '{}\n'.format(fielddisplay(self, 'tk', 'tidal load Love number (deg 2)'))
    75         s += '{}\n'.format(fielddisplay(self, 'tk2secular', 'secular fluid Love number'))
    76         return s
    77     #}}}
    78 
    79101    def marshall(self, prefix, md, fid):  #{{{
    80102        WriteData(fid, prefix, 'object', self, 'fieldname', 'h', 'name', 'md.solidearth.lovenumbers.h', 'format', 'DoubleMat', 'mattype', 1)
    81103        WriteData(fid, prefix, 'object', self, 'fieldname', 'k', 'name', 'md.solidearth.lovenumbers.k', 'format', 'DoubleMat', 'mattype', 1)
    … …  
    85107        WriteData(fid, prefix, 'object', self, 'fieldname', 'tk', 'name', 'md.solidearth.lovenumbers.tk', 'format', 'DoubleMat', 'mattype', 1)
    86108        WriteData(fid, prefix, 'object', self, 'fieldname', 'tl', 'name', 'md.solidearth.lovenumbers.tl', 'format', 'DoubleMat', 'mattype', 1)
    87109        WriteData(fid, prefix, 'object', self, 'data', self.tk2secular, 'fieldname', 'lovenumbers.tk2secular', 'format', 'Double')
     110
     111        if (self.istime):
     112            scale = md.constants.yts
     113        else:
     114            scale = 1.0 / md.constants.yts
     115        WriteData(fid, prefix, 'object', self, 'fieldname', 'istime', 'name', 'md.solidearth.lovenumbers.istime', 'format', 'Boolean')
     116        WriteData(fid, prefix, 'object', self, 'fieldname', 'timefreq', 'name', 'md.solidearth.lovenumbers.timefreq', 'format', 'DoubleMat', 'mattype', 1, 'scale', scale);
    88117    #}}}
    89118
    90119    def extrude(self, md):  #{{{

  • ../trunk-jpl/src/m/classes/materials.py

     
    2424            if not(self.nature[i] == 'litho' or self.nature[i] == 'ice' or self.nature[i] == 'hydro'):
    2525                raise RuntimeError("materials constructor error message: nature of the material not supported yet! ('ice' or 'litho' or 'hydro')")
    2626
     27        # Start filling in the dynamic fields (not truly dynamic under Python)
    2728        for i in range(len(self.nature)):
    2829            nat = self.nature[i]
    2930            if nat == 'ice':
    … …  
    6667                raise RuntimeError("materials constructor error message: nature of the material not supported yet! ('ice' or 'litho' or 'hydro')")
    6768        self.earth_density = 0
    6869
    69         # Set default parameters
    7070        self.setdefaultparameters()
    7171    #}}}
    7272
    … …  
    104104                s += '{}\n'.format(fielddisplay(self, 'ebm_alpha', 'array describing each layer\'s exponent parameter controlling the shape of shear modulus curve between taul and tauh, only for EBM rheology (numlayers)'))
    105105                s += '{}\n'.format(fielddisplay(self, 'ebm_delta', 'array describing each layer\'s amplitude of the transient relaxation (ratio between elastic rigity to pre-maxwell relaxation rigity), only for EBM rheology (numlayers)'))
    106106                s += '{}\n'.format(fielddisplay(self, 'ebm_taul', 'array describing each layer\'s starting period for transient relaxation, only for EBM rheology  (numlayers) [s]'))
    107                 s += '{}\n'.format(fielddisplay(self, 'ebm_tauh', 'array describing each layer''s array describing each layer\'s end period for transient relaxation, only for Burgers rheology  (numlayers) [s]'))
     107                s += '{}\n'.format(fielddisplay(self, 'ebm_tauh', 'array describing each layer''s array describing each layer\'s end period for transient relaxation, only for Burgers rheology (numlayers) [s]'))
    108108
    109109                s += '{}\n'.format(fielddisplay(self, 'rheologymodel', 'array describing whether we adopt a Maxwell (0), Burgers (1) or EBM (2) rheology (default: 0)'))
    110110                s += '{}\n'.format(fielddisplay(self, 'density', 'array describing each layer\'s density (numlayers) [kg/m^3]'))
    111                 s += '{}\n'.format(fielddisplay(self, 'issolid', 'array describing whether the layer is solid or liquid (default 1) (numlayers)'))
     111                s += '{}\n'.format(fielddisplay(self, 'issolid', 'array describing whether the layer is solid or liquid (default: 1) (numlayers)'))
    112112            elif nat == 'hydro':
    113113                s += 'Hydro:\n'
    114114                s += '{}\n'.format(fielddisplay(self, 'rho_ice', 'ice density [kg/m^3]'))
    … …  
    126126            nat = self.nature[i]
    127127            if nat == 'ice':
    128128                # Ice density (kg/m^3)
    129                 self.rho_ice = 917.
     129                self.rho_ice = 917
    130130
    131131                # Ocean water density (kg/m^3)
    132                 self.rho_water = 1023.
     132                self.rho_water = 1023
    133133
    134134                # Fresh water density (kg/m^3)
    135                 self.rho_freshwater = 1000.
     135                self.rho_freshwater = 1000
    136136
    137137                # Water viscosity (N.s/m^2)
    138138                self.mu_water = 0.001787
    139139
    140140                # Ice heat capacity cp (J/kg/K)
    141                 self.heatcapacity = 2093.
     141                self.heatcapacity = 2093
    142142
    143143                # Ice latent heat of fusion L (J/kg)
    144144                self.latentheat = 3.34 * 1e5
    … …  
    159159                self.beta = 9.8 * 1e-8
    160160
    161161                # Mixed layer (ice-water interface) heat capacity (J/kg/K)
    162                 self.mixed_layer_capacity = 3974.
     162                self.mixed_layer_capacity = 3974
    163163
    164164                # Thermal exchange velocity (ice-water interface) (m/s)
    165165                self.thermal_exchange_velocity = 1.00 * 1e-4
    … …  
    192192                self.ebm_taul = [np.nan, np.nan]
    193193                self.ebm_tauh = [np.nan, np.nan]
    194194                self.rheologymodel = [0, 0]
    195                 self.density = [5.51e3, 5.50e3]  # (Pa) # Mantle and lithosphere density [kg/m^3]
    196                 self.issolid = [1, 1]  # Is layer solid or liquid?
     195                self.density = [5.51e3, 5.50e3] # (Pa) # Mantle and lithosphere density [kg/m^3]
     196                self.issolid = [1, 1] # Is layer solid or liquid?
    197197            elif nat == 'hydro':
    198198                # Ice density (kg/m^3)
    199                 self.rho_ice = 917.
     199                self.rho_ice = 917
    200200
    201201                # Ocean water density (kg/m^3)
    202                 self.rho_water = 1023.
     202                self.rho_water = 1023
    203203
    204204                # Fresh water density (kg/m^3)
    205                 self.rho_freshwater = 1000.
     205                self.rho_freshwater = 1000
    206206            else:
    207207                raise RuntimeError("materials setdefaultparameters error message: nature of the material not supported yet! ('ice' or 'litho' or 'hydro')")
    208208
    … …  
    250250                    raise RuntimeError('First layer must be solid (issolid[0] > 0 AND lame_mu[0] > 0). Add a weak inner core if necessary.')
    251251                ind = np.where(md.materials.issolid == 0)[0]
    252252                if np.sum(np.in1d(np.diff(ind),1) >= 1): # If there are at least two consecutive indices that contain issolid = 0
    253                     raise RuntimeError('Fluid layers detected at layers #{} but having 2 or more adjacent fluid layers is not supported yet. Consider merging them.'.format(i))
     253                    raise RuntimeError('Fluid layers detected at layers #' + indices + ', but having 2 or more adjacent fluid layers is not supported yet. Consider merging them.')
    254254
    255255            elif nat == 'hydro':
    256256                md = checkfield(md, 'fieldname', 'materials.rho_ice', '>', 0)
    … …  
    266266    def marshall(self, prefix, md, fid): #{{{
    267267        #1: MatdamageiceEnum 2: MatestarEnum 3: MaticeEnum 4: MatenhancediceEnum 5: MaterialsEnum
    268268        WriteData(fid, prefix, 'name', 'md.materials.nature', 'data', naturetointeger(self.nature), 'format', 'IntMat', 'mattype', 3)
    269         WriteData(fid, prefix, 'name', 'md.materials.type', 'data', 5, 'format', 'Integer') #DANGER, this can evolve if you have classes
     269        WriteData(fid, prefix, 'name', 'md.materials.type', 'data', 5, 'format', 'Integer') #DANGER: this can evolve if you have classes
    270270        for i in range(len(self.nature)):
    271271            nat = self.nature[i]
    272272            if nat == 'ice':
    … …  
    306306                for i in range(self.numlayers):
    307307                    earth_density = earth_density + (pow(self.radius[i + 1], 3) - pow(self.radius[i], 3)) * self.density[i]
    308308                earth_density = earth_density / pow(self.radius[self.numlayers], 3)
     309                print(earth_density)
    309310                self.earth_density = earth_density
    310311            elif nat == 'hydro':
    311312                WriteData(fid, prefix, 'object', self, 'class', 'materials', 'fieldname', 'rho_ice', 'format', 'Double')

  • ../trunk-jpl/src/m/classes/solidearth.m

     
    4444                                        error('solidearth constructor error message: zero or one argument only!');
    4545                        end
    4646                end % }}}
     47                function disp(self) % {{{
     48                        disp(sprintf('   solidearth inputs, forcings and settings:'));
     49
     50                        fielddisplay(self,'planetradius','planet radius [m]');
     51                        fielddisplay(self,'transitions','indices into parts of the mesh that will be icecaps');
     52                        fielddisplay(self,'requested_outputs','additional outputs requested');
     53                        fielddisplay(self,'partitionice','ice partition vector for barystatic contribution');
     54                        fielddisplay(self,'partitionhydro','hydro partition vector for barystatic contribution');
     55                        fielddisplay(self,'partitionocean','ocean partition vector for barystatic contribution');
     56                        if isempty(self.external), fielddisplay(self,'external','external solution, of the type solidearthsolution'); end
     57                        self.settings.disp();
     58                        self.lovenumbers.disp();
     59                        self.rotational.disp();
     60                        if ~isempty(self.external),
     61                                self.external.disp();
     62                        end
     63
     64                end % }}}
    4765                function self = setdefaultparameters(self,planet) % {{{
    4866
    4967                        %output default:
    … …  
    86104                function list=defaultoutputs(self,md) % {{{
    87105                        list = {'Sealevel'};
    88106                end % }}}
    89                 function disp(self) % {{{
    90                         disp(sprintf('   solidearth inputs, forcings and settings:'));
    91 
    92                         fielddisplay(self,'planetradius','planet radius [m]');
    93                         fielddisplay(self,'transitions','indices into parts of the mesh that will be icecaps');
    94                         fielddisplay(self,'requested_outputs','additional outputs requested');
    95                         fielddisplay(self,'partitionice','ice partition vector for barystatic contribution');
    96                         fielddisplay(self,'partitionhydro','hydro partition vector for barystatic contribution');
    97                         fielddisplay(self,'partitionocean','ocean partition vector for barystatic contribution');
    98                         if isempty(self.external), fielddisplay(self,'external','external solution, of the type solidearthsolution'); end
    99                         self.settings.disp();
    100                         self.lovenumbers.disp();
    101                         self.rotational.disp();
    102                         if ~isempty(self.external),
    103                                 self.external.disp();
    104                         end
    105 
    106                 end % }}}
    107107                function marshall(self,prefix,md,fid) % {{{
    108108
    109109                        WriteData(fid,prefix,'object',self,'fieldname','planetradius','format','Double');

  • ../trunk-jpl/src/m/classes/solidearthsettings.py

     
    120120            if md.mesh.__class__.__name__ is 'mesh3dsurface':
    121121                if self.grdmodel == 2:
    122122                    raise RuntimeException('model requires a 2D mesh to run gia Ivins computations (change mesh from mesh3dsurface to mesh2d)')
    123                 else:
    124                     if self.grdmodel == 1:
    125                         raise RuntimeException('model requires a 3D surface mesh to run GRD computations (change mesh from mesh2d to mesh3dsurface)')
     123            else:
     124                if self.grdmodel == 1:
     125                    raise RuntimeException('model requires a 3D surface mesh to run GRD computations (change mesh from mesh2d to mesh3dsurface)')
    126126            if self.sealevelloading and not self.grdocean:
    127127                raise RuntimeException('solidearthsettings checkconsistency error message: need grdocean on if sealevelloading flag is set')
    128128
    … …  
    133133    #}}}
    134134
    135135    def marshall(self, prefix, md, fid): #{{{
    136         WriteData(fid, prefix, 'object', self, 'fieldname', 'reltol', 'name', 'md.solidearth.settings.reltol', 'format', 'Double')
    137         WriteData(fid, prefix, 'object', self, 'fieldname', 'abstol', 'name', 'md.solidearth.settings.abstol', 'format', 'Double')
    138         WriteData(fid, prefix, 'object', self, 'fieldname', 'maxiter', 'name', 'md.solidearth.settings.maxiter', 'format', 'Integer')
    139         WriteData(fid, prefix, 'object', self, 'fieldname', 'selfattraction', 'name', 'md.solidearth.settings.selfattraction', 'format', 'Boolean')
    140         WriteData(fid, prefix, 'object', self, 'fieldname', 'elastic', 'name', 'md.solidearth.settings.elastic', 'format', 'Boolean')
    141         WriteData(fid, prefix, 'object', self, 'fieldname', 'viscous', 'name', 'md.solidearth.settings.viscous', 'format', 'Boolean')
    142         WriteData(fid, prefix, 'object', self, 'fieldname', 'rotation', 'name', 'md.solidearth.settings.rotation', 'format', 'Boolean')
    143         WriteData(fid, prefix, 'object', self, 'fieldname', 'grdocean', 'name', 'md.solidearth.settings.grdocean', 'format', 'Boolean')
    144         WriteData(fid, prefix, 'object', self, 'fieldname', 'ocean_area_scaling', 'name', 'md.solidearth.settings.ocean_area_scaling', 'format', 'Integer')
    145         WriteData(fid, prefix, 'object', self, 'fieldname', 'runfrequency', 'name', 'md.solidearth.settings.runfrequency', 'format', 'Integer')
    146         WriteData(fid, prefix, 'object', self, 'fieldname', 'degacc', 'name', 'md.solidearth.settings.degacc', 'format', 'Double')
    147         WriteData(fid, prefix, 'object', self, 'fieldname', 'timeacc', 'name', 'md.solidearth.settings.timeacc', 'format', 'Double', 'scale', md.constants.yts)
    148         WriteData(fid, prefix, 'object', self, 'fieldname', 'horiz', 'name', 'md.solidearth.settings.horiz', 'format', 'Integer')
    149         WriteData(fid, prefix, 'object', self, 'fieldname', 'sealevelloading', 'name', 'md.solidearth.settings.sealevelloading', 'format', 'Integer')
    150         WriteData(fid, prefix, 'object', self, 'fieldname','isgrd', 'name', 'md.solidearth.settings.isgrd', 'format', 'Integer')
    151         WriteData(fid, prefix, 'object', self, 'fieldname', 'compute_bp_grd', 'name', 'md.solidearth.settings.compute_bp_grd', 'format', 'Integer')
    152         WriteData(fid, prefix, 'object', self, 'fieldname', 'grdmodel', 'name', 'md.solidearth.settings.grdmodel', 'format', 'Integer')
    153         WriteData(fid, prefix, 'object', self, 'fieldname', 'cross_section_shape', 'name', 'md.solidearth.settings.cross_section_shape', 'format', 'Integer')
     136        WriteData(fid, prefix, 'object', self, 'fieldname', 'reltol', 'name', 'md.solidearth.settings.reltol', 'format', 'Double');
     137        WriteData(fid, prefix, 'object', self, 'fieldname', 'abstol', 'name', 'md.solidearth.settings.abstol', 'format', 'Double');
     138        WriteData(fid, prefix, 'object', self, 'fieldname', 'maxiter', 'name', 'md.solidearth.settings.maxiter', 'format', 'Integer');
     139        WriteData(fid, prefix, 'object', self, 'fieldname', 'selfattraction', 'name', 'md.solidearth.settings.selfattraction', 'format', 'Boolean');
     140        WriteData(fid, prefix, 'object', self, 'fieldname', 'elastic', 'name', 'md.solidearth.settings.elastic', 'format', 'Boolean');
     141        WriteData(fid, prefix, 'object', self, 'fieldname', 'viscous', 'name', 'md.solidearth.settings.viscous', 'format', 'Boolean');
     142        WriteData(fid, prefix, 'object', self, 'fieldname', 'rotation', 'name', 'md.solidearth.settings.rotation', 'format', 'Boolean');
     143        WriteData(fid, prefix, 'object', self, 'fieldname', 'grdocean', 'name', 'md.solidearth.settings.grdocean', 'format', 'Boolean');
     144        WriteData(fid, prefix, 'object', self, 'fieldname', 'ocean_area_scaling', 'name', 'md.solidearth.settings.ocean_area_scaling', 'format', 'Boolean');
     145        WriteData(fid, prefix, 'object', self, 'fieldname', 'runfrequency', 'name', 'md.solidearth.settings.runfrequency', 'format', 'Integer');
     146        WriteData(fid, prefix, 'object', self, 'fieldname', 'degacc', 'name', 'md.solidearth.settings.degacc', 'format', 'Double');
     147        WriteData(fid, prefix, 'object', self, 'fieldname', 'timeacc', 'name', 'md.solidearth.settings.timeacc', 'format', 'Double', 'scale',md.constants.yts);
     148        WriteData(fid, prefix, 'object', self, 'fieldname', 'horiz', 'name', 'md.solidearth.settings.horiz', 'format', 'Integer');
     149        WriteData(fid, prefix, 'object', self, 'fieldname', 'sealevelloading', 'name', 'md.solidearth.settings.sealevelloading', 'format', 'Integer');
     150        WriteData(fid, prefix, 'object', self, 'fieldname', 'isgrd', 'name', 'md.solidearth.settings.isgrd', 'format', 'Integer');
     151        WriteData(fid, prefix, 'object', self, 'fieldname', 'compute_bp_grd', 'name', 'md.solidearth.settings.compute_bp_grd', 'format', 'Integer');
     152        WriteData(fid, prefix, 'object', self, 'fieldname', 'grdmodel', 'name', 'md.solidearth.settings.grdmodel', 'format', 'Integer');
     153        WriteData(fid, prefix, 'object', self, 'fieldname', 'cross_section_shape', 'name', 'md.solidearth.settings.cross_section_shape', 'format', 'Integer');
    154154    #}}}
    155155
    156156    def extrude(self, md): #{{{

  • ../trunk-jpl/src/m/classes/solidearthsolution.py

     
     1import numpy as np
     2
     3from checkfield import checkfield
     4from fielddisplay import fielddisplay
     5from WriteData import WriteData
     6
     7
     8class solidearthsolution(object):
     9    """SOLIDEARTHSOLUTION class definition
     10
     11    Usage:
     12        solidearthsolution = solidearthsolution()
     13    """
     14
     15    def __init__(self, *args): #{{{
     16        self.displacementeast = None
     17        self.displacementnorth = None
     18        self.displacementup = None
     19        self.geoid = None
     20
     21        if len(args) == 0:
     22            self.setdefaultparameters()
     23        else:
     24            raise RuntimeError('constructor not supported')
     25    #}}}
     26
     27    def __repr__(self): #{{{
     28        s = '         units for time series is (yr)\n'
     29        s += '{}\n'.format(fielddisplay(self, 'displacementeast', 'solid-Earth Eastwards bedrock displacement series (m)'))
     30        s += '{}\n'.format(fielddisplay(self, 'displacementnorth', 'solid-Earth Northwards bedrock displacement time series (m)'))
     31        s += '{}\n'.format(fielddisplay(self, 'displacementup', 'solid-Earth bedrock uplift time series (m)'))
     32        s += '{}\n'.format(fielddisplay(self, 'geoid', 'solid-Earth geoid time series (m)'))
     33
     34        return s
     35    #}}}
     36
     37    def setdefaultparameters(self): #{{{
     38        self.displacementeast = []
     39        self.displacementnorth = []
     40        self.displacementup = []
     41        self.geoid = []
     42    #}}}
     43
     44    def checkconsistency(self, md, solution, analyses): #{{{
     45        md = checkfield(md, 'fieldname', 'solidearth.external.displacementeast', 'Inf', 1, 'timeseries', 1)
     46        md = checkfield(md, 'fieldname', 'solidearth.external.displacementnorth', 'Inf', 1, 'timeseries', 1)
     47        md = checkfield(md, 'fieldname', 'solidearth.external.displacementup', 'Inf', 1, 'timeseries', 1)
     48        md = checkfield(md, 'fieldname', 'solidearth.external.geoid', 'Inf', 1, 'timeseries', 1)
     49    #}}}
     50
     51    def marshall(self, prefix, md, fid): #{{{
     52        yts = md.constants.yts
     53
     54        # Transform our time series into time series rates
     55        if np.shape(self.displacementeast, 1) == 1:
     56            print('External solidearthsolution warning: only one time step provided, assuming the values are rates per year')
     57            displacementeast_rate = np.append(np.array(self.displacementeast).reshape(-1, 1), 0)
     58            displacementnorth_rate = np.append(np.array(self.displacementnorth).reshape(-1, 1), 0)
     59            displacementup_rate = np.append(np.array(self.displacementup).reshape(-1, 1), 0)
     60            geoid_rate = np.append(np.array(self.geoid).reshape(-1, 1), 0)
     61        else:
     62            time = self.displacementeast[-1, :]
     63            dt = np.diff(time, 1, 1)
     64            displacementeast_rate = np.diff(self.displacementeast[0:-2, :], 1, 1) / dt
     65            displacementeast_rate.append(time[0:-2])
     66            displacementnorth_rate = np.diff(self.displacementnorth[0:-2, :], 1, 1) / dt
     67            displacementnorth_rate.append(time[0:-2])
     68            displacementup_rate = np.diff(self.displacementup[0:-2, :], 1, 1) / dt
     69            displacementup_rate.append(time[0:-2])
     70            geoid_rate = np.diff(self.geoid[0:-2, :], 1, 1) / dt
     71            geoid_rate.append(time[0:-2])
     72        WriteData(fid, prefix, 'object', self, 'fieldname', 'displacementeast', 'data', displacementeast_rate, 'format', 'DoubleMat', 'name', 'md.solidearth.external.displacementeast', 'mattype', 1, 'scale', 1 / yts,'timeserieslength', md.mesh.numberofvertices + 1, 'yts', yts);
     73        WriteData(fid, prefix, 'object', self, 'fieldname', 'displacementup', 'data', displacementup_rate,'format', 'DoubleMat', 'name', 'md.solidearth.external.displacementup', 'mattype', 1, 'scale', 1 / yts,'timeserieslength', md.mesh.numberofvertices + 1, 'yts', yts);
     74        WriteData(fid, prefix, 'object', self, 'fieldname', 'displacementnorth', 'data', displacementnorth_rate,'format', 'DoubleMat', 'name', 'md.solidearth.external.displacementnorth', 'mattype', 1, 'scale', 1 / yts,'timeserieslength', md.mesh.numberofvertices + 1, 'yts', yts);
     75        WriteData(fid, prefix, 'object', self, 'fieldname', 'geoid', 'data', geoid_rate,'format', 'DoubleMat', 'name', 'md.solidearth.external.geoid', 'mattype', 1, 'scale', 1 / yts,'timeserieslength', md.mesh.numberofvertices + 1, 'yts', yts);
     76    #}}}
     77
     78    def extrude(self, md): #{{{
     79        return self
     80    #}}}

  • ../trunk-jpl/test/NightlyRun/test2002.m

     
    1414%solidearth loading:  {{{
    1515md.masstransport.spcthickness=[md.geometry.thickness;0];
    1616md.smb.mass_balance=zeros(md.mesh.numberofvertices,1);
     17
    1718%antarctica
    1819xe=md.mesh.x(md.mesh.elements)*[1;1;1]/3;
    1920ye=md.mesh.y(md.mesh.elements)*[1;1;1]/3;
    … …  
    2526pos=find(late < -80);
    2627md.masstransport.spcthickness(md.mesh.elements(pos,:))= md.masstransport.spcthickness(md.mesh.elements(pos,:))-100;
    2728posant=pos;
     29
    2830%greenland
    2931pos=find(late>60 & late<90 & longe>-75 & longe<-15);
    3032md.masstransport.spcthickness(md.mesh.elements(pos,:))= md.masstransport.spcthickness(md.mesh.elements(pos,:))-100;
    … …  
    9597md=solve(md,'Transient');
    9698Seustatic=md.results.TransientSolution.Sealevel;
    9799Beustatic=md.results.TransientSolution.Bed;
     100pause
    98101
    99102%eustatic + selfattraction run:
    100103md.solidearth.settings.selfattraction=1;

  • ../trunk-jpl/test/NightlyRun/test2002.py

     
    11#Test Name: EarthSlc
    22import numpy as np
    33
     4from MatlabFuncs import *
     5
    46from gmshplanet import *
    57from gmtmask import *
    68from lovenumbers import *
    … …  
    1113from solve import *
    1214
    1315
    14 #mesh earth:
     16# Mesh earth
     17#
     18# NOTE: In MATLAB, we currently use cached mesh to account for differences in
     19# mesh generated under Linux versus under macOS
     20#
    1521md = model()
    16 md.mesh = gmshplanet('radius', 6.371012 * 1e3, 'resolution', 700.) #700 km resolution mesh
     22md.mesh = gmshplanet('radius', 6.371012 * 1e3, 'resolution', 700.) # 700 km resolution mesh
    1723
    18 #parameterize solidearth solution:
    19 #solidearth loading:
    20 md.solidearth.surfaceload.icethicknesschange = np.zeros((md.mesh.numberofelements, 1))
    21 md.solidearth.initialsealevel = np.zeros((md.mesh.numberofvertices, 1))
    22 md.dsl.global_average_thermosteric_sea_level_change = np.zeros((2, 1))
    23 md.dsl.sea_surface_height_change_above_geoid = np.zeros((md.mesh.numberofvertices + 1, 1))
    24 md.dsl.sea_water_pressure_change_at_sea_floor = np.zeros((md.mesh.numberofvertices + 1, 1))
     24# Geometry for the bed, arbitrary thickness of 100
     25md.geometry.bed = np.zeros((md.mesh.numberofvertices, 1))
     26md.geometry.base = md.geometry.bed
     27md.geometry.thickness = 100 * np.ones((md.mesh.numberofvertices, 1))
     28md.geometry.surface = md.geometry.bed + md.geometry.thickness
    2529
    26 #antarctica
    27 late = md.mesh.lat[md.mesh.elements - 1].sum(axis=1) / 3
    28 longe = md.mesh.long[md.mesh.elements - 1].sum(axis=1) / 3
     30# Solidearth loading #{{{
     31md.masstransport.spcthickness = np.append(md.geometry.thickness, 0)
     32md.smb.mass_balance = np.zeros((md.mesh.numberofvertices, 1))
     33
     34# Antarctica
     35xe = md.mesh.x[md.mesh.elements - 1].sum(axis=1) / 3
     36ye = md.mesh.y[md.mesh.elements - 1].sum(axis=1) / 3
     37ze = md.mesh.z[md.mesh.elements - 1].sum(axis=1) / 3
     38re = pow((pow(xe, 2) + pow(ye, 2) + pow(ze, 2)), 0.5)
     39
     40late = asind(ze / re)
     41longe = atan2d(ye, xe)
    2942pos = np.where(late < -80)[0]
    30 md.solidearth.surfaceload.icethicknesschange[pos] = -100
    31 #greenland
    32 pos = np.where(np.logical_and.reduce((late > 70, late < 80, longe > -60, longe < -30)))[0]
    33 md.solidearth.surfaceload.icethicknesschange[pos] = -100
     43md.masstransport.spcthickness[md.mesh.elements[pos]] = md.masstransport.spcthickness[md.mesh.elements[pos]] - 100
     44posant = pos
    3445
    35 #elastic loading from love numbers:
    36 md.solidearth.lovenumbers = lovenumbers('maxdeg', 100)
     46# Greenland
     47pos = np.where(np.logical_and.reduce((late > 60, late < 90, longe > -75, longe < -15)))[0]
     48md.masstransport.spcthickness[md.mesh.elements[pos]] = md.masstransport.spcthickness[md.mesh.elements[pos]] - 100
     49posgre = pos
     50
     51# Elastic loading from love numbers:
     52md.solidearth.lovenumbers = lovenumbers('maxdeg', 1000)
    3753#}}}
    3854
    39 #mask: {{{
     55# Mask: {{{
    4056mask = gmtmask(md.mesh.lat, md.mesh.long)
    41 icemask = np.ones((md.mesh.numberofvertices, 1))
     57oceanmask = -1 * np.ones((md.mesh.numberofvertices, 1))
    4258pos = np.where(mask == 0)[0]
    43 icemask[pos] = -1
    44 pos = np.where(mask[md.mesh.elements - 1].sum(axis=1) < 3)[0]
    45 icemask[md.mesh.elements[pos, :] - 1] = -1
    46 md.mask.ice_levelset = icemask
    47 md.mask.ocean_levelset = -icemask
     59oceanmask[pos] = 1
    4860
    49 #make sure that the elements that have loads are fully grounded
    50 pos = np.nonzero(md.solidearth.surfaceload.icethicknesschange)[0]
    51 md.mask.ocean_levelset[md.mesh.elements[pos, :] - 1] = 1
     61icemask = np.ones((md.mesh.numberofvertices, 1))
     62icemask[md.mesh.elements[posant]] = -1
     63icemask[md.mesh.elements[posgre]] = -1
    5264
    53 #make sure wherever there is an ice load, that the mask is set to ice:
    54 #pos = np.nonzero(md.solidearth.surfaceload.icethicknesschange)[0] # TODO: Do we need to do this twice?
    55 md.mask.ice_levelset[md.mesh.elements[pos, :] - 1] = -1
    56 # }}}
     65md.mask.ice_levelset = icemask
     66md.mask.ocean_levelset = oceanmask
     67#}}}
    5768
    58 md.solidearth.settings.ocean_area_scaling = 0
     69# Time stepping {{{
     70md.timestepping.start_time = 0
     71md.timestepping.time_step = 1
     72md.timestepping.final_time = 1
     73#}}}
    5974
    60 #geometry for the bed, arbitrary
    61 md.geometry.bed = -np.ones((md.mesh.numberofvertices, 1))
     75# Masstransport
     76md.basalforcings.groundedice_melting_rate = np.zeros((md.mesh.numberofvertices, 1))
     77md.basalforcings.floatingice_melting_rate = np.zeros((md.mesh.numberofvertices, 1))
     78md.initialization.vx = np.zeros((md.mesh.numberofvertices, 1))
     79md.initialization.vy = np.zeros((md.mesh.numberofvertices, 1))
     80md.initialization.sealevel = np.zeros((md.mesh.numberofvertices, 1))
     81md.initialization.str = 0
    6282
    63 #materials
    64 md.materials=materials('hydro')
     83# Materials
     84md.materials = materials('hydro')
    6585
    66 #Miscellaneous
     86# Miscellaneous
    6787md.miscellaneous.name = 'test2002'
    6888
    69 #Solution parameters
     89# Solution parameters
     90md.cluster.np = 3
    7091md.solidearth.settings.reltol = np.nan
    7192md.solidearth.settings.abstol = 1e-3
    72 md.solidearth.settings.computesealevelchange = 1
     93md.solidearth.settings.sealevelloading = 1
     94md.solidearth.settings.isgrd = 1
     95md.solidearth.settings.ocean_area_scaling = 0
     96md.solidearth.settings.grdmodel = 1
    7397
    74 #max number of iterations reverted back to 10 (i.e., the original default value)
     98# Physics
     99md.transient.issmb = 0
     100md.transient.isstressbalance = 0
     101md.transient.isthermal = 0
     102md.transient.ismasstransport = 1
     103md.transient.isslc = 1
     104md.solidearth.requested_outputs = ['Sealevel', 'Bed']
     105
     106# Max number of iterations reverted back to 10 (i.e., the original default value)
    75107md.solidearth.settings.maxiter = 10
    76108
    77 #eustatic run:
    78 md.solidearth.settings.rigid = 0
     109# Eustatic run
     110md.solidearth.settings.selfattraction = 0
    79111md.solidearth.settings.elastic = 0
    80112md.solidearth.settings.rotation = 0
    81 md = solve(md, 'Sealevelrise')
    82 Seustatic = md.results.SealevelriseSolution.Sealevel
     113md.solidearth.settings.viscous = 0
    83114
    84 #eustatic + rigid run:
    85 md.solidearth.settings.rigid = 1
     115md = solve(md, 'Transient')
     116Seustatic = md.results.TransientSolution.Sealevel
     117Beustatic = md.results.TransientSolution.Bed
     118
     119# Eustatic + selfattraction run
     120md.solidearth.settings.selfattraction = 1
    86121md.solidearth.settings.elastic = 0
    87122md.solidearth.settings.rotation = 0
    88 md = solve(md, 'Sealevelrise')
    89 Srigid = md.results.SealevelriseSolution.Sealevel
     123md.solidearth.settings.viscous = 0
     124md = solve(md, 'tr')
     125Sselfattraction = md.results.TransientSolution.Sealevel
     126Bselfattraction = md.results.TransientSolution.Bed
    90127
    91 #eustatic + rigid + elastic run:
    92 md.solidearth.settings.rigid = 1
     128# Eustatic + selfattraction + elastic run
     129md.solidearth.settings.selfattraction = 1
    93130md.solidearth.settings.elastic = 1
    94131md.solidearth.settings.rotation = 0
    95 md = solve(md, 'Sealevelrise')
    96 Selastic = md.results.SealevelriseSolution.Sealevel
     132md.solidearth.settings.viscous = 0
     133md = solve(md, 'tr')
     134Selastic = md.results.TransientSolution.Sealevel
     135Belastic = md.results.TransientSolution.Bed
    97136
    98 #eustatic + rigid + elastic + rotation run:
    99 md.solidearth.settings.rigid = 1
     137# Eustatic + selfattraction + elastic + rotation run
     138md.solidearth.settings.selfattraction = 1
    100139md.solidearth.settings.elastic = 1
    101140md.solidearth.settings.rotation = 1
    102 md = solve(md, 'Sealevelrise')
    103 Srotation = md.results.SealevelriseSolution.Sealevel
     141md.solidearth.settings.viscous = 0
     142md = solve(md, 'tr')
     143Srotation = md.results.TransientSolution.Sealevel
     144Brotation = md.results.TransientSolution.Bed
    104145
    105 #Fields and tolerances to track changes
    106 field_names = ['Eustatic', 'Rigid', 'Elastic', 'Rotation']
    107 field_tolerances = [1e-13, 1e-13, 1e-13, 1e-13]
    108 field_values = [Seustatic, Srigid, Selastic, Srotation]
     146# Fields and tolerances to track changes
     147field_names = ['Seustatic', 'Sselfattraction', 'Selastic', 'Srotation', 'Beustatic', 'Bselfattraction', 'Belastic', 'Brotation']
     148field_tolerances = [1e-13, 1e-13, 1e-13, 1e-13, 1e-13, 1e-13, 1e-13, 1e-13]
     149field_values = [Seustatic, Sselfattraction, Selastic, Srotation, Beustatic, Bselfattraction, Belastic, Brotation]
Note: See TracBrowser for help on using the repository browser.