source: issm/oecreview/Archive/25834-26739/ISSM-26058-26059.diff@ 28275

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

CHG: added 25834-26739
File size: 119.2 KB

  • ../trunk-jpl/src/m/boundaryconditions/SetIceSheetBC.m

     
    3030
    3131%Initialize surface and basal forcings
    3232md.smb = initialize(md.smb,md);
    33 md.basalforcings   = initialize(md.basalforcings,md);
     33md.basalforcings = initialize(md.basalforcings,md);
    3434
    3535%Initialize ocean forcings and sealevel
    36 md.dsl= initialize(md.dsl,md);
     36md.dsl = initialize(md.dsl,md);
    3737
    3838%Deal with other boundary conditions
    3939if isnan(md.balancethickness.thickening_rate),

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

     
    22
    33
    44def SetIceSheetBC(md):
    5     """
    6     SETICESHEETBC - Create the boundary conditions for stressbalance and thermal models for an IceSheet with no Ice Front
     5    """SETICESHEETBC - Create the boundary conditions for stressbalance and thermal models for an IceSheet with no Ice Front
    76
    8        Usage:
    9           md = SetIceSheetBC(md)
     7    Usage:
     8        md = SetIceSheetBC(md)
    109
    11        See also: SETICESHELFBC, SETMARINEICESHEETBC
     10    See also: SETICESHELFBC, SETMARINEICESHEETBC
    1211    """
    1312
    1413    #node on Dirichlet
    … …  
    3231
    3332    #No ice front -> do nothing
    3433
    35     #Create zeros basalforcings and smb
     34    #Initialize surface and basal forcings
    3635    md.smb.initialize(md)
    3736    md.basalforcings.initialize(md)
    3837
     38    #Initialize ocean forcings and sealevel
     39    md.dsl.initialize(md)
     40
    3941    #Deal with other boundary conditions
    4042    if np.all(np.isnan(md.balancethickness.thickening_rate)):
    4143        md.balancethickness.thickening_rate = np.zeros((md.mesh.numberofvertices))

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

     
    66
    77
    88class hydrologydc(object):
    9     """
    10     Hydrologydc class definition
     9    """HYDROLOGYDC class definition
    1110
    1211    Usage:
    13             hydrologydc = hydrologydc()
     12        hydrologydc = hydrologydc()
    1413    """
    1514
    1615    def __init__(self):  # {{{
    … …  
    4847        self.epl_conductivity = 0
    4948        self.eplflip_lock = 0
    5049
    51     #set defaults
    5250        self.setdefaultparameters()
    5351    # }}}
    5452
    5553    def __repr__(self):  # {{{
     54        # TODO:
     55        # - Convert all formatting to calls to <string>.format (see any
     56        #   already converted <class>.__repr__ method for examples)
     57        #
    5658        string = '   hydrology Dual Porous Continuum Equivalent parameters:'
    5759        string = ' - general parameters'
    5860        string = "%s\n%s" % (string, fielddisplay(self, 'water_compressibility', 'compressibility of water [Pa^ - 1]'))

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

     
    2323                end % }}}
    2424                function list = defaultoutputs(self,md) % {{{
    2525                        list = {''};
    26                 end % }}}   
     26                end % }}}
    2727
    2828                function self = setdefaultparameters(self) % {{{
    2929                        self.requested_outputs = {'default'};
    … …  
    4545                function marshall(self,prefix,md,fid) % {{{
    4646                        WriteData(fid,prefix,'name','md.hydrology.model','data',2,'format','Integer');
    4747                        WriteData(fid,prefix,'object',self,'fieldname','spcwatercolumn','format','DoubleMat','mattype',1,'timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts);
    48       outputs = self.requested_outputs;
    49       pos  = find(ismember(outputs,'default'));
    50       if ~isempty(pos),
    51         outputs(pos) = [];  %remove 'default' from outputs
    52         outputs      = [outputs defaultoutputs(self,md)]; %add defaults
    53       end
    54       WriteData(fid,prefix,'data',outputs,'name','md.hydrology.requested_outputs','format','StringArray');
     48                        outputs = self.requested_outputs;
     49                        pos  = find(ismember(outputs,'default'));
     50                        if ~isempty(pos),
     51                                outputs(pos) = [];  %remove 'default' from outputs
     52                                outputs      = [outputs defaultoutputs(self,md)]; %add defaults
     53                        end
     54                        WriteData(fid,prefix,'data',outputs,'name','md.hydrology.requested_outputs','format','StringArray');
    5555                end % }}}
    5656                function savemodeljs(self,fid,modelname) % {{{
    57                
     57
    5858                        writejs1Darray(fid,[modelname '.hydrology.spcwatercolumn'],self.spcwatercolumn);
    5959
    6060                end % }}}

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

     
     1import numpy as np
     2
     3from structtoobj import structtoobj
     4
     5class hydrologytws(object):
     6    """HYDROLOGYTWS class definition
     7
     8    Usage:
     9        hydrologytws = hydrologytws()
     10    """
     11
     12    def __init__(self):  # {{{
     13        self.spcwatercolumn = np.nan
     14        self.requested_outputs = np.nan
     15
     16        nargs = len(args)
     17        if nargs == 0:
     18            self.setdefaultparameters()
     19        elif nargs == 1:
     20            self = structtoobj(self, args[0])
     21        else:
     22            raise RuntimeError('constructor not supported')
     23    #}}}
     24
     25    def __repr__(self):  # {{{
     26        s = '   hydrologytws solution parameters:\n'
     27        s += '{}\n'.format(fielddisplay(self, 'spcwatercolumn', 'water thickness constraints (NaN means no constraint) [m]'))
     28        s += '{}\n'.format(fielddisplay(self, 'requested_outputs', 'additional outputs requested'))
     29        return s
     30    #}}}
     31
     32    def defaultoutputs(self, md): # {{{
     33        return ['']
     34    #}}}
     35
     36    def setdefaultparameters(self):  # {{{
     37        self.requested_outputs = ['defualt']
     38        return self
     39    #}}}
     40
     41    def extrude(self, md):  # {{{
     42        return self
     43    #}}}
     44
     45    def checkconsistency(self, md, solution, analyses):  # {{{
     46        # Early return
     47        if 'HydrologyTwsAnalysis' not in analyses:
     48            return
     49        md = checkfield(md, 'fieldname', 'hydrology.spcwatercolumn', 'Inf', 1, 'timeseries', 1)
     50    #}}}
     51
     52    def marshall(self, prefix, md, fid):  # {{{
     53        WriteData(fid, prefix, 'name', 'md.hydrology.model', 'data', 2, 'format', 'Integer')
     54        WriteData(fid, prefix, 'object', self, 'fieldname', 'spcwatercolumn', 'format', 'DoubleMat', 'mattype', 1, 'timeserieslength', md.mesh.numberofvertices + 1, 'yts', md.constants.yts)
     55        outputs = self.requested_outputs
     56        pos  = find(ismember(outputs,'default'))
     57        if not len(pos),
     58            outputs[pos] = [];  # remove 'default' from outputs
     59            outputs.extend(defaultoutputs(self, md)) # add defaults
     60        end
     61        WriteData(fid, prefix, 'data', outputs, 'name', 'md.hydrology.requested_outputs', 'format', 'StringArray')
     62    # }}}
     63
     64

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

     
    55
    66classdef matenhancedice
    77        properties (SetAccess=public)
    8                 rho_ice                    = 0.;
    9                 rho_water                  = 0.;
    10                 rho_freshwater             = 0.;
    11                 mu_water                   = 0.;
    12                 heatcapacity               = 0.;
    13                 latentheat                 = 0.;
    14                 thermalconductivity        = 0.;
    15                 temperateiceconductivity   = 0.;
     8                rho_ice                         = 0.;
     9                rho_water                       = 0.;
     10                rho_freshwater                  = 0.;
     11                mu_water                        = 0.;
     12                heatcapacity                    = 0.;
     13                latentheat                      = 0.;
     14                thermalconductivity             = 0.;
     15                temperateiceconductivity        = 0.;
    1616                effectiveconductivity_averaging = 0.;
    17                 meltingpoint               = 0.;
    18                 beta                       = 0.;
    19                 mixed_layer_capacity       = 0.;
    20                 thermal_exchange_velocity  = 0.;
    21                 rheology_E   = NaN;
    22                 rheology_B   = NaN;
    23                 rheology_n   = NaN;
    24                 rheology_law = '';
     17                meltingpoint                    = 0.;
     18                beta                            = 0.;
     19                mixed_layer_capacity            = 0.;
     20                thermal_exchange_velocity       = 0.;
     21                rheology_E                      = NaN;
     22                rheology_B                      = NaN;
     23                rheology_n                      = NaN;
     24                rheology_law                    = '';
    2525
    26                 %slr
    27                 earth_density              = 0;
     26                %SLC
     27                earth_density                   = 0;
    2828
    2929        end
    3030        methods

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

     
    55
    66classdef matestar
    77        properties (SetAccess=public)
    8                 rho_ice                    = 0.;
    9                 rho_water                  = 0.;
    10                 rho_freshwater             = 0.;
    11                 mu_water                   = 0.;
    12                 heatcapacity               = 0.;
    13                 latentheat                 = 0.;
    14                 thermalconductivity        = 0.;
    15                 temperateiceconductivity   = 0.;
     8                rho_ice                         = 0.;
     9                rho_water                       = 0.;
     10                rho_freshwater                  = 0.;
     11                mu_water                        = 0.;
     12                heatcapacity                    = 0.;
     13                latentheat                      = 0.;
     14                thermalconductivity             = 0.;
     15                temperateiceconductivity        = 0.;
    1616                effectiveconductivity_averaging = 0;
    17                 meltingpoint               = 0.;
    18                 beta                       = 0.;
    19                 mixed_layer_capacity       = 0.;
    20                 thermal_exchange_velocity  = 0.;
    21                 rheology_B    = NaN;
    22                 rheology_Ec   = NaN;
    23                 rheology_Es   = NaN;
    24                 rheology_law = '';
     17                meltingpoint                    = 0.;
     18                beta                            = 0.;
     19                mixed_layer_capacity            = 0.;
     20                thermal_exchange_velocity       = 0.;
     21                rheology_B                      = NaN;
     22                rheology_Ec                     = NaN;
     23                rheology_Es                     = NaN;
     24                rheology_law                    = '';
    2525
    2626                %slc
    27                 earth_density              = 0;
     27                earth_density                   = 0;
    2828
    2929        end
    3030        methods

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

     
    77
    88
    99class matice(object):
    10     '''
    11     MATICE class definition
     10    """MATICE class definition
    1211
    1312    Usage:
    1413            matice = matice()
    15     '''
     14    """
    1615
    1716    def __init__(self): #{{{
    1817        self.rho_ice = 0.
    … …  
    3231        self.rheology_n = np.nan
    3332        self.rheology_law = ''
    3433
    35         #giaivins
    36         self.lithosphere_shear_modulus = 0.
    37         self.lithosphere_density = 0.
    38         self.mantle_shear_modulus = 0.
    39         self.mantle_density = 0.
    40 
    41         #slr
     34        #slc
    4235        self.earth_density = 0
    4336        self.setdefaultparameters()
    4437    #}}}
    … …  
    6154        s = "%s\n%s" % (s, fielddisplay(self, "rheology_B", "flow law parameter [Pa s^(1/n)]"))
    6255        s = "%s\n%s" % (s, fielddisplay(self, "rheology_n", "Glen's flow law exponent"))
    6356        s = "%s\n%s" % (s, fielddisplay(self, "rheology_law", "law for the temperature dependance of the rheology: 'None', 'BuddJacka', 'Cuffey', 'CuffeyTemperate', 'Paterson', 'Arrhenius', 'LliboutryDuval', 'NyeCO2', or 'NyeH2O'"))
    64         s = "%s\n%s" % (s, fielddisplay(self, "lithosphere_shear_modulus", "Lithosphere shear modulus [Pa]"))
    65         s = "%s\n%s" % (s, fielddisplay(self, "lithosphere_density", "Lithosphere density [g/cm^-3]"))
    66         s = "%s\n%s" % (s, fielddisplay(self, "mantle_shear_modulus", "Mantle shear modulus [Pa]"))
    67         s = "%s\n%s" % (s, fielddisplay(self, "mantle_density", "Mantle density [g/cm^-3]"))
    6857        s = "%s\n%s" % (s, fielddisplay(self, "earth_density", "Mantle density [kg/m^-3]"))
    6958
    7059        return s
    … …  
    10796        #available: none, paterson and arrhenius
    10897        self.rheology_law = 'Paterson'
    10998
    110         # GIA:
    111         self.lithosphere_shear_modulus = 6.7e10 # (Pa)
    112         self.lithosphere_density = 3.32  # (g/cm^-3)
    113         self.mantle_shear_modulus = 1.45e11  # (Pa)
    114         self.mantle_density = 3.34  # (g/cm^-3)
     99        # Rheology for ice
     100        self.rheology_B = 2.1 * 1e8
     101        self.rheology_n = 3
    115102
    116103        # SLR
    117104        self.earth_density = 5512  # average density of the Earth, (kg/m^3)
    … …  
    118105    #}}}
    119106
    120107    def checkconsistency(self, md, solution, analyses): #{{{
    121         if solution != 'SealevelriseSolution':
     108        if solution == 'TransientSolution' and md.transient.isslc:
     109            md = checkfield(md, 'fieldname', 'materials.earth_density', '>', 0, 'numel', [1])
     110        else:
    122111            md = checkfield(md, 'fieldname', 'materials.rho_ice', '>', 0)
    123112            md = checkfield(md, 'fieldname', 'materials.rho_water', '>', 0)
    124113            md = checkfield(md, 'fieldname', 'materials.rho_freshwater', '>', 0)
    125114            md = checkfield(md, 'fieldname', 'materials.mu_water', '>', 0)
    126115            md = checkfield(md, 'fieldname', 'materials.rheology_B', '>', 0, 'universal', 1, 'NaN', 1, 'Inf', 1)
    127             md = checkfield(md, 'fieldname', 'materials.rheology_n', '>', 0, 'size', [md.mesh.numberofelements])
     116            md = checkfield(md, 'fieldname', 'materials.rheology_n', '>', 0, 'universal',1, 'NaN', 1, 'Inf', 1)
    128117            md = checkfield(md, 'fieldname', 'materials.rheology_law', 'values', ['None', 'BuddJacka', 'Cuffey', 'CuffeyTemperate', 'Paterson', 'Arrhenius', 'LliboutryDuval', 'NyeCO2', 'NyeH2O'])
    129             md = checkfield(md, 'fieldname', 'materials.effectiveconductivity_averaging', 'numel', [1], 'values', [0, 1, 2])
     118            md = checkfield(md,'fieldname','materials.effectiveconductivity_averaging', 'numel', [1], 'values', [0, 1, 2])
    130119
    131         if 'GiaAnalysis' in analyses:
    132             md = checkfield(md, 'fieldname', 'materials.lithosphere_shear_modulus', '>', 0, 'numel', [1])
    133             md = checkfield(md, 'fieldname', 'materials.lithosphere_density', '>', 0, 'numel', [1])
    134             md = checkfield(md, 'fieldname', 'materials.mantle_shear_modulus', '>', 0, 'numel', [1])
    135             md = checkfield(md, 'fieldname', 'materials.mantle_density', '>', 0, 'numel', [1])
    136 
    137         if 'SealevelriseAnalysis' in analyses:
    138             md = checkfield(md, 'fieldname', 'materials.earth_density', '>', 0, 'numel', [1])
    139 
    140120        return md
    141121    #}}}
    142122
    … …  
    158138        WriteData(fid, prefix, 'object', self, 'class', 'materials', 'fieldname', 'rheology_B', 'format', 'DoubleMat', 'mattype', 1, 'timeserieslength', md.mesh.numberofvertices + 1, 'yts', md.constants.yts)
    159139        WriteData(fid, prefix, 'object', self, 'class', 'materials', 'fieldname', 'rheology_n', 'format', 'DoubleMat', 'mattype', 2)
    160140        WriteData(fid, prefix, 'data', self.rheology_law, 'name', 'md.materials.rheology_law', 'format', 'String')
    161         WriteData(fid, prefix, 'object', self, 'class', 'materials', 'fieldname', 'lithosphere_shear_modulus', 'format', 'Double')
    162         WriteData(fid, prefix, 'object', self, 'class', 'materials', 'fieldname', 'lithosphere_density', 'format', 'Double', 'scale', 10.**3.)
    163         WriteData(fid, prefix, 'object', self, 'class', 'materials', 'fieldname', 'mantle_shear_modulus', 'format', 'Double')
    164         WriteData(fid, prefix, 'object', self, 'class', 'materials', 'fieldname', 'mantle_density', 'format', 'Double', 'scale', 10.**3.)
    165141        WriteData(fid, prefix, 'object', self, 'class', 'materials', 'fieldname', 'earth_density', 'format', 'Double')
    166142    #}}}

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

     
    1919                        end
    2020                end % }}}
    2121                function self = setdefaultparameters(self) % {{{
    22                
     22
    2323                        icethicknesschange=[];
    2424                        waterheightchange=[];
    2525                        otherchange=[];
    26                
     26
    2727                end % }}}
    28                         function md = checkconsistency(self,md,solution,analyses) % {{{
     28                function md = checkconsistency(self,md,solution,analyses) % {{{
    2929
    3030                        if ~ismember('SealevelchangeAnalysis',analyses) | (strcmp(solution,'TransientSolution') & md.transient.isslc==0),
    3131                                return;

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

     
    77        properties (SetAccess=public)
    88                issmb             = 0;
    99                ismasstransport   = 0;
    10                 isoceantransport   = 0;
     10                isoceantransport  = 0;
    1111                isstressbalance   = 0;
    1212                isthermal         = 0;
    1313                isgroundingline   = 0;
    … …  
    1515                isdamageevolution = 0;
    1616                ismovingfront     = 0;
    1717                ishydrology       = 0;
    18         issampling        = 0;
     18                issampling        = 0;
    1919                isslc             = 0;
    2020                amr_frequency     = 0;
    2121                isoceancoupling   = 0;
    … …  
    4343                        self.isdamageevolution = 0;
    4444                        self.ismovingfront   =0;
    4545                        self.ishydrology     = 0;
    46             self.issampling      = 0;
     46                        self.issampling      = 0;
    4747                        self.isslc           = 0;
    4848                        self.isoceancoupling = 0;
    4949                        self.amr_frequency      = 0;
    … …  
    6464                        self.isdamageevolution = 0;
    6565                        self.ismovingfront   = 0;
    6666                        self.ishydrology     = 0;
    67             self.issampling      = 0;
     67                        self.issampling      = 0;
    6868                        self.isslc           = 0;
    6969                        self.isoceancoupling = 0;
    7070                        self.amr_frequency      = 0;
    … …  
    9797                        md = checkfield(md,'fieldname','transient.requested_outputs','stringrow',1);
    9898                        md = checkfield(md,'fieldname','transient.isslc','numel',[1],'values',[0 1]);
    9999                        md = checkfield(md,'fieldname','transient.isoceancoupling','numel',[1],'values',[0 1]);
    100             md = checkfield(md,'fieldname','transient.issampling','numel',[1],'values',[0 1]); 
     100                        md = checkfield(md,'fieldname','transient.issampling','numel',[1],'values',[0 1]); 
    101101                        md = checkfield(md,'fieldname','transient.amr_frequency','numel',[1],'>=',0,'NaN',1,'Inf',1);
    102102
    103103                        if (~strcmp(solution,'TransientSolution') & md.transient.iscoupling==1),
    … …  
    120120                        fielddisplay(self,'isdamageevolution','indicates whether damage evolution is used in the transient');
    121121                        fielddisplay(self,'ismovingfront','indicates whether a moving front capability is used in the transient');
    122122                        fielddisplay(self,'ishydrology','indicates whether an hydrology model is used');
    123             fielddisplay(self,'issampling','indicates whether sampling is used in the transient')
     123                        fielddisplay(self,'issampling','indicates whether sampling is used in the transient')
    124124                        fielddisplay(self,'isslc','indicates whether a sea-level change solution is used in the transient');
    125125                        fielddisplay(self,'isoceancoupling','indicates whether a coupling with an ocean model is used in the transient');
    126126                        fielddisplay(self,'amr_frequency','frequency at which mesh is refined in simulations with multiple time_steps');
    … …  
    138138                        WriteData(fid,prefix,'object',self,'fieldname','isdamageevolution','format','Boolean');
    139139                        WriteData(fid,prefix,'object',self,'fieldname','ishydrology','format','Boolean');
    140140                        WriteData(fid,prefix,'object',self,'fieldname','ismovingfront','format','Boolean');
    141             WriteData(fid,prefix,'object',self,'fieldname','issampling','format','Boolean');
     141                        WriteData(fid,prefix,'object',self,'fieldname','issampling','format','Boolean');
    142142                        WriteData(fid,prefix,'object',self,'fieldname','isslc','format','Boolean');
    143143                        WriteData(fid,prefix,'object',self,'fieldname','isoceancoupling','format','Boolean');
    144144                        WriteData(fid,prefix,'object',self,'fieldname','amr_frequency','format','Integer');
    … …  
    164164                        writejsdouble(fid,[modelname '.trans.isdamageevolution'],self.isdamageevolution);
    165165                        writejsdouble(fid,[modelname '.trans.ismovingfront'],self.ismovingfront);
    166166                        writejsdouble(fid,[modelname '.trans.ishydrology'],self.ishydrology);
    167             writejsdouble(fid,[modelname '.trans.issampling'],self.issampling);
     167                        writejsdouble(fid,[modelname '.trans.issampling'],self.issampling);
    168168                        writejsdouble(fid,[modelname '.trans.isslc'],self.isslc);
    169169                        writejsdouble(fid,[modelname '.trans.isoceancoupling'],self.isoceancoupling);
    170170                        writejsdouble(fid,[modelname '.trans.amr_frequency'],self.amr_frequency);

  • ../trunk-jpl/src/m/consistency/ismodelselfconsistent.m

     
    22%ISMODELSELFCONSISTENT - check that model forms a closed form solvable problem.
    33%
    44%   Usage:
    5 %      ismodelselfconsistent(md),
     5%      ismodelselfconsistent(md);
    66
    77%initialize consistency as true
    88md.private.isconsistent=true;

  • ../trunk-jpl/src/m/plot/radarpower.py

     
     1def solveslm(slm, solutionstringi, *args):
     2    """RADARPOWER - overlay a power radar image on an existing mesh
     3
     4    This routine will overlay a power radar image on an existing mesh.
     5    The power amplitude will be output to vel for now.
     6    In the future, think about a field to hold this value.
     7
     8    Usage:
     9        md=radarpower(md,options);
     10        md=radarpower(md)
     11
     12    TODO:
     13    - Translate from MATLAB API as we bring Python plotting capabilities online
     14    """
     15
     16    return md

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

     
    11function slm=solveslm(slm,solutionstringi,varargin)
    2 %SOLVESLR - apply solution sequence for this sealevel model
     2%SOLVESLM - apply solution sequence for this sealevel model
    33%
    44%   Usage:
    55%      slm=solve(slm,solutionstring,varargin)
    … …  
    4242slm.private.solution=solutionstring;
    4343cluster=slm.cluster;
    4444batch=0;
    45 %now, go through icecaps, glacies and earth, and upload all the data independently:
     45%now, go through icecaps, glaciers and earth, and upload all the data independently:
    4646disp('solving ice caps first');
    4747for i=1:length(slm.icecaps),
    4848        slm.icecaps{i}=solve(slm.icecaps{i},solutionstringi,'batch','yes');
    … …  
    5050disp('solving earth now');
    5151slm.earth=solve(slm.earth,solutionstringi,'batch','yes');
    5252
    53 %Firs, build a runtime name that is unique
     53%First, build a runtime name that is unique
    5454c=clock;
    5555slm.private.runtimename=sprintf('%s-%02i-%02i-%04i-%02i-%02i-%02i-%i',slm.miscellaneous.name,c(2),c(3),c(1),c(4),c(5),floor(c(6)),feature('GetPid'));
    5656

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

     
     1from datetime import datetime
     2import os
     3
     4import numpy as np
     5
     6from loadresultsfromcluster import loadresultsfromcluster
     7from pairoptions import pairoptions
     8from waitonlock import waitonlock
     9
     10
     11def solveslm(slm, solutionstringi, *args):
     12    """SOLVESLM - apply solution sequence for this sealevel model
     13
     14    Usage:
     15        slm=solve(slm,solutionstring,varargin)
     16        where varargin is a lit of paired arguments of string OR enums
     17
     18    solution types available comprise:
     19        - 'Transient'
     20
     21    extra options:
     22
     23    Examples:
     24        slm=solve(slm,'Transient');
     25    """
     26
     27    # Recover and process solve options
     28    if solutionstringi.lower() == 'tr' or solutionstringi.lower() == 'transient':
     29        solutionstring = 'TransientSolution'
     30    else:
     31        raise RuntimeError('solutionstring {} not supported!'.format(solutionstringi))
     32
     33    # Default settings for debugging
     34    valgrind = 0
     35    #slm.cluster.interactive = 0
     36    #valgrind = 1
     37
     38    # Check consistency
     39    slm.checkconsistency(solutionstring)
     40
     41    # Process options
     42    options = pairoptions('solutionstring', solutionstring, *args)
     43
     44    # Make sure we request sum of cluster processors
     45    totalnp = 0
     46    for i in range(len(slm.icecaps)):
     47        totalnp = totalnp + slm.icecaps[i].cluster.np
     48    totalnp = totalnp + slm.earth.cluster.np
     49    if totalnp != slm.cluster.np:
     50        raise RuntimeError('sum of all icecaps and earch cluster processors requestes should be equal to slm.cluster.np')
     51
     52    # Recover some fields
     53    slm.private.solution = solutionstring
     54    cluster = slm.cluster
     55    batch = 0
     56    # Now, go through icecaps, glaciers and earth, and upload all the data independently
     57    print('solving ice caps first')
     58    for i in range(len(slm.icecaps)):
     59        slm.icecaps[i] = solve(slm.icecaps[i], solutionastringi,'batch','yes')
     60    print('solving earth now')
     61    slm.earth = solve(slm.earth, solutionstringi, 'batch', 'yes')
     62
     63    # First, build a runtime name that is unique
     64    c = datetime.now()
     65    md.private.runtimename = "%s-%02i-%02i-%04i-%02i-%02i-%02i-%i" % (md.miscellaneous.name, c.month, c.day, c.year, c.hour, c.minute, c.second, os.getpid())
     66
     67    # Write all input files
     68    privateruntimenames = []
     69    miscellaneousnames = []
     70    nps = []
     71    for i in range(len(slm.icecaps)):
     72        privateruntimenames.append(slm.icecaps[i],private.runtimename)
     73        miscellaneousnames.append(slm.earth.miscellaneous.name)
     74        nps.append(slm.earth.cluster.np)
     75
     76    BuildQueueScriptMultipleModels(cluster, slm.private.runtimename, slm.miscellaneous.name, slm.private.solution, valgrind, privateruntimenames, miscellaneousnames, nps)
     77
     78    # Upload all required files, given that each individual solution for icecaps and earth model already did
     79    filelist = [slm.miscellaneous.name + '.queue']
     80    UploadQueueJob(cluster, slm.miscellaneous.name, slm.private.runtimename, filelist)
     81
     82    # Launch queue job
     83    LaunchQueueJob(cluster, slm.miscellaneous.name, slm.private.runtimename, filelist, '', batch)
     84
     85    # Wait on lock
     86    if slm.settings.waitonlock > 0:
     87        islock = waitonlock(slm)
     88        if islock == 0:  # no results to be loaded
     89            print('The results must be loaded manually with md = loadresultsfromcluster(md).')
     90        else: # load results
     91            if slm.verbose.solution:
     92                print('loading results from cluster')
     93            slm = loadresultsfromcluster(slm)
     94
     95    return slm

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

     
    66classdef dsl
    77        properties (SetAccess=public)
    88
    9                 global_average_thermosteric_sea_level; %corresponds to zostoga field in CMIP5 archives. Specified as a temporally variable quantity (in m)
    10                 sea_surface_height_above_geoid; %corresponds to zos field in CMIP5 archives. Spatial average is 0. Specified as a spatio-temporally variable (in m)
    11                 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!)
     9                global_average_thermosteric_sea_level; %Corresponds to zostoga field in CMIP5 archives. Specified as a temporally variable quantity (in m).
     10                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).
     11                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!).
    1212
    1313        end
    1414        methods
    … …  
    5050                function disp(self) % {{{
    5151
    5252                        disp(sprintf('   dsl parameters:'));
    53                         fielddisplay(self,'global_average_thermosteric_sea_level','corresponds to zostoga field in CMIP5 archives. Specified as a temporally variable quantity (in m)');
    54                         fielddisplay(self,'sea_surface_height_above_geoid','corresponds to zos field in CMIP5 archives. Spatial average is 0. Specified as a spatio-temporally quantity (in m)');
    55                         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)');
     53                        fielddisplay(self,'global_average_thermosteric_sea_level','Corresponds to zostoga field in CMIP5 archives. Specified as a temporally variable quantity (in m).');
     54                        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).');
     55                        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!).');
    5656
    5757                end % }}}
    5858                function marshall(self,prefix,md,fid) % {{{

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

     
    1414    """
    1515
    1616    def __init__(self): #{{{
    17         self.global_average_thermosteric_sea_level_change = 0 #corresponds to zostoga field in CMIP5 archives. Specified as a temporally variable global rate (mm/yr)
    18         self.sea_surface_height_change_above_geoid = float('NaN') #corresponds to zos field in CMIP5 archives. Spatial average is 0. Specified as a spatio-temporally variable rate (mm/yr)
    19         self.sea_water_pressure_change_at_sea_floor = float('NaN') #corresponds to bpo field in CMIP5 archives. Specified as a spatio-temporally variable rate (in Pa/yr)
    20         self.compute_fingerprints = 0; #do we use the sea water pressure change to compute fingerprints and correct sea_surface_height_change_above_geoid
     17        self.global_average_thermosteric_sea_level = np.nan # Corresponds to zostoga field in CMIP5 archives. Specified as a temporally variable quantity (in m).
     18        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).
     19        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!).
    2120    #}}}
    2221
    2322    def __repr__(self): #{{{
    2423        s = '   dsl parameters:\n'
    25         s += '{}\n'.format(fielddisplay(self, 'global_average_thermosteric_sea_level_change', 'corresponds to zostoga field in CMIP5 archives. Specified as a temporally variable global rate (mm/yr)'))
    26         s += '{}\n'.format(fielddisplay(self, 'sea_surface_height_change_above_geoid', 'corresponds to zos field in CMIP5 archives. Spatial average is 0. Specified as a spatio-temporally variable rate (mm/yr)'))
    27         s += '{}\n'.format(fielddisplay(self, 'sea_water_pressure_change_at_sea_floor', 'corresponds to bpo field in CMIP5 archives. Specified as a spatio-temporally variable rate (in Pa/yr)'))
    28         s += '{}\n'.format(fielddisplay(self, 'compute_fingerprints', 'do we use the sea water pressure change to compute fingerprints and correct sea_surface_height_change_above_geoid'))
     24        s += '{}\n'.format(fielddisplay(self, 'global_average_thermosteric_sea_level', 'Corresponds to zostoga field in CMIP5 archives. Specified as a temporally variable quantity (in m).'))
     25        s += '{}\n'.format(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).'))
     26        s += '{}\n'.format(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!).'))
    2927        return s
    3028    #}}}
    3129
    3230    def extrude(self, md): #{{{
    33         self.sea_surface_height_change_above_geoid = project3d(md, 'vector', self.sea_surface_height_change_above_geoid, 'type', 'node')
    34         self.sea_water_pressure_change_at_sea_floor = project3d(md, 'vector', self.sea_water_pressure_change_at_sea_floor, 'type', 'node')
     31        self.sea_surface_height_above_geoid = project3d(md, 'vector', self.sea_surface_height_above_geoid, 'type', 'node')
     32        self.sea_water_pressure_at_sea_floor = project3d(md, 'vector', self.sea_water_pressure_at_sea_floor, 'type', 'node')
    3533        return self
    3634    #}}}
    3735
    … …  
    4139
    4240    def checkconsistency(self, md, solution, analyses): #{{{
    4341        # Early return
    44         if 'SealevelriseAnalysis' not in analyses:
     42        if ('SealevelriseAnalysis' not in analyses) or (solution == 'TransientSolution' and not md.transient.isslc):
    4543            return md
    46         if solution == 'TransientSolution' and md.transient.isslc == 0:
    47             return md
    48         md = checkfield(md, 'fieldname', 'dsl.global_average_thermosteric_sea_level_change', 'NaN', 1, 'Inf', 1)
    49         md = checkfield(md, 'fieldname', 'dsl.sea_surface_height_change_above_geoid', 'NaN', 1, 'Inf', 1, 'timeseries', 1)
    50         md = checkfield(md, 'fieldname', 'dsl.sea_water_pressure_change_at_sea_floor', 'NaN', 1, 'Inf', 1, 'timeseries', 1)
    51         md = checkfield(md, 'fieldname', 'dsl.compute_fingerprints', 'NaN', 1, 'Inf', 1, 'values', [0, 1])
    52         if self.compute_fingerprints:
    53             # Check if geodetic flag of slr is on
    54             if not md.slr.geodetic:
    55                 raise RuntimeError('DSL checkconsistency error message: if bottom pressure fingerprints computations are requested, slr class should have geodetic flag on')
     44
     45        md = checkfield(md, 'fieldname', 'dsl.global_average_thermosteric_sea_level', 'NaN', 1, 'Inf', 1)
     46        md = checkfield(md, 'fieldname', 'dsl.sea_surface_height_above_geoid', 'NaN', 1, 'Inf', 1, 'timeseries', 1)
     47        md = checkfield(md, 'fieldname', 'dsl.sea_water_pressure_at_sea_floor', 'NaN', 1, 'Inf', 1, 'timeseries', 1)
     48
     49        if md.solidearth.settings.compute_bp_grd:
     50            md = checkfield(md, 'fieldname', dsl.sea_water_pressure_at_sea_floor, 'empty', 1)
     51
    5652        return md
    5753    # }}}
    5854
    … …  
    5955    def marshall(self, prefix, md, fid):   #{{{
    6056        yts = md.constants.yts
    6157        WriteData(fid, prefix, 'name', 'md.dsl.model', 'data', 1, 'format', 'Integer')
    62         WriteData(fid, prefix, 'object', self, 'class', 'dsl', 'fieldname', 'compute_fingerprints', 'format', 'Integer')
    63         WriteData(fid, prefix, 'object', self, 'class', 'dsl', 'fieldname', 'global_average_thermosteric_sea_level_change', 'format', 'DoubleMat', 'mattype', 1, 'timeserieslength', 1+1, 'yts', md.constants.yts, 'scale', 1e-3/md.constants.yts)
    64         WriteData(fid, prefix, 'object', self, 'class', 'dsl', 'fieldname', 'sea_water_pressure_change_at_sea_floor', 'format', 'DoubleMat', 'mattype', 1, 'timeserieslength', md.mesh.numberofvertices+1, 'yts', md.constants.yts, 'scale', 1e-3/md.constants.yts)
    65         WriteData(fid, prefix, 'object', self, 'class', 'dsl', 'fieldname', 'sea_surface_height_change_above_geoid', 'format', 'DoubleMat', 'mattype', 1, 'timeserieslength', md.mesh.numberofvertices+1, 'yts', md.constants.yts)
     58        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.
     59        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.
     60        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.
    6661    # }}}

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

     
    77        properties (SetAccess=public)
    88
    99                modelid; %index into the multi-model ensemble, determine which field will be used.
    10                 global_average_thermosteric_sea_level; %corresponds to zostoga fields in CMIP5 archives. Specified as a temporally quantity (in m) for each ensemble.
    11                 sea_surface_height_above_geoid; %corresponds to zos fields in CMIP5 archives. Spatial average is 0. Specified as a spatio-temporally quantity (in m) for each ensemble
    12                 sea_water_pressure_at_sea_floor; %corresponds to bpo field in CMIP5 archives. Specified as a spatio-temporally variable quantity (in m equivalent, not in Pa!)  for each ensemble
     10                global_average_thermosteric_sea_level; %Corresponds to zostoga field in CMIP5 archives. Specified as a temporally variable quantity (in m) for each ensemble.
     11                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) for each ensemble.
     12                sea_water_pressure_at_sea_floor; %Corresponds to bpo field in CMIP5 archives. Specified as a spatio-temporally variable quantity (in m equivalent, not in Pa!) for each ensemble.
    1313
    1414        end
    1515        methods
    … …  
    5252
    5353                        disp(sprintf('   dsl mme parameters:'));
    5454                        fielddisplay(self,'modelid','index into the multi-model ensemble, determine which field will be used.');
    55                         fielddisplay(self,'global_average_thermosteric_sea_level','corresponds to zostoga fields in CMIP5 archives. Specified as a temporally variable quantity (in m) for each ensemble.');
    56                         fielddisplay(self,'sea_surface_height_above_geoid','corresponds to zos fields in CMIP5 archives. Spatial average is 0. Specified as a spatio-temporally quantity (in m) for each ensemble.');
    57                         fielddisplay(self,'sea_water_pressure_at_sea_floor','corresponds to bpo fields in CMIP5 archives. Specified as a spatio-temporally quantity (in m) for each ensemble.');
     55                        fielddisplay(self,'global_average_thermosteric_sea_level','Corresponds to zostoga field in CMIP5 archives. Specified as a temporally variable quantity (in m) for each ensemble.');
     56                        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) for each ensemble.');
     57                        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!) for each ensemble.');
    5858                end % }}}
    5959                function marshall(self,prefix,md,fid) % {{{
    6060

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

     
    1414    """
    1515
    1616    def __init__(self, *args): #{{{
    17         self.modelid                                        = 0 #index into the multi-model ensemble
    18         self.global_average_thermosteric_sea_level_change   = [] #corresponds to zostoga fields in CMIP5 archives. Specified as a temporally variable global rate (mm/yr) for each ensemble.
    19         self.sea_surface_height_change_above_geoid          = [] #corresponds to zos fields in CMIP5 archives. Spatial average is 0. Specified as a spatio-temporally variable rate (mm/yr) for each ensemble
    20         self.sea_water_pressure_change_at_sea_floor         = [] #corresponds to bpo fields in CMIP5 archives. Specified as a spatio-temporally variable rate (in mm/yr equivalent, not in Pa/yr!) for each ensemble
    21         self.compute_fingerprints                           = 0 #corresponds to zos fields in CMIP5 archives. Spatial average is 0. Specified as a spatio-temporally variable rate (mm/yr) for each ensemble
    22        
    23         nargin = len(args)
     17        self.modelid                               = 0 # Index into the multi-model ensemble
     18        self.global_average_thermosteric_sea_level = [] # Corresponds to zostoga field in CMIP5 archives. Specified as a temporally variable quantity (in m) for each ensemble.
     19        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) for each ensemble.
     20        self.sea_water_pressure_at_sea_floor       = [] #Corresponds to bpo field in CMIP5 archives. Specified as a spatio-temporally variable quantity (in m equivalent, not in Pa!) for each ensemble.
    2421
    25         if nargin == 0:
     22        nargs = len(args)
     23
     24        if nargs == 0:
    2625            self.setdefaultparameters()
    2726        else:
    2827            raise Exception('constructor not supported')
    … …  
    3130    def __repr__(self): # {{{
    3231        s = '   dsl mme parameters:\n'
    3332        s += '{}\n'.format(fielddisplay(self, 'modelid', 'index into the multi-model ensemble, determines which field will be used.'))
    34         s += '{}\n'.format(fielddisplay(self, 'global_average_thermosteric_sea_level_change', 'corresponds to zostoga fields in CMIP5 archives. Specified as a temporally variable global rate (mm/yr) for each ensemble.'))
    35         s += '{}\n'.format(fielddisplay(self, 'sea_surface_height_change_above_geoid', 'corresponds to zos fields in CMIP5 archives. Spatial average is 0. Specified as a spatio-temporally variable rate (mm/yr) for each ensemble.'))
    36         s += '{}\n'.format(fielddisplay(self, 'sea_water_pressure_change_at_sea_floor', 'corresponds to bpo fields in CMIP5 archives. Specified as a spatio-temporally variable rate (in mm/yr) for each ensemble.'))
    37         s += '{}\n'.format(fielddisplay(self, 'compute_fingerprints', 'do we use the sea water pressure change to compute fingerprints and correct sea_surface_height_change_above_geoid'))
     33        s += '{}\n'.format(fielddisplay(self, 'global_average_thermosteric_sea_level', 'Corresponds to zostoga field in CMIP5 archives. Specified as a temporally variable quantity (in m) for each ensemble.'))
     34        s += '{}\n'.format(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) for each ensemble.'))
     35        s += '{}\n'.format(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!) for each ensemble.'))
    3836        return s
    3937    #}}}
    4038
    … …  
    4341    #}}}
    4442
    4543    def checkconsistency(self, md, solution, analyses): # {{{
    46         if ('SealevelriseAnalysis' not in analyses) or (solution == 'TransientSolution' and not md.transient.isslr):
     44        # Early return
     45        if ('SealevelriseAnalysis' not in analyses) or (solution == 'TransientSolution' and not md.transient.isslc):
    4746            return md
    48         for i in range(len(self.global_average_thermosteric_sea_level_change)):
    49             md = checkfield(md, 'field', self.global_average_thermosteric_sea_level_change[i], 'NaN', 1, 'Inf', 1)
    50             md = checkfield(md, 'field', self.sea_surface_height_change_above_geoid[i], 'NaN', 1, 'Inf', 1, 'timeseries', 1)
    51             md = checkfield(md, 'field', self.sea_water_pressure_change_at_sea_floor[i], 'NaN', 1, 'Inf', 1, 'timeseries', 1)
    52         md = checkfield(md, 'field', self.modelid, 'NaN', 1, 'Inf', 1, '>=', 1, '<=',len(self.global_average_thermosteric_sea_level_change))
    53         if self.compute_fingerprints:
    54             #check geodetic flag of slr is on
    55             if not md.solidearth.settings.computesealevelchange:
    56                 raise Exception('DSL checkconsistency error message: if bottom pressure fingerprints computations are requested, slr class should have geodetic flag on')
     47
     48        for i in range(len(self.global_average_thermosteric_sea_level)):
     49            md = checkfield(md, 'field', self.global_average_thermosteric_sea_level[i], 'NaN', 1, 'Inf', 1)
     50            md = checkfield(md, 'field', self.sea_surface_height_above_geoid[i], 'NaN', 1, 'Inf', 1, 'timeseries', 1)
     51            md = checkfield(md, 'field', self.sea_water_pressure_at_sea_floor[i], 'NaN', 1, 'Inf', 1, 'timeseries', 1)
     52        md = checkfield(md, 'field', self.modelid, 'NaN', 1, 'Inf', 1, '>=', 1, '<=',len(self.global_average_thermosteric_sea_level))
     53
     54        if self.solidearth.settings.compute_bp_grd:
     55            md = checkfield(md, 'field', self.sea_water_pressure_at_sea_floor, 'empty', 1)
     56
    5757        return md
    5858    #}}}
    5959
    6060    def marshall(self, prefix, md, fid): #{{{
    6161        WriteData(fid, prefix, 'name', 'md.dsl.model', 'data', 2, 'format', 'Integer')
    62         WriteData(fid, prefix, 'object', self, 'fieldname', 'compute_fingerprints', 'format', 'Integer')
    6362        WriteData(fid, prefix, 'object', self, 'fieldname', 'modelid', 'format', 'Double')
    64         WriteData(fid, prefix, 'name', 'md.dsl.nummodels', 'data', len(self.global_average_thermosteric_sea_level_change), 'format', 'Integer')
    65         WriteData(fid, prefix, 'object', self, 'fieldname', 'global_average_thermosteric_sea_level_change', 'format', 'MatArray', 'timeseries', 1, 'timeserieslength', 2, 'yts', md.constants.yts, 'scale', 1e-3 / md.constants.yts)
    66         WriteData(fid, prefix, 'object', self, 'fieldname', 'sea_water_pressure_change_at_sea_floor', 'format', 'MatArray', 'timeserieslength', md.mesh.numberofvertices + 1, 'yts', md.constants.yts, 'scale', 1e-3)
    67         WriteData(fid, prefix, 'object', self, 'fieldname', 'sea_surface_height_change_above_geoid', 'format', 'MatArray', 'timeserieslength', md.mesh.numberofvertices + 1, 'yts', md.constants.yts, 'scale', 1e-3 / md.constants.yts)
     63        WriteData(fid, prefix, 'name', 'md.dsl.nummodels', 'data', len(self.global_average_thermosteric_sea_level), 'format', 'Integer')
     64        WriteData(fid, prefix, 'object', self, 'fieldname', 'global_average_thermosteric_sea_level', 'format', 'MatArray', 'timeseries', 1, 'timeserieslength', 2, 'yts', md.constants.yts, 'scale', 1e-3 / md.constants.yts)
     65        WriteData(fid, prefix, 'object', self, 'fieldname', 'sea_water_pressure_at_sea_floor', 'format', 'MatArray', 'timeserieslength', md.mesh.numberofvertices + 1, 'yts', md.constants.yts, 'scale', 1e-3)
     66        WriteData(fid, prefix, 'object', self, 'fieldname', 'sea_surface_height_above_geoid', 'format', 'MatArray', 'timeserieslength', md.mesh.numberofvertices + 1, 'yts', md.constants.yts, 'scale', 1e-3 / md.constants.yts)
    6867    #}}}
    6968
    7069    def extrude(self, md): #{{{
    71         for i in range(len(self.global_average_thermosteric_sea_level_change)):
    72             self.sea_surface_height_change_above_geoid[i] = project3d(md, 'vector', self.self.sea_surface_height_change_above_geoid[i], 'type', 'node', 'layer', 1)
    73             self.sea_water_pressure_change_at_sea_floor[i] = project3d(md, 'vector', self.sea_water_pressure_change_at_sea_floor[i], 'type', 'node', 'layer', 1)
     70        for i in range(len(self.global_average_thermosteric_sea_level)):
     71            self.sea_surface_height_above_geoid[i] = project3d(md, 'vector', self.self.sea_surface_height_above_geoid[i], 'type', 'node', 'layer', 1)
     72            self.sea_water_pressure_at_sea_floor[i] = project3d(md, 'vector', self.sea_water_pressure_at_sea_floor[i], 'type', 'node', 'layer', 1)
    7473
    7574        return self
    7675    #}}}

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

     
    7979                        end
    8080
    8181                        %need 'litho' material:
    82                         if ~isa(md.materials,'materials')
     82                        if ~isa(md.materials,'materials') | ~sum(strcmpi(md.materials.nature,'litho'))
    8383                                error('Need a ''litho'' material to run a Fourier Love number analysis');
    84                         else
    85                                 if ~sum(strcmpi(md.materials.nature,'litho')),
    86                                         error('Need a ''litho'' material to run a Fourier Love number analysis');
    87                                 end
    8884                        end
    8985
    9086                end % }}}

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

     
    44
    55
    66class fourierlove(object):
    7     """
    8     Fourier Love Number class definition
     7    """FOURIERLOVE - Fourier Love Number class definition
    98
    10        Usage:
    11           fourierlove = fourierlove()
     9    Usage:
     10        fourierlove = fourierlove()
    1211    """
    1312    def __init__(self):  # {{{
    1413        self.nfreq = 0
    … …  
    2221        self.love_kernels = 0
    2322        self.forcing_type = 0
    2423
    25     #set defaults
    2624        self.setdefaultparameters()
    2725    #}}}
    2826
    2927    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        #
    3032        string = '   Fourier Love class:'
    3133        string = "%s\n%s" % (string, fielddisplay(self, 'nfreq', 'number of frequencies sampled (default 1, elastic) [Hz]'))
    3234        string = "%s\n%s" % (string, fielddisplay(self, 'frequencies', 'frequencies sampled (convention defaults to 0 for the elastic case) [Hz]'))
    … …  
    7577    #}}}
    7678
    7779    def checkconsistency(self, md, solution, analyses):  # {{{
     80        if 'LoveAnalysis' not in analyses:
     81            return md
    7882        md = checkfield(md, 'fieldname', 'love.nfreq', 'NaN', 1, 'Inf', 1, 'numel', [1], '>', 0)
    7983        md = checkfield(md, 'fieldname', 'love.frequencies', 'NaN', 1, 'Inf', 1, 'numel', [md.love.nfreq])
    8084        md = checkfield(md, 'fieldname', 'love.sh_nmax', 'NaN', 1, 'Inf', 1, 'numel', [1], '>', 0)
    … …  
    8892        if md.love.sh_nmin <= 1 and md.love.forcing_type == 9:
    8993            raise RuntimeError("Degree 1 not supported for Volumetric Potential forcing. Use sh_min >= 2 for this kind of calculation.")
    9094
     95        # need 'litho' material
     96        if not hasattr(md.materials, 'materials') or 'litho' not in md.materials.nature:
     97            raise RuntimeError('Need a \'litho\' material to run a Fourier Love number analysis')
    9198        return md
    9299    # }}}
    93100

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

     
    8484
    8585                end % }}}
    8686                function marshall(self,prefix,md,fid) % {{{
    87                        
    88                         if (size(self.thickness==md.mesh.numberofvertices) | (self.thickness==md.mesh.numberofvertices+1)),
     87
     88                        if (size(self.thickness)==md.mesh.numberofvertices) | (size(self.thickness)==md.mesh.numberofvertices+1)),
    8989                                WriteData(fid,prefix,'object',self,'fieldname','thickness','format','DoubleMat','mattype',1,'timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts);
    90                         elseif (size(self.thickness==md.mesh.numberofelements) | (self.thickness==md.mesh.numberofelements+1)),
     90                        elseif (size(self.thickness)==md.mesh.numberofelements) | (size(self.thickness)==md.mesh.numberofelements+1)),
    9191                                WriteData(fid,prefix,'object',self,'fieldname','thickness','format','DoubleMat','mattype',1,'timeserieslength',md.mesh.numberofelements+1,'yts',md.constants.yts);
    9292                        else
    9393                                error('geometry thickness time series should be a vertex or element time series');

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

     
    5050    #}}}
    5151
    5252    def checkconsistency(self, md, solution, analyses): #{{{
    53         if (solution == 'TransientSolution' and md.transient.isgia) or (solution == 'GiaSolution'):
    54             md = checkfield(md, 'fieldname', 'geometry.thickness', 'timeseries', 1, 'NaN', 1, 'Inf', 1)
    55         elif solution == 'SealevelriseSolution':
    56             md = checkfield(md, 'fieldname', 'geometry.bed', 'NaN', 1, 'Inf', 1, 'size', [md.mesh.numberofvertices])
    57         elif solution == 'LoveSolution':
    58             return
     53        if solution == 'LoveSolution':
     54            return md
    5955        else:
    6056            md = checkfield(md, 'fieldname', 'geometry.surface', 'NaN', 1, 'Inf', 1, 'size', [md.mesh.numberofvertices])
    6157            md = checkfield(md, 'fieldname', 'geometry.base', 'NaN', 1, 'Inf', 1, 'size', [md.mesh.numberofvertices])
    … …  
    7066                if np.any(np.abs(self.bed[pos] - self.base[pos]) > 10**-9):
    7167                    md.checkmessage('equality base = bed on grounded ice violated')
    7268                md = checkfield(md, 'fieldname', 'geometry.bed', 'NaN', 1, 'Inf', 1, 'size', [md.mesh.numberofvertices])
    73 
    7469        return md
    7570    # }}}
    7671
    7772    def marshall(self, prefix, md, fid): #{{{
     73        if (len(self.thickness) == md.mesh.numberofvertices) or (len(self.thickness) == md.mesh.numberofvertices + 1):
     74            WriteData(fid, prefix, 'object', self, 'fieldname', 'thickness', 'format', 'DoubleMat', 'mattype', 1, 'timeserieslength', md.mesh.numberofvertices + 1, 'yts', md.constants.yts)
     75        elif (len(self.thickness) == md.mesh.numberofelements) or (len(self.thickness) == md.mesh.numberofelements + 1):
     76            WriteData(fid, prefix, 'object', self, 'fieldname', 'thickness', 'format', 'DoubleMat', 'mattype', 1, 'timeserieslength', md.mesh.numberofvertices + 1, 'yts', md.constants.yts)
     77        else:
     78            raise RuntimeError('geometry thickness time series should be a vertex or element time series')
     79
    7880        WriteData(fid, prefix, 'object', self, 'fieldname', 'surface', 'format', 'DoubleMat', 'mattype', 1)
    79         WriteData(fid, prefix, 'object', self, 'fieldname', 'thickness', 'format', 'DoubleMat', 'mattype', 1, 'timeserieslength', md.mesh.numberofvertices + 1, 'yts', md.constants.yts)
    8081        WriteData(fid, prefix, 'object', self, 'fieldname', 'base', 'format', 'DoubleMat', 'mattype', 1)
    8182        WriteData(fid, prefix, 'object', self, 'fieldname', 'bed', 'format', 'DoubleMat', 'mattype', 1)
    8283        WriteData(fid, prefix, 'object', self, 'fieldname', 'hydrostatic_ratio', 'format', 'DoubleMat', 'mattype', 1)

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

     
    3030
    3131                        %Type of stabilization to use 0:nothing 1:artificial_diffusivity
    3232                        self.stabilization     = 1;
    33       self.requested_outputs = {'default'};
     33                        self.requested_outputs = {'default'};
    3434                end % }}}
    3535                function md = checkconsistency(self,md,solution,analyses) % {{{
    3636                       

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

     
    44
    55
    66class hydrologyshreve(object):
    7     """
    8     HYDROLOGYSHREVE class definition
     7    """HYDROLOGYSHREVE class definition
    98
    10        Usage:
    11           hydrologyshreve = hydrologyshreve()
     9    Usage:
     10        hydrologyshreve = hydrologyshreve()
    1211    """
    1312
    1413    def __init__(self):  # {{{
    … …  
    1514        self.spcwatercolumn = float('NaN')
    1615        self.stabilization = 0
    1716        self.requested_outputs = []
    18     #set defaults
     17
    1918        self.setdefaultparameters()
    20 
    2119    #}}}
    2220
    2321    def __repr__(self):  # {{{
    24 
     22        # TODO:
     23        # - Convert all formatting to calls to <string>.format (see any
     24        #   already converted <class>.__repr__ method for examples)
     25        #
    2526        string = '   hydrologyshreve solution parameters:'
    2627        string = "%s\n%s" % (string, fielddisplay(self, 'spcwatercolumn', 'water thickness constraints (NaN means no constraint) [m]'))
    2728        string = "%s\n%s" % (string, fielddisplay(self, 'stabilization', 'artificial diffusivity (default is 1). can be more than 1 to increase diffusivity.'))
    … …  
    4748
    4849    def checkconsistency(self, md, solution, analyses):  # {{{
    4950        #Early return
    50         if 'HydrologyShreveAnalysis' not in analyses:
     51        if 'HydrologyShreveAnalysis' not in analyses or (solution == 'TransientSolution' and not md.transient.ishydrology):
    5152            return md
    5253
    5354        md = checkfield(md, 'fieldname', 'hydrology.spcwatercolumn', 'Inf', 1, 'timeseries', 1)
    … …  
    6061        WriteData(fid, prefix, 'name', 'md.hydrology.model', 'data', 2, 'format', 'Integer')
    6162        WriteData(fid, prefix, 'object', self, 'fieldname', 'spcwatercolumn', 'format', 'DoubleMat', 'mattype', 1, 'timeserieslength', md.mesh.numberofvertices + 1, 'yts', md.constants.yts)
    6263        WriteData(fid, prefix, 'object', self, 'fieldname', 'stabilization', 'format', 'Double')
    63     #process requested outputs
     64        #process requested outputs
    6465        outputs = self.requested_outputs
    6566        indices = [i for i, x in enumerate(outputs) if x == 'default']
    6667        if len(indices) > 0:

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

     
    2121                channelarea         = NaN;
    2222                sealevel            = NaN;
    2323                bottompressure      = NaN;
    24         sample              = NaN;
     24                sample              = NaN;
    2525        end
    2626        methods
    2727                function self = extrude(self,md) % {{{
    … …  
    5151                        end
    5252                end % }}}
    5353                function self = setdefaultparameters(self) % {{{
    54 
    5554                end % }}}
    5655                function md = checkconsistency(self,md,solution,analyses) % {{{
    5756                        if ismember('StressbalanceAnalysis',analyses) & ~(strcmp(solution,'TransientSolution') & md.transient.isstressbalance == 0),
    … …  
    9796                        end
    9897                        if ismember('HydrologyShreveAnalysis',analyses),
    9998                                if isa(md.hydrology,'hydrologyshreve'),
    100                                         if strcmp(solution,'TransientSolution') & md.transient.ishydrology,
     99                                        if (strcmp(solution,'TransientSolution') & md.transient.ishydrology) | strcmp(solution,'HydrologySolution'),
    101100                                                md = checkfield(md,'fieldname','initialization.watercolumn','NaN',1,'Inf',1,'size',[md.mesh.numberofvertices 1]);
    102101                                        end
    103                                         if strcmp(solution,'HydrologySolution'),
    104                                                 md = checkfield(md,'fieldname','initialization.watercolumn','NaN',1,'Inf',1,'size',[md.mesh.numberofvertices 1]);
    105                                         end
    106102                                end
    107103                        end
    108104                        if ismember('HydrologyTwsAnalysis',analyses),
    … …  
    110106                                        md = checkfield(md,'fieldname','initialization.watercolumn','NaN',1,'Inf',1,'size',[md.mesh.numberofvertices 1]);
    111107                                end
    112108                        end
    113                         if ismember('SealevelchangeAnalysis',analyses) & ~(strcmp(solution,'TransientSolution') & md.transient.isslc == 0),
     109                        if ismember('SealevelchangeAnalysis',analyses),
    114110                                if strcmp(solution,'TransientSolution') & md.transient.isslc,
    115111                                        md = checkfield(md,'fieldname','initialization.sealevel','NaN',1,'Inf',1,'size',[md.mesh.numberofvertices 1]);
    116112                                end
    … …  
    133129                                                md = checkfield(md,'fieldname','initialization.epl_head','NaN',1,'Inf',1,'size',[md.mesh.numberofvertices 1]);
    134130                                                md = checkfield(md,'fieldname','initialization.epl_thickness','NaN',1,'Inf',1,'size',[md.mesh.numberofvertices 1]);
    135131                                        end
    136                 end
    137             end
    138             if ismember('SamplingAnalysis',analyses) & ~(strcmp(solution,'TransientSolution') & md.transient.issampling == 0),
    139                 if ~isnan(md.initialization.sample)
    140                     md = checkfield(md,'fieldname','initialization.sample','NaN',1,'Inf',1,'size',[md.mesh.numberofvertices 1]);
    141                 end
     132                                end
    142133                        end
     134                        if ismember('SamplingAnalysis',analyses) & ~(strcmp(solution,'TransientSolution') & md.transient.issampling == 0),
     135                                if ~isnan(md.initialization.sample)
     136                                        md = checkfield(md,'fieldname','initialization.sample','NaN',1,'Inf',1,'size',[md.mesh.numberofvertices 1]);
     137                                end
     138                        end
    143139                end % }}}
    144140                function disp(self) % {{{
    145141                        disp(sprintf('   initial field values:'));
    … …  
    158154                        fielddisplay(self,'watercolumn','subglacial water sheet thickness (for Shreve and GlaDS) [m]');
    159155                        fielddisplay(self,'hydraulic_potential','Hydraulic potential (for GlaDS) [Pa]');
    160156                        fielddisplay(self,'channelarea','subglacial water channel area (for GlaDS) [m2]');
    161             fielddisplay(self,'sample','Realization of a Gaussian random field');
     157                        fielddisplay(self,'sample','Realization of a Gaussian random field');
    162158                end % }}}
    163159                function marshall(self,prefix,md,fid) % {{{
    164160
    … …  
    178174                        WriteData(fid,prefix,'object',self,'fieldname','watercolumn','format','DoubleMat','mattype',1);
    179175                        WriteData(fid,prefix,'object',self,'fieldname','channelarea','format','DoubleMat','mattype',1);
    180176                        WriteData(fid,prefix,'object',self,'fieldname','hydraulic_potential','format','DoubleMat','mattype',1);
    181             WriteData(fid,prefix,'object',self,'fieldname','sample','format','DoubleMat','mattype',1);
    182            
     177                        WriteData(fid,prefix,'object',self,'fieldname','sample','format','DoubleMat','mattype',1);
     178
    183179                        if md.thermal.isenthalpy,
    184180                                if numel(self.enthalpy) <= 1,
    185181                                        %reconstruct enthalpy
    … …  
    194190                        end
    195191                end % }}}
    196192                function savemodeljs(self,fid,modelname) % {{{
    197                
     193
    198194                        writejs1Darray(fid,[modelname '.initialization.vx'],self.vx);
    199195                        writejs1Darray(fid,[modelname '.initialization.vy'],self.vy);
    200196                        writejs1Darray(fid,[modelname '.initialization.vz'],self.vz);
    … …  
    209205                        writejs1Darray(fid,[modelname '.initialization.watercolumn'],self.watercolumn);
    210206                        writejs1Darray(fid,[modelname '.initialization.hydraulic_potential'],self.hydraulic_potential);
    211207                        writejs1Darray(fid,[modelname '.initialization.channel'],self.channelarea);
    212             writejs1Darray(fid,[modelname '.initialization.sample'],self.sample);
    213            
     208                        writejs1Darray(fid,[modelname '.initialization.sample'],self.sample);
     209
    214210                end % }}}
    215211        end
    216212end

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

     
    2828        self.epl_thickness = np.nan
    2929        self.hydraulic_potential = np.nan
    3030        self.channelarea = np.nan
     31        self.sealevel = np.nan
     32        self.bottompressure = np.nan
    3133        self.sample = np.nan
    3234
    3335        #set defaults
    … …  
    6668        self.sediment_head = project3d(md, 'vector', self.sediment_head, 'type', 'node', 'layer', 1)
    6769        self.epl_head = project3d(md, 'vector', self.epl_head, 'type', 'node', 'layer', 1)
    6870        self.epl_thickness = project3d(md, 'vector', self.epl_thickness, 'type', 'node', 'layer', 1)
     71        self.sealevel = project3d(md, 'vector', self.sealevel, 'type', 'node', 'layer', 1)
     72        self.bottompressure = project3d(md, 'vector', self.bottompressure, 'type', 'node', 'layer', 1)
    6973
    7074        #Lithostatic pressure by default
    7175        if np.ndim(md.geometry.surface) == 2:
    … …  
    8993        if 'MasstransportAnalysis' in analyses and not solution == 'TransientSolution' and not md.transient.ismasstransport:
    9094            md = checkfield(md, 'fieldname', 'initialization.vx', 'NaN', 1, 'Inf', 1, 'size', [md.mesh.numberofvertices])
    9195            md = checkfield(md, 'fieldname', 'initialization.vy', 'NaN', 1, 'Inf', 1, 'size', [md.mesh.numberofvertices])
     96        if 'OceantransportAnalysis' in analyses:
     97            if solution == 'TransientSolution' and md.transient.isslc and md.transient.isoceantransport:
     98                md = checkfield(md, 'fieldname', 'initialization.bottompressure', 'NaN', 1, 'Inf', 1, 'size', [md.mesh.numberofvertices])
    9299        if 'BalancethicknessAnalysis' in analyses and solution == 'BalancethicknessSolution':
    93100            md = checkfield(md, 'fieldname', 'initialization.vx', 'NaN', 1, 'Inf', 1, 'size', [md.mesh.numberofvertices])
    94101            md = checkfield(md, 'fieldname', 'initialization.vy', 'NaN', 1, 'Inf', 1, 'size', [md.mesh.numberofvertices])
    … …  
    111118                md = checkfield(md, 'fieldname', 'delta Tpmp', 'field', np.absolute(md.initialization.temperature[pos] - (md.materials.meltingpoint - md.materials.beta * md.initialization.pressure[pos])), '<', 1e-11, 'message', 'set temperature to pressure melting point at locations with waterfraction > 0')
    112119        if 'HydrologyShreveAnalysis' in analyses:
    113120            if hasattr(md.hydrology, 'hydrologyshreve'):
     121                if (solution == 'TransientSolution' and md.transient.ishydrology) or solution == 'HydrologySolution':
     122                    md = checkfield(md, 'fieldname', 'initialization.watercolumn', 'NaN', 1, 'Inf', 1, 'size', [md.mesh.numberofvertices])
     123        if 'HydrologyTwsAnalysis' in analyses:
     124            if hasattr(md.hydrology, 'hydrologytws'):
    114125                md = checkfield(md, 'fieldname', 'initialization.watercolumn', 'NaN', 1, 'Inf', 1, 'size', [md.mesh.numberofvertices])
    115         if 'HydrologyDCInefficientAnalysis' in analyses:
    116             if hasattr(md.hydrology, 'hydrologydc'):
    117                 md = checkfield(md, 'fieldname', 'initialization.sediment_head', 'NaN', 1, 'Inf', 1, 'size', [md.mesh.numberofvertices])
    118         if 'HydrologyDCEfficientAnalysis' in analyses:
    119             if hasattr(md.hydrology, 'hydrologydc'):
    120                 if md.hydrology.isefficientlayer == 1:
    121                     md = checkfield(md, 'fieldname', 'initialization.epl_head', 'NaN', 1, 'Inf', 1, 'size', [md.mesh.numberofvertices])
    122                     md = checkfield(md, 'fieldname', 'initialization.epl_thickness', 'NaN', 1, 'Inf', 1, 'size', [md.mesh.numberofvertices])
     126        if 'SealevelchangeAnalysis' in analyses:
     127            if solution == 'TransientSolution' and md.transient.isslc:
     128                md = checkfield(md, 'fieldname', 'initialization.sealevel', 'NaN', 1, 'Inf', 1, 'size', [md.mesh.numberofvertices])
    123129        if 'HydrologyGlaDSAnalysis' in analyses:
    124130            if hasattr(md.hydrology, 'hydrologyglads'):
    125131                md = checkfield(md, 'fieldname', 'initialization.watercolumn', 'NaN', 1, 'Inf', 1, 'size', [md.mesh.numberofvertices])
    … …  
    137143        WriteData(fid, prefix, 'object', self, 'fieldname', 'vy', 'format', 'DoubleMat', 'mattype', 1, 'scale', 1. / yts)
    138144        WriteData(fid, prefix, 'object', self, 'fieldname', 'vz', 'format', 'DoubleMat', 'mattype', 1, 'scale', 1. / yts)
    139145        WriteData(fid, prefix, 'object', self, 'fieldname', 'pressure', 'format', 'DoubleMat', 'mattype', 1)
     146        WriteData(fid, prefix, 'object', self, 'fieldname', 'sealevel', 'format', 'DoubleMat', 'mattype', 1)
     147        WriteData(fid, prefix, 'object', self, 'fieldname', 'bottompressure', 'format', 'DoubleMat', 'mattype', 1)
    140148        WriteData(fid, prefix, 'object', self, 'fieldname', 'temperature', 'format', 'DoubleMat', 'mattype', 1)
    141149        WriteData(fid, prefix, 'object', self, 'fieldname', 'waterfraction', 'format', 'DoubleMat', 'mattype', 1)
    142150        WriteData(fid, prefix, 'object', self, 'fieldname', 'sediment_head', 'format', 'DoubleMat', 'mattype', 1)

  • ../trunk-jpl/src/m/consistency/ismodelselfconsistent.py

     
    11def ismodelselfconsistent(md):  #{{{
    2     '''
    3     ISMODELSELFCONSISTENT - check that model forms a closed form solvable problem.
     2    """ISMODELSELFCONSISTENT - check that model forms a closed form solvable problem.
    43
    5        Usage:
    6           ismodelselfconsistent(md),
    7     '''
     4    Usage:
     5        ismodelselfconsistent(md)
     6    """
    87
    98    #initialize consistency as true
    109    md.private.isconsistent = True
    … …  
    5251        analyses = ['EnthalpyAnalysis', 'ThermalAnalysis', 'MeltingAnalysis']
    5352    elif solutiontype == 'MasstransportSolution':
    5453        analyses = ['MasstransportAnalysis']
     54    elif solutiontype == 'OceantransportSolution':
     55        analyses = ['OceantransportAnalysis']
    5556    elif solutiontype == 'BalancethicknessSolution':
    5657        analyses = ['BalancethicknessAnalysis']
    5758    elif solutiontype == 'Balancethickness2Solution':
    … …  
    7172    elif solutiontype == 'EsaSolution':
    7273        analyses = ['EsaAnalysis']
    7374    elif solutiontype == 'TransientSolution':
    74         analyses = ['StressbalanceAnalysis', 'StressbalanceVerticalAnalysis', 'StressbalanceSIAAnalysis', 'L2ProjectionBaseAnalysis', 'ThermalAnalysis', 'MeltingAnalysis', 'EnthalpyAnalysis', 'MasstransportAnalysis', 'HydrologyShaktiAnalysis', 'HydrologyGladsAnalysis', 'HydrologyDCInefficientAnalysis', 'HydrologyDCEfficientAnalysis', 'SealevelriseAnalysis']
     75        analyses = ['StressbalanceAnalysis', 'StressbalanceVerticalAnalysis', 'StressbalanceSIAAnalysis', 'L2ProjectionBaseAnalysis', 'ThermalAnalysis', 'MeltingAnalysis', 'EnthalpyAnalysis', 'MasstransportAnalysis', 'OceantransportAnalysis', 'HydrologyShaktiAnalysis', 'HydrologyGladsAnalysis', 'HydrologyShreveAnalysis', 'HydrologyTwsAnalysis', 'HydrologyDCInefficientAnalysis', 'HydrologyDCEfficientAnalysis', 'SealevelriseAnalysis']
    7576    elif solutiontype == 'SealevelriseSolution':
    7677        analyses = ['SealevelriseAnalysis']
    7778    elif solutiontype == 'HydrologySolution':
    78         analyses = ['L2ProjectionBaseAnalysis', 'HydrologyShreveAnalysis', 'HydrologyDCInefficientAnalysis', 'HydrologyDCEfficientAnalysis']
     79        analyses = ['L2ProjectionBaseAnalysis', 'HydrologyShreveAnalysis', 'HydrologyDCInefficientAnalysis', 'HydrologyDCEfficientAnalysis', 'HydrologyGladsAnalysis', 'HydrologyShaktiAnalysis', 'HydrologyTwsAnalysis']
    7980    elif 'DamageEvolutionSolution':
    8081        analyses = ['DamageEvolutionAnalysis']
    8182    elif 'SamplingSolution':

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

     
    99%   Solution types available comprise:
    1010%   - 'Stressbalance'      or 'sb'
    1111%   - 'Masstransport'      or 'mt'
    12 %   - 'Oceantransport' or 'oceant'
     12%   - 'Oceantransport'     or 'oceant'
    1313%   - 'Thermal'            or 'th'
    1414%   - 'Steadystate'        or 'ss'
    1515%   - 'Transient'          or 'tr'

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

     
    2121%Loading history
    2222md.timestepping.start_time=-2400000; %4,800 kyr :: EVALUATION TIME
    2323md.timestepping.time_step= 2400000; %2,400 kyr :: EVALUATION TIME
    24 % to get rid of default final_time: make sure final_time>start_time
    25 md.timestepping.final_time=2400000; %2,500 kyr
     24% to get rid of default final_time, make sure final_time > start_time
     25md.timestepping.final_time=2400000; %2,400 kyr
    2626md.masstransport.spcthickness=[...
    2727        [md.geometry.thickness; 0],...
    2828        [md.geometry.thickness; 2400000]...
    … …  
    4545%md.cluster=generic('name',oshostname(),'np',3);
    4646md.verbose=verbose('11111111111');
    4747md.verbose.solver=0;
    48 md=solve(md,'tr');
     48md=solve(md,'Transient');
    4949
    5050%Fields and tolerances to track changes
    5151field_names     ={'UGrd'};

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

     
    55
    66classdef matdamageice
    77        properties (SetAccess=public)
    8                 rho_ice                    = 0.;
    9                 rho_water                  = 0.;
    10                 rho_freshwater             = 0.;
    11                 mu_water                   = 0.;
    12                 heatcapacity               = 0.;
    13                 latentheat                 = 0.;
    14                 thermalconductivity        = 0.;
    15                 temperateiceconductivity   = 0.;
     8                rho_ice                         = 0.;
     9                rho_water                       = 0.;
     10                rho_freshwater                  = 0.;
     11                mu_water                        = 0.;
     12                heatcapacity                    = 0.;
     13                latentheat                      = 0.;
     14                thermalconductivity             = 0.;
     15                temperateiceconductivity        = 0.;
    1616                effectiveconductivity_averaging = 0.;
    17                 meltingpoint               = 0.;
    18                 beta                       = 0.;
    19                 mixed_layer_capacity       = 0.;
    20                 thermal_exchange_velocity  = 0.;
    21                 rheology_B   = NaN;
    22                 rheology_n   = NaN;
    23                 rheology_law = '';
     17                meltingpoint                    = 0.;
     18                beta                            = 0.;
     19                mixed_layer_capacity            = 0.;
     20                thermal_exchange_velocity       = 0.;
     21                rheology_B                      = NaN;
     22                rheology_n                      = NaN;
     23                rheology_law                    = '';
    2424
    2525                %slc
    26                 earth_density              = 0;
     26                earth_density                   = 0;
    2727
    2828        end
    2929        methods

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

     
    3030        self.rheology_n = float('NaN')
    3131        self.rheology_law = ''
    3232
    33     #giaivins:
    34         self.lithosphere_shear_modulus = 0.
    35         self.lithosphere_density = 0.
    36         self.mantle_shear_modulus = 0.
    37         self.mantle_density = 0.
     33        #SLC
     34        self.earth_density = 5512  # average density of the Earth, (kg / m^3)
    3835
    39     #SLR
    40         self.earth_density = 5512  # average density of the Earth, (kg / m^3)
    4136        self.setdefaultparameters()
    4237    #}}}
    4338
    4439    def __repr__(self):  # {{{
     40        # TODO:
     41        # - Convert all formatting to calls to <string>.format (see any
     42        #   already converted <class>.__repr__ method for examples)
     43        #
    4544        string = "   Materials:"
    4645        string = "%s\n%s" % (string, fielddisplay(self, "rho_ice", "ice density [kg / m^3]"))
    4746        string = "%s\n%s" % (string, fielddisplay(self, "rho_water", "water density [kg / m^3]"))
    … …  
    5958        string = "%s\n%s" % (string, fielddisplay(self, "rheology_B", "flow law parameter [Pa s^(1 / n)]"))
    6059        string = "%s\n%s" % (string, fielddisplay(self, "rheology_n", "Glen's flow law exponent"))
    6160        string = "%s\n%s" % (string, fielddisplay(self, "rheology_law", "law for the temperature dependance of the rheology: 'None', 'BuddJacka', 'Cuffey', 'CuffeyTemperate', 'Paterson', 'Arrhenius' or 'LliboutryDuval'"))
    62         string = "%s\n%s" % (string, fielddisplay(self, "lithosphere_shear_modulus", "Lithosphere shear modulus [Pa]"))
    63         string = "%s\n%s" % (string, fielddisplay(self, "lithosphere_density", "Lithosphere density [g / cm^ - 3]"))
    64         string = "%s\n%s" % (string, fielddisplay(self, "mantle_shear_modulus", "Mantle shear modulus [Pa]"))
    65         string = "%s\n%s" % (string, fielddisplay(self, "mantle_density", "Mantle density [g / cm^ - 3]"))
    6661        string = "%s\n%s" % (string, fielddisplay(self, "earth_density", "Mantle density [kg / m^ - 3]"))
    6762        return string
    6863    #}}}
    … …  
    10499        #available: none, paterson and arrhenius
    105100        self.rheology_law = 'Paterson'
    106101
    107         # GIA:
    108         self.lithosphere_shear_modulus = 6.7e10  # (Pa)
    109         self.lithosphere_density = 3.32  # (g / cm^ - 3)
    110         self.mantle_shear_modulus = 1.45e11  # (Pa)
    111         self.mantle_density = 3.34  # (g / cm^ - 3)
    112 
    113         #SLR
     102        #SLC
    114103        self.earth_density = 5512  #average density of the Earth, (kg / m^3)
    115104        return self
    116105    #}}}
    … …  
    130119        md = checkfield(md, 'fieldname', 'materials.mantle_density', '>', 0, 'numel', [1])
    131120        md = checkfield(md, 'fieldname', 'materials.earth_density', '>', 0, 'numel', [1])
    132121
    133         if 'GiaAnalysis' in analyses:
    134             md = checkfield(md, 'fieldname', 'materials.lithosphere_shear_modulus', '>', 0, 'numel', 1)
    135             md = checkfield(md, 'fieldname', 'materials.lithosphere_density', '>', 0, 'numel', 1)
    136             md = checkfield(md,'fieldname', 'materials.mantle_shear_modulus', '>', 0, 'numel', 1)
    137             md = checkfield(md,'fieldname', 'materials.mantle_density', '>', 0, 'numel', 1)
    138 
    139122        if 'SealevelriseAnalysis' in analyses:
    140123                md = checkfield(md, 'fieldname', 'materials.earth_density', '>', 0, 'numel', 1)
    141124
    … …  
    160143        WriteData(fid, prefix, 'object', self, 'class', 'materials', 'fieldname', 'rheology_B', 'format', 'DoubleMat', 'mattype', 1)
    161144        WriteData(fid, prefix, 'object', self, 'class', 'materials', 'fieldname', 'rheology_n', 'format', 'DoubleMat', 'mattype', 2)
    162145        WriteData(fid, prefix, 'data', self.rheology_law, 'name', 'md.materials.rheology_law', 'format', 'String')
    163         WriteData(fid, prefix, 'object', self, 'class', 'materials', 'fieldname', 'lithosphere_shear_modulus', 'format', 'Double')
    164         WriteData(fid, prefix, 'object', self, 'class', 'materials', 'fieldname', 'lithosphere_density', 'format', 'Double', 'scale', 10.**3.)
    165         WriteData(fid, prefix, 'object', self, 'class', 'materials', 'fieldname', 'mantle_shear_modulus', 'format', 'Double')
    166         WriteData(fid, prefix, 'object', self, 'class', 'materials', 'fieldname', 'mantle_density', 'format', 'Double', 'scale', 10.**3.)
    167146        WriteData(fid, prefix, 'object', self, 'class', 'materials', 'fieldname', 'earth_density', 'format', 'Double')
    168147
    169148    # }}}

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

     
     1from checkfield import checkfield
    12from fielddisplay import fielddisplay
    23from project3d import project3d
    3 from checkfield import checkfield
    44from WriteData import WriteData
    55
    66
    77class matenhancedice(object):
    8     """
    9     MATICE class definition
     8    """MATICE class definition
    109
    11         Usage:
    12             matenhancedice = matenhancedice()
     10    Usage:
     11        matenhancedice = matenhancedice()
    1312    """
    1413
    1514    def __init__(self): #{{{
    … …  
    3130        self.rheology_n = float('NaN')
    3231        self.rheology_law = ''
    3332
    34         #GIA
    35         self.lithosphere_shear_modulus = 0.
    36         self.lithosphere_density = 0.
    37         self.mantle_shear_modulus = 0.
    38         self.mantle_density = 0.
     33        #SLC
     34        self.earth_density = 0  # average density of the Earth, (kg/m^3)
    3935
    40         #SLR
    41         self.earth_density = 0  # average density of the Earth, (kg/m^3)
    4236        self.setdefaultparameters()
    4337    #}}}
    4438
    4539    def __repr__(self): #{{{
     40        # TODO:
     41        # - Convert all formatting to calls to <string>.format (see any
     42        #   already converted <class>.__repr__ method for examples)
     43        #
    4644        s = "   Materials:"
    4745        s = "%s\n%s" % (s, fielddisplay(self, "rho_ice", "ice density [kg/m^3]"))
    4846        s = "%s\n%s" % (s, fielddisplay(self, "rho_water", "water density [kg/m^3]"))
    … …  
    6159        s = "%s\n%s" % (s, fielddisplay(self, "rheology_B", "flow law parameter [Pa s^(1/n)]"))
    6260        s = "%s\n%s" % (s, fielddisplay(self, "rheology_n", "Glen's flow law exponent"))
    6361        s = "%s\n%s" % (s, fielddisplay(self, "rheology_law", "law for the temperature dependance of the rheology: 'None', 'BuddJacka', 'Cuffey', 'CuffeyTemperate', 'Paterson', 'Arrhenius' or 'LliboutryDuval'"))
    64         s = "%s\n%s" % (s, fielddisplay(self, "lithosphere_shear_modulus", "Lithosphere shear modulus [Pa]"))
    65         s = "%s\n%s" % (s, fielddisplay(self, "lithosphere_density", "Lithosphere density [g/cm^-3]"))
    66         s = "%s\n%s" % (s, fielddisplay(self, "mantle_shear_modulus", "Mantle shear modulus [Pa]"))
    67         s = "%s\n%s" % (s, fielddisplay(self, "mantle_density", "Mantle density [g/cm^-3]"))
    6862        s = "%s\n%s" % (s, fielddisplay(self, "earth_density", "Mantle density [kg/m^-3]"))
    6963
    7064        return s
    … …  
    108102        #available: none, paterson and arrhenius
    109103        self.rheology_law = 'Paterson'
    110104
    111     # GIA:
     105        #GIA
    112106        self.lithosphere_shear_modulus = 6.7 * 10**10  # (Pa)
    113107        self.lithosphere_density = 3.32  # (g / cm^ - 3)
    114108        self.mantle_shear_modulus = 1.45 * 10**11  # (Pa)
    115109        self.mantle_density = 3.34  # (g / cm^ - 3)
    116110
    117     #SLR
     111        #SLC
    118112        self.earth_density = 5512  #average density of the Earth, (kg / m^3)
    119113
    120114        return self
    … …  
    131125        md = checkfield(md, 'fieldname', 'materials.rheology_law', 'values', ['None', 'BuddJacka', 'Cuffey', 'CuffeyTemperate', 'Paterson', 'Arrhenius', 'LliboutryDuval'])
    132126        md = checkfield(md, 'fieldname', 'materials.effectiveconductivity_averaging', 'numel', [1], 'values', [0, 1, 2])
    133127
    134         if 'GiaAnalysis' in analyses:
    135             md = checkfield(md, 'fieldname', 'materials.lithosphere_shear_modulus', '>', 0, 'numel', 1)
    136             md = checkfield(md, 'fieldname', 'materials.lithosphere_density', '>', 0, 'numel', 1)
    137             md = checkfield(md, 'fieldname', 'materials.mantle_shear_modulus', '>', 0, 'numel', 1)
    138             md = checkfield(md, 'fieldname', 'materials.mantle_density', '>', 0, 'numel', 1)
    139            
    140128        if 'SealevelriseAnalysis' in analyses:
    141129            md = checkfield(md, 'fieldname', 'materials.earth_density', '>', 0, 'numel', 1)
    142130        return md
    … …  
    161149        WriteData(fid, prefix, 'object', self, 'class', 'materials', 'fieldname', 'rheology_B', 'format', 'DoubleMat', 'mattype', 1, 'timeserieslength', md.mesh.numberofvertices + 1, 'yts', md.constants.yts)
    162150        WriteData(fid, prefix, 'object', self, 'class', 'materials', 'fieldname', 'rheology_n', 'format', 'DoubleMat', 'mattype', 2)
    163151        WriteData(fid, prefix, 'data', self.rheology_law, 'name', 'md.materials.rheology_law', 'format', 'String')
    164         WriteData(fid, prefix, 'object', self, 'class', 'materials', 'fieldname', 'lithosphere_shear_modulus', 'format', 'Double')
    165         WriteData(fid, prefix, 'object', self, 'class', 'materials', 'fieldname', 'lithosphere_density', 'format', 'Double', 'scale', 10**3)
    166         WriteData(fid, prefix, 'object', self, 'class', 'materials', 'fieldname', 'mantle_shear_modulus', 'format', 'Double')
    167         WriteData(fid, prefix, 'object', self, 'class', 'materials', 'fieldname', 'mantle_density', 'format', 'Double', 'scale', 10**3)
    168152        WriteData(fid, prefix, 'object', self, 'class', 'materials', 'fieldname', 'earth_density', 'format', 'Double')
    169153    # }}}

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

     
    122122                                        self.rheology_B   = 1*1e8;
    123123                                        self.rheology_n   = 3;
    124124
    125 
    126125                                case 'litho'
    127126                                        %we default to a configuration that enables running GIA solutions using giacaron and/or giaivins.
    128127                                        self.numlayers=2;

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

     
    77
    88
    99class materials(object):
    10     '''
    11     MATERIALS class definition
     10    """MATERIALS class definition
    1211
    13        Usage:
    14           materials = materials()
    15     '''
     12    Usage:
     13        materials = materials()
     14    """
    1615
    1716    def __init__(self, *args): #{{{
    1817        self.nature = []
    … …  
    3736                setattr(self, 'latentheat', 0)
    3837                setattr(self, 'thermalconductivity', 0)
    3938                setattr(self, 'temperateiceconductivity', 0)
     39                setattr(self, 'effectiveconductivity_averaging', 0)
    4040                setattr(self, 'meltingpoint', 0)
    4141                setattr(self, 'beta', 0)
    4242                setattr(self, 'mixed_layer_capacity', 0)
    … …  
    5959                setattr(self, 'rho_ice', 0)
    6060                setattr(self, 'rho_water', 0)
    6161                setattr(self, 'rho_freshwater', 0)
    62                 setattr(self, 'earth_density', 0)
    6362            else:
    6463                raise RuntimeError("materials constructor error message: nature of the material not supported yet! ('ice' or 'litho' or 'hydro')")
     64        setattr(self, 'earth_density', 0)
    6565
    6666        #set default parameters:
    6767        self.setdefaultparameters()
    … …  
    7373            if nat == 'ice':
    7474                #ice density (kg/m^3)
    7575                self.rho_ice = 917.
     76
    7677                #ocean water density (kg/m^3)
    7778                self.rho_water = 1023.
     79
    7880                #fresh water density (kg/m^3)
    7981                self.rho_freshwater = 1000.
     82
    8083                #water viscosity (N.s/m^2)
    8184                self.mu_water = 0.001787
     85
    8286                #ice heat capacity cp (J/kg/K)
    8387                self.heatcapacity = 2093.
     88
    8489                #ice latent heat of fusion L (J / kg)
    8590                self.latentheat = 3.34e5
     91
    8692                #ice thermal conductivity (W/m/K)
    8793                self.thermalconductivity = 2.4
     94
    8895                #wet ice thermal conductivity (W/m/K)
    8996                self.temperateiceconductivity = 0.24
     97
     98                #computation of effective conductivity
     99                self.effectiveconductivity_averaging = 1
     100
    90101                #the melting point of ice at 1 atmosphere of pressure in K
    91102                self.meltingpoint = 273.15
     103
    92104                #rate of change of melting point with pressure (K/Pa)
    93105                self.beta = 9.8e-8
     106
    94107                #mixed layer (ice-water interface) heat capacity (J/kg/K)
    95108                self.mixed_layer_capacity = 3974.
     109
    96110                #thermal exchange velocity (ice-water interface) (m/s)
    97111                self.thermal_exchange_velocity = 1.00e-4
     112
    98113                #Rheology law: what is the temperature dependence of B with T
    99114                #available: none, paterson and arrhenius
    100115                self.rheology_law = 'Paterson'
     116
     117                #Rheology fields default
     118                self.rheology_B = 1e8
     119                self.rheology_n = 3
    101120            elif nat == 'litho':
    102                 #we default to a configuration that enables running GIA solutions using giacaron and / or giaivins.
     121                #we default to a configuration that enables running GIA solutions using giacaron and/or giaivins.
    103122                self.numlayers = 2
     123
    104124                #center of the earth (approximation, must not be 0), then the lab (lithosphere / asthenosphere boundary) then the surface
    105125                #(with 1d3 to avoid numerical singularities)
    106126                self.radius = [1e3, 6278e3, 6378e3]
     127
    107128                self.viscosity = [1e21, 1e40]  #mantle and lithosphere viscosity (respectively) [Pa.s]
    108129                self.lame_mu = [1.45e11, 6.7e10]  # (Pa)  #lithosphere and mantle shear modulus (respectively) [Pa]
    109130                self.lame_lambda = self.lame_mu  # (Pa)  #mantle and lithosphere lamba parameter (respectively) [Pa]
    … …  
    115136            elif nat == 'hydro':
    116137                #ice density (kg/m^3)
    117138                self.rho_ice = 917.
     139
    118140                #ocean water density (kg/m^3)
    119141                self.rho_water = 1023.
    120                 #average density of the Earth (kg/m^3)
    121                 self.earth_density = 5512
     142
    122143                #fresh water density (kg/m^3)
    123144                self.rho_freshwater = 1000.
    124145            else:
    125146                raise RuntimeError("materials setdefaultparameters error message: nature of the material not supported yet! ('ice' or 'litho' or 'hydro')")
     147
     148            #average density of the Earth (kg/m^3)
     149            self.earth_density = 5512
    126150    #}}}
    127151
    128152    def __repr__(self): #{{{
    … …  
    222246        #1: MatdamageiceEnum 2: MatestarEnum 3: MaticeEnum 4: MatenhancediceEnum 5: MaterialsEnum
    223247        WriteData(fid, prefix, 'name', 'md.materials.nature', 'data', naturetointeger(self.nature), 'format', 'IntMat', 'mattype', 3)
    224248        WriteData(fid, prefix, 'name', 'md.materials.type', 'data', 5, 'format', 'Integer') #DANGER, this can evolve if you have classes
    225 
    226249        for i in range(len(self.nature)):
    227250            nat = self.nature[i]
    228251            if nat == 'ice':
    … …  
    235258                WriteData(fid, prefix, 'object', self, 'class', 'materials', 'fieldname', 'latentheat', 'format', 'Double')
    236259                WriteData(fid, prefix, 'object', self, 'class', 'materials', 'fieldname', 'thermalconductivity', 'format', 'Double')
    237260                WriteData(fid, prefix, 'object', self, 'class', 'materials', 'fieldname', 'temperateiceconductivity', 'format', 'Double')
     261                WriteData(fid, prefix, 'object', self, 'class', 'materials', 'fieldname', 'effectiveconductivity_averaging', 'format', 'Integer')
    238262                WriteData(fid, prefix, 'object', self, 'class', 'materials', 'fieldname', 'meltingpoint', 'format', 'Double')
    239263                WriteData(fid, prefix, 'object', self, 'class', 'materials', 'fieldname', 'beta', 'format', 'Double')
    240264                WriteData(fid, prefix, 'object', self, 'class', 'materials', 'fieldname', 'mixed_layer_capacity', 'format', 'Double')
    … …  
    253277                WriteData(fid, prefix, 'object', self, 'class', 'materials', 'fieldname', 'isburgers', 'format', 'DoubleMat', 'mattype', 3)
    254278                WriteData(fid, prefix, 'object', self, 'class', 'materials', 'fieldname', 'burgers_viscosity', 'format', 'DoubleMat', 'mattype', 3)
    255279                WriteData(fid, prefix, 'object', self, 'class', 'materials', 'fieldname', 'burgers_mu', 'format', 'DoubleMat', 'mattype', 3)
     280                # Compute earth density compatible with our layer density distribution
     281                earth_density = 0
     282                for i in range(len(self.numlayers)):
     283                    earth_density = earth_density + (self.radius[i + 1] ** 3 - self.radius[i] ** 3) * self.density[i]
     284                earth_density = earth_density / self.radius[self.numlayers + 1] ** 3
     285                self.earth_density = earth_density
    256286            elif nat == 'hydro':
    257287                WriteData(fid, prefix, 'object', self, 'class', 'materials', 'fieldname', 'rho_ice', 'format', 'Double')
    258288                WriteData(fid, prefix, 'object', self, 'class', 'materials', 'fieldname', 'rho_water', 'format', 'Double')
    259                 WriteData(fid, prefix, 'object', self, 'class', 'materials', 'fieldname', 'earth_density', 'format', 'Double')
    260289                WriteData(fid, prefix, 'object', self, 'class', 'materials', 'fieldname', 'rho_freshwater', 'format', 'Double')
    261290            else:
    262291                raise RuntimeError("materials constructor error message: nature of the material not supported yet! ('ice' or 'litho' or 'hydro')")
     292        WriteData(fid, prefix, 'data', self.earth_density, 'name', 'md.materials.earth_density', 'format', 'Double')
    263293    #}}}
    264294
    265295    def extrude(self, md): #{{{

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

     
    3333        self.rheology_Es = np.nan
    3434        self.rheology_law = ''
    3535
    36         #GIA
    37         self.lithosphere_shear_modulus = 0.
    38         self.lithosphere_density = 0.
    39         self.mantle_shear_modulus = 0.
    40         self.mantle_density = 0.
    41 
    42         #SLR
     36        #slc
    4337        self.earth_density = 0
    4438
    4539        #set default parameters
    … …  
    6559        s = "%s\n%s" % (s, fielddisplay(self, 'rheology_Ec', 'compressive enhancement factor'))
    6660        s = "%s\n%s" % (s, fielddisplay(self, 'rheology_Es', 'shear enhancement factor'))
    6761        s = "%s\n%s" % (s, fielddisplay(self, 'rheology_law', ['law for the temperature dependance of the rheology: ''None'', ''BuddJacka'', ''Cuffey'', ''CuffeyTemperate'', ''Paterson'', ''Arrhenius'' or ''LliboutryDuval''']))
    68         s = "%s\n%s" % (s, fielddisplay(self, 'lithosphere_shear_modulus', 'Lithosphere shear modulus [Pa]'))
    69         s = "%s\n%s" % (s, fielddisplay(self, 'lithosphere_density', 'Lithosphere density [g/cm^-3]'))
    70         s = "%s\n%s" % (s, fielddisplay(self, 'mantle_shear_modulus', 'Mantle shear modulus [Pa]'))
    71         s = "%s\n%s" % (s, fielddisplay(self, 'mantle_density', 'Mantle density [g/cm^-3]'))
    7262        s = "%s\n%s" % (s, fielddisplay(self, 'earth_density', 'Mantle density [kg/m^-3]'))
    7363
    7464        return s
    … …  
    112102        #Rheology law: what is the temperature dependence of B with T
    113103        #available: none, paterson and arrhenius
    114104        self.rheology_law = 'Paterson'
    115         # GIA
    116         self.lithosphere_shear_modulus = 6.7e10  # (Pa)
    117         self.lithosphere_density = 3.32  # (g / cm^ - 3)
    118         self.mantle_shear_modulus = 1.45e11  # (Pa)
    119         self.mantle_density = 3.34  # (g / cm^ - 3)
    120         #SLR
     105        #slc
    121106        self.earth_density = 5512  # average density of the Earth, (kg / m^3)
    122107
    123108        return self
    … …  
    134119        md = checkfield(md, 'fieldname', 'materials.rheology_law', 'values', ['None', 'BuddJacka', 'Cuffey', 'CuffeyTemperate', 'Paterson', 'Arrhenius', 'LliboutryDuval'])
    135120        md = checkfield(md, 'fieldname', 'materials.effectiveconductivity_averaging', 'numel', [1], 'values', [0, 1, 2])
    136121
    137         if 'GiaAnalysis' in analyses:
    138             md = checkfield(md, 'fieldname', 'materials.lithosphere_shear_modulus', '>', 0, 'numel', 1)
    139             md = checkfield(md, 'fieldname', 'materials.lithosphere_density', '>', 0, 'numel', 1)
    140             md = checkfield(md, 'fieldname', 'materials.mantle_shear_modulus', '>', 0, 'numel', 1)
    141             md = checkfield(md, 'fieldname', 'materials.mantle_density', '>', 0, 'numel', 1)
    142 
    143122        if 'SealevelriseAnalysis' in analyses:
    144123            md = checkfield(md, 'fieldname', 'materials.earth_density', '>', 0, 'numel', 1)
    145124
    … …  
    165144        WriteData(fid, prefix, 'object', self, 'class', 'materials', 'fieldname', 'rheology_Ec', 'format', 'DoubleMat', 'mattype', 1)
    166145        WriteData(fid, prefix, 'object', self, 'class', 'materials', 'fieldname', 'rheology_Es', 'format', 'DoubleMat', 'mattype', 1)
    167146        WriteData(fid, prefix, 'data', self.rheology_law, 'name', 'md.materials.rheology_law', 'format', 'String')
    168         WriteData(fid, prefix, 'object', self, 'class', 'materials', 'fieldname', 'lithosphere_shear_modulus', 'format', 'Double')
    169         WriteData(fid, prefix, 'object', self, 'class', 'materials', 'fieldname', 'lithosphere_density', 'format', 'Double', 'scale', 1.0e3)
    170         WriteData(fid, prefix, 'object', self, 'class', 'materials', 'fieldname', 'mantle_shear_modulus', 'format', 'Double')
    171         WriteData(fid, prefix, 'object', self, 'class', 'materials', 'fieldname', 'mantle_density', 'format', 'Double', 'scale', 1.0e3)
    172147        WriteData(fid, prefix, 'object', self, 'class', 'materials', 'fieldname', 'earth_density', 'format', 'Double')
    173148    # }}}

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

     
    103103
    104104                end % }}}
    105105                function md = checkconsistency(self,md,solution,analyses) % {{{
    106                        
     106
    107107                        if strcmpi(solution,'TransientSolution') & md.transient.isslc,
    108108                                md = checkfield(md,'fieldname','materials.earth_density','>',0,'numel',1);
    109                                 return;
    110                         end
    111        
    112                         md = checkfield(md,'fieldname','materials.rho_ice','>',0);
    113                         md = checkfield(md,'fieldname','materials.rho_water','>',0);
    114                         md = checkfield(md,'fieldname','materials.rho_freshwater','>',0);
    115                         md = checkfield(md,'fieldname','materials.mu_water','>',0);
    116                         md = checkfield(md,'fieldname','materials.rheology_B','>',0,'universal',1,'NaN',1,'Inf',1);
    117                         md = checkfield(md,'fieldname','materials.rheology_n','>',0,'universal',1,'NaN',1,'Inf',1);
    118                         md = checkfield(md,'fieldname','materials.rheology_law','values',{'None' 'BuddJacka' 'Cuffey' 'CuffeyTemperate' 'Paterson' 'Arrhenius' 'LliboutryDuval' 'NyeCO2' 'NyeH2O'});
    119                         md = checkfield(md,'fieldname','materials.effectiveconductivity_averaging','numel',[1],'values',[0 1 2]);
     109                        else
     110                                md = checkfield(md,'fieldname','materials.rho_ice','>',0);
     111                                md = checkfield(md,'fieldname','materials.rho_water','>',0);
     112                                md = checkfield(md,'fieldname','materials.rho_freshwater','>',0);
     113                                md = checkfield(md,'fieldname','materials.mu_water','>',0);
     114                                md = checkfield(md,'fieldname','materials.rheology_B','>',0,'universal',1,'NaN',1,'Inf',1);
     115                                md = checkfield(md,'fieldname','materials.rheology_n','>',0,'universal',1,'NaN',1,'Inf',1);
     116                                md = checkfield(md,'fieldname','materials.rheology_law','values',{'None' 'BuddJacka' 'Cuffey' 'CuffeyTemperate' 'Paterson' 'Arrhenius' 'LliboutryDuval' 'NyeCO2' 'NyeH2O'});
     117                                md = checkfield(md,'fieldname','materials.effectiveconductivity_averaging','numel',[1],'values',[0 1 2]);
    120118
    121119                end % }}}
    122120                function disp(self) % {{{

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

     
    4343                frontalforcings  = 0;
    4444                love                     = 0;
    4545                esa              = 0;
    46         sampling         = 0;
     46                sampling         = 0;
    4747
    4848                autodiff         = 0;
    4949                inversion        = 0;
    … …  
    155155                        %2019 Jan..
    156156                        if isa(md.frontalforcings,'double');
    157157                                if(isprop('meltingrate',md.calving) & ~isnan(md.calving.meltingrate))
    158                                         disp('Warning: md.calving.meltingrate is now in md.frontalforcings');
     158                gia                     disp('Warning: md.calving.meltingrate is now in md.frontalforcings');
    159159                                end
    160160                                md.frontalforcings=frontalforcings(md.calving);
    161161                        end
    … …  
    186186                                                md.smb.isconstrainsurfaceT = 0;
    187187                                        end
    188188                                end
    189             end
    190             %2021 February 17
     189                        end
     190                        %2021 February 17
    191191                        if isa(md.sampling,'double'); md.sampling=sampling(); end
    192192                end% }}}
    193193        end
    … …  
    241241                        md.frontalforcings  = frontalforcings();
    242242                        md.love             = fourierlove();
    243243                        md.esa              = esa();
    244             md.sampling         = sampling();
     244                        md.sampling         = sampling();
    245245                        md.autodiff         = autodiff();
    246246                        md.inversion        = inversion();
    247247                        md.qmu              = qmu();
    … …  
    614614                                                %size = number of elements * n
    615615                                                elseif fieldsize(1)==numberofelements1
    616616                                                        md2.(model_fields{i}).(object_fields{j})=field(pos_elem,:);
    617                         elseif (fieldsize(1)==numberofelements1+1)
    618                             md2.(model_fields{i}).(object_fields{j})=[field(pos_elem,:); field(end,:)];
     617                                                elseif (fieldsize(1)==numberofelements1+1)
     618                                                        md2.(model_fields{i}).(object_fields{j})=[field(pos_elem,:); field(end,:)];
    619619                                                end
    620620                                        end
    621621                                else
    … …  
    627627                                        %size = number of elements * n
    628628                                        elseif fieldsize(1)==numberofelements1
    629629                                                md2.(model_fields{i})=field(pos_elem,:);
    630                     elseif (fieldsize(1)==numberofelements1+1)
     630                                        elseif (fieldsize(1)==numberofelements1+1)
    631631                                                md2.(model_fields{i})=[field(pos_elem,:); field(end,:)];
    632632                                        end
    633633                                end
    … …  
    774774                                                %get subfields
    775775                                                % loop over time steps
    776776                                                for p=1:length(md1.results.(solutionfields{i}))
    777                                                     current = md1.results.(solutionfields{i})(p);
    778                                                     solutionsubfields=fields(current);
    779                                                     for j=1:length(solutionsubfields),
     777                                                        current = md1.results.(solutionfields{i})(p);
     778                                                        solutionsubfields=fields(current);
     779                                                        for j=1:length(solutionsubfields),
    780780                                                        field=md1.results.(solutionfields{i})(p).(solutionsubfields{j});
    781781                                                        if length(field)==numberofvertices1,
    782                                                             md2.results.(solutionfields{i})(p).(solutionsubfields{j})=field(pos_node);
     782                                                                md2.results.(solutionfields{i})(p).(solutionsubfields{j})=field(pos_node);
    783783                                                        elseif length(field)==numberofelements1,
    784                                                             md2.results.(solutionfields{i})(p).(solutionsubfields{j})=field(pos_elem);
     784                                                                md2.results.(solutionfields{i})(p).(solutionsubfields{j})=field(pos_elem);
    785785                                                        else
    786                                                             md2.results.(solutionfields{i})(p).(solutionsubfields{j})=field;
     786                                                                md2.results.(solutionfields{i})(p).(solutionsubfields{j})=field;
    787787                                                        end
    788                                                     end
     788                                                        end
    789789                                                end
    790790                                        else
    791791                                                field=md1.results.(solutionfields{i});
    … …  
    15001500                        disp(sprintf('%19s: %-22s -- %s','frontalforcings' ,['[1x1 ' class(self.frontalforcings) ']'],'parameters for frontalforcings'));
    15011501                        disp(sprintf('%19s: %-22s -- %s','esa'             ,['[1x1 ' class(self.esa) ']'],'parameters for elastic adjustment solution'));
    15021502                        disp(sprintf('%19s: %-22s -- %s','love'            ,['[1x1 ' class(self.love) ']'],'parameters for love solution'));
    1503             disp(sprintf('%19s: %-22s -- %s','sampling'        ,['[1x1 ' class(self.sampling) ']'],'parameters for stochastic sampler'));
     1503                        disp(sprintf('%19s: %-22s -- %s','sampling'        ,['[1x1 ' class(self.sampling) ']'],'parameters for stochastic sampler'));
    15041504                        disp(sprintf('%19s: %-22s -- %s','autodiff'        ,['[1x1 ' class(self.autodiff) ']'],'automatic differentiation parameters'));
    15051505                        disp(sprintf('%19s: %-22s -- %s','inversion'       ,['[1x1 ' class(self.inversion) ']'],'parameters for inverse methods'));
    15061506                        disp(sprintf('%19s: %-22s -- %s','qmu'             ,['[1x1 ' class(self.qmu) ']'],'Dakota properties'));

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

     
    5454from thermal import thermal
    5555from steadystate import steadystate
    5656from transient import transient
    57 from giaivins import giaivins
    58 from giamme import giamme
    5957from esa import esa
    6058from autodiff import autodiff
    6159from inversion import inversion
    … …  
    113111        self.levelset = None
    114112        self.calving = None
    115113        self.frontalforcings = None
    116         self.gia = None
    117114        self.love = None
    118115        self.esa = None
    119116        self.sampling = None
    … …  
    171168                'levelset',
    172169                'calving',
    173170                'frontalforcings',
    174                 'gia',
    175171                'love',
    176172                'esa',
    177173                'sampling',
    … …  
    224220        s = "%s\n%s" % (s, "%19s: %-22s -- %s" % ("levelset", "[%s %s]" % ("1x1", obj.levelset.__class__.__name__), "parameters for moving boundaries (level - set method)"))
    225221        s = "%s\n%s" % (s, "%19s: %-22s -- %s" % ("calving", "[%s %s]" % ("1x1", obj.calving.__class__.__name__), "parameters for calving"))
    226222        s = "%s\n%s" % (s, "%19s: %-22s -- %s" % ("frontalforcings", "[%s %s]" % ("1x1", obj.frontalforcings.__class__.__name__), "parameters for frontalforcings"))
    227         s = "%s\n%s" % (s, "%19s: %-22s -- %s" % ("gia", "[%s %s]" % ("1x1", obj.gia.__class__.__name__), "parameters for gia solution"))
    228223        s = "%s\n%s" % (s, '%19s: %-22s -- %s' % ("esa", "[%s %s]" % ("1x1", obj.esa.__class__.__name__), "parameters for elastic adjustment solution"))
    229224        s = "%s\n%s" % (s, '%19s: %-22s -- %s' % ("sampling", "[%s %s]" % ("1x1", obj.sampling.__class__.__name__), "parameters for stochastic sampler"))
    230225        s = "%s\n%s" % (s, '%19s: %-22s -- %s' % ("love", "[%s %s]" % ("1x1", obj.love.__class__.__name__), "parameters for love solution"))
    … …  
    271266        self.levelset = levelset()
    272267        self.calving = calving()
    273268        self.frontalforcings = frontalforcings()
    274         self.gia = giamme()
    275269        self.love = fourierlove()
    276270        self.esa = esa()
    277271        self.sampling = sampling()
    … …  
    851845        if not np.isnan(md.initialization.watercolumn).all():
    852846            md.initialization.watercolumn = project2d(md, md.initialization.watercolumn, 1)
    853847
    854         # giaivins
    855         if md.gia.__class__.__name__ == 'giaivins':
    856             if not np.isnan(md.gia.mantle_viscosity).all():
    857                 md.gia.mantle_viscosity = project2d(md, md.gia.mantle_viscosity, 1)
    858             if not np.isnan(md.gia.lithosphere_thickness).all():
    859                 md.gia.lithosphere_thickness = project2d(md, md.gia.lithosphere_thickness, 1)
    860 
    861848        # elementstype
    862849        if not np.isnan(md.flowequation.element_equation).all():
    863850            md.flowequation.element_equation = project2d(md, md.flowequation.element_equation, 1)

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

     
    9090                        for i=1:length(slm.icecaps),
    9191                                md= slm.icecaps{i};
    9292                                if ~isempty(md.solidearth.external),
    93                                         error(sprintf('cannot run external forcings on an ice sheet when runing a coupling earth/ice sheet model');
     93                                        error('cannot run external forcings on an ice sheet when running a coupling earth/ice sheet model');
    9494                                end
    9595
    9696                        end
    97                         %make sure that we have the rigth grd model for computing ouf sealevel patterns:
     97                        %make sure that we have the right grd model for computing out sealevel patterns:
    9898                        for i=1:length(slm.icecaps),
    9999                                md= slm.icecaps{i};
    100100                                if md.solidearth.settings.grdmodel~=0
    101                                         error(sprintf('sealevelmodel checkconsistenty error message: ice sheets do not run GRD module, specify solidearth.settings.grdmodel=0 on ice cap %i',i));
     101                                        error(sprintf('sealevelmodel checkconsistency error message: ice sheets do not run GRD module, specify solidearth.settings.grdmodel=0 on ice cap %i',i));
    102102                                end
    103103                        end
    104104

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

     
    110110            md = slm.icecaps[i]
    111111            if np.nonzero(md.slr.steric_rate - slm.earth.slr.steric_rate[slm.earth.slr.transitions[i]]) != []:
    112112                raise Exception('steric rate on ice cap {} is not the same as for the earth'.format(md.miscellaneous.name))
     113
     114        # Make sure grd is the same everywhere
     115        for i in range(len(slm.icecaps)):
     116            md = slm.icecaps[i]
     117            if md.solidearthsettings.isgrd != slm.earth.solidearthsettings.isgrd:
     118                raise RuntimeError('isgrd on ice cap {} is not the same as for the earth\n'.format(md.miscellaneous.name))
     119        # Make sure that there is no solid earth external forcing on the basins
     120        for i in range(len(slm.icecaps)):
     121            md = slm.icecaps[i]
     122            if md.solidearth.external:
     123                raise RuntimeError('cannot run external forcings on an ice sheet when running a coupling earth/ice sheet model')
     124        # Make sure that we have the right grd model for computing our sealevel patterns
     125        for i in range(len(slm.icecaps)):
     126            md = slm.icecaps[i]
     127            if md.solidearth.settings.grdmodel != 0:
     128                raise RuntimeError('sealevelmodel checkconsistency error message: ice sheets do not run GRD module, specify solidearth.settings.grdmodel=0 on ice cap {}'.format(i))
    113129    #}}}
    114130
    115131    def mergeresults(self):  # {{{

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

     
    1313    def __init__(self):  # {{{
    1414        self.issmb = 0
    1515        self.ismasstransport = 0
     16        self.isoceantransport = 0
    1617        self.isstressbalance = 0
    1718        self.isthermal = 0
    1819        self.isgroundingline = 0
    19         self.isgia = 0
    2020        self.isesa = 0
    2121        self.isdamageevolution = 0
    2222        self.ismovingfront = 0
    … …  
    2323        self.ishydrology = 0
    2424        self.issampling = 0
    2525        self.isslc = 0
    26         self.iscoupler = 0
    2726        self.amr_frequency = 0
    2827        self.isoceancoupling = 0
    2928        self.requested_outputs = []
    … …  
    3534        s = '   transient solution parameters:\n'
    3635        s += '{}\n'.format(fielddisplay(self, 'issmb', 'indicates if a surface mass balance solution is used in the transient'))
    3736        s += '{}\n'.format(fielddisplay(self, 'ismasstransport', 'indicates if a masstransport solution is used in the transient'))
     37        s += '{}\n'.format(fielddisplay(self, 'isoceantransport', 'indicates whether an ocean masstransport solution is used in the transient'))
    3838        s += '{}\n'.format(fielddisplay(self, 'isstressbalance', 'indicates if a stressbalance solution is used in the transient'))
    3939        s += '{}\n'.format(fielddisplay(self, 'isthermal', 'indicates if a thermal solution is used in the transient'))
    4040        s += '{}\n'.format(fielddisplay(self, 'isgroundingline', 'indicates if a groundingline migration is used in the transient'))
    41         s += '{}\n'.format(fielddisplay(self, 'isgia', 'indicates if a postglacial rebound is used in the transient'))
    4241        s += '{}\n'.format(fielddisplay(self, 'isesa', 'indicates whether an elastic adjustment model is used in the transient'))
    4342        s += '{}\n'.format(fielddisplay(self, 'isdamageevolution', 'indicates whether damage evolution is used in the transient'))
    4443        s += '{}\n'.format(fielddisplay(self, 'ismovingfront', 'indicates whether a moving front capability is used in the transient'))
    … …  
    4645        s += '{}\n'.format(fielddisplay(self, 'issampling', 'indicates whether sampling is used in the transient'))
    4746        s += '{}\n'.format(fielddisplay(self, 'isslc', 'indicates if a sea level change solution is used in the transient'))
    4847        s += '{}\n'.format(fielddisplay(self, 'isoceancoupling', 'indicates whether coupling with an ocean model is used in the transient'))
    49         s += '{}\n'.format(fielddisplay(self, 'iscoupler', 'indicates whether different models are being run with need for coupling'))
    5048        s += '{}\n'.format(fielddisplay(self, 'amr_frequency', 'frequency at which mesh is refined in simulations with multiple time_steps'))
    5149        s += '{}\n'.format(fielddisplay(self, 'requested_outputs', 'list of additional outputs requested'))
    5250        return s
    … …  
    6260    def deactivateall(self):  #{{{
    6361        self.issmb = 0
    6462        self.ismasstransport = 0
     63        self.isoceantransport = 0
    6564        self.isstressbalance = 0
    6665        self.isthermal = 0
    6766        self.isgroundingline = 0
    68         self.isgia = 0
    6967        self.isesa = 0
    7068        self.isdamageevolution = 0
    7169        self.ismovingfront = 0
    … …  
    7371        self.issampling = 0
    7472        self.isslc = 0
    7573        self.isoceancoupling = 0
    76         self.iscoupler = 0
    7774        self.amr_frequency = 0
    7875
    7976        # Default output
    … …  
    8582        #full analysis: Stressbalance, Masstransport and Thermal but no groundingline migration for now
    8683        self.issmb = 1
    8784        self.ismasstransport = 1
     85        self.isoceantransport = 0
    8886        self.isstressbalance = 1
    8987        self.isthermal = 1
    9088        self.isgroundingline = 0
    91         self.isgia = 0
    9289        self.isesa = 0
    9390        self.isdamageevolution = 0
    9491        self.ismovingfront = 0
    … …  
    9693        self.issampling = 0
    9794        self.isslc = 0
    9895        self.isoceancoupling = 0
    99         self.iscoupler = 0
    10096        self.amr_frequency = 0
    10197
    10298        # Default output
    … …  
    111107
    112108        md = checkfield(md, 'fieldname', 'transient.issmb', 'numel', [1], 'values', [0, 1])
    113109        md = checkfield(md, 'fieldname', 'transient.ismasstransport', 'numel', [1], 'values', [0, 1])
     110        md = checkfield(md, 'fieldname', 'transient.isocealtransport', 'numel', [1], 'values', [0, 1])
    114111        md = checkfield(md, 'fieldname', 'transient.isstressbalance', 'numel', [1], 'values', [0, 1])
    115112        md = checkfield(md, 'fieldname', 'transient.isthermal', 'numel', [1], 'values', [0, 1])
    116113        md = checkfield(md, 'fieldname', 'transient.isgroundingline', 'numel', [1], 'values', [0, 1])
    117         md = checkfield(md, 'fieldname', 'transient.isgia', 'numel', [1], 'values', [0, 1])
    118114        md = checkfield(md, 'fieldname', 'transient.isesa', 'numel', [1], 'values', [0, 1])
    119115        md = checkfield(md, 'fieldname', 'transient.isdamageevolution', 'numel', [1], 'values', [0, 1])
    120116        md = checkfield(md, 'fieldname', 'transient.ishydrology', 'numel', [1], 'values', [0, 1])
    … …  
    122118        md = checkfield(md, 'fieldname', 'transient.ismovingfront', 'numel', [1], 'values', [0, 1])
    123119        md = checkfield(md, 'fieldname', 'transient.isslc', 'numel', [1], 'values', [0, 1])
    124120        md = checkfield(md, 'fieldname', 'transient.isoceancoupling', 'numel', [1], 'values', [0, 1])
    125         md = checkfield(md, 'fieldname', 'transient.iscoupler', 'numel', [1], 'values', [0, 1])
    126121        md = checkfield(md, 'fieldname', 'transient.amr_frequency', 'numel', [1], '>=', 0, 'NaN', 1, 'Inf', 1)
    127122
    128123        if solution != 'TransientSolution' and md.transient.iscoupling:
    … …  
    135130    def marshall(self, prefix, md, fid):  # {{{
    136131        WriteData(fid, prefix, 'object', self, 'fieldname', 'issmb', 'format', 'Boolean')
    137132        WriteData(fid, prefix, 'object', self, 'fieldname', 'ismasstransport', 'format', 'Boolean')
     133        WriteData(fid, prefix, 'object', self, 'fieldname', 'isoceantransport', 'format', 'Boolean')
    138134        WriteData(fid, prefix, 'object', self, 'fieldname', 'isstressbalance', 'format', 'Boolean')
    139135        WriteData(fid, prefix, 'object', self, 'fieldname', 'isthermal', 'format', 'Boolean')
    140136        WriteData(fid, prefix, 'object', self, 'fieldname', 'isgroundingline', 'format', 'Boolean')
    141         WriteData(fid, prefix, 'object', self, 'fieldname', 'isgia', 'format', 'Boolean')
    142137        WriteData(fid, prefix, 'object', self, 'fieldname', 'isesa', 'format', 'Boolean')
    143138        WriteData(fid, prefix, 'object', self, 'fieldname', 'isdamageevolution', 'format', 'Boolean')
    144139        WriteData(fid, prefix, 'object', self, 'fieldname', 'ishydrology', 'format', 'Boolean')
    … …  
    146141        WriteData(fid, prefix, 'object', self, 'fieldname', 'issampling', 'format', 'Boolean')
    147142        WriteData(fid, prefix, 'object', self, 'fieldname', 'isslc', 'format', 'Boolean')
    148143        WriteData(fid, prefix, 'object', self, 'fieldname', 'isoceancoupling', 'format', 'Boolean')
    149         WriteData(fid, prefix, 'object', self, 'fieldname', 'iscoupler', 'format', 'Boolean')
    150144        WriteData(fid, prefix, 'object', self, 'fieldname', 'amr_frequency', 'format', 'Integer')
    151145
    152146        # Process requested outputs

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

     
    11from datetime import datetime
    22import os
    3 import shutil
    43
    54from ismodelselfconsistent import ismodelselfconsistent
    65from loadresultsfromcluster import loadresultsfromcluster
    … …  
    2120    Solution types available comprise:
    2221    - 'Stressbalance'      or 'sb'
    2322    - 'Masstransport'      or 'mt'
     23    - 'Oceantransport'     or 'oceant'
    2424    - 'Thermal'            or 'th'
    2525    - 'Steadystate'        or 'ss'
    2626    - 'Transient'          or 'tr'
    … …  
    3131    - 'Hydrology'          or 'hy'
    3232    - 'DamageEvolution'    or 'da'
    3333    - 'Gia'                or 'gia'
     34    - 'Love'               or 'lv'
    3435    - 'Esa'                or 'esa'
    35     - 'Sealevelchange'     or 'slc'
    36     - 'Love'               or 'lv'
    3736    - 'Sampling'           or 'smp'
    3837
    3938    Extra options:
    … …  
    5453        solutionstring = 'StressbalanceSolution'
    5554    elif solutionstring.lower() == 'mt' or solutionstring.lower() == 'masstransport':
    5655        solutionstring = 'MasstransportSolution'
     56    elif solutionstring.lower() == 'oceant' or solutionstring.lower() == 'oceantransport':
     57        solutionstring = 'OceantransportSolution'
    5758    elif solutionstring.lower() == 'th' or solutionstring.lower() == 'thermal':
    5859        solutionstring = 'ThermalSolution'
    5960    elif solutionstring.lower() == 'st' or solutionstring.lower() == 'steadystate':
    … …  
    7879        solutionstring = 'LoveSolution'
    7980    elif solutionstring.lower() == 'esa':
    8081        solutionstring = 'EsaSolution'
    81     elif solutionstring.lower() == 'slc' or solutionstring.lower() == 'sealevelchange':
    82         solutionstring = 'SealevelchangeSolution'
    8382    elif solutionstring.lower() == 'smp' or solutionstring.lower() == 'sampling':
    8483        solutionstring = 'SamplingSolution'
    8584    else:
    … …  
    150149    # Wait on lock
    151150    if md.settings.waitonlock > 0:
    152151        islock = waitonlock(md)
    153         if islock == 0:  # no results to be loaded
     152        if islock == 0: # no results to be loaded
    154153            print('The results must be loaded manually with md = loadresultsfromcluster(md).')
    155         else:  # load results
     154        else: # load results
    156155            if md.verbose.solution:
    157156                print('loading results from cluster')
    158157            md = loadresultsfromcluster(md)

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

     
    33
    44import numpy as np
    55
    6 from giaivins import giaivins
     6from materials import *
    77from model import *
    88from parameterize import *
    99from setflowequation import *
    … …  
    1616md = setmask(md, '', '')
    1717md = parameterize(md, '../Par/SquareSheetConstrained.py')
    1818
    19 #GIA
    20 md.gia = giaivins()
    21 md.gia.lithosphere_thickness = 100. * np.ones(md.mesh.numberofvertices) # in km
    22 md.gia.mantle_viscosity = 1.0e21 * np.ones(md.mesh.numberofvertices) # in Pa.s
    23 md.materials.lithosphere_shear_modulus = 6.7e10 # in Pa
    24 md.materials.lithosphere_density = 3.32 # in g/cm^3
    25 md.materials.mantle_shear_modulus = 1.45e11 # in Pa
    26 md.materials.mantle_density = 3.34 # in g/cm^3
     19# GIA Ivins, 2 layer model
     20md.solidearth.settings.grdmodel = 2
    2721
    28 #Indicate what you want to compute
    29 md.gia.cross_section_shape = 1 # for square-edged x-section
     22md.materials = materials('litho','ice')
     23md.materials.numlayers = 2;
     24md.materials.radius = [10, 6271e3, 6371e3]
     25md.materials.density = [3.34e3, 3.32e3]
     26md.materials.lame_mu = [1.45e11, 6.7e10]
     27md.materials.viscosity = [1e21, 0]
     28md.initialization.sealevel = np.zeros(md.mesh.numberofvertices)
     29md.solidearth.settings.cross_section_shape = 1 # for square-edged x-section
     30md.solidearth.settings.computesealevelchange = 0 # do not compute sea level, only deformation
     31md.solidearth.requested_outputs = ['Sealevel', 'SealevelUGrd']
    3032
    31 #Define loading history (see test2001.m for the description)
    32 md.timestepping.start_time = 2400000  # 2, 400 kyr
    33 md.timestepping.final_time = 2500000  # 2, 500 kyr
    34 md.geometry.thickness = np.array([
    35     np.append(md.geometry.thickness * 0.0, 0.0),
    36     np.append(md.geometry.thickness / 2.0, 0.1),
    37     np.append(md.geometry.thickness, 0.2),
    38     np.append(md.geometry.thickness, 1.0),
    39     np.append(md.geometry.thickness, md.timestepping.start_time)
     33# Loading history
     34md.timestepping.start_time = -2400000  # 4,800 kyr :: EVALUATION TIME
     35md.timestepping.time_step = 2400000  # 2,400 kyr :: EVALUATION TIME
     36# To get rid of default final_time, make sure final_time > start_time
     37md.timestepping.final_time = 2400000  # 2,400 kyr
     38md.masstransport.spcthickness np.array([
     39    np.append(md.geometry.thickness, 0),
     40    np.append(md.geometry.thickness, 2400000)
    4041    ]).T
    4142
     43# Geometry at 0 initially
     44md.geometry.thickness = np.zeros(md.mesh.numberofvertices)
     45md.geometry.surface = np.zeros(md.mesh.numberofvertices)
     46md.geometry.base = np.zeros(md.mesh.numberofvertices)
     47
     48# Physics
     49md.transient.issmb = 0
     50md.transient.isstressbalance = 0
     51md.transient.isthermal = 0
     52md.transient.ismasstransport = 0
     53md.transient.isslc = 0
     54
    4255#Solve for GIA deflection
     56md.cluster = generic('name', gethostname(), 'np', 1)
    4357md.cluster = generic('name', gethostname(), 'np', 3)
    44 md.verbose = verbose('1111111')
    45 md = solve(md, 'Gia')
     58md.verbose = verbose('11111111111')
     59md.verbose.solver = 0
     60md = solve(md, 'Transient')
    4661
    4762#Fields and tolerances to track changes
    48 field_names = ['UGia', 'UGiaRate']
    49 field_tolerances = [1e-13, 1e-13]
    50 field_values = [md.results.GiaSolution.UGia, md.results.GiaSolution.UGiaRate]
     63field_names = ['UGrd']
     64field_tolerances = [1e-13]
     65field_values = [md.results.TransientSolution[0].SealevelUGrd]
Note: See TracBrowser for help on using the repository browser.