Changeset 25688


Ignore:
Timestamp:
10/19/20 13:31:26 (4 years ago)
Author:
jdquinn
Message:

CHG: MATLAB -> Python translations in support of tests 2005 and 2006 (still have errors to work out here); clean up; various minor improvements

Location:
issm/trunk-jpl
Files:
3 added
86 edited

Legend:

Unmodified
Added
Removed

  • issm/trunk-jpl/src

  • issm/trunk-jpl/src/c

    • Property svn:ignore

      •  

        old new  
        2323issm_ocean
        2424issm_dakota
         25issm_post

  • issm/trunk-jpl/src/c/classes/Dakota

    • Property svn:ignore

      •  

        old new  
        11.deps
         2.dirstamp

  • issm/trunk-jpl/src/c/classes/Inputs

    • Property svn:ignore set to
      .deps
      .dirstamp

  • issm/trunk-jpl/src/c/modules/FrontalForcingsx

    • Property svn:ignore set to
      .deps
      .dirstamp

  • issm/trunk-jpl/src/c/modules/KillIcebergsx

    • Property svn:ignore set to
      .dirstamp
      .deps

  • issm/trunk-jpl/src/c/modules/OceanExchangeDatax

    • Property svn:ignore set to
      .deps

  • issm/trunk-jpl/src/c/modules/QmuStatisticsx

    • Property svn:ignore set to
      .deps
      .dirstamp

  • issm/trunk-jpl/src/c/modules/QmuStatisticsx/QmuStatisticsx.cpp

    r25639 r25688  
    10011001        }
    10021002
    1003         _printf0_("Done with MeanVariance :\n");
     1003        _printf0_("Done with MeanVariance:\n");
    10041004        IssmComm::SetComm(ISSM_MPI_COMM_WORLD);
    10051005
    … …  
    11611161                }
    11621162        }
    1163         _printf0_("Done with SampleSeries :\n");
     1163        _printf0_("Done with SampleSeries:\n");
    11641164        IssmComm::SetComm(ISSM_MPI_COMM_WORLD);
    11651165

  • issm/trunk-jpl/src/c/toolkits/codipack

    • Property svn:ignore set to
      .deps

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

    r24793 r25688  
     1import numpy as np
     2
     3from checkfield import *
    14from fielddisplay import fielddisplay
    2 from checkfield import *
    35from project3d import *
    46from WriteData import *
    … …  
    68
    79class SMBcomponents(object):
    8     """
    9     SMBcomponents Class definition
     10    """SMBCOMPONENTS class definition
    1011
    11        Usage:
    12           SMBcomponents = SMBcomponents()
     12    Usage:
     13        SMBcomponents = SMBcomponents()
    1314    """
    1415
    15     def __init__(self):  # {{{
    16         self.accumulation = float('NaN')
    17         self.runoff = float('NaN')
    18         self.evaporation = float('NaN')
     16    def __init__(self, *args):  # {{{
    1917        self.isclimatology = 0
     18        self.accumulation = np.nan
     19        self.runoff = np.nan
     20        self.evaporation = np.nan
    2021        self.steps_per_step = 1
    2122        self.averaging = 0
    2223        self.requested_outputs = []
     24
     25        nargs = len(args)
     26        if nargs == 0:
     27            pass
     28        else:
     29            raise Exception('constructor not supported')
    2330    # }}}
    2431
    2532    def __repr__(self):  # {{{
    26         string = "   surface forcings parameters (SMB = accumulation-runoff-evaporation) :"
    27         string = "%s\n%s" % (string, fielddisplay(self, 'accumulation', 'accumulated snow [m/yr ice eq]'))
    28         string = "%s\n%s" % (string, fielddisplay(self, 'runoff', 'amount of ice melt lost from the ice column [m/yr ice eq]'))
    29         string = "%s\n%s" % (string, fielddisplay(self, 'evaporation', 'mount of ice lost to evaporative processes [m/yr ice eq]'))
    30         string = "%s\n%s" % (string, fielddisplay(self, 'isclimatology', 'repeat all forcings when past last forcing time (default false)'))
    31         string = "%s\n%s" % (string, fielddisplay(self, 'steps_per_step', 'number of smb steps per time step'))
    32         string = "%s\n%s" % (string, fielddisplay(self, 'averaging', 'averaging methods from short to long steps'))
    33         string = "%s\n\t\t%s" % (string, '0: Arithmetic (default)')
    34         string = "%s\n\t\t%s" % (string, '1: Geometric')
    35         string = "%s\n\t\t%s" % (string, '2: Harmonic')
    36         string = "%s\n%s" % (string, fielddisplay(self, 'requested_outputs', 'additional outputs requested'))
    37         return string
     33        s = '   surface forcings parameters (SMB = accumulation-runoff-evaporation):\n'
     34        s += '{}\n'.format(fielddisplay(self, 'accumulation', 'accumulated snow [m/yr ice eq]'))
     35        s += '{}\n'.format(fielddisplay(self, 'runoff', 'amount of ice melt lost from the ice column [m/yr ice eq]'))
     36        s += '{}\n'.format(fielddisplay(self, 'evaporation', 'mount of ice lost to evaporative processes [m/yr ice eq]'))
     37        s += '{}\n'.format(fielddisplay(self, 'isclimatology', 'repeat all forcings when past last forcing time (default false)'))
     38        s += '{}\n'.format(fielddisplay(self, 'steps_per_step', 'number of smb steps per time step'))
     39        s += '{}\n'.format(fielddisplay(self, 'averaging', 'averaging methods from short to long steps'))
     40        s += '\t\t{}\n'.format('0: Arithmetic (default)')
     41        s += '\t\t{}\n'.format('1: Geometric')
     42        s += '\t\t{}\n'.format('2: Harmonic')
     43        s += '{}\n'.format(fielddisplay(self, 'requested_outputs', 'additional outputs requested'))
     44        return s
    3845    # }}}
    3946
    … …  
    5360            self.accumulation = np.zeros((md.mesh.numberofvertices))
    5461            print("      no SMB.accumulation specified: values set as zero")
    55 
    5662        if np.all(np.isnan(self.runoff)):
    5763            self.runoff = np.zeros((md.mesh.numberofvertices))
    … …  
    6167            self.evaporation = np.zeros((md.mesh.numberofvertices))
    6268            print("      no SMB.evaporation specified: values set as zero")
    63 
    6469        return self
    6570    # }}}
    … …  
    7479            md = checkfield(md, 'fieldname', 'smb.runoff', 'size', [md.mesh.numberofvertices], 'NaN', 1, 'Inf', 1)
    7580            md = checkfield(md, 'fieldname', 'smb.evaporation', 'size', [md.mesh.numberofvertices], 'NaN', 1, 'Inf', 1)
    76 
    7781        md = checkfield(md, 'fieldname', 'smb.steps_per_step', '>=', 1, 'numel', [1])
    7882        md = checkfield(md, 'fieldname', 'smb.averaging', 'numel', [1], 'values', [0, 1, 2])
    7983        md = checkfield(md, 'fieldname', 'masstransport.requested_outputs', 'stringrow', 1)
    8084        md = checkfield(md, 'fieldname', 'smb.isclimatology', 'values', [0, 1])
    81 
     85        if self.isclimatology:
     86            md = checkfield(md, 'fieldname', 'smb.accumulation', 'size', [md.mesh.numberofvertices + 1],
     87                'message', 'accumulation must have md.mesh.numberofvertices+1 rows in order to force a climatology')
     88            md = checkfield(md, 'fieldname', 'smb.runoff', 'size', [md.mesh.numberofvertices + 1], 'message', 'runoff must have md.mesh.numberofvertices+1 rows in order to force a climatology')
     89            md = checkfield(md, 'fieldname', 'smb.evaporation', 'size', [md.mesh.numberofvertices + 1], 'message', 'evaporation must have md.mesh.numberofvertices+1 rows in order to force a climatology')
    8290        return md
    8391    # }}}

  • issm/trunk-jpl/src/m/classes/SMBforcing.m

    r24806 r25688  
    3232                end % }}}
    3333                function list = defaultoutputs(self,md) % {{{
    34 
    3534                        list = {''};
    36 
    3735                end % }}}
    3836                function self = extrude(self,md) % {{{
    39 
    4037                        self.mass_balance=project3d(md,'vector',self.mass_balance,'type','node');
    41 
    4238                end % }}}
    4339                function self = initialize(self,md) % {{{
    … …  
    5046                end % }}}
    5147                function md = checkconsistency(self,md,solution,analyses) % {{{
    52 
    5348                        if (strcmp(solution,'TransientSolution') & md.transient.issmb == 0), return; end
    54 
    5549                        if ismember('MasstransportAnalysis',analyses),
    5650                                md = checkfield(md,'fieldname','smb.mass_balance','timeseries',1,'NaN',1,'Inf',1);

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

    r24793 r25688  
    11import numpy as np
     2
     3from checkfield import checkfield
    24from fielddisplay import fielddisplay
    3 from checkfield import checkfield
     5from project3d import project3d
    46from WriteData import WriteData
    5 from project3d import project3d
    67
    78
    89class SMBforcing(object):
    9     """
    10     SMBforcing Class definition
     10    """SMBFORCING class definition
    1111
    12        Usage:
    13           SMB = SMBforcing()
     12    Usage:
     13        SMB = SMBforcing()
    1414    """
    1515
    16     def __init__(self):  # {{{
    17         self.mass_balance = float('NaN')
     16    def __init__(self, *args):  # {{{
     17        self.isclimatology = 0
     18        self.mass_balance = np.nan
     19        self.steps_per_step = 1
    1820        self.requested_outputs = []
    19         self.isclimatology = 0
    20         self.steps_per_step = 1
    2121        self.averaging = 0
     22
     23        nargs = len(args)
     24        if nargs == 0:
     25            pass
     26        else:
     27            raise Exception('constructor not supported')
    2228    #}}}
    2329
    2430    def __repr__(self):  # {{{
    25         string = "   surface forcings parameters:"
    26         string = "%s\n%s" % (string, fielddisplay(self, 'mass_balance', 'surface mass balance [m/yr ice eq]'))
    27         string = "%s\n%s" % (string, fielddisplay(self, 'isclimatology', 'repeat all forcings when past last forcing time (default false)'))
    28         string = "%s\n%s" % (string, fielddisplay(self, 'steps_per_step', 'number of smb steps per time step'))
    29         string = "%s\n%s" % (string, fielddisplay(self, 'averaging', 'averaging methods from short to long steps'))
    30         string = "%s\n\t\t%s" % (string, '0: Arithmetic (default)')
    31         string = "%s\n\t\t%s" % (string, '1: Geometric')
    32         string = "%s\n\t\t%s" % (string, '2: Harmonic')
    33         string = "%s\n%s" % (string, fielddisplay(self, 'requested_outputs', 'additional outputs requested'))
    34         return string
     31        s = '   surface forcings parameters:\n'
     32        s += '{}\n'.format(fielddisplay(self, 'mass_balance', 'surface mass balance [m/yr ice eq]'))
     33        s += '{}\n'.format(fielddisplay(self, 'isclimatology', 'repeat all forcings when past last forcing time (default false)'))
     34        s += '{}\n'.format(fielddisplay(self, 'steps_per_step', 'number of smb steps per time step'))
     35        s += '{}\n'.format(fielddisplay(self, 'averaging', 'averaging methods from short to long steps'))
     36        s += '\t\t{}\n'.format('0: Arithmetic (default)')
     37        s += '\t\t{}\n'.format('1: Geometric')
     38        s += '\t\t{}\n'.format('2: Harmonic')
     39        s += '{}\n'.format(fielddisplay(self, 'requested_outputs', 'additional outputs requested'))
     40        return s
    3541    #}}}
    3642
    3743    def extrude(self, md):  # {{{
    38 
    3944        self.mass_balance = project3d(md, 'vector', self.mass_balance, 'type', 'node')
    4045        return self
    … …  
    4954            self.mass_balance = np.zeros((md.mesh.numberofvertices))
    5055            print("      no SMBforcing.mass_balance specified: values set as zero")
    51 
    5256        return self
    5357    #}}}
    5458
    5559    def checkconsistency(self, md, solution, analyses):  # {{{
     60        if solution == 'TransientSolution' and not md.transient.issmb:
     61            return
    5662        if 'MasstransportAnalysis' in analyses:
    5763            md = checkfield(md, 'fieldname', 'smb.mass_balance', 'timeseries', 1, 'NaN', 1, 'Inf', 1)
    58 
    5964        if 'BalancethicknessAnalysis' in analyses:
    6065            md = checkfield(md, 'fieldname', 'smb.mass_balance', 'size', [md.mesh.numberofvertices], 'NaN', 1, 'Inf', 1)
    61 
    6266        md = checkfield(md, 'fieldname', 'smb.steps_per_step', '>=', 1, 'numel', [1])
    6367        md = checkfield(md, 'fieldname', 'smb.averaging', 'numel', [1], 'values', [0, 1, 2])
    6468        md = checkfield(md, 'fieldname', 'masstransport.requested_outputs', 'stringrow', 1)
    6569        md = checkfield(md, 'fieldname', 'smb.isclimatology', 'values', [0, 1])
    66         if (self.isclimatology > 0):
     70        if self.isclimatology:
    6771            md = checkfield(md, 'fieldname', 'smb.mass_balance', 'size', [md.mesh.numberofvertices + 1], 'message', 'mass_balance must have md.mesh.numberofvertices + 1 rows in order to force a climatology')
    68 
    6972        return md
    7073    # }}}
    … …  
    8588        WriteData(fid, prefix, 'data', outputs, 'name', 'md.smb.requested_outputs', 'format', 'StringArray')
    8689        WriteData(fid, prefix, 'object', self, 'class', 'smb', 'fieldname', 'isclimatology', 'format', 'Boolean')
    87 
    8890    # }}}

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

    r24793 r25688  
    11import numpy as np
     2
     3from checkfield import checkfield
    24from fielddisplay import fielddisplay
    3 from checkfield import checkfield
     5from helpers import *
     6from MatlabFuncs import *
     7from project3d import project3d
    48from WriteData import WriteData
    5 from project3d import project3d
    6 from MatlabFuncs import *
    7 from helpers import *
    89
    910
    1011class SMBpddSicopolis(object):
    11     """
    12     SMBpddSicopolis Class definition
     12    """SMBPDDSICOPOLIS class definition
    1313
    1414    Usage:
    1515        SMBpddSicopolis = SMBpddSicopolis()
    16 """
     16    """
    1717
    18     def __init__(self):  # {{{
    19         self.precipitation = float('NaN')
    20         self.monthlytemperatures = float('NaN')
    21         self.temperature_anomaly = float('NaN')
    22         self.precipitation_anomaly = float('NaN')
    23         self.smb_corr = float('NaN')
     18    def __init__(self, *args):  # {{{
     19        self.precipitation = np.nan
     20        self.monthlytemperatures = np.nan
     21        self.temperature_anomaly = np.nan
     22        self.precipitation_anomaly = np.nan
     23        self.smb_corr = np.nan
    2424        self.desfac = 0
    25         self.s0p = float('NaN')
    26         self.s0t = float('NaN')
     25        self.s0p = np.nan
     26        self.s0t = np.nan
    2727        self.rlaps = 0
    2828        self.isfirnwarming = 0
    … …  
    3131        self.requested_outputs = []
    3232
    33         self.setdefaultparameters()
     33        nargs = len(args)
     34        if nargs == 0:
     35            self.setdefaultparameters()
     36        else:
     37            raise Exception('constructor not supported')
    3438    # }}}
    3539
    3640    def __repr__(self):  # {{{
    37         string = '   surface forcings parameters:'
    38         string += '\n   SICOPOLIS PDD scheme (Calov & Greve, 2005) :'
    39 
    40         string = "%s\n%s" % (string, fielddisplay(self, 'monthlytemperatures', 'monthly surface temperatures [K]'))
    41         string = "%s\n%s" % (string, fielddisplay(self, 'precipitation', 'monthly surface precipitation [m/yr water eq]'))
    42         string = "%s\n%s" % (string, fielddisplay(self, 'temperature_anomaly', 'anomaly to monthly reference temperature (additive [K])'))
    43         string = "%s\n%s" % (string, fielddisplay(self, 'precipitation_anomaly', 'anomaly to monthly precipitation (multiplicative, e.g. q = q0*exp(0.070458*DeltaT) after Huybrechts (2002)) [no unit])'))
    44         string = "%s\n%s" % (string, fielddisplay(self, 'smb_corr', 'correction of smb after PDD call [m/a]'))
    45         string = "%s\n%s" % (string, fielddisplay(self, 's0p', 'should be set to elevation from precip source (between 0 and a few 1000s m, default is 0) [m]'))
    46         string = "%s\n%s" % (string, fielddisplay(self, 's0t', 'should be set to elevation from temperature source (between 0 and a few 1000s m, default is 0) [m]'))
    47         string = "%s\n%s" % (string, fielddisplay(self, 'rlaps', 'present day lapse rate (default is 7.4 degree/km)'))
    48         string = "%s\n%s" % (string, fielddisplay(self, 'desfac', 'desertification elevation factor (default is -log(2.0)/1000)'))
    49         string = "%s\n%s" % (string, fielddisplay(self, 'isfirnwarming', 'is firnwarming (Reeh 1991) activated (0 or 1, default is 1)'))
    50         string = "%s\n%s" % (string, fielddisplay(self, 'steps_per_step', 'number of smb steps per time step'))
    51         string = "%s\n%s" % (string, fielddisplay(self, 'averaging', 'averaging methods from short to long steps'))
    52         string = "%s\n\t\t%s" % (string, '0: Arithmetic (default)')
    53         string = "%s\n\t\t%s" % (string, '1: Geometric')
    54         string = "%s\n\t\t%s" % (string, '2: Harmonic')
    55 
    56         string = "%s\n%s" % (string, fielddisplay(self, 'requested_outputs', 'additional outputs requested (TemperaturePDD, SmbAccumulation, SmbMelt)'))
    57     # }}}
    58 
    59     @staticmethod
    60     def SMBpddSicopolis(*args):  # {{{
    61         nargin = len(args)
    62 
    63         if nargin == 0:
    64             return SMBpddSicopolis()
    65         else:
    66             raise RuntimeError('SMBpddSicopolis: constructor not supported')
     41        s = '   surface forcings parameters:\n'
     42        s += '   SICOPOLIS PDD scheme (Calov & Greve, 2005):\n'
     43        s += '{}\n'.format(fielddisplay(self, 'monthlytemperatures', 'monthly surface temperatures [K]'))
     44        s += '{}\n'.format(fielddisplay(self, 'precipitation', 'monthly surface precipitation [m/yr water eq]'))
     45        s += '{}\n'.format(fielddisplay(self, 'temperature_anomaly', 'anomaly to monthly reference temperature (additive [K])'))
     46        s += '{}\n'.format(fielddisplay(self, 'precipitation_anomaly', 'anomaly to monthly precipitation (multiplicative, e.g. q = q0*exp(0.070458*DeltaT) after Huybrechts (2002)) [no unit])'))
     47        s += '{}\n'.format(fielddisplay(self, 'smb_corr', 'correction of smb after PDD call [m/a]'))
     48        s += '{}\n'.format(fielddisplay(self, 's0p', 'should be set to elevation from precip source (between 0 and a few 1000s m, default is 0) [m]'))
     49        s += '{}\n'.format(fielddisplay(self, 's0t', 'should be set to elevation from temperature source (between 0 and a few 1000s m, default is 0) [m]'))
     50        s += '{}\n'.format(fielddisplay(self, 'rlaps', 'present day lapse rate (default is 7.4 degree/km)'))
     51        s += '{}\n'.format(fielddisplay(self, 'desfac', 'desertification elevation factor (default is -log(2.0)/1000)'))
     52        s += '{}\n'.format(fielddisplay(self, 'isfirnwarming', 'is firnwarming (Reeh 1991) activated (0 or 1, default is 1)'))
     53        s += '{}\n'.format(fielddisplay(self, 'steps_per_step', 'number of smb steps per time step'))
     54        s += '{}\n'.format(fielddisplay(self, 'averaging', 'averaging methods from short to long steps'))
     55        s += '\t\t{}\n'.format('0: Arithmetic (default)')
     56        s += '\t\t{}\n'.format('1: Geometric')
     57        s += '\t\t{}\n'.format('2: Harmonic')
     58        s += '{}\n'.format(fielddisplay(self, 'requested_outputs', 'additional outputs requested (TemperaturePDD, SmbAccumulation, SmbMelt)'))
     59        return s
    6760    # }}}
    6861
    … …  
    8376
    8477    def initialize(self, md):  # {{{
    85 
    8678        if np.isnan(self.s0p):
    8779            self.s0p = np.zeros((md.mesh.numberofvertices, ))
    … …  
    10395            self.smb_corr = np.zeros((md.mesh.numberofvertices, ))
    10496            print('      no SMBpddSicopolis.smb_corr specified: values set as zero')
     97        return self
    10598    # }}}
    10699
    … …  
    109102        self.desfac = -np.log(2.0) / 1000
    110103        self.rlaps = 7.4
    111 
     104        return self
    112105    # }}}
    113106
    114107    def checkconsistency(self, md, solution, analyses):  # {{{
    115         if (strcmp(solution, 'TransientSolution') and md.transient.issmb == 0):
     108        if solution == 'TransientSolution' and not md.transient.issmb:
    116109            return
    117 
    118110        if 'MasstransportAnalysis' in analyses:
    119111            md = checkfield(md, 'fieldname', 'smb.desfac', '<=', 1, 'numel', 1)
    … …  
    123115            md = checkfield(md, 'fieldname', 'smb.monthlytemperatures', 'NaN', 1, 'Inf', 1, 'size', [md.mesh.numberofvertices, 12])
    124116            md = checkfield(md, 'fieldname', 'smb.precipitation', 'NaN', 1, 'Inf', 1, 'size', [md.mesh.numberofvertices, 12])
    125 
    126117        md = checkfield(md, 'fieldname', 'smb.steps_per_step', '>=', 1, 'numel', [1])
    127118        md = checkfield(md, 'fieldname', 'smb.averaging', 'numel', [1], 'values', [0, 1, 2])
    128119        md = checkfield(md, 'fieldname', 'smb.requested_outputs', 'stringrow', 1)
    129 
    130120        return md
    131121    # }}}
    132122
    133123    def marshall(self, prefix, md, fid):  # {{{
    134 
    135124        yts = md.constants.yts
    136 
    137125        WriteData(fid, prefix, 'name', 'md.smb.model', 'data', 10, 'format', 'Integer')
    138 
    139126        WriteData(fid, prefix, 'object', self, 'class', 'smb', 'fieldname', 'isfirnwarming', 'format', 'Boolean')
    140127        WriteData(fid, prefix, 'object', self, 'class', 'smb', 'fieldname', 'desfac', 'format', 'Double')
    … …  
    142129        WriteData(fid, prefix, 'object', self, 'class', 'smb', 'fieldname', 's0t', 'format', 'DoubleMat', 'mattype', 1)
    143130        WriteData(fid, prefix, 'object', self, 'class', 'smb', 'fieldname', 'rlaps', 'format', 'Double')
    144 
    145131        WriteData(fid, prefix, 'object', self, 'class', 'smb', 'fieldname', 'monthlytemperatures', 'format', 'DoubleMat', 'mattype', 1, 'timeserieslength', md.mesh.numberofvertices + 1, 'yts', md.constants.yts)
    146132        WriteData(fid, prefix, 'object', self, 'class', 'smb', 'fieldname', 'precipitation', 'format', 'DoubleMat', 'mattype', 1, 'scale', 1. / yts, 'timeserieslength', md.mesh.numberofvertices + 1, 'yts', yts)
    … …  
    159145
    160146        WriteData(fid, prefix, 'data', outputs, 'name', 'md.smb.requested_outputs', 'format', 'StringArray')
    161 
    162147    # }}}

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

    r24213 r25688  
     1import numpy as np
     2
     3from checkfield import checkfield
    14from fielddisplay import fielddisplay
    25from project3d import project3d
    3 from checkfield import checkfield
    46from WriteData import WriteData
    5 import numpy as np
    67
    78
    89class basalforcings(object):
    9     """
    10     BASAL FORCINGS class definition
     10    """BASAL FORCINGS class definition
    1111
    12        Usage:
    13           basalforcings = basalforcings()
     12    Usage:
     13        basalforcings = basalforcings()
    1414    """
    1515
    1616    def __init__(self):  # {{{
    17         self.groundedice_melting_rate = float('NaN')
    18         self.floatingice_melting_rate = float('NaN')
    19         self.geothermalflux = float('NaN')
     17        self.groundedice_melting_rate = np.nan
     18        self.floatingice_melting_rate = np.nan
     19        self.geothermalflux = np.nan
    2020
    21     #set defaults
    2221        self.setdefaultparameters()
    23 
    2422    #}}}
    2523
    2624    def __repr__(self):  # {{{
    27         string = "   basal forcings parameters:"
    28 
    29         string = "%s\n%s" % (string, fielddisplay(self, "groundedice_melting_rate", "basal melting rate (positive if melting) [m / yr]"))
    30         string = "%s\n%s" % (string, fielddisplay(self, "floatingice_melting_rate", "basal melting rate (positive if melting) [m / yr]"))
    31         string = "%s\n%s" % (string, fielddisplay(self, "geothermalflux", "geothermal heat flux [W / m^2]"))
    32         return string
     25        s = '   basal forcings parameters:\n'
     26        s += '{}\n'.format(fielddisplay(self, "groundedice_melting_rate", "basal melting rate (positive if melting) [m / yr]"))
     27        s += '{}\n'.format(fielddisplay(self, "floatingice_melting_rate", "basal melting rate (positive if melting) [m / yr]"))
     28        s += '{}\n'.format(fielddisplay(self, "geothermalflux", "geothermal heat flux [W / m^2]"))
     29        return s
    3330    #}}}
    3431
    … …  
    3633        self.groundedice_melting_rate = project3d(md, 'vector', self.groundedice_melting_rate, 'type', 'node', 'layer', 1)
    3734        self.floatingice_melting_rate = project3d(md, 'vector', self.floatingice_melting_rate, 'type', 'node', 'layer', 1)
    38         self.geothermalflux = project3d(md, 'vector', self.geothermalflux, 'type', 'node', 'layer', 1)  #bedrock only gets geothermal flux
     35        self.geothermalflux = project3d(md, 'vector', self.geothermalflux, 'type', 'node', 'layer', 1) # Bedrock only gets geothermal flux
    3936        return self
    4037    #}}}
    … …  
    4239    def initialize(self, md):  # {{{
    4340        if np.all(np.isnan(self.groundedice_melting_rate)):
    44             self.groundedice_melting_rate = np.zeros((md.mesh.numberofvertices))
     41            self.groundedice_melting_rate = np.zeros((md.mesh.numberofvertices, ))
    4542            print("      no basalforcings.groundedice_melting_rate specified: values set as zero")
    46 
    4743        if np.all(np.isnan(self.floatingice_melting_rate)):
    48             self.floatingice_melting_rate = np.zeros((md.mesh.numberofvertices))
     44            self.floatingice_melting_rate = np.zeros((md.mesh.numberofvertices, ))
    4945            print("      no basalforcings.floatingice_melting_rate specified: values set as zero")
    50     #if np.all(np.isnan(self.geothermalflux)):
    51     #self.geothermalflux = np.zeros((md.mesh.numberofvertices))
    52     #print "      no basalforcings.geothermalflux specified: values set as zero"
    53 
    5446        return self
    5547    #}}}
    … …  
    6052
    6153    def checkconsistency(self, md, solution, analyses):  # {{{
    62 
    63         if 'MasstransportAnalysis' in analyses and not (solution == 'TransientSolution' and not md.transient.ismasstransport):
     54        if 'MasstransportAnalysis' in analyses and not solution == 'TransientSolution' and not md.transient.ismasstransport:
    6455            md = checkfield(md, 'fieldname', 'basalforcings.groundedice_melting_rate', 'NaN', 1, 'Inf', 1, 'timeseries', 1)
    6556            md = checkfield(md, 'fieldname', 'basalforcings.floatingice_melting_rate', 'NaN', 1, 'Inf', 1, 'timeseries', 1)
    66 
    6757        if 'BalancethicknessAnalysis' in analyses:
    6858            md = checkfield(md, 'fieldname', 'basalforcings.groundedice_melting_rate', 'NaN', 1, 'Inf', 1, 'size', [md.mesh.numberofvertices])
    6959            md = checkfield(md, 'fieldname', 'basalforcings.floatingice_melting_rate', 'NaN', 1, 'Inf', 1, 'size', [md.mesh.numberofvertices])
    70 
    71         if 'ThermalAnalysis' in analyses and not (solution == 'TransientSolution' and not md.transient.isthermal):
     60        if 'ThermalAnalysis' in analyses and not solution == 'TransientSolution' and not md.transient.isthermal:
    7261            md = checkfield(md, 'fieldname', 'basalforcings.groundedice_melting_rate', 'NaN', 1, 'Inf', 1, 'timeseries', 1)
    7362            md = checkfield(md, 'fieldname', 'basalforcings.floatingice_melting_rate', 'NaN', 1, 'Inf', 1, 'timeseries', 1)
    7463            md = checkfield(md, 'fieldname', 'basalforcings.geothermalflux', 'NaN', 1, 'Inf', 1, 'timeseries', 1, '>=', 0)
    75 
    7664        return md
    7765    # }}}
    7866
    7967    def marshall(self, prefix, md, fid):  # {{{
    80 
    8168        yts = md.constants.yts
    82 
    8369        WriteData(fid, prefix, 'name', 'md.basalforcings.model', 'data', 1, 'format', 'Integer')
    8470        WriteData(fid, prefix, 'object', self, 'fieldname', 'groundedice_melting_rate', 'format', 'DoubleMat', 'mattype', 1, 'scale', 1. / yts, 'timeserieslength', md.mesh.numberofvertices + 1, 'yts', md.constants.yts)

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

    r25169 r25688  
    88
    99class dsl(object):
    10     '''
    11     DSL - slass definition
     10    """DSL - slass definition
    1211
    13         Usage:
    14             dsl = dsl()
    15     '''
     12    Usage:
     13        dsl = dsl()
     14    """
    1615
    1716    def __init__(self): #{{{
    … …  
    2322
    2423    def __repr__(self): #{{{
    25         s = "   dsl parameters:"
    26         s = "%s\n%s" % (s, fielddisplay(self, 'global_average_thermosteric_sea_level_change', 'corresponds to zostoga field in CMIP5 archives. Specified as a temporally variable global rate (mm/yr)'))
    27         s = "%s\n%s" % (s, 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)'))
    28         s = "%s\n%s" % (s, 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)'))
    29         s = "%s\n%s" % (s, fielddisplay(self, 'compute_fingerprints', 'do we use the sea water pressure change to compute fingerprints and correct sea_surface_height_change_above_geoid'))
    30 
     24        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'))
    3129        return s
    3230    #}}}
    … …  
    4341
    4442    def checkconsistency(self, md, solution, analyses): #{{{
    45         # Early retun
     43        # Early return
    4644        if 'SealevelriseAnalysis' not in analyses:
    4745            return md
    48 
    4946        if solution == 'TransientSolution' and md.transient.isslr == 0:
    5047            return md
    51 
    5248        md = checkfield(md, 'fieldname', 'dsl.global_average_thermosteric_sea_level_change', 'NaN', 1, 'Inf', 1)
    5349        md = checkfield(md, 'fieldname', 'dsl.sea_surface_height_change_above_geoid', 'NaN', 1, 'Inf', 1, 'timeseries', 1)
    5450        md = checkfield(md, 'fieldname', 'dsl.sea_water_pressure_change_at_sea_floor', 'NaN', 1, 'Inf', 1, 'timeseries', 1)
    5551        md = checkfield(md, 'fieldname', 'dsl.compute_fingerprints', 'NaN', 1, 'Inf', 1, 'values', [0, 1])
    56 
    5752        if self.compute_fingerprints:
    58             #check geodetic flag of slr is on:
    59             if md.slr.geodetic == 0:
     53            # Check if geodetic flag of slr is on
     54            if not md.slr.geodetic:
    6055                raise RuntimeError('DSL checkconsistency error message: if bottom pressure fingerprints computations are requested, slr class should have geodetic flag on')
    6156        return md

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

    r25343 r25688  
    88
    99class dslmme(object):
    10     '''
    11     DSLMME class definition
     10    """DSLMME class definition
    1211
    13         Usage:
    14             dsl = dslmme() #dynamic sea level class based on a multi-model ensemble of CMIP5 outputs
    15     '''
     12    Usage:
     13        dsl = dslmme() #dynamic sea level class based on a multi-model ensemble of CMIP5 outputs
     14    """
    1615
    1716    def __init__(self, *args): #{{{
    … …  
    3736        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.'))
    3837        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'))
    39 
    4038        return s
    4139    #}}}
    4240
    43     def setdefaultparameters(self): # {{{
    44         return
     41    def setdefaultparameters(self): #{{{
     42        return self
    4543    #}}}
    4644
    4745    def checkconsistency(self, md, solution, analyses): # {{{
    48         if ('SealevelriseAnalysis' not in analyses) or (solution == 'TransientSolution' and md.transient.isslr == 0):
     46        if ('SealevelriseAnalysis' not in analyses) or (solution == 'TransientSolution' and not md.transient.isslr):
    4947            return md
    50 
    5148        for i in range(len(self.global_average_thermosteric_sea_level_change)):
    5249            md = checkfield(md, 'field', self.global_average_thermosteric_sea_level_change[i], 'NaN', 1, 'Inf', 1)
    5350            md = checkfield(md, 'field', self.sea_surface_height_change_above_geoid[i], 'NaN', 1, 'Inf', 1, 'timeseries', 1)
    5451            md = checkfield(md, 'field', self.sea_water_pressure_change_at_sea_floor[i], 'NaN', 1, 'Inf', 1, 'timeseries', 1)
    55 
    5652        md = checkfield(md, 'field', self.modelid, 'NaN', 1, 'Inf', 1, '>=', 1, '<=',len(self.global_average_thermosteric_sea_level_change))
    57        
    5853        if self.compute_fingerprints:
    5954            #check geodetic flag of slr is on
    60             if md.solidearth.settings.computesealevelchange == 0,
     55            if not md.solidearth.settings.computesealevelchange:
    6156                raise Exception('DSL checkconsistency error message: if bottom pressure fingerprints computations are requested, slr class should have geodetic flag on')
    62 
    6357        return md
    6458    #}}}

  • issm/trunk-jpl/src/m/classes/flowequation.m

    r25661 r25688  
    1313                isFS                           = 0;
    1414                isNitscheBC                    = 0;
    15       FSNitscheGamma                 = 1e6;
     15                FSNitscheGamma                 = 1e6;
    1616                fe_SSA                         = '';
    1717                fe_HO                          = '';

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

    r25626 r25688  
    2323        self.isFS = 0
    2424        self.isNitscheBC = 0
    25         self.FSNitscheGamma = 1e-6
     25        self.FSNitscheGamma = 1e6
    2626        self.fe_SSA = ''
    2727        self.fe_HO = ''
    … …  
    3232        self.augmented_lagrangian_rholambda = 1.
    3333        self.XTH_theta = 0.
    34         self.vertex_equation = float('NaN')
    35         self.element_equation = float('NaN')
    36         self.borderSSA = float('NaN')
    37         self.borderHO = float('NaN')
    38         self.borderFS = float('NaN')
    39 
    40     #set defaults
     34        self.vertex_equation = np.nan
     35        self.element_equation = np.nan
     36        self.borderSSA = np.nan
     37        self.borderHO = np.nan
     38        self.borderFS = np.nan
     39       
    4140        self.setdefaultparameters()
    42 
    4341    #}}}
    4442
    4543    def __repr__(self):  # {{{
    4644        s = '   flow equation parameters:\n'
    47         s += "{}\n".format(fielddisplay(self, 'isSIA', "is the Shallow Ice Approximation (SIA) used?"))
    48         s += "{}\n".format(fielddisplay(self, 'isSSA', "is the Shelfy-Stream Approximation (SSA) used?"))
    49         s += "{}\n".format(fielddisplay(self, 'isL1L2', "are L1L2 equations used?"))
    50         s += "{}\n".format(fielddisplay(self, 'isMLHO', "are Mono-layer Higher-Order equations used?"))
    51         s += "{}\n".format(fielddisplay(self, 'isHO', "is the Higher-Order (HO) approximation used?"))
    52         s += "{}\n".format(fielddisplay(self, 'isFS', "are the Full-FS (FS) equations used?"))
    53         s += "{}\n".format(fielddisplay(self, 'isNitscheBC', "is weakly imposed condition used?"))
    54         s += "{}\n".format(fielddisplay(self, 'FSNitscheGamma', "Gamma value for the Nitsche term (default: 1e6)"))
    55         s += "{}\n".format(fielddisplay(self, 'fe_SSA', "Finite Element for SSA: 'P1', 'P1bubble' 'P1bubblecondensed' 'P2'"))
    56         s += "{}\n".format(fielddisplay(self, 'fe_HO', "Finite Element for HO:  'P1', 'P1bubble', 'P1bubblecondensed', 'P1xP2', 'P2xP1', 'P2', 'P2bubble', 'P1xP3', 'P2xP4'"))
    57         s += "{}\n".format(fielddisplay(self, 'fe_FS', "Finite Element for FS:  'P1P1' (debugging only) 'P1P1GLS' 'MINIcondensed' 'MINI' 'TaylorHood' 'LATaylorHood' 'XTaylorHood'"))
    58         s += "{}\n".format(fielddisplay(self, 'vertex_equation', "flow equation for each vertex"))
    59         s += "{}\n".format(fielddisplay(self, 'element_equation', "flow equation for each element"))
    60         s += "{}\n".format(fielddisplay(self, 'borderSSA', "vertices on SSA's border (for tiling)"))
    61         s += "{}\n".format(fielddisplay(self, 'borderHO', "vertices on HO's border (for tiling)"))
    62         s += "{}".format(fielddisplay(self, 'borderFS', "vertices on FS' border (for tiling)"))
     45        s += '{}\n'.format(fielddisplay(self, 'isSIA', "is the Shallow Ice Approximation (SIA) used?"))
     46        s += '{}\n'.format(fielddisplay(self, 'isSSA', "is the Shelfy-Stream Approximation (SSA) used?"))
     47        s += '{}\n'.format(fielddisplay(self, 'isL1L2', "are L1L2 equations used?"))
     48        s += '{}\n'.format(fielddisplay(self, 'isMLHO', "are Mono-layer Higher-Order equations used?"))
     49        s += '{}\n'.format(fielddisplay(self, 'isHO', "is the Higher-Order (HO) approximation used?"))
     50        s += '{}\n'.format(fielddisplay(self, 'isFS', "are the Full-FS (FS) equations used?"))
     51        s += '{}\n'.format(fielddisplay(self, 'isNitscheBC', "is weakly imposed condition used?"))
     52        s += '{}\n'.format(fielddisplay(self, 'FSNitscheGamma', "Gamma value for the Nitsche term (default: 1e6)"))
     53        s += '{}\n'.format(fielddisplay(self, 'fe_SSA', "Finite Element for SSA: 'P1', 'P1bubble' 'P1bubblecondensed' 'P2'"))
     54        s += '{}\n'.format(fielddisplay(self, 'fe_HO', "Finite Element for HO:  'P1', 'P1bubble', 'P1bubblecondensed', 'P1xP2', 'P2xP1', 'P2', 'P2bubble', 'P1xP3', 'P2xP4'"))
     55        s += '{}\n'.format(fielddisplay(self, 'fe_FS', "Finite Element for FS:  'P1P1' (debugging only) 'P1P1GLS' 'MINIcondensed' 'MINI' 'TaylorHood' 'LATaylorHood' 'XTaylorHood'"))
     56        s += '{}\n'.format(fielddisplay(self, 'vertex_equation', "flow equation for each vertex"))
     57        s += '{}\n'.format(fielddisplay(self, 'element_equation', "flow equation for each element"))
     58        s += '{}\n'.format(fielddisplay(self, 'borderSSA', "vertices on SSA's border (for tiling)"))
     59        s += '{}\n'.format(fielddisplay(self, 'borderHO', "vertices on HO's border (for tiling)"))
     60        s += '{}\n'.format(fielddisplay(self, 'borderFS', "vertices on FS' border (for tiling)"))
     61        return s
     62    #}}}
    6363
    64         return s
     64    def setdefaultparameters(self):  # {{{
     65        # P1 for SSA
     66        self.fe_SSA = 'P1'
     67
     68        # P1 for HO
     69        self.fe_HO = 'P1'
     70
     71        # MINI condensed element for FS by default
     72        self.fe_FS = 'MINIcondensed'
     73        return self
    6574    #}}}
    6675
    … …  
    7483    #}}}
    7584
    76     def setdefaultparameters(self):  # {{{
    77         #P1 for SSA
    78         self.fe_SSA = 'P1'
    79 
    80         #P1 for HO
    81         self.fe_HO = 'P1'
    82 
    83         #MINI condensed element for FS by default
    84         self.fe_FS = 'MINIcondensed'
    85 
    86         return self
    87     #}}}
    88 
    8985    def checkconsistency(self, md, solution, analyses):  # {{{
    90 
    91         #Early return
     86        # Early return
    9287        if ('StressbalanceAnalysis' not in analyses and 'StressbalanceSIAAnalysis' not in analyses) or (solution == 'TransientSolution' and not md.transient.isstressbalance):
    9388            return md
    94 
    9589        md = checkfield(md, 'fieldname', 'flowequation.isSIA', 'numel', [1], 'values', [0, 1])
    9690        md = checkfield(md, 'fieldname', 'flowequation.isSSA', 'numel', [1], 'values', [0, 1])
    … …  
    122116            md = checkfield(md, 'fieldname', 'flowequation.element_equation', 'size', [md.mesh.numberofelements], 'values', np.arange(0, 9 + 1))
    123117        else:
    124             raise RuntimeError('mesh type not supported yet')
     118            raise RuntimeError('Case not supported yet')
    125119        if not (self.isSIA or self.isSSA or self.isL1L2 or self.isMLHO or self.isHO or self.isFS):
    126120            md.checkmessage("no element types set for this model")
    127 
    128121        if 'StressbalanceSIAAnalysis' in analyses:
    129122            if any(self.element_equation == 1):
    130123                if np.any(np.logical_and(self.vertex_equation, md.mask.ocean_levelset)):
    131124                    print("\n !!! Warning: SIA's model is not consistent on ice shelves !!!\n")
    132 
    133125        return md
    134126    # }}}
    … …  
    154146        WriteData(fid, prefix, 'object', self, 'fieldname', 'borderHO', 'format', 'DoubleMat', 'mattype', 1)
    155147        WriteData(fid, prefix, 'object', self, 'fieldname', 'borderFS', 'format', 'DoubleMat', 'mattype', 1)
    156         #convert approximations to enums
     148        # Convert approximations to enums
    157149        WriteData(fid, prefix, 'data', self.vertex_equation, 'name', 'md.flowequation.vertex_equation', 'format', 'DoubleMat', 'mattype', 1)
    158150        WriteData(fid, prefix, 'data', self.element_equation, 'name', 'md.flowequation.element_equation', 'format', 'DoubleMat', 'mattype', 2)
    159 
    160151    # }}}

  • issm/trunk-jpl/src/m/classes/friction.m

    r25499 r25688  
    6666                               
    6767                                otherwise
    68                                         error('not supported yet');             
     68                                        error('not supported yet');
    6969                        end
    7070                end % }}}

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

    r25499 r25688  
    2121        self.effective_pressure = np.nan
    2222        self.effective_pressure_limit = 0
    23 
    24         #set defaults
    2523        self.setdefaultparameters()
    2624    #}}}
    2725
    2826    def __repr__(self):  # {{{
    29         s = "Basal shear stress parameters: Sigma_b = coefficient^2 * Neff ^r * |u_b|^(s - 1) * u_b, \n(effective stress Neff = rho_ice * g * thickness + rho_water * g * base, r = q / p and s = 1 / p)\n"
    30         s = "{}\n".format(fielddisplay(self, "coefficient", "friction coefficient [SI]"))
    31         s = "{}\n".format(fielddisplay(self, "p", "p exponent"))
    32         s = "{}\n".format(fielddisplay(self, "q", "q exponent"))
    33         s = "{}\n".format(fielddisplay(self, 'coupling', 'Coupling flag 0: uniform sheet (negative pressure ok, default), 1: ice pressure only, 2: water pressure assuming uniform sheet (no negative pressure), 3: use provided effective_pressure, 4: used coupled model (not implemented yet)'))
    34         s = "{}\n".format(fielddisplay(self, 'effective_pressure', 'Effective Pressure for the forcing if not coupled [Pa]'))
    35         s = "{}\n".format(fielddisplay(self, 'effective_pressure_limit', 'Neff do not allow to fall below a certain limit: effective_pressure_limit * rho_ice * g * thickness (default 0)'))
    36        
     27        s = 'Basal shear stress parameters: Sigma_b = coefficient^2 * Neff ^r * |u_b|^(s - 1) * u_b,\n'
     28        s += '(effective stress Neff = rho_ice * g * thickness + rho_water * g * base, r = q / p and s = 1 / p)\n'
     29        s += '{}\n'.format(fielddisplay(self, "coefficient", "friction coefficient [SI]"))
     30        s += '{}\n'.format(fielddisplay(self, "p", "p exponent"))
     31        s += '{}\n'.format(fielddisplay(self, "q", "q exponent"))
     32        s += '{}\n'.format(fielddisplay(self, 'coupling', 'Coupling flag 0: uniform sheet (negative pressure ok, default), 1: ice pressure only, 2: water pressure assuming uniform sheet (no negative pressure), 3: use provided effective_pressure, 4: used coupled model (not implemented yet)'))
     33        s += '{}\n'.format(fielddisplay(self, 'effective_pressure', 'Effective Pressure for the forcing if not coupled [Pa]'))
     34        s += '{}\n'.format(fielddisplay(self, 'effective_pressure_limit', 'Neff do not allow to fall below a certain limit: effective_pressure_limit * rho_ice * g * thickness (default 0)'))
    3735        return s
    3836    #}}}
    … …  
    4038    def setdefaultparameters(self):  # {{{
    4139        self.effective_pressure_limit = 0
    42 
    4340        return self
    4441    #}}}
    … …  
    4845        self.p = project3d(md, 'vector', self.p, 'type', 'element')
    4946        self.q = project3d(md, 'vector', self.q, 'type', 'element')
    50     #if self.coupling == 0:  #doesnt work with empty loop, so just skip it?
     47        # If self.coupling == 0:  #doesnt work with empty loop, so just skip it?
    5148        if self.coupling in[3, 4]:
    5249            self.effective_pressure = project3d(md, 'vector', self.effective_pressure, 'type', 'node', 'layer', 1)
    … …  
    6259
    6360    def checkconsistency(self, md, solution, analyses):  # {{{
    64         #Early return
     61        # Early return
    6562        if 'StressbalanceAnalysis' not in analyses and 'ThermalAnalysis' not in analyses:
    6663            return md
    67         if solution == 'TransientSoluition' and not md.transient.isstressbalance and not md.transient.isthermal:
     64        if solution == 'TransientSolution' and not md.transient.isstressbalance and not md.transient.isthermal:
    6865            return md
    69 
    7066        md = checkfield(md, 'fieldname', 'friction.coefficient', 'timeseries', 1, 'NaN', 1, 'Inf', 1)
    7167        md = checkfield(md, 'fieldname', 'friction.q', 'NaN', 1, 'Inf', 1, 'size', [md.mesh.numberofelements])
    … …  
    7773        elif self.coupling > 4:
    7874            raise ValueError('not supported yet')
    79 
    8075        return md
    8176    # }}}
    … …  
    8378    def marshall(self, prefix, md, fid):  # {{{
    8479        WriteData(fid, prefix, 'name', 'md.friction.law', 'data', 1, 'format', 'Integer')
    85 
    8680        if type(self.coefficient) in [np.ndarray] and (self.coefficient.shape[0] == md.mesh.numberofvertices or self.coefficient.shape[0] == (md.mesh.numberofvertices + 1)):
    8781            mattype = 1
    … …  
    9084            mattype = 2
    9185            tsl = md.mesh.numberofelements
    92 
    9386        WriteData(fid, prefix, 'object', self, 'fieldname', 'coefficient', 'format', 'DoubleMat', 'mattype', mattype, 'timeserieslength', tsl + 1, 'yts', md.constants.yts)
    9487        WriteData(fid, prefix, 'object', self, 'fieldname', 'p', 'format', 'DoubleMat', 'mattype', 2)

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

    r25499 r25688  
    2323        self.effective_pressure_limit = 0
    2424
    25         #set defaults
    2625        self.setdefaultparameters()
    2726    #}}}
    2827
    2928    def __repr__(self):  # {{{
    30         s = "Basal shear stress parameters: Sigma_b = min(coefficient^2 * Neff ^r * |u_b|^(s - 1) * u_b, \n coefficientcoulomb^2 * rho_i * g * (h - h_f)), (effective stress Neff = rho_ice * g * thickness + rho_water * g * bed, r = q / p and s = 1 / p).\n"
    31         s = "{}\n".format(fielddisplay(self, "coefficient", "power law (Weertman) friction coefficient [SI]"))
    32         s = "{}\n".format(fielddisplay(self, "coefficientcoulomb", "Coulomb friction coefficient [SI]"))
    33         s = "{}\n".format(fielddisplay(self, "p", "p exponent"))
    34         s = "{}\n".format(fielddisplay(self, "q", "q exponent"))
    35         s = "{}\n".format(fielddisplay(self, 'coupling', 'Coupling flag: 0 for default, 1 for forcing(provide md.friction.effective_pressure)  and 2 for coupled(not implemented yet)'))
    36         s = "{}\n".format(fielddisplay(self, 'effective_pressure', 'Effective Pressure for the forcing if not coupled [Pa]'))
    37         s = "{}\n".format(fielddisplay(self, 'effective_pressure_limit', 'Neff do not allow to fall below a certain limit: effective_pressure_limit * rho_ice * g * thickness (default 0)'))
    38        
     29        s = 'Basal shear stress parameters: Sigma_b = min(coefficient^2 * Neff ^r * |u_b|^(s - 1) * u_b,\n'
     30        s += 'coefficientcoulomb^2 * rho_i * g * (h - h_f)), (effective stress Neff = rho_ice * g * thickness + rho_water * g * bed, r = q / p and s = 1 / p).\n'
     31        s += '{}\n'.format(fielddisplay(self, "coefficient", "power law (Weertman) friction coefficient [SI]"))
     32        s += '{}\n'.format(fielddisplay(self, "coefficientcoulomb", "Coulomb friction coefficient [SI]"))
     33        s += '{}\n'.format(fielddisplay(self, "p", "p exponent"))
     34        s += '{}\n'.format(fielddisplay(self, "q", "q exponent"))
     35        s += '{}\n'.format(fielddisplay(self, 'coupling', 'Coupling flag: 0 for default, 1 for forcing(provide md.friction.effective_pressure)  and 2 for coupled(not implemented yet)'))
     36        s += '{}\n'.format(fielddisplay(self, 'effective_pressure', 'Effective Pressure for the forcing if not coupled [Pa]'))
     37        s += '{}\n'.format(fielddisplay(self, 'effective_pressure_limit', 'Neff do not allow to fall below a certain limit: effective_pressure_limit * rho_ice * g * thickness (default 0)'))
    3938        return s
    4039    #}}}
    … …  
    4241    def setdefaultparameters(self):  # {{{
    4342        self.effective_pressure_limit = 0
    44 
    4543        return self
    4644    #}}}
    … …  
    5755        elif self.coupling > 2:
    5856            raise ValueError('not supported yet')
    59 
    6057        return self
    6158    #}}}
    6259    def checkconsistency(self, md, solution, analyses):  # {{{
    63         #Early return
     60        # Early return
    6461        if 'StressbalanceAnalysis' not in analyses and 'ThermalAnalysis' not in analyses:
    6562            return md
    66 
    6763        md = checkfield(md, 'fieldname', 'friction.coefficient', 'timeseries', 1, 'NaN', 1, 'Inf', 1)
    6864        md = checkfield(md, 'fieldname', 'friction.coefficientcoulomb', 'timeseries', 1, 'NaN', 1, 'Inf', 1)
    … …  
    7773        elif self.coupling > 2:
    7874            raise ValueError('not supported yet')
    79 
    8075        return md
    8176    # }}}

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

    r25385 r25688  
     1import numpy as np
     2
    13from checkfield import checkfield
    24from fielddisplay import fielddisplay
    … …  
    68
    79class frictionjosh(object):
    8     '''
    9     FRICTION class definition
     10    """FRICTIONJOSH class definition
    1011
    11         Usage:
    12             frictionjosh = frictionjosh()
    13     '''
     12    Usage:
     13        frictionjosh = frictionjosh()
     14    """
    1415
    1516    def __init__(self):  # {{{
    16         self.coefficient = float('NaN')
    17         self.pressure_adjusted_temperature = float('NaN')
     17        self.coefficient = np.nan
     18        self.pressure_adjusted_temperature = np.nan
    1819        self.gamma = 0
    1920        self.effective_pressure_limit = 0
    20         #set defaults
     21
    2122        self.setdefaultparameters()
    2223        #self.requested_outputs = []
    … …  
    2425
    2526    def __repr__(self):  # {{{
    26         string = "Basal shear stress parameters: Sigma_b = coefficient^2 * Neff ^r * |u_b|^(s - 1) * u_b, \n(effective stress Neff = rho_ice * g * thickness + rho_water * g * base, r = q / p and s = 1 / p)"
    27 
    28         string = "%s\n%s" % (string, fielddisplay(self, "coefficient", "friction coefficient [SI]"))
    29         string = "%s\n%s" % (string, fielddisplay(self, 'pressure_adjusted_temperature', 'friction pressure_adjusted_temperature (T - Tpmp) [K]'))
    30         string = "%s\n%s" % (string, fielddisplay(self, 'gamma', '(T - Tpmp)/gamma [K]'))
    31         string = "%s\n%s" % (string, fielddisplay(self, 'effective_pressure_limit', 'Neff do not allow to fall below a certain limit: effective_pressure_limit * rho_ice * g * thickness (default 0)'))
    32         #string = "%s\n%s" % (string, fielddisplay(self, 'requested_outputs', 'additional outputs requested'))
    33         return string
     27        s = 'Basal shear stress parameters: Sigma_b = coefficient^2 * Neff ^r * |u_b|^(s - 1) * u_b,\n'
     28        s += '(effective stress Neff = rho_ice * g * thickness + rho_water * g * base, r = q / p and s = 1 / p)\n'
     29        s += '{}\n'.format(fielddisplay(self, "coefficient", "friction coefficient [SI]"))
     30        s += '{}\n'.format(fielddisplay(self, 'pressure_adjusted_temperature', 'friction pressure_adjusted_temperature (T - Tpmp) [K]'))
     31        s += '{}\n'.format(fielddisplay(self, 'gamma', '(T - Tpmp)/gamma [K]'))
     32        s += '{}\n'.format(fielddisplay(self, 'effective_pressure_limit', 'Neff do not allow to fall below a certain limit: effective_pressure_limit * rho_ice * g * thickness (default 0)'))
     33        #s += '{}\n'.format(fielddisplay(self, 'requested_outputs', 'additional outputs requested'))
     34        return s
    3435    #}}}
    3536
    … …  
    5354
    5455    def checkconsistency(self, md, solution, analyses):  # {{{
    55         #Early return
     56        # Early return
    5657        if 'StressbalanceAnalysis' not in analyses and 'ThermalAnalysis' not in analyses:
    5758            return md
    58 
    5959        md = checkfield(md, 'fieldname', 'friction.coefficient', 'timeseries', 1, 'NaN', 1, 'Inf', 1)
    6060        md = checkfield(md, 'fieldname', 'friction.pressure_adjusted_temperature','NaN',1,'Inf',1)
    6161        md = checkfield(md, 'fieldname', 'friction.gamma','numel',1,'NaN',1,'Inf',1,'>',0.)
    6262        md = checkfield(md, 'fieldname', 'friction.effective_pressure_limit', 'numel', [1], '>=', 0)
    63 
    64         #Check that temperature is provided
     63        # Check that temperature is provided
    6564        md = checkfield(md,'fieldname','initialization.temperature','NaN',1,'Inf',1,'size','universal')
    6665        return md
    … …  
    7372        WriteData(fid,prefix,'class','friction','object',self,'fieldname','gamma','format','Double')
    7473        WriteData(fid,prefix,'object',self,'class','friction','fieldname','effective_pressure_limit','format','Double')
    75 
    7674    # }}}

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

    r25023 r25688  
    11import numpy as np
     2
     3from checkfield import checkfield
    24from fielddisplay import fielddisplay
    35from project3d import project3d
    4 from checkfield import checkfield
    56from WriteData import WriteData
    67
    78
    89class frictionpism(object):
    9     """
    10     frictionpism class definition
     10    """FRICTIONPISM class definition
    1111
    1212    Usage:
    13       frictionpism = frictionpism()
     13        frictionpism = frictionpism()
    1414    """
     15
    1516    def __init__(self):
    1617        self.pseudoplasticity_exponent = 0.
    … …  
    2930        self.sediment_compressibility_coefficient = project3d(md, 'vector', self.sediment_compressibility_coefficient, 'type', 'node', 'layer', 1)
    3031        return self
    31         # }}}
     32    # }}}
    3233
    3334    def setdefaultparameters(self):  # {{{
    … …  
    4041
    4142    def checkconsistency(self, md, solution, analyses):  # {{{
    42 
    43         #Early return
     43        # Early return
    4444        if 'StressbalanceAnalysis' not in analyses and 'ThermalAnalysis' not in analyses:
    4545            return md
    46         if solution == 'TransientSolution' and md.transient.isstressbalance == 0 and md.transient.isthermal == 0:
     46        if solution == 'TransientSolution' and not md.transient.isstressbalance and not md.transient.isthermal:
    4747            return md
    48 
    4948        md = checkfield(md, 'fieldname', 'friction.pseudoplasticity_exponent', 'numel', [1], '>', 0, 'NaN', 1, 'Inf', 1)
    5049        md = checkfield(md, 'fieldname', 'friction.threshold_speed', 'numel', [1], '>', 0, 'NaN', 1, 'Inf', 1)
    … …  
    5352        md = checkfield(md, 'fieldname', 'friction.till_friction_angle', 'NaN', 1, 'Inf', 1, '<', 360., '>', 0., 'size', [md.mesh.numberofvertices])  #User should give angle in degrees, Matlab calculates in rad
    5453        md = checkfield(md, 'fieldname', 'friction.sediment_compressibility_coefficient', 'NaN', 1, 'Inf', 1, '<', 1., '>', 0., 'size', [md.mesh.numberofvertices])
     54        return md
    5555    # }}}
    5656
    5757    def __repr__(self):  # {{{
    58         string = 'Basal shear stress parameters for the PISM friction law (See Aschwanden et al. 2016 for more details)'
    59         string = "%s\n%s" % (string, fielddisplay(self, 'pseudoplasticity_exponent', 'pseudoplasticity exponent [dimensionless]'))
    60         string = "%s\n%s" % (string, fielddisplay(self, 'threshold_speed', 'threshold speed [m / yr]'))
    61         string = "%s\n%s" % (string, fielddisplay(self, 'delta', 'lower limit of the effective pressure, expressed as a fraction of overburden pressure [dimensionless]'))
    62         string = "%s\n%s" % (string, fielddisplay(self, 'void_ratio', 'void ratio at a reference effective pressure [dimensionless]'))
    63         string = "%s\n%s" % (string, fielddisplay(self, 'till_friction_angle', 'till friction angle [deg], recommended default: 30 deg'))
    64         string = "%s\n%s" % (string, fielddisplay(self, 'sediment_compressibility_coefficient', 'coefficient of compressibility of the sediment [dimensionless], recommended default: 0.12'))
    65         # }}}
     58        s = 'Basal shear stress parameters for the PISM friction law (See Aschwanden et al. 2016 for more details)\n'
     59        s += "{}\n".format(fielddisplay(self, 'pseudoplasticity_exponent', 'pseudoplasticity exponent [dimensionless]'))
     60        s += "{}\n".format(fielddisplay(self, 'threshold_speed', 'threshold speed [m / yr]'))
     61        s += "{}\n".format(fielddisplay(self, 'delta', 'lower limit of the effective pressure, expressed as a fraction of overburden pressure [dimensionless]'))
     62        s += "{}\n".format(fielddisplay(self, 'void_ratio', 'void ratio at a reference effective pressure [dimensionless]'))
     63        s += "{}\n".format(fielddisplay(self, 'till_friction_angle', 'till friction angle [deg], recommended default: 30 deg'))
     64        s += "{}\n".format(fielddisplay(self, 'sediment_compressibility_coefficient', 'coefficient of compressibility of the sediment [dimensionless], recommended default: 0.12'))
     65    # }}}
    6666
    6767    def marshall(self, prefix, md, fid):  # {{{
    6868        yts = md.constants.yts
    69 
    7069        WriteData(fid, prefix, 'name', 'md.friction.law', 'data', 10, 'format', 'Integer')
    7170        WriteData(fid, prefix, 'class', 'friction', 'object', self, 'fieldname', 'pseudoplasticity_exponent', 'format', 'Double')
    … …  
    7574        WriteData(fid, prefix, 'class', 'friction', 'object', self, 'fieldname', 'till_friction_angle', 'format', 'DoubleMat', 'mattype', 1)
    7675        WriteData(fid, prefix, 'class', 'friction', 'object', self, 'fieldname', 'sediment_compressibility_coefficient', 'format', 'DoubleMat', 'mattype', 1)
    77         # }}}
     76    # }}}

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

    r25502 r25688  
    2626
    2727    def __repr__(self):  # {{{
    28         s = "Basal shear stress parameters: Sigma_b = coefficient^2 * Neff * u_b\n(effective stress Neff = rho_ice * g * thickness + rho_water * g * (head - b))\n"
     28        s = 'Basal shear stress parameters: Sigma_b = coefficient^2 * Neff * u_b\n'
     29        s += '(effective stress Neff = rho_ice * g * thickness + rho_water * g * (head - b))\n'
    2930        s += "{}\n".format(fielddisplay(self, "coefficient", "friction coefficient [SI]"))
    30 
    3131        return s
    3232    #}}}
    … …  
    4545        if 'StressbalanceAnalysis' not in analyses and 'ThermalAnalysis' not in analyses:
    4646            return md
    47 
    4847        md = checkfield(md, 'fieldname', 'friction.coefficient', 'timeseries', 1, 'NaN', 1, 'Inf', 1)
    49 
    5048        return md
    5149    # }}}

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

    r25503 r25688  
    3636        s = "{}\n".format(fielddisplay(self, 'q', 'q exponent'))
    3737        s = "{}\n".format(fielddisplay(self, 'water_layer', 'water thickness at the base of the ice (m)'))
    38 
    3938        return s
    4039    #}}}
    … …  
    4544
    4645    def checkconsistency(self, md, solution, analyses):  #{{{
    47         #Early return
    48         if ('StressbalanceAnalysis' not in analyses) and ('ThermalAnalysis' not in analyses):
     46        # Early return
     47        if 'StressbalanceAnalysis' not in analyses and 'ThermalAnalysis' not in analyses:
    4948            return
    50 
    5149        md = checkfield(md, 'fieldname', 'friction.coefficient', 'timeseries', 1, 'NaN', 1, 'Inf', 1)
    5250        md = checkfield(md, 'fieldname', 'friction.f', 'size', [1], 'NaN', 1, 'Inf', 1)
    … …  
    5452        md = checkfield(md, 'fieldname', 'friction.p', 'NaN', 1, 'Inf', 1, 'size', [md.mesh.numberofelements])
    5553        md = checkfield(md, 'fieldname', 'thermal.spctemperature', 'Inf', 1, 'timeseries', 1, '>=', 0.)
     54        return md
    5655    # }}}
    5756
    … …  
    6160        self.q = project3d(md, 'vector', self.q, 'type', 'element')
    6261        self.water_layer = project3d(md, 'vector', self.water_layer, 'type', 'node', 'layer', 1)
    63 
    6462        return self
    6563    # }}}

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

    r25169 r25688  
    88
    99class geometry(object):
    10     '''
    11     GEOMETRY class definition
     10    """GEOMETRY class definition
    1211
    13         Usage:
    14             geometry = geometry()
    15     '''
     12    Usage:
     13        geometry = geometry()
     14    """
    1615
    17     def __init__(self): #{{{
     16    def __init__(self, *args): #{{{
    1817        self.surface = np.nan
    1918        self.thickness = np.nan
    … …  
    2221        self.hydrostatic_ratio = np.nan
    2322
    24         #set defaults
    25         self.setdefaultparameters()
     23        nargs = len(args)
     24        if nargs == 0:
     25            self.setdefaultparameters()
     26        else:
     27            raise Exception('constructor not supported')
    2628    #}}}
    2729
    2830    def __repr__(self): #{{{
    29         s = "   geometry parameters:"
    30         s = "%s\n%s" % (s, fielddisplay(self, 'surface', 'ice upper surface elevation [m]'))
    31         s = "%s\n%s" % (s, fielddisplay(self, 'thickness', 'ice thickness [m]'))
    32         s = "%s\n%s" % (s, fielddisplay(self, 'base', 'ice base elevation [m]'))
    33         s = "%s\n%s" % (s, fielddisplay(self, 'bed', 'bed elevation [m]'))
    34        
     31        s = '   geometry parameters:\n'
     32        s += '{}\n'.format(fielddisplay(self, 'surface', 'ice upper surface elevation [m]'))
     33        s += '{}\n'.format(fielddisplay(self, 'thickness', 'ice thickness [m]'))
     34        s += '{}\n'.format(fielddisplay(self, 'base', 'ice base elevation [m]'))
     35        s += '{}\n'.format(fielddisplay(self, 'bed', 'bed elevation [m]'))
    3536        return s
    3637    #}}}

  • issm/trunk-jpl/src/m/classes/giaivins.m

    r25129 r25688  
    3535                                if size(md.geometry.thickness,1)==md.mesh.numberofvertices+1,
    3636                                        %recover the furthest time "in time":
    37                                         if(thickness(end,end)~=md.timestepping.start_time),
     37                                        if(md.geometry.thickness(end,end)~=md.timestepping.start_time),
    3838                                                md = checkmessage(md,['if ismasstransport is on, transient thickness forcing'...
    3939                                                        ' for the giaivins model should not be provided in the future.'...
    … …  
    5151                        fielddisplay(self,'lithosphere_thickness','lithosphere thickness (km)');
    5252                        fielddisplay(self,'cross_section_shape','1: square-edged (default). 2: elliptical.  See iedge in GiaDeflectionCore');
    53 
     53True
    5454                end % }}}
    5555                function marshall(self,prefix,md,fid) % {{{

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

    r25158 r25688  
    88
    99class giaivins(object):
    10     """
    11     GIA class definition for Ivins and James model
     10    """GIA class definition for Ivins and James model
    1211
    13        Usage:
    14           giaivins = giaivins()
     12    Usage:
     13        giaivins = giaivins()
    1514    """
    1615
    … …  
    2120
    2221        nargin = len(args)
    23 
    2422        if nargin == 0:
    2523            self.setdefaultparameters()
    … …  
    3028    def __repr__(self): #{{{
    3129        s = '   giaivins parameters:\n'
    32         s += "{}\n".format(fielddisplay(self, 'mantle_viscosity', 'mantle viscosity [Pa s]'))
    33         s += "{}\n".format(fielddisplay(self, 'lithosphere_thickness', 'lithosphere thickness (km)'))
    34         s += "{}\n".format(fielddisplay(self, 'cross_section_shape', "1: square-edged (default). 2: elliptical. See iedge in GiaDeflectionCore"))
    35 
     30        s += '{}\n'.format(fielddisplay(self, 'mantle_viscosity', 'mantle viscosity [Pa s]'))
     31        s += '{}\n'.format(fielddisplay(self, 'lithosphere_thickness', 'lithosphere thickness (km)'))
     32        s += '{}\n'.format(fielddisplay(self, 'cross_section_shape', "1: square-edged (default). 2: elliptical. See iedge in GiaDeflectionCore"))
    3633        return s
    3734    #}}}
    … …  
    3936    def setdefaultparameters(self): #{{{
    4037        self.cross_section_shape = 1
     38        return self
    4139    #}}}
    4240
    … …  
    4947        md = checkfield(md, 'fieldname', 'gia.cross_section_shape', 'numel', [1], 'values', [1, 2])
    5048
    51         #be sure that if we are running a masstransport ice flow model coupled with giaivins, that thickness forcings
    52         #are not provided into the future.
    53 
     49        # Be sure that if we are running a masstransport ice flow model coupled
     50        # with giaivins, that thickness forcings are not provided into the
     51        # future.
     52        if solution == 'TransientSolution' and md.transient.ismasstransport and md.transient.isgia:
     53            # Figure out if thickness is a transient forcing
     54            if md.geometry.thickness.shape[0] == md.mesh.numberofvertices + 1:
     55                # Recover the furthest time "in time"
     56                if md.geometry.thickness[-1, -1] != md.timestepping.start_time:
     57                    md.checkmessage('if masstransport is on, transient thickness forcing for the giaivins model should not be provided in the future. Synchronize your start_time to correspond to the most recent transient thickness forcing timestep')
    5458        return md
    5559    #}}}

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

    r25169 r25688  
    88
    99class giamme(object):
    10     '''
    11     GIAMME class definition
     10    """GIAMME class definition
    1211
    13         Usage:
    14             gia = giamme() #gia class based on a multi-model ensemble (ex: Caron et al 2017 statistics)
    15     '''
     12    Usage:
     13        gia = giamme() #gia class based on a multi-model ensemble (ex: Caron et al 2017 statistics)
     14    """
    1615
    1716    def __init__(self, *args): #{{{
    … …  
    3332        s += '{}\n'.format(fielddisplay(self, 'Ngia', 'geoid (mm/yr).'))
    3433        s += '{}\n'.format(fielddisplay(self, 'Ugia', 'uplift (mm/yr).'))
    35 
    3634        return s
    3735    #}}}
    3836
    3937    def setdefaultparameters(self): # {{{
    40         return
     38        return self
    4139    #}}}
    4240
    4341    def checkconsistency(self, md, solution, analyses): # {{{
    44         if ('SealevelriseAnalysis' not in analyses) or (solution == 'TransientSolution' and md.transient.isslr == 0):
     42        if 'SealevelriseAnalysis' not in analyses:
    4543            return md
    46 
     44        if solution == 'TransientSolution' and not md.transient.isslr:
     45            return md
    4746        md = checkfield(md, 'field', self.Ngia, 'NaN', 1, 'Inf', 1)
    4847        md = checkfield(md, 'field', self.Ugia, 'NaN', 1, 'Inf', 1)
    4948        md = checkfield(md, 'field', self.modelid, 'NaN', 1, 'Inf', 1, '>=', 0, '<=', len(self.Ngia))
    50 
    5149        return md
    5250    #}}}
    … …  
    7169            self.Ngia[i] = project3d(md, 'vector', self.Ngia[i], 'type', 'node', 'layer', 1)
    7270            self.Ugia[i] = project3d(md, 'vector', self.Ugia[i], 'type', 'node', 'layer', 1)
    73 
    7471        return self
    7572    #}}}

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

    r25385 r25688  
    11import numpy as np
     2
     3from checkfield import checkfield
     4from fielddisplay import fielddisplay
    25from project3d import project3d
    3 from fielddisplay import fielddisplay
    4 from checkfield import checkfield
    56from WriteData import WriteData
    67
    78
    89class initialization(object):
    9     """
    10     INITIALIZATION class definition
     10    """INITIALIZATION class definition
    1111
    1212    Usage:
    … …  
    1515
    1616    def __init__(self):  # {{{
     17        self.vx = np.nan
     18        self.vy = np.nan
     19        self.vz = np.nan
     20        self.vel = np.nan
     21        self.enthalpy = np.nan
     22        self.pressure = np.nan
     23        self.temperature = np.nan
     24        self.waterfraction = np.nan
     25        self.watercolumn = np.nan
     26        self.sediment_head = np.nan
     27        self.epl_head = np.nan
     28        self.epl_thickness = np.nan
     29        self.hydraulic_potential = np.nan
     30        self.channelarea = np.nan
    1731
    18         self.vx = float('NaN')
    19         self.vy = float('NaN')
    20         self.vz = float('NaN')
    21         self.vel = float('NaN')
    22         self.enthalpy = float('NaN')
    23         self.pressure = float('NaN')
    24         self.temperature = float('NaN')
    25         self.waterfraction = float('NaN')
    26         self.watercolumn = float('NaN')
    27         self.sediment_head = float('NaN')
    28         self.epl_head = float('NaN')
    29         self.epl_thickness = float('NaN')
    30         self.hydraulic_potential = float('NaN')
    31         self.channelarea = float('NaN')
    32 
    33     #set defaults
     32        #set defaults
    3433        self.setdefaultparameters()
    35 
    3634    #}}}
    3735
    3836    def __repr__(self):  # {{{
    39         string = '   initial field values:'
    40         string = "%s\n%s" % (string, fielddisplay(self, 'vx', 'x component of velocity [m / yr]'))
    41         string = "%s\n%s" % (string, fielddisplay(self, 'vy', 'y component of velocity [m / yr]'))
    42         string = "%s\n%s" % (string, fielddisplay(self, 'vz', 'z component of velocity [m / yr]'))
    43         string = "%s\n%s" % (string, fielddisplay(self, 'vel', 'velocity norm [m / yr]'))
    44         string = "%s\n%s" % (string, fielddisplay(self, 'pressure', 'pressure [Pa]'))
    45         string = "%s\n%s" % (string, fielddisplay(self, 'temperature', 'temperature [K]'))
    46         string = "%s\n%s" % (string, fielddisplay(self, 'enthalpy', 'enthalpy [J]'))
    47         string = "%s\n%s" % (string, fielddisplay(self, 'waterfraction', 'fraction of water in the ice'))
    48         string = "%s\n%s" % (string, fielddisplay(self, 'watercolumn', 'thickness of subglacial water [m]'))
    49         string = "%s\n%s" % (string, fielddisplay(self, 'sediment_head', 'sediment water head of subglacial system [m]'))
    50         string = "%s\n%s" % (string, fielddisplay(self, 'epl_head', 'epl water head of subglacial system [m]'))
    51         string = "%s\n%s" % (string, fielddisplay(self, 'epl_thickness', 'thickness of the epl [m]'))
    52         string = "%s\n%s" % (string, fielddisplay(self, 'hydraulic_potential', 'Hydraulic potential (for GlaDS) [Pa]'))
    53         string = "%s\n%s" % (string, fielddisplay(self, 'channelarea', 'subglaciale water channel area (for GlaDS) [m2]'))
    54 
    55         return string
     37        s = '   initial field values:\n'
     38        s += '{}\n'.format(fielddisplay(self, 'vx', 'x component of velocity [m / yr]'))
     39        s += '{}\n'.format(fielddisplay(self, 'vy', 'y component of velocity [m / yr]'))
     40        s += '{}\n'.format(fielddisplay(self, 'vz', 'z component of velocity [m / yr]'))
     41        s += '{}\n'.format(fielddisplay(self, 'vel', 'velocity norm [m / yr]'))
     42        s += '{}\n'.format(fielddisplay(self, 'pressure', 'pressure [Pa]'))
     43        s += '{}\n'.format(fielddisplay(self, 'temperature', 'temperature [K]'))
     44        s += '{}\n'.format(fielddisplay(self, 'enthalpy', 'enthalpy [J]'))
     45        s += '{}\n'.format(fielddisplay(self, 'waterfraction', 'fraction of water in the ice'))
     46        s += '{}\n'.format(fielddisplay(self, 'watercolumn', 'thickness of subglacial water [m]'))
     47        s += '{}\n'.format(fielddisplay(self, 'sediment_head', 'sediment water head of subglacial system [m]'))
     48        s += '{}\n'.format(fielddisplay(self, 'epl_head', 'epl water head of subglacial system [m]'))
     49        s += '{}\n'.format(fielddisplay(self, 'epl_thickness', 'thickness of the epl [m]'))
     50        s += '{}\n'.format(fielddisplay(self, 'hydraulic_potential', 'Hydraulic potential (for GlaDS) [Pa]'))
     51        s += '{}\n'.format(fielddisplay(self, 'channelarea', 'subglaciale water channel area (for GlaDS) [m2]'))
     52        return s
    5653    #}}}
    5754
    … …  
    7067
    7168        #Lithostatic pressure by default
    72         #        self.pressure = md.constants.g * md.materials.rho_ice * (md.geometry.surface[:, 0] - md.mesh.z)
    73         #self.pressure = md.constants.g * md.materials.rho_ice * (md.geometry.surface-md.mesh.z.reshape(-1, ))
    74 
    7569        if np.ndim(md.geometry.surface) == 2:
    76             print('Reshaping md.geometry.surface for you convenience but you should fix it in you files')
     70            print('Reshaping md.geometry.surface for your convenience but you should fix it in your model set up')
    7771            self.pressure = md.constants.g * md.materials.rho_ice * (md.geometry.surface.reshape(-1, ) - md.mesh.z)
    7872        else:
    … …  
    8781
    8882    def checkconsistency(self, md, solution, analyses):  # {{{
    89         if 'StressbalanceAnalysis' in analyses:
     83        if 'StressbalanceAnalysis' in analyses and not solution == 'TransientSolution' and not md.transient.isstressbalance:
    9084            if not np.any(np.logical_or(np.isnan(md.initialization.vx), np.isnan(md.initialization.vy))):
    9185                md = checkfield(md, 'fieldname', 'initialization.vx', 'NaN', 1, 'Inf', 1, 'size', [md.mesh.numberofvertices])
    9286                md = checkfield(md, 'fieldname', 'initialization.vy', 'NaN', 1, 'Inf', 1, 'size', [md.mesh.numberofvertices])
    93         if 'MasstransportAnalysis' in analyses:
     87        if 'MasstransportAnalysis' in analyses and not solution == 'TransientSolution' and not md.transient.ismasstransport:
    9488            md = checkfield(md, 'fieldname', 'initialization.vx', 'NaN', 1, 'Inf', 1, 'size', [md.mesh.numberofvertices])
    9589            md = checkfield(md, 'fieldname', 'initialization.vy', 'NaN', 1, 'Inf', 1, 'size', [md.mesh.numberofvertices])
    96         if 'BalancethicknessAnalysis' in analyses:
     90        if 'BalancethicknessAnalysis' in analyses and solution == 'BalancethicknessSolution':
    9791            md = checkfield(md, 'fieldname', 'initialization.vx', 'NaN', 1, 'Inf', 1, 'size', [md.mesh.numberofvertices])
    9892            md = checkfield(md, 'fieldname', 'initialization.vy', 'NaN', 1, 'Inf', 1, 'size', [md.mesh.numberofvertices])
    99             #Triangle with zero velocity
     93            # Triangle with zero velocity
    10094            if np.any(np.logical_and(np.sum(np.abs(md.initialization.vx[md.mesh.elements - 1]), axis=1) == 0,
    10195                                     np.sum(np.abs(md.initialization.vy[md.mesh.elements - 1]), axis=1) == 0)):
    10296                md.checkmessage("at least one triangle has all its vertices with a zero velocity")
    103         if 'ThermalAnalysis' in analyses and not (solution == 'TransientSolution' and md.transient.isthermal == 0):
     97        if 'ThermalAnalysis' in analyses and not solution == 'TransientSolution' and not md.transient.isthermal:
    10498            md = checkfield(md, 'fieldname', 'initialization.vx', 'NaN', 1, 'Inf', 1, 'size', [md.mesh.numberofvertices])
    10599            md = checkfield(md, 'fieldname', 'initialization.vy', 'NaN', 1, 'Inf', 1, 'size', [md.mesh.numberofvertices])
    … …  
    108102                md = checkfield(md, 'fieldname', 'initialization.vz', 'NaN', 1, 'Inf', 1, 'size', [md.mesh.numberofvertices])
    109103            md = checkfield(md, 'fieldname', 'initialization.pressure', 'NaN', 1, 'Inf', 1, 'size', [md.mesh.numberofvertices])
    110             if ('EnthalpyAnalysis' in analyses and md.thermal.isenthalpy):
    111                 md = checkfield(md, 'fieldname', 'initialization.waterfraction', '>=', 0, 'size', [md.mesh.numberofvertices])
    112                 md = checkfield(md, 'fieldname', 'initialization.watercolumn', '>=', 0, 'size', [md.mesh.numberofvertices])
    113                 pos = np.nonzero(md.initialization.waterfraction > 0.)[0]
    114                 if(pos.size):
    115                     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')
     104        if 'EnthalpyAnalysis' in analyses and md.thermal.isenthalpy:
     105            md = checkfield(md, 'fieldname', 'initialization.waterfraction', '>=', 0, 'size', [md.mesh.numberofvertices])
     106            md = checkfield(md, 'fieldname', 'initialization.watercolumn', '>=', 0, 'size', [md.mesh.numberofvertices])
     107            pos = np.nonzero(md.initialization.waterfraction > 0.)[0]
     108            if(pos.size):
     109                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')
    116110        if 'HydrologyShreveAnalysis' in analyses:
    117111            if hasattr(md.hydrology, 'hydrologyshreve'):
    … …  
    134128
    135129    def marshall(self, prefix, md, fid):  # {{{
    136 
    137130        yts = md.constants.yts
    138 
    139131        WriteData(fid, prefix, 'object', self, 'fieldname', 'vx', 'format', 'DoubleMat', 'mattype', 1, 'scale', 1. / yts)
    140132        WriteData(fid, prefix, 'object', self, 'fieldname', 'vy', 'format', 'DoubleMat', 'mattype', 1, 'scale', 1. / yts)

  • issm/trunk-jpl/src/m/classes/inversion.m

    r22303 r25688  
    108108                                end
    109109                        end
    110 
    111110                        if strcmp(solution,'BalancethicknessSolution')
    112111                                md = checkfield(md,'fieldname','inversion.thickness_obs','size',[md.mesh.numberofvertices 1],'NaN',1,'Inf',1);

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

    r24261 r25688  
    11import numpy as np
     2
     3from checkfield import checkfield
     4from fielddisplay import fielddisplay
     5from marshallcostfunctions import marshallcostfunctions
    26from project3d import project3d
    3 from fielddisplay import fielddisplay
    4 from checkfield import checkfield
    5 from WriteData import WriteData
    67from supportedcontrols import supportedcontrols
    78from supportedcostfunctions import supportedcostfunctions
    8 from marshallcostfunctions import marshallcostfunctions
     9from WriteData import WriteData
    910
    1011
    1112class inversion(object):
    12     """
    13     INVERSION class definition
     13    """INVERSION class definition
    1414
    15        Usage:
    16           inversion = inversion()
     15    Usage:
     16        inversion = inversion()
    1717    """
    1818
    … …  
    2020        self.iscontrol = 0
    2121        self.incomplete_adjoint = 0
    22         self.control_parameters = float('NaN')
     22        self.control_parameters = np.nan
    2323        self.nsteps = 0
    24         self.maxiter_per_step = float('NaN')
     24        self.maxiter_per_step = np.nan
    2525        self.cost_functions = ''
    26         self.cost_functions_coefficients = float('NaN')
    27         self.gradient_scaling = float('NaN')
     26        self.cost_functions_coefficients = np.nan
     27        self.gradient_scaling = np.nan
    2828        self.cost_function_threshold = 0
    29         self.min_parameters = float('NaN')
    30         self.max_parameters = float('NaN')
    31         self.step_threshold = float('NaN')
    32         self.vx_obs = float('NaN')
    33         self.vy_obs = float('NaN')
    34         self.vz_obs = float('NaN')
    35         self.vel_obs = float('NaN')
    36         self.thickness_obs = float('NaN')
    37         self.surface_obs = float('NaN')
     29        self.min_parameters = np.nan
     30        self.max_parameters = np.nan
     31        self.step_threshold = np.nan
     32        self.vx_obs = np.nan
     33        self.vy_obs = np.nan
     34        self.vz_obs = np.nan
     35        self.vel_obs = np.nan
     36        self.thickness_obs = np.nan
     37        self.surface_obs = np.nan
    3838
    39     #set defaults
    4039        self.setdefaultparameters()
    41 
    4240    #}}}
    4341
    4442    def __repr__(self):  # {{{
    45         string = '   inversion parameters:'
    46         string = "%s\n%s" % (string, fielddisplay(self, 'iscontrol', 'is inversion activated?'))
    47         string = "%s\n%s" % (string, fielddisplay(self, 'incomplete_adjoint', '1: linear viscosity, 0: non - linear viscosity'))
    48         string = "%s\n%s" % (string, fielddisplay(self, 'control_parameters', 'ex: {''FrictionCoefficient''}, or {''MaterialsRheologyBbar''}'))
    49         string = "%s\n%s" % (string, fielddisplay(self, 'nsteps', 'number of optimization searches'))
    50         string = "%s\n%s" % (string, fielddisplay(self, 'cost_functions', 'indicate the type of response for each optimization step'))
    51         string = "%s\n%s" % (string, fielddisplay(self, 'cost_functions_coefficients', 'cost_functions_coefficients applied to the misfit of each vertex and for each control_parameter'))
    52         string = "%s\n%s" % (string, fielddisplay(self, 'cost_function_threshold', 'misfit convergence criterion. Default is 1%, NaN if not applied'))
    53         string = "%s\n%s" % (string, fielddisplay(self, 'maxiter_per_step', 'maximum iterations during each optimization step'))
    54         string = "%s\n%s" % (string, fielddisplay(self, 'gradient_scaling', 'scaling factor on gradient direction during optimization, for each optimization step'))
    55         string = "%s\n%s" % (string, fielddisplay(self, 'step_threshold', 'decrease threshold for misfit, default is 30%'))
    56         string = "%s\n%s" % (string, fielddisplay(self, 'min_parameters', 'absolute minimum acceptable value of the inversed parameter on each vertex'))
    57         string = "%s\n%s" % (string, fielddisplay(self, 'max_parameters', 'absolute maximum acceptable value of the inversed parameter on each vertex'))
    58         string = "%s\n%s" % (string, fielddisplay(self, 'vx_obs', 'observed velocity x component [m / yr]'))
    59         string = "%s\n%s" % (string, fielddisplay(self, 'vy_obs', 'observed velocity y component [m / yr]'))
    60         string = "%s\n%s" % (string, fielddisplay(self, 'vel_obs', 'observed velocity magnitude [m / yr]'))
    61         string = "%s\n%s" % (string, fielddisplay(self, 'thickness_obs', 'observed thickness [m]'))
    62         string = "%s\n%s" % (string, fielddisplay(self, 'surface_obs', 'observed surface elevation [m]'))
    63         string = "%s\n%s" % (string, 'Available cost functions:')
    64         string = "%s\n%s" % (string, '   101: SurfaceAbsVelMisfit')
    65         string = "%s\n%s" % (string, '   102: SurfaceRelVelMisfit')
    66         string = "%s\n%s" % (string, '   103: SurfaceLogVelMisfit')
    67         string = "%s\n%s" % (string, '   104: SurfaceLogVxVyMisfit')
    68         string = "%s\n%s" % (string, '   105: SurfaceAverageVelMisfit')
    69         string = "%s\n%s" % (string, '   201: ThicknessAbsMisfit')
    70         string = "%s\n%s" % (string, '   501: DragCoefficientAbsGradient')
    71         string = "%s\n%s" % (string, '   502: RheologyBbarAbsGradient')
    72         string = "%s\n%s" % (string, '   503: ThicknessAbsGradient')
    73         return string
    74     #}}}
    75 
    76     def extrude(self, md):  # {{{
    77         self.vx_obs = project3d(md, 'vector', self.vx_obs, 'type', 'node')
    78         self.vy_obs = project3d(md, 'vector', self.vy_obs, 'type', 'node')
    79         self.vel_obs = project3d(md, 'vector', self.vel_obs, 'type', 'node')
    80         self.thickness_obs = project3d(md, 'vector', self.thickness_obs, 'type', 'node')
    81         if not np.any(np.isnan(self.cost_functions_coefficients)):
    82             self.cost_functions_coefficients = project3d(md, 'vector', self.cost_functions_coefficients, 'type', 'node')
    83         if not np.any(np.isnan(self.min_parameters)):
    84             self.min_parameters = project3d(md, 'vector', self.min_parameters, 'type', 'node')
    85         if not np.any(np.isnan(self.max_parameters)):
    86             self.max_parameters = project3d(md, 'vector', self.max_parameters, 'type', 'node')
    87         return self
     43        s = '   inversion parameters:\n'
     44        s += '{}\n'.format(fielddisplay(self, 'iscontrol', 'is inversion activated?'))
     45        s += '{}\n'.format(fielddisplay(self, 'incomplete_adjoint', '1: linear viscosity, 0: non - linear viscosity'))
     46        s += '{}\n'.format(fielddisplay(self, 'control_parameters', 'ex: {''FrictionCoefficient''}, or {''MaterialsRheologyBbar''}'))
     47        s += '{}\n'.format(fielddisplay(self, 'nsteps', 'number of optimization searches'))
     48        s += '{}\n'.format(fielddisplay(self, 'cost_functions', 'indicate the type of response for each optimization step'))
     49        s += '{}\n'.format(fielddisplay(self, 'cost_functions_coefficients', 'cost_functions_coefficients applied to the misfit of each vertex and for each control_parameter'))
     50        s += '{}\n'.format(fielddisplay(self, 'cost_function_threshold', 'misfit convergence criterion. Default is 1%, NaN if not applied'))
     51        s += '{}\n'.format(fielddisplay(self, 'maxiter_per_step', 'maximum iterations during each optimization step'))
     52        s += '{}\n'.format(fielddisplay(self, 'gradient_scaling', 'scaling factor on gradient direction during optimization, for each optimization step'))
     53        s += '{}\n'.format(fielddisplay(self, 'step_threshold', 'decrease threshold for misfit, default is 30%'))
     54        s += '{}\n'.format(fielddisplay(self, 'min_parameters', 'absolute minimum acceptable value of the inversed parameter on each vertex'))
     55        s += '{}\n'.format(fielddisplay(self, 'max_parameters', 'absolute maximum acceptable value of the inversed parameter on each vertex'))
     56        s += '{}\n'.format(fielddisplay(self, 'vx_obs', 'observed velocity x component [m / yr]'))
     57        s += '{}\n'.format(fielddisplay(self, 'vy_obs', 'observed velocity y component [m / yr]'))
     58        s += '{}\n'.format(fielddisplay(self, 'vel_obs', 'observed velocity magnitude [m / yr]'))
     59        s += '{}\n'.format(fielddisplay(self, 'thickness_obs', 'observed thickness [m]'))
     60        s += '{}\n'.format(fielddisplay(self, 'surface_obs', 'observed surface elevation [m]'))
     61        s += '{}\n'.format('Available cost functions:')
     62        s += '{}\n'.format('   101: SurfaceAbsVelMisfit')
     63        s += '{}\n'.format('   102: SurfaceRelVelMisfit')
     64        s += '{}\n'.format('   103: SurfaceLogVelMisfit')
     65        s += '{}\n'.format('   104: SurfaceLogVxVyMisfit')
     66        s += '{}\n'.format('   105: SurfaceAverageVelMisfit')
     67        s += '{}\n'.format('   201: ThicknessAbsMisfit')
     68        s += '{}\n'.format('   501: DragCoefficientAbsGradient')
     69        s += '{}\n'.format('   502: RheologyBbarAbsGradient')
     70        s += '{}\n'.format('   503: ThicknessAbsGradient')
     71        return s
    8872    #}}}
    8973
    9074    def setdefaultparameters(self):  # {{{
    91 
    9275        #default is incomplete adjoint for now
    9376        self.incomplete_adjoint = 1
    … …  
    11598        #if J[n] - J[n - 1] / J[n] < criteria, the control run stops
    11699        #NaN if not applied
    117         self.cost_function_threshold = float('NaN')  #not activated
     100        self.cost_function_threshold = np.nan  #not activated
     101        return self
     102    #}}}
    118103
     104    def extrude(self, md):  # {{{
     105        self.vx_obs = project3d(md, 'vector', self.vx_obs, 'type', 'node')
     106        self.vy_obs = project3d(md, 'vector', self.vy_obs, 'type', 'node')
     107        self.vel_obs = project3d(md, 'vector', self.vel_obs, 'type', 'node')
     108        self.thickness_obs = project3d(md, 'vector', self.thickness_obs, 'type', 'node')
     109        if not np.any(np.isnan(self.cost_functions_coefficients)):
     110            self.cost_functions_coefficients = project3d(md, 'vector', self.cost_functions_coefficients, 'type', 'node')
     111        if not np.any(np.isnan(self.min_parameters)):
     112            self.min_parameters = project3d(md, 'vector', self.min_parameters, 'type', 'node')
     113        if not np.any(np.isnan(self.max_parameters)):
     114            self.max_parameters = project3d(md, 'vector', self.max_parameters, 'type', 'node')
    119115        return self
    120116    #}}}
    121117
    122118    def checkconsistency(self, md, solution, analyses):  # {{{
    123         #Early return
     119        # Early return
    124120        if not self.iscontrol:
    125121            return md
    … …  
    140136        md = checkfield(md, 'fieldname', 'inversion.max_parameters', 'size', [md.mesh.numberofvertices, num_controls])
    141137
    142     #Only SSA, HO and FS are supported right now
     138        # Only SSA, HO and FS are supported right now
    143139        if solution == 'StressbalanceSolution':
    144140            if not (md.flowequation.isSSA or md.flowequation.isHO or md.flowequation.isFS or md.flowequation.isL1L2):
    145141                md.checkmessage("'inversion can only be performed for SSA, HO or FS ice flow models")
    146 
    147142        if solution == 'BalancethicknessSolution':
     143            md = checkfield(md, 'fieldname', 'inversion.thickness_obs', 'size', [md.mesh.numberofvertices], 'NaN', 1, 'Inf', 1)
     144        elif solution == 'BalancethicknessSoftSolution':
    148145            md = checkfield(md, 'fieldname', 'inversion.thickness_obs', 'size', [md.mesh.numberofvertices], 'NaN', 1, 'Inf', 1)
    149146        else:
    150147            md = checkfield(md, 'fieldname', 'inversion.vx_obs', 'size', [md.mesh.numberofvertices], 'NaN', 1, 'Inf', 1)
    151148            md = checkfield(md, 'fieldname', 'inversion.vy_obs', 'size', [md.mesh.numberofvertices], 'NaN', 1, 'Inf', 1)
    152 
    153149        return md
    154150    # }}}
    … …  
    177173        WriteData(fid, prefix, 'object', self, 'fieldname', 'surface_obs', 'format', 'DoubleMat', 'mattype', 1)
    178174
    179     #process control parameters
     175        # Process control parameters
    180176        num_control_parameters = len(self.control_parameters)
    181177        WriteData(fid, prefix, 'object', self, 'fieldname', 'control_parameters', 'format', 'StringArray')
    182178        WriteData(fid, prefix, 'data', num_control_parameters, 'name', 'md.inversion.num_control_parameters', 'format', 'Integer')
    183179
    184     #process cost functions
     180        # Process cost functions
    185181        num_cost_functions = np.size(self.cost_functions)
    186182        data = marshallcostfunctions(self.cost_functions)

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

    r25158 r25688  
    33from checkfield import checkfield
    44from fielddisplay import fielddisplay
     5from project3d import project3d
    56from WriteData import WriteData
    67
    78
    89class linearbasalforcings(object):
    9     """
    10     LINEAR BASAL FORCINGS class definition
     10    """LINEAR BASAL FORCINGS class definition
    1111
    12        Usage:
    13           basalforcings = linearbasalforcings()
     12    Usage:
     13        basalforcings = linearbasalforcings()
    1414    """
    1515
    16     def __init__(self, *args):  # {{{
    17 
    18         if not len(args):
     16    def __init__(self, *args): #{{{
     17        nargs = len(args)
     18        if nargs == 0:
    1919            print('empty init')
    2020            self.deepwater_melting_rate = 0.
    … …  
    2626            self.geothermalflux = float('NaN')
    2727
    28     #set defaults
     28            # set defaults
    2929            self.setdefaultparameters()
    30         elif len(args) == 1 and args[0].__module__ == 'basalforcings':
     30        elif nargs == 1 and args[0].__module__ == 'basalforcings':
    3131            print('converting basalforings to linearbasalforcings')
    3232            inv = args[0]
    … …  
    3939            self.upperwater_elevation = 0.
    4040
    41     #set defaults
     41            # Set defaults
    4242            self.setdefaultparameters()
    4343        else:
    4444            raise Exception('constructor not supported')
    4545    #}}}
    46     def __repr__(self):  # {{{
    47         string = "   linear basal forcings parameters:"
    48         string = "%s\n%s" % (string, fielddisplay(self, "deepwater_melting_rate", "basal melting rate (positive if melting applied for floating ice whith base < deepwater_elevation) [m/yr]"))
    49         string = "%s\n%s" % (string, fielddisplay(self, "deepwater_elevation", "elevation of ocean deepwater [m]"))
    50         string = "%s\n%s" % (string, fielddisplay(self, "upperwater_melting_rate", "upper melting rate (positive if melting applied for floating ice whith base >= upperwater_elevation) [m/yr]"))
    51         string = "%s\n%s" % (string, fielddisplay(self, "upperwater_elevation", "elevation of ocean upper water [m]"))
    52         string = "%s\n%s" % (string, fielddisplay(self, "groundedice_melting_rate", "basal melting rate (positive if melting) [m/yr]"))
    53         string = "%s\n%s" % (string, fielddisplay(self, "perturbation_melting_rate", "perturbation applied to computed melting rate (positive if melting) [m/yr]"))
    54         string = "%s\n%s" % (string, fielddisplay(self, "geothermalflux", "geothermal heat flux [W/m^2]"))
    55         return string
     46
     47    def __repr__(self): #{{{
     48        s = '   linear basal forcings parameters:\n'
     49        s += '{}\n'.format(fielddisplay(self, "deepwater_melting_rate", "basal melting rate (positive if melting applied for floating ice whith base < deepwater_elevation) [m/yr]"))
     50        s += '{}\n'.format(fielddisplay(self, "deepwater_elevation", "elevation of ocean deepwater [m]"))
     51        s += '{}\n'.format(fielddisplay(self, "upperwater_melting_rate", "upper melting rate (positive if melting applied for floating ice whith base >= upperwater_elevation) [m/yr]"))
     52        s += '{}\n'.format(fielddisplay(self, "upperwater_elevation", "elevation of ocean upper water [m]"))
     53        s += '{}\n'.format(fielddisplay(self, "groundedice_melting_rate", "basal melting rate (positive if melting) [m/yr]"))
     54        s += '{}\n'.format(fielddisplay(self, "perturbation_melting_rate", "perturbation applied to computed melting rate (positive if melting) [m/yr]"))
     55        s += '{}\n'.format(fielddisplay(self, "geothermalflux", "geothermal heat flux [W/m^2]"))
     56        return s
    5657    #}}}
    57     def initialize(self, md):  # {{{
     58
     59    def extrude(self, md): #{{{
     60        self.perturbation_melting_rate = project3d(md, 'vector', self.perturbation_melting_rate, 'type', 'node', 'layer', 1)
     61        self.groundedice_melting_rate = project3d(md, 'vector', self.groundedice_melting_rate, 'type', 'node', 'layer', 1)
     62        self.geothermalflux = project3d(md, 'vector', self.geothermalflux, 'type', 'node', 'layer', 1) # Bedrock only gets geothermal flux
     63        return self
     64    #}}}
     65
     66    def initialize(self, md): #{{{
    5867        if np.all(np.isnan(self.groundedice_melting_rate)):
    5968            self.groundedice_melting_rate = np.zeros((md.mesh.numberofvertices))
    … …  
    6170        return self
    6271    #}}}
     72
    6373    def setdefaultparameters(self):  # {{{
    6474        self.deepwater_melting_rate = 50.0
    … …  
    6676        self.upperwater_melting_rate = 0.0
    6777        self.upperwater_elevation = -400.0
     78        return self
    6879    #}}}
    69     def checkconsistency(self, md, solution, analyses):  # {{{
    7080
     81    def checkconsistency(self, md, solution, analyses): #{{{
    7182        if not np.all(np.isnan(self.perturbation_melting_rate)):
    7283            md = checkfield(md, 'fieldname', 'basalforcings.perturbation_melting_rate', 'NaN', 1, 'Inf', 1, 'timeseries', 1)
    73 
    74         if 'MasstransportAnalysis' in analyses and not (solution == 'TransientSolution' and not md.transient.ismasstransport):
     84        if 'MasstransportAnalysis' in analyses and not solution == 'TransientSolution' and not md.transient.ismasstransport:
    7585            md = checkfield(md, 'fieldname', 'basalforcings.groundedice_melting_rate', 'NaN', 1, 'Inf', 1, 'timeseries', 1)
    7686            md = checkfield(md, 'fieldname', 'basalforcings.deepwater_melting_rate', '>=', 0, 'singletimeseries', 1)
    … …  
    7888            md = checkfield(md, 'fieldname', 'basalforcings.deepwater_elevation', '<', 'basalforcings.upperwater_elevation', 'singletimeseries', 1)
    7989            md = checkfield(md, 'fieldname', 'basalforcings.upperwater_elevation', '<=', 0, 'singletimeseries', 1)
    80 
    8190        if 'BalancethicknessAnalysis' in analyses:
    8291            md = checkfield(md, 'fieldname', 'basalforcings.groundedice_melting_rate', 'NaN', 1, 'Inf', 1, 'size', [md.mesh.numberofvertices])
    … …  
    8594            md = checkfield(md, 'fieldname', 'basalforcings.deepwater_elevation', '<', 'basalforcings.upperwater_elevation', 'singletimeseries', 1)
    8695            md = checkfield(md, 'fieldname', 'basalforcings.upperwater_elevation', '<=', 0, 'singletimeseries', 1)
    87 
    88         if 'ThermalAnalysis' in analyses and not (solution == 'TransientSolution' and not md.transient.isthermal):
     96        if 'ThermalAnalysis' in analyses and not solution == 'TransientSolution' and not md.transient.isthermal:
    8997            md = checkfield(md, 'fieldname', 'basalforcings.groundedice_melting_rate', 'NaN', 1, 'Inf', 1, 'timeseries', 1)
    9098            md = checkfield(md, 'fieldname', 'basalforcings.deepwater_melting_rate', '>=', 0, 'singletimeseries', 1)
    … …  
    95103
    96104        return md
    97     # }}}
    98     def marshall(self, prefix, md, fid):  # {{{
     105    #}}}
     106
     107    def marshall(self, prefix, md, fid): #{{{
    99108        yts = md.constants.yts
    100109
    … …  
    107116        WriteData(fid, prefix, 'object', self, 'fieldname', 'upperwater_melting_rate', 'format', 'DoubleMat', 'mattype', 3, 'timeserieslength', 2, 'name', 'md.basalforcings.upperwater_melting_rate', 'scale', 1. / yts, 'yts', yts)
    108117        WriteData(fid, prefix, 'object', self, 'fieldname', 'upperwater_elevation', 'format', 'DoubleMat', 'mattype', 3, 'name', 'md.basalforcings.upperwater_elevation', 'yts', yts)
    109     # }}}
     118    #}}}

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

    r25183 r25688  
    99
    1010class lovenumbers(object): #{{{
    11     '''
    12     LOVENUMBERS numbers class definition
     11    """LOVENUMBERS class definition
    1312
    14         Usage:
    15             lovenumbers = lovenumbers() #will setup love numbers deg 1001 by default
    16             lovenumbers = lovenumbers('maxdeg', 10001);  #supply numbers of degrees required (here, 10001)
    17     '''
     13    Usage:
     14        lovenumbers = lovenumbers() #will setup love numbers deg 1001 by default
     15        lovenumbers = lovenumbers('maxdeg', 10001);  #supply numbers of degrees required (here, 10001)
     16    """
    1817
    1918    def __init__(self, *args): #{{{
    20         #regular love numbers:
    21         self.h = []  #provided by PREM model
    22         self.k = []  #idem
    23         self.l = []  #idem
     19        # Regular love numbers
     20        self.h = []  # Provided by PREM model
     21        self.k = []  # idem
     22        self.l = []  # idem
    2423
    25         #tidal love numbers for computing rotational feedback:
     24        # Tidal love numbers for computing rotational feedback
    2625        self.th = []
    2726        self.tk = []
    2827        self.tl = []
    29         self.tk2secular = 0 #deg 2 secular number.
     28        self.tk2secular = 0 # deg 2 secular number
    3029
    3130        options = pairoptions(*args)
    … …  
    3635
    3736    def setdefaultparameters(self, maxdeg, referenceframe): #{{{
    38         #initialize love numbers:
     37        # Initialize love numbers
    3938        self.h = getlovenumbers('type', 'loadingverticaldisplacement', 'referenceframe', referenceframe, 'maxdeg', maxdeg)
    4039        self.k = getlovenumbers('type', 'loadinggravitationalpotential', 'referenceframe', referenceframe, 'maxdeg', maxdeg)
    … …  
    4443        self.tl = getlovenumbers('type', 'tidalhorizontaldisplacement', 'referenceframe', referenceframe, 'maxdeg', maxdeg)
    4544
    46         #secular fluid love number:
     45        # Secular fluid love number
    4746        self.tk2secular = 0.942
     47        return self
    4848    #}}}
    4949
    5050    def checkconsistency(self, md, solution, analyses): #{{{
    51         if ('SealevelriseAnalysis' not in analyses) or (solution == 'TransientSolution' and md.transient.isslr == 0):
     51        if ('SealevelriseAnalysis' not in analyses) or (solution == 'TransientSolution' and not md.transient.isslr):
    5252            return
    53 
    5453        md = checkfield(md, 'fieldname', 'solidearth.lovenumbers.h', 'NaN', 1, 'Inf', 1)
    5554        md = checkfield(md, 'fieldname', 'solidearth.lovenumbers.k', 'NaN', 1, 'Inf', 1)
    … …  
    5958        md = checkfield(md, 'fieldname', 'solidearth.lovenumbers.tk', 'NaN', 1, 'Inf', 1)
    6059        md = checkfield(md, 'fieldname', 'solidearth.lovenumbers.tk2secular', 'NaN', 1, 'Inf', 1)
    61 
    62         #check that love numbers are provided at the same level of accuracy:
     60        # Check that love numbers are provided at the same level of accuracy
    6361        if (self.h.shape[0] != self.k.shape[0]) or (self.h.shape[0] != self.l.shape[0]):
    6462            raise ValueError('lovenumbers error message: love numbers should be provided at the same level of accuracy')
     63        return md
    6564    #}}}
    6665
    … …  
    7069
    7170    def __repr__(self): #{{{
    72         s = '   lovenumbers parameters:'
    73 
     71        s = '   lovenumbers parameters:\n'
    7472        s += '{}\n'.format(fielddisplay(self, 'h', 'load Love number for radial displacement'))
    7573        s += '{}\n'.format(fielddisplay(self, 'k', 'load Love number for gravitational potential perturbation'))
    … …  
    7876        s += '{}\n'.format(fielddisplay(self, 'tk', 'tidal load Love number (deg 2)'))
    7977        s += '{}\n'.format(fielddisplay(self, 'tk2secular', 'secular fluid Love number'))
    80 
    8178        return s
    8279    #}}}

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

    r24213 r25688  
    11import numpy as np
     2
     3from checkfield import checkfield
     4from fielddisplay import fielddisplay
     5from marshallcostfunctions import marshallcostfunctions
    26from project3d import project3d
    3 from fielddisplay import fielddisplay
    4 from checkfield import checkfield
    5 from WriteData import WriteData
    67from supportedcontrols import supportedcontrols
    78from supportedcostfunctions import supportedcostfunctions
    8 from marshallcostfunctions import marshallcostfunctions
     9from WriteData import WriteData
    910
    1011
    1112class m1qn3inversion(object):
    12     '''
    13     M1QN3 class definition
     13    """M1QN3 class definition
    1414
    15    Usage:
    16       m1qn3inversion = m1qn3inversion()
    17     '''
     15    Usage:
     16        m1qn3inversion = m1qn3inversion()
     17    """
    1818
    1919    def __init__(self, *args):  # {{{
    20 
    2120        if not len(args):
    2221            print('empty init')
    2322            self.iscontrol = 0
    2423            self.incomplete_adjoint = 0
    25             self.control_parameters = float('NaN')
    26             self.control_scaling_factors = float('NaN')
     24            self.control_parameters = np.nan
     25            self.control_scaling_factors = np.nan
    2726            self.maxsteps = 0
    2827            self.maxiter = 0
    2928            self.dxmin = 0.
    3029            self.gttol = 0.
    31             self.cost_functions = float('NaN')
    32             self.cost_functions_coefficients = float('NaN')
    33             self.min_parameters = float('NaN')
    34             self.max_parameters = float('NaN')
    35             self.vx_obs = float('NaN')
    36             self.vy_obs = float('NaN')
    37             self.vz_obs = float('NaN')
    38             self.vel_obs = float('NaN')
    39             self.thickness_obs = float('NaN')
     30            self.cost_functions = np.nan
     31            self.cost_functions_coefficients = np.nan
     32            self.min_parameters = np.nan
     33            self.max_parameters = np.nan
     34            self.vx_obs = np.nan
     35            self.vy_obs = np.nan
     36            self.vz_obs = np.nan
     37            self.vel_obs = np.nan
     38            self.thickness_obs = np.nan
    4039
    41             #set defaults
    4240            self.setdefaultparameters()
    4341        elif len(args) == 1 and args[0].__module__ == 'inversion':
    … …  
    6664
    6765    def __repr__(self):  # {{{
    68         string = '   m1qn3inversion parameters:'
    69         string = "%s\n%s" % (string, fielddisplay(self, 'iscontrol', 'is inversion activated?'))
    70         string = "%s\n%s" % (string, fielddisplay(self, 'incomplete_adjoint', '1: linear viscosity, 0: non - linear viscosity'))
    71         string = "%s\n%s" % (string, fielddisplay(self, 'control_parameters', 'ex: [''FrictionCoefficient''], or [''MaterialsRheologyBbar'']'))
    72         string = "%s\n%s" % (string, fielddisplay(self, 'control_scaling_factors', 'order of magnitude of each control (useful for multi - parameter optimization)'))
    73         string = "%s\n%s" % (string, fielddisplay(self, 'maxsteps', 'maximum number of iterations (gradient computation)'))
    74         string = "%s\n%s" % (string, fielddisplay(self, 'maxiter', 'maximum number of Function evaluation (forward run)'))
    75         string = "%s\n%s" % (string, fielddisplay(self, 'dxmin', 'convergence criterion: two points less than dxmin from eachother (sup - norm) are considered identical'))
    76         string = "%s\n%s" % (string, fielddisplay(self, 'gttol', '||g(X)||/||g(X0)|| (g(X0): gradient at initial guess X0)'))
    77         string = "%s\n%s" % (string, fielddisplay(self, 'cost_functions', 'indicate the type of response for each optimization step'))
    78         string = "%s\n%s" % (string, fielddisplay(self, 'cost_functions_coefficients', 'cost_functions_coefficients applied to the misfit of each vertex and for each control_parameter'))
    79         string = "%s\n%s" % (string, fielddisplay(self, 'min_parameters', 'absolute minimum acceptable value of the inversed parameter on each vertex'))
    80         string = "%s\n%s" % (string, fielddisplay(self, 'max_parameters', 'absolute maximum acceptable value of the inversed parameter on each vertex'))
    81         string = "%s\n%s" % (string, fielddisplay(self, 'vx_obs', 'observed velocity x component [m / yr]'))
    82         string = "%s\n%s" % (string, fielddisplay(self, 'vy_obs', 'observed velocity y component [m / yr]'))
    83         string = "%s\n%s" % (string, fielddisplay(self, 'vel_obs', 'observed velocity magnitude [m / yr]'))
    84         string = "%s\n%s" % (string, fielddisplay(self, 'thickness_obs', 'observed thickness [m]'))
    85         string = "%s\n%s" % (string, 'Available cost functions:')
    86         string = "%s\n%s" % (string, '   101: SurfaceAbsVelMisfit')
    87         string = "%s\n%s" % (string, '   102: SurfaceRelVelMisfit')
    88         string = "%s\n%s" % (string, '   103: SurfaceLogVelMisfit')
    89         string = "%s\n%s" % (string, '   104: SurfaceLogVxVyMisfit')
    90         string = "%s\n%s" % (string, '   105: SurfaceAverageVelMisfit')
    91         string = "%s\n%s" % (string, '   201: ThicknessAbsMisfit')
    92         string = "%s\n%s" % (string, '   501: DragCoefficientAbsGradient')
    93         string = "%s\n%s" % (string, '   502: RheologyBbarAbsGradient')
    94         string = "%s\n%s" % (string, '   503: ThicknessAbsGradient')
    95         return string
    96     #}}}
    97 
    98     def extrude(self, md):  # {{{
    99         self.vx_obs = project3d(md, 'vector', self.vx_obs, 'type', 'node')
    100         self.vy_obs = project3d(md, 'vector', self.vy_obs, 'type', 'node')
    101         self.vel_obs = project3d(md, 'vector', self.vel_obs, 'type', 'node')
    102         self.thickness_obs = project3d(md, 'vector', self.thickness_obs, 'type', 'node')
    103         if not np.any(np.isnan(self.cost_functions_coefficients)):
    104             self.cost_functions_coefficients = project3d(md, 'vector', self.cost_functions_coefficients, 'type', 'node')
    105         if not np.any(np.isnan(self.min_parameters)):
    106             self.min_parameters = project3d(md, 'vector', self.min_parameters, 'type', 'node')
    107         if not np.any(np.isnan(self.max_parameters)):
    108             self.max_parameters = project3d(md, 'vector', self.max_parameters, 'type', 'node')
    109         return self
     66        s = '   m1qn3inversion parameters:\n'
     67        s += '{}\n'.format(fielddisplay(self, 'iscontrol', 'is inversion activated?'))
     68        s += '{}\n'.format(fielddisplay(self, 'incomplete_adjoint', '1: linear viscosity, 0: non - linear viscosity'))
     69        s += '{}\n'.format(fielddisplay(self, 'control_parameters', 'ex: [''FrictionCoefficient''], or [''MaterialsRheologyBbar'']'))
     70        s += '{}\n'.format(fielddisplay(self, 'control_scaling_factors', 'order of magnitude of each control (useful for multi - parameter optimization)'))
     71        s += '{}\n'.format(fielddisplay(self, 'maxsteps', 'maximum number of iterations (gradient computation)'))
     72        s += '{}\n'.format(fielddisplay(self, 'maxiter', 'maximum number of Function evaluation (forward run)'))
     73        s += '{}\n'.format(fielddisplay(self, 'dxmin', 'convergence criterion: two points less than dxmin from eachother (sup - norm) are considered identical'))
     74        s += '{}\n'.format(fielddisplay(self, 'gttol', '||g(X)||/||g(X0)|| (g(X0): gradient at initial guess X0)'))
     75        s += '{}\n'.format(fielddisplay(self, 'cost_functions', 'indicate the type of response for each optimization step'))
     76        s += '{}\n'.format(fielddisplay(self, 'cost_functions_coefficients', 'cost_functions_coefficients applied to the misfit of each vertex and for each control_parameter'))
     77        s += '{}\n'.format(fielddisplay(self, 'min_parameters', 'absolute minimum acceptable value of the inversed parameter on each vertex'))
     78        s += '{}\n'.format(fielddisplay(self, 'max_parameters', 'absolute maximum acceptable value of the inversed parameter on each vertex'))
     79        s += '{}\n'.format(fielddisplay(self, 'vx_obs', 'observed velocity x component [m / yr]'))
     80        s += '{}\n'.format(fielddisplay(self, 'vy_obs', 'observed velocity y component [m / yr]'))
     81        s += '{}\n'.format(fielddisplay(self, 'vel_obs', 'observed velocity magnitude [m / yr]'))
     82        s += '{}\n'.format(fielddisplay(self, 'thickness_obs', 'observed thickness [m]'))
     83        s += '{}\n'.format('Available cost functions:')
     84        s += '{}\n'.format('   101: SurfaceAbsVelMisfit')
     85        s += '{}\n'.format('   102: SurfaceRelVelMisfit')
     86        s += '{}\n'.format('   103: SurfaceLogVelMisfit')
     87        s += '{}\n'.format('   104: SurfaceLogVxVyMisfit')
     88        s += '{}\n'.format('   105: SurfaceAverageVelMisfit')
     89        s += '{}\n'.format('   201: ThicknessAbsMisfit')
     90        s += '{}\n'.format('   501: DragCoefficientAbsGradient')
     91        s += '{}\n'.format('   502: RheologyBbarAbsGradient')
     92        s += '{}\n'.format('   503: ThicknessAbsGradient')
     93        return s
    11094    #}}}
    11195
    … …  
    130114    #}}}
    131115
     116    def extrude(self, md):  # {{{
     117        self.vx_obs = project3d(md, 'vector', self.vx_obs, 'type', 'node')
     118        self.vy_obs = project3d(md, 'vector', self.vy_obs, 'type', 'node')
     119        self.vel_obs = project3d(md, 'vector', self.vel_obs, 'type', 'node')
     120        self.thickness_obs = project3d(md, 'vector', self.thickness_obs, 'type', 'node')
     121        if not np.any(np.isnan(self.cost_functions_coefficients)):
     122            self.cost_functions_coefficients = project3d(md, 'vector', self.cost_functions_coefficients, 'type', 'node')
     123        if not np.any(np.isnan(self.min_parameters)):
     124            self.min_parameters = project3d(md, 'vector', self.min_parameters, 'type', 'node')
     125        if not np.any(np.isnan(self.max_parameters)):
     126            self.max_parameters = project3d(md, 'vector', self.max_parameters, 'type', 'node')
     127        return self
     128    #}}}
     129
    132130    def checkconsistency(self, md, solution, analyses):  # {{{
    133         #Early return
     131        # Early return
    134132        if not self.iscontrol:
    135133            return md
    … …  
    153151        if solution == 'BalancethicknessSolution':
    154152            md = checkfield(md, 'fieldname', 'inversion.thickness_obs', 'size', [md.mesh.numberofvertices], 'NaN', 1, 'Inf', 1)
     153        elif solution == 'BalancethicknessSoftSolution':
     154            md = checkfield(md, 'fieldname', 'inversion.thickness_obs', 'size', [md.mesh.numberofvertices], 'NaN', 1, 'Inf', 1)
    155155        else:
    156156            md = checkfield(md, 'fieldname', 'inversion.vx_obs', 'size', [md.mesh.numberofvertices], 'NaN', 1, 'Inf', 1)
    157157            md = checkfield(md, 'fieldname', 'inversion.vy_obs', 'size', [md.mesh.numberofvertices], 'NaN', 1, 'Inf', 1)
    158 
    159158        return md
    160159    # }}}
    … …  
    180179        WriteData(fid, prefix, 'object', self, 'class', 'inversion', 'fieldname', 'thickness_obs', 'format', 'DoubleMat', 'mattype', 1)
    181180
    182     #process control parameters
     181        # Process control parameters
    183182        num_control_parameters = len(self.control_parameters)
    184183        WriteData(fid, prefix, 'object', self, 'fieldname', 'control_parameters', 'format', 'StringArray')
    185184        WriteData(fid, prefix, 'data', num_control_parameters, 'name', 'md.inversion.num_control_parameters', 'format', 'Integer')
    186185
    187     #process cost functions
     186        # Process cost functions
    188187        num_cost_functions = np.size(self.cost_functions)
    189188        data = marshallcostfunctions(self.cost_functions)

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

    r24292 r25688  
    11import numpy as np
     2
     3from checkfield import checkfield
    24from fielddisplay import fielddisplay
    35from project3d import project3d
    4 from checkfield import checkfield
    56from WriteData import WriteData
    67
    78
    89class masstransport(object):
    9     """
    10     MASSTRANSPORT class definition
     10    """MASSTRANSPORT class definition
    1111
    12        Usage:
    13           masstransport = masstransport()
     12    Usage:
     13        masstransport = masstransport()
    1414    """
    1515
    … …  
    2424        self.requested_outputs = []
    2525
    26     #set defaults
     26        # Set defaults
    2727        self.setdefaultparameters()
    2828
    … …  
    3030
    3131    def __repr__(self):  # {{{
    32         string = '   Masstransport solution parameters:'
    33         string = "%s\n%s" % (string, fielddisplay(self, 'spcthickness', 'thickness constraints (NaN means no constraint) [m]'))
    34         string = "%s\n%s" % (string, fielddisplay(self, 'isfreesurface', 'do we use free surfaces (FS only) or mass conservation'))
    35         string = "%s\n%s" % (string, fielddisplay(self, 'min_thickness', 'minimum ice thickness allowed [m]'))
    36         string = "%s\n%s" % (string, fielddisplay(self, 'hydrostatic_adjustment', 'adjustment of ice shelves surface and bed elevations: ''Incremental'' or ''Absolute'' '))
    37         string = "%s\n%s" % (string, fielddisplay(self, 'stabilization', '0: no, 1: artificial_diffusivity, 2: streamline upwinding, 3: discontinuous Galerkin, 4: Flux Correction Transport'))
    38         string = "%s\n%s" % (string, fielddisplay(self, 'requested_outputs', 'additional outputs requested'))
    39 
    40         return string
     32        s = '   Masstransport solution parameters:\n'
     33        s += '{}\n'.format(fielddisplay(self, 'spcthickness', 'thickness constraints (NaN means no constraint) [m]'))
     34        s += '{}\n'.format(fielddisplay(self, 'isfreesurface', 'do we use free surfaces (FS only) or mass conservation'))
     35        s += '{}\n'.format(fielddisplay(self, 'min_thickness', 'minimum ice thickness allowed [m]'))
     36        s += '{}\n'.format(fielddisplay(self, 'hydrostatic_adjustment', 'adjustment of ice shelves surface and bed elevations: ''Incremental'' or ''Absolute'' '))
     37        s += '{}\n'.format(fielddisplay(self, 'stabilization', '0: no, 1: artificial_diffusivity, 2: streamline upwinding, 3: discontinuous Galerkin, 4: Flux Correction Transport'))
     38        s += '{}\n'.format(fielddisplay(self, 'requested_outputs', 'additional outputs requested'))
     39        return s
    4140    #}}}
    4241
    … …  
    5150    #}}}
    5251    def setdefaultparameters(self):  # {{{
    53         #Type of stabilization to use 0:nothing 1:artificial_diffusivity 3:Discontinuous Galerkin
     52        # Type of stabilization to use 0:nothing 1:artificial_diffusivity 3:Discontinuous Galerkin
    5453        self.stabilization = 1
    55         #Factor applied to compute the penalties kappa = max(stiffness matrix) * 1.0**penalty_factor
     54        # Factor applied to compute the penalties kappa = max(stiffness matrix) * 1.0**penalty_factor
    5655        self.penalty_factor = 3
    57         #Minimum ice thickness that can be used
     56        # Minimum ice thickness that can be used
    5857        self.min_thickness = 1
    59         #Hydrostatic adjustment
     58        # Hydrostatic adjustment
    6059        self.hydrostatic_adjustment = 'Absolute'
    61         #default output
     60        # Default output
    6261        self.requested_outputs = ['default']
    6362        return self
    … …  
    6564
    6665    def checkconsistency(self, md, solution, analyses):  # {{{
    67         #Early return
     66        # Early return
    6867        if ('MasstransportAnalysis' not in analyses) or (solution == 'TransientSolution' and not md.transient.ismasstransport):
    6968            return md
    … …  
    9089        WriteData(fid, prefix, 'object', self, 'fieldname', 'penalty_factor', 'format', 'Double')
    9190
    92     #process requested outputs
     91        # Process requested outputs
    9392        outputs = self.requested_outputs
    9493        indices = [i for i, x in enumerate(outputs) if x == 'default']

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

    r24261 r25688  
     1import numpy as np
     2
     3from checkfield import checkfield
    14from fielddisplay import fielddisplay
    25from project3d import project3d
    3 from checkfield import checkfield
    46from WriteData import WriteData
    5 import numpy as np
    67
    78
    89class mismipbasalforcings(object):
    9     """
    10     MISMIP Basal Forcings class definition
     10    """MISMIP Basal Forcings class definition
    1111
    1212    Usage:
    13     mismipbasalforcings = mismipbasalforcings()
     13        mismipbasalforcings = mismipbasalforcings()
    1414    """
    1515
    1616    def __init__(self):  # {{{
    17         self.groundedice_melting_rate = float('NaN')
    18         self.meltrate_factor = float('NaN')
    19         self.threshold_thickness = float('NaN')
    20         self.upperdepth_melt = float('NaN')
    21         self.geothermalflux = float('NaN')
     17        self.groundedice_melting_rate = np.nan
     18        self.meltrate_factor = np.nan
     19        self.threshold_thickness = np.nan
     20        self.upperdepth_melt = np.nan
     21        self.geothermalflux = np.nan
    2222        self.setdefaultparameters()
    23 
    2423    #}}}
    2524
    2625    def __repr__(self):  # {{{
    27         string = " MISMIP + basal melt parameterization\n"
    28         string = "%s\n%s" % (string, fielddisplay(self, "groundedice_melting_rate", "basal melting rate (positive if melting) [m / yr]"))
    29         string = "%s\n%s" % (string, fielddisplay(self, "meltrate_factor", "Melt - rate rate factor [1 / yr] (sign is opposite to MISMIP + benchmark to remain consistent with ISSM convention of positive values for melting)"))
    30         string = "%s\n%s" % (string, fielddisplay(self, "threshold_thickness", "Threshold thickness for saturation of basal melting [m]"))
    31         string = "%s\n%s" % (string, fielddisplay(self, "upperdepth_melt", "Depth above which melt rate is zero [m]"))
    32         string = "%s\n%s" % (string, fielddisplay(self, "geothermalflux", "Geothermal heat flux [W / m^2]"))
    33         return string
     26        s = '   MISMIP + basal melt parameterization\n'
     27        s += '{}\n'.format(fielddisplay(self, "groundedice_melting_rate", "basal melting rate (positive if melting) [m / yr]"))
     28        s += '{}\n'.format(fielddisplay(self, "meltrate_factor", "Melt - rate rate factor [1 / yr] (sign is opposite to MISMIP + benchmark to remain consistent with ISSM convention of positive values for melting)"))
     29        s += '{}\n'.format(fielddisplay(self, "threshold_thickness", "Threshold thickness for saturation of basal melting [m]"))
     30        s += '{}\n'.format(fielddisplay(self, "upperdepth_melt", "Depth above which melt rate is zero [m]"))
     31        s += '{}\n'.format(fielddisplay(self, "geothermalflux", "Geothermal heat flux [W / m^2]"))
     32        return s
    3433    #}}}
    3534
    … …  
    4342        if np.all(np.isnan(self.groundedice_melting_rate)):
    4443            self.groundedice_melting_rate = np.zeros((md.mesh.numberofvertices))
    45             print(' no basalforcings.groundedice_melting_rate specified: values set as zero')
     44            print('no basalforcings.groundedice_melting_rate specified: values set as zero')
    4645        if np.all(np.isnan(self.geothermalflux)):
    4746            self.geothermalflux = np.zeros((md.mesh.numberofvertices))
    … …  
    5958
    6059    def checkconsistency(self, md, solution, analyses):  # {{{
    61         #Early return
    62         if 'MasstransportAnalysis' in analyses and not (solution == 'TransientSolution' and md.transient.ismasstransport == 0):
     60        # Early return
     61        if 'MasstransportAnalysis' in analyses and not solution == 'TransientSolution' and not md.transient.ismasstransport:
    6362            md = checkfield(md, 'fieldname', 'basalforcings.groundedice_melting_rate', 'NaN', 1, 'Inf', 1, 'timeseries', 1)
    6463            md = checkfield(md, 'fieldname', 'basalforcings.meltrate_factor', '>=', 0, 'numel', [1])
    6564            md = checkfield(md, 'fieldname', 'basalforcings.threshold_thickness', '>=', 0, 'numel', [1])
    6665            md = checkfield(md, 'fieldname', 'basalforcings.upperdepth_melt', '<=', 0, 'numel', [1])
    67 
    6866        if 'BalancethicknessAnalysis' in analyses:
    6967            md = checkfield(md, 'fieldname', 'basalforcings.groundedice_melting_rate', 'NaN', 1, 'Inf', 1, 'size', [md.mesh.numberofvertices])
    … …  
    7169            md = checkfield(md, 'fieldname', 'basalforcings.threshold_thickness', '>=', 0, 'numel', [1])
    7270            md = checkfield(md, 'fieldname', 'basalforcings.upperdepth_melt', '<=', 0, 'numel', [1])
    73 
    74         if 'ThermalAnalysis' in analyses and not (solution == 'TransientSolution' and md.transient.isthermal == 0):
     71        if 'ThermalAnalysis' in analyses and not (solution == 'TransientSolution' and not md.transient.isthermal):
    7572            md = checkfield(md, 'fieldname', 'basalforcings.groundedice_melting_rate', 'NaN', 1, 'Inf', 1, 'timeseries', 1)
    7673            md = checkfield(md, 'fieldname', 'basalforcings.meltrate_factor', '>=', 0, 'numel', [1])
    … …  
    8582        if yts != 365.2422 * 24. * 3600.:
    8683            print('WARNING: value of yts for MISMIP + runs different from ISSM default!')
    87 
    8884        WriteData(fid, prefix, 'name', 'md.basalforcings.model', 'data', 3, 'format', 'Integer')
    8985        WriteData(fid, prefix, 'object', self, 'fieldname', 'groundedice_melting_rate', 'format', 'DoubleMat', 'name', 'md.basalforcings.groundedice_melting_rate', 'mattype', 1, 'scale', 1. / yts, 'timeserieslength', md.mesh.numberofvertices + 1, 'yts', yts)

  • issm/trunk-jpl/src/m/classes/nodalvalue.m

    r21445 r25688  
    6060                function md = marshall(self,prefix,md,fid) % {{{
    6161
    62                 WriteData(fid,prefix,'data',self.name,'name','md.nodalvalue.name','format','String');
    63                 WriteData(fid,prefix,'data',self.definitionstring,'name','md.nodalvalue.definitionenum','format','String');
    64                 WriteData(fid,prefix,'data',self.model_string,'name','md.nodalvalue.model_enum','format','String');
    65                 WriteData(fid,prefix,'data',self.node,'name','md.nodalvalue.node','format','Integer');
     62                        WriteData(fid,prefix,'data',self.name,'name','md.nodalvalue.name','format','String');
     63                        WriteData(fid,prefix,'data',self.definitionstring,'name','md.nodalvalue.definitionenum','format','String');
     64                        WriteData(fid,prefix,'data',self.model_string,'name','md.nodalvalue.model_enum','format','String');
     65                        WriteData(fid,prefix,'data',self.node,'name','md.nodalvalue.node','format','Integer');
    6666
    6767                end % }}}

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

    r25169 r25688  
    33from checkfield import checkfield
    44from fielddisplay import fielddisplay
     5from pairoptions import pairoptions
    56from WriteData import WriteData
    67
    78
    8 class dslmme(object):
    9     '''
    10     NODALVALUE class definition
     9class nodalvalue(object):
     10    """NODALVALUE class definition
    1111
    12         Usage:
    13             nodalvalue=nodalvalue()
    14             nodalvalue=nodalvalue(
    15                 'name', 'SealevelriseSNodalValue',
    16                 'definitionstring', 'Outputdefinition1',
    17                 'model_string', 'SealevelriseS',
    18                 'node', 1
    19                 )
    20     '''
     12    Usage:
     13        nodalvalue=nodalvalue()
     14        nodalvalue=nodalvalue(
     15            'name', 'SealevelriseSNodalValue',
     16            'definitionstring', 'Outputdefinition1',
     17            'model_string', 'SealevelriseS',
     18            'node', 1
     19        )
     20    """
    2121
    2222    def __init__(self, *args): #{{{
    2323        self.name               = ''
    24         self.definitionstring   = '' #string that identifies this output definition uniquely, from 'Outputdefinition[1-10]'
    25         self.model_string       = '' #string for field that is being retrieved
     24        self.definitionstring   = '' # string that identifies this output definition uniquely, from 'Outputdefinition[1-10]'
     25        self.model_string       = '' # string for field that is being retrieved
    2626        self.node               = np.nan #for which node are we retrieving the value?
    2727       
    … …  
    2929        options = pairoptions(*args)
    3030
    31         #get name
    32         self.name               = options.getfieldvalue(options, 'name', '')
    33         self.definitionstring   = options.getfieldvalue(options, 'definitionstring', '')
    34         self.model_string       = options.getfieldvalue(options, 'model_string', '')
    35         self.node               = options.getfieldvalue(options, 'node', '')
     31        # Get name
     32        self.name               = options.getfieldvalue('name', '')
     33        self.definitionstring   = options.getfieldvalue('definitionstring', '')
     34        self.model_string       = options.getfieldvalue('model_string', '')
     35        self.node               = options.getfieldvalue('node', '')
    3636    #}}}
    3737
    … …  
    4242        s += '{}\n'.format(fielddisplay(self, 'model_string', 'string for field that is being retrieved'))
    4343        s += '{}\n'.format(fielddisplay(self, 'node', 'vertex index at which we retrieve the value'))
    44 
    4544        return s
    4645    #}}}
    4746
    4847    def setdefaultparameters(self): # {{{
    49         return
     48        return self
    5049    #}}}
    5150
    … …  
    5554        OutputdefinitionStringArray = []
    5655        for i in range(100):
    57             OutputdefinitionStringArray[i] = 'Outputdefinition{}'.format(i)
     56            OutputdefinitionStringArray.append('Outputdefinition{}'.format(i))
    5857        md = checkfield(md, 'fieldname', 'self.definitionstring', 'field', self.definitionstring, 'values', OutputdefinitionStringArray)
    5958        md = checkfield(md, 'fieldname', 'self.node', 'field', self.node, 'values', range(md.mesh.numberofvertices))
    60 
    6159        return md
    6260    #}}}

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

    r24261 r25688  
    11import numpy as np
     2
     3from checkfield import checkfield
    24from fielddisplay import fielddisplay
    3 from checkfield import checkfield
     5from project3d import project3d
     6from structtoobj import structtoobj
    47from WriteData import WriteData
    5 from project3d import project3d
    68
    79
    810class plumebasalforcings(object):
    9     '''
    10     PLUME BASAL FORCINGS class definition
     11    """PLUME BASAL FORCINGS class definition
    1112
    12         Usage:
    13             plumebasalforcings = plumebasalforcings()
    14     '''
     13    Usage:
     14        plumebasalforcings = plumebasalforcings()
     15    """
    1516
    16     def __init__(self):  # {{{
    17         self.floatingice_melting_rate = float('NaN')
    18         self.groundedice_melting_rate = float('NaN')
    19         self.mantleconductivity = float('NaN')
    20         self.nusselt = float('NaN')
    21         self.dtbg = float('NaN')
    22         self.plumeradius = float('NaN')
    23         self.topplumedepth = float('NaN')
    24         self.bottomplumedepth = float('NaN')
    25         self.plumex = float('NaN')
    26         self.plumey = float('NaN')
    27         self.crustthickness = float('NaN')
    28         self.uppercrustthickness = float('NaN')
    29         self.uppercrustheat = float('NaN')
    30         self.lowercrustheat = float('NaN')
     17    def __init__(self, *args):  # {{{
     18        self.floatingice_melting_rate   = np.nan
     19        self.groundedice_melting_rate   = np.nan
     20        self.mantleconductivity         = np.nan
     21        self.nusselt                    = np.nan
     22        self.dtbg                       = np.nan
     23        self.plumeradius                = np.nan
     24        self.topplumedepth              = np.nan
     25        self.bottomplumedepth           = np.nan
     26        self.plumex                     = np.nan
     27        self.plumey                     = np.nan
     28        self.crustthickness             = np.nan
     29        self.uppercrustthickness        = np.nan
     30        self.uppercrustheat             = np.nan
     31        self.lowercrustheat             = np.nan
    3132
    32         self.setdefaultparameters()
     33        nargs = len(args)
     34        if nargs == 0:
     35            self.setdefaultparameters()
     36        elif nargs == 1:
     37            # TODO: This has not been tested
     38            self = structtoobj(self, args[0])
     39        else:
     40            error('constuctor not supported')
    3341    #}}}
    3442
    3543    def __repr__(self):  # {{{
    36         string = '   mantle plume basal melt parameterization:'
    37 
    38         string = "%s\n%s" % (string, fielddisplay(self, 'groundedice_melting_rate', 'basal melting rate (positive if melting) [m / yr]'))
    39         string = "%s\n%s" % (string, fielddisplay(self, 'floatingice_melting_rate', 'basal melting rate (positive if melting) [m / yr]'))
    40         string = "%s\n%s" % (string, fielddisplay(self, 'mantleconductivity', 'mantle heat conductivity [W / m^3]'))
    41         string = "%s\n%s" % (string, fielddisplay(self, 'nusselt', 'nusselt number, ratio of mantle to plume [1]'))
    42         string = "%s\n%s" % (string, fielddisplay(self, 'dtbg', 'background temperature gradient [degree / m]'))
    43         string = "%s\n%s" % (string, fielddisplay(self, 'plumeradius', 'radius of the mantle plume [m]'))
    44         string = "%s\n%s" % (string, fielddisplay(self, 'topplumedepth', 'depth of the mantle plume top below the crust [m]'))
    45         string = "%s\n%s" % (string, fielddisplay(self, 'bottomplumedepth', 'depth of the mantle plume base below the crust [m]'))
    46         string = "%s\n%s" % (string, fielddisplay(self, 'plumex', 'x coordinate of the center of the plume [m]'))
    47         string = "%s\n%s" % (string, fielddisplay(self, 'plumey', 'y coordinate of the center of the plume [m]'))
    48         string = "%s\n%s" % (string, fielddisplay(self, 'crustthickness', 'thickness of the crust [m]'))
    49         string = "%s\n%s" % (string, fielddisplay(self, 'uppercrustthickness', 'thickness of the upper crust [m]'))
    50         string = "%s\n%s" % (string, fielddisplay(self, 'uppercrustheat', 'volumic heat of the upper crust [w / m^3]'))
    51         string = "%s\n%s" % (string, fielddisplay(self, 'lowercrustheat', 'volumic heat of the lowercrust [w / m^3]'))
    52 
    53         return string
     44        s = '   mantle plume basal melt parameterization:\n'
     45        s += '{}\n'.format(fielddisplay(self, 'groundedice_melting_rate', 'basal melting rate (positive if melting) [m / yr]'))
     46        s += '{}\n'.format(fielddisplay(self, 'floatingice_melting_rate', 'basal melting rate (positive if melting) [m / yr]'))
     47        s += '{}\n'.format(fielddisplay(self, 'mantleconductivity', 'mantle heat conductivity [W / m^3]'))
     48        s += '{}\n'.format(fielddisplay(self, 'nusselt', 'nusselt number, ratio of mantle to plume [1]'))
     49        s += '{}\n'.format(fielddisplay(self, 'dtbg', 'background temperature gradient [degree / m]'))
     50        s += '{}\n'.format(fielddisplay(self, 'plumeradius', 'radius of the mantle plume [m]'))
     51        s += '{}\n'.format(fielddisplay(self, 'topplumedepth', 'depth of the mantle plume top below the crust [m]'))
     52        s += '{}\n'.format(fielddisplay(self, 'bottomplumedepth', 'depth of the mantle plume base below the crust [m]'))
     53        s += '{}\n'.format(fielddisplay(self, 'plumex', 'x coordinate of the center of the plume [m]'))
     54        s += '{}\n'.format(fielddisplay(self, 'plumey', 'y coordinate of the center of the plume [m]'))
     55        s += '{}\n'.format(fielddisplay(self, 'crustthickness', 'thickness of the crust [m]'))
     56        s += '{}\n'.format(fielddisplay(self, 'uppercrustthickness', 'thickness of the upper crust [m]'))
     57        s += '{}\n'.format(fielddisplay(self, 'uppercrustheat', 'volumic heat of the upper crust [w / m^3]'))
     58        s += '{}\n'.format(fielddisplay(self, 'lowercrustheat', 'volumic heat of the lowercrust [w / m^3]'))
     59        return s
    5460    #}}}
    5561
    … …  
    6167            self.floatingice_melting_rate = np.zeros((md.mesh.numberofvertices, ))
    6268            print('      no basalforcings.floatingice_melting_rate specified: values set as zero')
    63         return
     69        return self
    6470    #}}}
    6571
    … …  
    7177
    7278    def setdefaultparameters(self):  # {{{
    73         #default values for melting parameterization
     79        # Default values for melting parameterization
    7480        self.mantleconductivity = 2.2
    7581        self.nusselt = 300

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

    r25685 r25688  
    1212
    1313class qmu(object):
    14     """
    15     QMU class definition
    16 
    17        Usage:
    18           qmu = qmu()
     14    """QMU class definition
     15
     16    Usage:
     17        qmu = qmu()
    1918    """
    2019
    … …  
    4443        self.vertex_weight = float('NaN')
    4544
    46     #set defaults
    4745        self.setdefaultparameters()
    48 
    4946    #}}}
    5047    def __repr__(self):  # {{{
    5148        s = '   qmu parameters:\n'
    52 
    53         s += "%s\n" % fielddisplay(self, 'isdakota', 'is qmu analysis activated?')
    54         s += "%s\n" % fielddisplay(self, 'output', 'are we outputting ISSM results, default is 0')
     49        s += '{}\n'.format(fielddisplay(self, 'isdakota', 'is QMU analysis activated?'))
     50        s += '{}\n'.format(fielddisplay(self, 'output', 'are we outputting ISSM results, default is 0'))
    5551        maxlen = 0
    56         s += "         variables:  (arrays of each variable class)\n"
     52        s += '         variables:  (arrays of each variable class)\n'
    5753
    5854    # OrderedStruct's iterator returns individual name / array - of - functions pairs
    … …  
    7470            s += "            %-*s:    [%ix%i]    '%s'\n" % (maxlen + 1, fname, a, b, type(response[1][0]))
    7571
    76         s += "%s\n" % fielddisplay(self, 'numberofresponses', 'number of responses')
     72        s += '{}\n'.format(fielddisplay(self, 'numberofresponses', 'number of responses'))
    7773
    7874        if type(self.method) != OrderedDict:
    … …  
    9187            print(type(param))
    9288            print(param)
    93             s += "         params:  (array of method - independent parameters)\n"
     89            s += "         params:  (array of method-independent parameters)\n"
    9490            fnames = vars(param)
    9591            maxlen = 0
    … …  
    116112                s += "            %-*s:    [%ix%i]    '%s'\n" % (maxlen + 1, fname, a, b, type(getattr(result, fname)))
    117113
    118         s += "%s\n" % fielddisplay(self, 'variablepartitions', '')
    119         s += "%s\n" % fielddisplay(self, 'variablepartitions_npart', '')
    120         s += "%s\n" % fielddisplay(self, 'variablepartitions_nt', '')
    121         s += "%s\n" % fielddisplay(self, 'variabledescriptors', '')
    122         s += "%s\n" % fielddisplay(self, 'responsedescriptors', '')
    123         s += "%s\n" % fielddisplay(self, 'method', 'array of dakota_method class')
    124         s += "%s\n" % fielddisplay(self, 'mass_flux_profile_directory', 'directory for mass flux profiles')
    125         s += "%s\n" % fielddisplay(self, 'mass_flux_profiles', 'list of mass_flux profiles')
    126         s += "%s\n" % fielddisplay(self, 'mass_flux_segments', '')
    127         s += "%s\n" % fielddisplay(self, 'adjacency', '')
    128         s += "%s\n" % fielddisplay(self, 'vertex_weight', 'weight applied to each mesh vertex')
    129 
     114        s += '{}\n'.format(fielddisplay(self, 'variablepartitions', ''))
     115        s += '{}\n'.format(fielddisplay(self, 'variablepartitions_npart', ''))
     116        s += '{}\n'.format(fielddisplay(self, 'variablepartitions_nt', ''))
     117        s += '{}\n'.format(fielddisplay(self, 'variabledescriptors', ''))
     118        s += '{}\n'.format(fielddisplay(self, 'responsedescriptors', ''))
     119        s += '{}\n'.format(fielddisplay(self, 'method', 'array of dakota_method class'))
     120        s += '{}\n'.format(fielddisplay(self, 'mass_flux_profile_directory', 'directory for mass flux profiles'))
     121        s += '{}\n'.format(fielddisplay(self, 'mass_flux_profiles', 'list of mass_flux profiles'))
     122        s += '{}\n'.format(fielddisplay(self, 'mass_flux_segments', ''))
     123        s += '{}\n'.format(fielddisplay(self, 'adjacency', ''))
     124        s += '{}\n'.format(fielddisplay(self, 'vertex_weight', 'weight applied to each mesh vertex'))
    130125        return s
    131126    # }}}

  • issm/trunk-jpl/src/m/classes/qmu/normal_uncertain.m

    r25222 r25688  
    161161                end % }}}
    162162                %new methods:
    163                         function distributed=isdistributed(self) % {{{
     163                function distributed=isdistributed(self) % {{{
    164164                        if strncmp(self.descriptor,'distributed_',12),
    165165                                distributed=1;
    … …  
    168168                        end
    169169                end % }}}
    170         function scaled=isscaled(self) % {{{
     170                function scaled=isscaled(self) % {{{
    171171                        if strncmp(self.descriptor,'scaled_',7),
    172172                                scaled=1;

  • issm/trunk-jpl/src/m/classes/qmu/normal_uncertain.py

    r25097 r25688  
    1010
    1111class normal_uncertain(object):
    12     '''
    13     NORMAL_UNCERTAIN class definition
     12    """NORMAL_UNCERTAIN class definition
    1413
    1514    Usage:
    … …  
    3837            'partition',vpartition
    3938            )
    40     '''
     39    """
     40
    4141    def __init__(self): #{{{
    4242        self.descriptor = ''
    … …  
    231231
    232232    #new methods:
     233    def isdistributed(self): #{{{
     234        if strncmp(self.descriptor, 'distributed_', 12):
     235            return True
     236        else:
     237            return False
     238    #}}}
     239   
    233240    def isscaled(self): #{{{
    234241        if strncmp(self.descriptor, 'scaled_', 7):

  • issm/trunk-jpl/src/m/classes/qmu/response_function.m

    r24998 r25688  
    116116        end % }}}
    117117                %new methods:
    118         function scaled =isscaled(self) % {{{
    119                 if strncmp(self.descriptor,'scaled_',7),
    120                         scaled=1;
    121                 else
    122                         scaled=0;
    123                 end
    124         end % }}}
     118        function scaled =isscaled(self) % {{{
     119                if strncmp(self.descriptor,'scaled_',7),
     120                        scaled=1;
     121                else
     122                        scaled=0;
     123                end
     124        end % }}}
    125125    end
    126126    methods (Static)

  • issm/trunk-jpl/src/m/classes/qmustatistics.m

    r25649 r25688  
    117117                                for i=1:length(self.method),
    118118                                        m=self.method(i);
    119                                         WriteData(fid,prefix,'name',sprintf('md.qmu.statistics.method(%i).name',i),'data',m.name,'Format','String');
     119                                        WriteData(fid,prefix,'name',sprintf('md.qmu.statistics.method(%i).name',i),'data',m.name,'format','String');
    120120                                        WriteData(fid,prefix,'data',m.fields,'name',sprintf('md.qmu.statistics.method(%i).fields',i),'format','StringArray');
    121121                                        WriteData(fid,prefix,'data',m.steps,'name',sprintf('md.qmu.statistics.method(%i).steps',i),'format','IntMat','mattype',3);

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

    r25685 r25688  
    2626        # fields : fields for the statistics being requested, ex: 'Sealevel', 'BslrIce', 'BsrlHydro'
    2727        # steps : time steps at which each field statistic is computed, ex: [1, 2, 5, 20] or [range(1:100)]
    28         # nbins : number of bins for 'Histgogram' statistics
     28        # nbins : number of bins for 'Histogram' statistics
    2929        # indices : vertex indices at which to retrieve samples
    3030
    3131        nargs = len(args)
    32 
    3332        if nargs == 0:
    3433            # Create a default object
    3534            self.setdefaultparameters()
    3635        elif nargs == 1:
    37             # NOTE: The following has not been tested. Remove this note when it has.
     36            # NOTE: The following has not been tested. Remove this note when it has
    3837            inputstruct = args[0]
    3938            list1 = properties('qmustatistics')
    … …  
    4847
    4948    def __repr__(self):  # {{{
    50         s = '{}\n'.format('qmustatistics: post-Dakota run processing of QMU statistics:')
    51 
     49        s = 'qmustatistics: post-Dakota run processing of QMU statistics:\n'
    5250        if self.method[0]['name'] == 'None':
    5351            return s
    54 
    5552        s += '{}\n'.format(fielddisplay(self, 'nfiles_per_directory', 'number of files per output directory'))
    5653        s += '{}\n'.format(fielddisplay(self, 'ndirectories', 'number of output directories; should be < numcpus'))
    57 
    5854        for i in range(len(self.method)):
    5955            s += '{}\n'.format('   method # {}'.format(i))
    6056            s += '{}\n'.format(self.method[i])
    61 
    6257        return s
    6358    #}}}
    … …  
    6762        self.nfiles_per_directory   = 5 # Number of files per output directory
    6863        self.ndirectories           = 50 # Number of output directories; should be < numcpus
     64        return self
    6965    #}}}
    7066
    … …  
    115111            WriteData(fid, prefix, 'name', 'md.qmu.statistics.ndirectories', 'data', self.ndirectories, 'format', 'Integer')
    116112            WriteData(fid, prefix, 'name', 'md.qmu.statistics.numstatistics', 'data', len(self.method), 'format', 'Integer')
    117             for i in range(len(self.method)):
    118                 m = self.method[i]
    119                 WriteData(fid, prefix, 'name', 'md.qmu.statistics.method[{}][\'name\']'.format(i), 'data', m['name'], 'Format', 'String')
    120                 WriteData(fid, prefix, 'data', m['fields'], 'name', 'md.qmu.statistics.method[{}][\'fields\']'.format(i), 'format', 'StringArray')
    121                 WriteData(fid, prefix, 'data', m['steps'], 'name', 'md.qmu.statistics.method[{}][\'steps\']'.format(i), 'format', 'IntMat', 'mattype', 3)
    122 
     113            for i in range(1, len(self.method) + 1):
     114                m = self.method[i - 1]
     115                WriteData(fid, prefix, 'name', 'md.qmu.statistics.method({}).name'.format(i), 'data', m['name'], 'format', 'String')
     116                WriteData(fid, prefix, 'data', m['fields'], 'name', 'md.qmu.statistics.method({}).fields'.format(i), 'format', 'StringArray')
     117                WriteData(fid, prefix, 'data', m['steps'], 'name', 'md.qmu.statistics.method({}).steps'.format(i), 'format', 'IntMat', 'mattype', 3)
    123118                if m['name'] == 'Histogram':
    124                     WriteData(fid, prefix, 'name', 'md.qmu.statistics.method[{}][\'nbins\']'.format(i), 'data', m['nbins'], 'format', 'Integer')
     119                    WriteData(fid, prefix, 'name', 'md.qmu.statistics.method({}).nbins'.format(i), 'data', m['nbins'], 'format', 'Integer')
    125120                elif m['name'] == 'MeanVariance':
    126121                    pass # do nothing
    127122                elif m['name'] == 'SampleSeries':
    128                     WriteData(fid, prefix, 'data', m['indices'], 'name', 'md.qmu.statistics.method[{}][\'indices\']'.format(i), 'format', 'IntMat', 'mattype', 3)
     123                    WriteData(fid, prefix, 'data', m['indices'], 'name', 'md.qmu.statistics.method({}).indices'.format(i), 'format', 'IntMat', 'mattype', 3)
    129124                else:
    130125                    raise Exception('qmustatistics marshall error message: unknown type ''{}'' for qmu.statistics.method[{}]'.format(m['name'], i))
    … …  
    134129        return self
    135130    #}}}
     131
     132    def addmethod(self, *args): #{{{
     133        """ADDMETHOD - Add new, empty method or passed dict to self.method
     134        """
     135        nargs = len(args)
     136        if nargs == 0:
     137            self.method.append({})
     138        elif nargs == 1:
     139            self.method.append(args[0])
     140        else:
     141            raise Exception('Number of args should be 0 (appends empty dict to methods member) or 1 (appends passed dict to methods member)')
     142    #}}}

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

    r25168 r25688  
    33
    44class results(object):
    5     '''
    6     RESULTS class definition
     5    """RESULTS class definition
    76
    8        Usage:
    9           results = results()
     7    Usage:
     8        results = results()
    109
    11         TODO:
    12         - Modify output so that it matches that of
     10    TODO:
     11    - Modify output so that it matches that of
    1312
    14             disp(md.results.<<solutionstring>>)
     13        disp(md.results.<<solutionstring>>)
    1514
    16         where <<solutionstring>> is one of the values from solve.m
    17     '''
     15    where <<solutionstring>> is one of the values from solve.m
     16    """
    1817
    1918    def __init__(self, *args):  # {{{

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

    r25171 r25688  
    77
    88class rotational(object):
    9     '''
    10     ROTATIONAL class definition
     9    """ROTATIONAL class definition
    1110
    12         Usage:
    13             rotational = rotational()
    14     '''
     11    Usage:
     12        rotational = rotational()
     13    """
    1514
    1615    def __init__(self, *args): #{{{
    … …  
    2019
    2120        nargin = len(args)
    22 
    2321        if nargin == 0:
    2422            self.setdefaultparameters()
    … …  
    3230        s += '{}\n'.format(fielddisplay(self, 'polarmoi', 'polar moment of inertia [kg m^2]'))
    3331        s += '{}\n'.format(fielddisplay(self, 'angularvelocity', 'mean rotational velocity of earth [per second]'))
    34 
    3532        return s
    3633    #}}}
    3734
    3835    def setdefaultparameters(self): # {{{
    39         #moment of inertia
     36        # Moment of inertia
    4037        self.equatorialmoi  = 8.0077e37 # [kg m^2]
    4138        self.polarmoi       = 8.0345e37 # [kg m^2]
    4239
    43         #mean rotational velocity of earth
     40        # Mean rotational velocity of earth
    4441        self.angularvelocity = 7.2921e-5 # [s^-1]
     42        return self
    4543    #}}}
    4644
    4745    def checkconsistency(self, md, solution, analyses): # {{{
    48         if ('SealevelriseAnalysis' not in analyses) or (solution == 'TransientSolution' and md.transient.isslr == 0):
     46        if ('SealevelriseAnalysis' not in analyses) or (solution == 'TransientSolution' and not md.transient.isslr):
    4947            return md
    50 
    5148        md = checkfield(md, 'fieldname', 'solidearth.rotational.equatorialmoi', 'NaN', 1, 'Inf', 1)
    5249        md = checkfield(md, 'fieldname', 'solidearth.rotational.polarmoi', 'NaN', 1, 'Inf', 1)
    5350        md = checkfield(md, 'fieldname', 'solidearth.rotational.angularvelocity', 'NaN', 1, 'Inf', 1)
    54 
    5551        return md
    5652    #}}}

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

    r25677 r25688  
    2727        methods
    2828                function slm = sealevelmodel(varargin) % {{{
    29 
    30                         if nargin==0,
    31                                 slm=setdefaultparameters(slm);
    32                         else
    33                                 slm=setdefaultparameters(slm);
    34 
    35                                 options=pairoptions(varargin{:});
    36                        
     29                        slm=setdefaultparameters(slm);
     30
     31                        if nargin==1,
     32
     33                                options=pairoptions(varargin{:});
     34
    3735                                %recover all the icecap models:
    3836                                slm.icecaps=getfieldvalues(options,'ice_cap',{});
    … …  
    419417
    420418                end % }}}
    421 function viscousiterations(self) % {{{
    422         for  i=1:length(self.icecaps),
    423                 ic=self.icecaps{i};
    424                 mvi=ic.results.TransientSolution(1).StressbalanceConvergenceNumSteps;
    425                 for j=2:length(ic.results.TransientSolution)-1,
    426                         mvi=max(mvi,ic.results.TransientSolution(j).StressbalanceConvergenceNumSteps);
    427                 end
    428                 disp(sprintf('%i, %s: %i',i,self.icecaps{i}.miscellaneous.name,mvi));
    429         end
    430 end % }}}
     419                function viscousiterations(self) % {{{
     420                        for  i=1:length(self.icecaps),
     421                                ic=self.icecaps{i};
     422                                mvi=ic.results.TransientSolution(1).StressbalanceConvergenceNumSteps;
     423                                for j=2:length(ic.results.TransientSolution)-1,
     424                                        mvi=max(mvi,ic.results.TransientSolution(j).StressbalanceConvergenceNumSteps);
     425                                end
     426                                disp(sprintf('%i, %s: %i',i,self.icecaps{i}.miscellaneous.name,mvi));
     427                        end
     428                end % }}}
    431429                function maxtimestep(self) % {{{
    432430                        for  i=1:length(self.icecaps),
    433                                 ic=self.icecaps{i}; 
     431                                ic=self.icecaps{i};
    434432                                mvi=length(ic.results.TransientSolution);
    435433                                timei=ic.results.TransientSolution(end).time;

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

    r25679 r25688  
    11from copy import deepcopy
    22import warnings
     3
    34import numpy as np
     5
    46from fielddisplay import fielddisplay
    57from generic import generic
    … …  
    4749        self.planet = ''
    4850
    49         # Create a default object
    5051        self.setdefaultparameters()
    5152
    5253        if len(args):
    53             # Use provided options to set fields
    5454            options = pairoptions(*args)
    5555
    … …  
    9090        # Is the coupler turned on?
    9191        for i in range(len(slm.icecaps)):
    92             if slm.icecaps[i].transient.iscoupler == 0:
     92            if not slm.icecaps[i].transient.iscoupler:
    9393                warnings.warn('sealevelmodel checkconsistency error: icecap model {} should have the transient coupler option turned on!'.format(slm.icecaps[i].miscellaneous.name))
    9494
    95         if slm.earth.transient.iscoupler == 0:
     95        if not slm.earth.transient.iscoupler:
    9696            warnings.warn('sealevelmodel checkconsistency error: earth model should have the transient coupler option turned on!')
    9797
    … …  
    9999        for i in range(len(slm.icecaps)):
    100100            if slm.icecaps[i].mesh.numberofvertices != len(slm.earth.slr.transitions[i]):
    101                 raise RuntimeError('sealevelmodel checkconsistency issue with size of transition vector for ice cap: {} name: {}'.format(i, slm.icecaps[i].miscellaneous.name))
     101                raise Exception('sealevelmodel checkconsistency issue with size of transition vector for ice cap: {} name: {}'.format(i, slm.icecaps[i].miscellaneous.name))
    102102
    103103        # Check that run frequency is the same everywhere
    104104        for i in range(len(slm.icecaps)):
    105105            if slm.icecaps[i].slr.geodetic_run_frequency != slm.earth.geodetic_run_frequency:
    106                 raise RuntimeError('sealevelmodel checkconsistency error: icecap model {} should have the same run frequency as earth!'.format(slm.icecaps[i].miscellaneous.name))
     106                raise Exception('sealevelmodel checkconsistency error: icecap model {} should have the same run frequency as earth!'.format(slm.icecaps[i].miscellaneous.name))
    107107
    108108        # Make sure steric_rate is the same everywhere
    … …  
    110110            md = slm.icecaps[i]
    111111            if np.nonzero(md.slr.steric_rate - slm.earth.slr.steric_rate[slm.earth.slr.transitions[i]]) != []:
    112                 raise RuntimeError('steric rate on ice cap {} is not the same as for the earth'.format(md.miscellaneous.name))
     112                raise Exception('steric rate on ice cap {} is not the same as for the earth'.format(md.miscellaneous.name))
    113113    #}}}
    114114
    … …  
    168168    def addbasin(self, bas):  # {{{
    169169        if bas.__class__.__name__ != 'basin':
    170             raise RuntimeError('addbasin method only takes a \'basin\' class object as input')
     170            raise Exception('addbasin method only takes a \'basin\' class object as input')
    171171        self.basins.append(bas)
    172172    #}}}
    … …  
    387387        for i in range(len(self.icecaps)):
    388388            ic = self.icecaps[i]
    389             mvi = ic.resutls.TransientSolution[0].StressbalanceConvergenceNumSteps
     389            mvi = ic.results.TransientSolution[0].StressbalanceConvergenceNumSteps
    390390            for j in range(1, len(ic.results.TransientSolution) - 1):
    391391                mvi = np.max(mvi, ic.results.TransientSolution[j].StressbalanceConvergenceNumSteps)
    … …  
    465465
    466466        if not noearth:
    467             ic.results.TransientSolution = ic.resutls.TransientSolution[:mintimestep]
     467            ic.results.TransientSolution = ic.results.TransientSolution[:mintimestep]
    468468
    469469        self.earth = ic

  • issm/trunk-jpl/src/m/classes/slr.m

    r25168 r25688  
    136136                        %end
    137137
    138                         %check that  if geodetic is requested, we are a mesh3dsurface model (planet), or if we are not,
     138                        %check that if geodetic is requested, we are a mesh3dsurface model (planet), or if we are not,
    139139                        %a coupler to a planet model is provided.
    140140                        if self.geodetic,

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

    r25171 r25688  
    5050
    5151    def __repr__(self):  # {{{
    52         s = '   slr parameters:'
    53         s = "%s\n%s" % (s, fielddisplay(self, 'deltathickness', 'thickness change: ice height equivalent [m]'))
    54         s = "%s\n%s" % (s, fielddisplay(self, 'sealevel', 'current sea level (prior to computation) [m]'))
    55         s = "%s\n%s" % (s, fielddisplay(self, 'spcthickness', 'thickness constraints (NaN means no constraint) [m]'))
    56         s = "%s\n%s" % (s, fielddisplay(self, 'reltol', 'sea level rise relative convergence criterion, (NaN: not applied)'))
    57         s = "%s\n%s" % (s, fielddisplay(self, 'abstol', 'sea level rise absolute convergence criterion, (default, NaN: not applied)'))
    58         s = "%s\n%s" % (s, fielddisplay(self, 'maxiter', 'maximum number of nonlinear iterations'))
    59         s = "%s\n%s" % (s, fielddisplay(self, 'love_h', 'load Love number for radial displacement'))
    60         s = "%s\n%s" % (s, fielddisplay(self, 'love_k', 'load Love number for gravitational potential perturbation'))
    61         s = "%s\n%s" % (s, fielddisplay(self, 'love_l', 'load Love number for horizontal displaements'))
    62         s = "%s\n%s" % (s, fielddisplay(self, 'tide_love_k', 'tidal load Love number (degree 2)'))
    63         s = "%s\n%s" % (s, fielddisplay(self, 'tide_love_h', 'tidal load Love number (degree 2)'))
    64         s = "%s\n%s" % (s, fielddisplay(self, 'fluid_love', 'secular fluid Love number'))
    65         s = "%s\n%s" % (s, fielddisplay(self, 'equatorial_moi', 'mean equatorial moment of inertia [kg m^2]'))
    66         s = "%s\n%s" % (s, fielddisplay(self, 'polar_moi', 'polar moment of inertia [kg m^2]'))
    67         s = "%s\n%s" % (s, fielddisplay(self, 'angular_velocity', 'mean rotational velocity of earth [per second]'))
    68         s = "%s\n%s" % (s, fielddisplay(self, 'ocean_area_scaling', 'correction for model representation of ocean area [default: No correction]'))
    69         s = "%s\n%s" % (s, fielddisplay(self, 'hydro_rate', 'rate of hydrological expansion [mm / yr]'))
    70         s = "%s\n%s" % (s, fielddisplay(self, 'geodetic', 'compute geodetic SLR? (in addition to steric?) default 0'))
    71         s = "%s\n%s" % (s, fielddisplay(self, 'geodetic_run_frequency', 'how many time steps we skip before we run SLR solver during transient (default: 1)'))
    72         s = "%s\n%s" % (s, fielddisplay(self, 'rigid', 'rigid earth graviational potential perturbation'))
    73         s = "%s\n%s" % (s, fielddisplay(self, 'elastic', 'elastic earth graviational potential perturbation'))
    74         s = "%s\n%s" % (s, fielddisplay(self, 'rotation', 'earth rotational potential perturbation'))
    75         s = "%s\n%s" % (s, fielddisplay(self, 'degacc', 'accuracy (default .01 deg) for numerical discretization of the Green''s functions'))
    76         s = "%s\n%s" % (s, fielddisplay(self, 'transitions', 'indices into parts of the mesh that will be icecaps'))
    77         s = "%s\n%s" % (s, fielddisplay(self, 'requested_outputs', 'additional outputs requested'))
    78 
     52        s = '   slr parameters:\n'
     53        s += '{}\n'.format(fielddisplay(self, 'deltathickness', 'thickness change: ice height equivalent [m]'))
     54        s += '{}\n'.format(fielddisplay(self, 'sealevel', 'current sea level (prior to computation) [m]'))
     55        s += '{}\n'.format(fielddisplay(self, 'spcthickness', 'thickness constraints (NaN means no constraint) [m]'))
     56        s += '{}\n'.format(fielddisplay(self, 'reltol', 'sea level rise relative convergence criterion, (NaN: not applied)'))
     57        s += '{}\n'.format(fielddisplay(self, 'abstol', 'sea level rise absolute convergence criterion, (default, NaN: not applied)'))
     58        s += '{}\n'.format(fielddisplay(self, 'maxiter', 'maximum number of nonlinear iterations'))
     59        s += '{}\n'.format(fielddisplay(self, 'love_h', 'load Love number for radial displacement'))
     60        s += '{}\n'.format(fielddisplay(self, 'love_k', 'load Love number for gravitational potential perturbation'))
     61        s += '{}\n'.format(fielddisplay(self, 'love_l', 'load Love number for horizontal displaements'))
     62        s += '{}\n'.format(fielddisplay(self, 'tide_love_k', 'tidal load Love number (degree 2)'))
     63        s += '{}\n'.format(fielddisplay(self, 'tide_love_h', 'tidal load Love number (degree 2)'))
     64        s += '{}\n'.format(fielddisplay(self, 'fluid_love', 'secular fluid Love number'))
     65        s += '{}\n'.format(fielddisplay(self, 'equatorial_moi', 'mean equatorial moment of inertia [kg m^2]'))
     66        s += '{}\n'.format(fielddisplay(self, 'polar_moi', 'polar moment of inertia [kg m^2]'))
     67        s += '{}\n'.format(fielddisplay(self, 'angular_velocity', 'mean rotational velocity of earth [per second]'))
     68        s += '{}\n'.format(fielddisplay(self, 'ocean_area_scaling', 'correction for model representation of ocean area [default: No correction]'))
     69        s += '{}\n'.format(fielddisplay(self, 'hydro_rate', 'rate of hydrological expansion [mm / yr]'))
     70        s += '{}\n'.format(fielddisplay(self, 'geodetic', 'compute geodetic SLR? (in addition to steric?) default 0'))
     71        s += '{}\n'.format(fielddisplay(self, 'geodetic_run_frequency', 'how many time steps we skip before we run SLR solver during transient (default: 1)'))
     72        s += '{}\n'.format(fielddisplay(self, 'rigid', 'rigid earth graviational potential perturbation'))
     73        s += '{}\n'.format(fielddisplay(self, 'elastic', 'elastic earth graviational potential perturbation'))
     74        s += '{}\n'.format(fielddisplay(self, 'rotation', 'earth rotational potential perturbation'))
     75        s += '{}\n'.format(fielddisplay(self, 'degacc', 'accuracy (default .01 deg) for numerical discretization of the Green''s functions'))
     76        s += '{}\n'.format(fielddisplay(self, 'transitions', 'indices into parts of the mesh that will be icecaps'))
     77        s += '{}\n'.format(fielddisplay(self, 'requested_outputs', 'additional outputs requested'))
    7978        return s
    8079    # }}}
    8180
    8281    def setdefaultparameters(self):  # {{{
    83         #Convergence criterion: absolute, relative and residual
    84         self.reltol = 0.01 #default
    85         self.abstol = np.nan #1 mm of sea level rise
    86         #maximum of non - linear iterations.
     82        # Convergence criterion: absolute, relative and residual
     83        self.reltol = 0.01 # default
     84        self.abstol = np.nan # 1 mm of sea level rise
     85        # Maximum number of non-linear iterations
    8786        self.maxiter = 5
    88         #computational flags:
     87        # Computational flags
    8988        self.geodetic = 0
    9089        self.rigid = 1
    … …  
    9291        self.ocean_area_scaling = 0
    9392        self.rotation = 1
    94         #tidal love numbers:
    95         self.tide_love_h = 0.6149  #degree 2
    96         self.tide_love_k = 0.3055  #degree 2
    97         #secular fluid love number:
     93        # Tidal love numbers
     94        self.tide_love_h = 0.6149 # degree 2
     95        self.tide_love_k = 0.3055 # degree 2
     96        # Secular fluid love number
    9897        self.fluid_love = 0.942
    99         #moment of inertia:
     98        # Moment of inertia
    10099        self.equatorial_moi = 8.0077e37  # [kg m^2]
    101100        self.polar_moi = 8.0345e37  # [kg m^2]
    102         #mean rotational velocity of earth
    103         self.angular_velocity = 7.2921e-5  # [s^ - 1]
    104         #numerical discretization accuracy
     101        # Mean rotational velocity of earth
     102        self.angular_velocity = 7.2921e-5  # [s^-1]
     103        # Numerical discretization accuracy
    105104        self.degacc = 0.01
    106         #hydro
     105        # Hydro
    107106        self.hydro_rate = 0
    108         #how many time steps we skip before we run SLR solver during transient
     107        # How many time steps we skip before we run SLR solver during transient
    109108        self.geodetic_run_frequency = 1
    110         #output default:
     109        # Output default
    111110        self.requested_outputs = ['default']
    112         #transitions should be a cell array of vectors:
     111        # Transitions should be a list of vectors
    113112        self.transitions = []
    114         #horizontal displacement?  (not by default)
     113        # Horizontal displacement? (not on by default)
    115114        self.horiz = 0
    116         #earth area
     115        # Earth area
    117116        self.planetradius = planetradius('earth')
    118 
    119117        return self
    120118    #}}}
    121119
    122120    def checkconsistency(self, md, solution, analyses):  # {{{
    123         #Early return
    124         if (solution != 'SealevelriseAnalysis') or (solution == 'TransientSolution' and md.transient.isslr == 0):
     121        # Early return
     122        if 'SealevelriseAnalysis' not in analyses or (solution == 'TransientSolution' and not md.transient.isslr):
    125123            return md
    126124
    … …  
    146144        md = checkfield(md, 'fieldname', 'slr.horiz', 'NaN', 1, 'Inf', 1, 'values', [0, 1])
    147145
    148         #check that love numbers are provided at the same level of accuracy:
     146        # Check that love numbers are provided at the same level of accuracy:
    149147        if (size(self.love_h, 0) != size(self.love_k, 0)) or (size(self.love_h, 0) != size(self.love_l, 0)):
    150148            raise Exception('slr error message: love numbers should be provided at the same level of accuracy')
    151149
    152         #cross check that whereever we have an ice load, the mask is < 0 on each vertex:
     150        # Cross check that whereever we have an ice load, the mask is < 0 on each vertex:
    153151        pos = np.where(self.deltathickness)
    154152        maskpos = md.mask.ice_levelset[md.mesh.elements[pos, :]]
    … …  
    157155            warnings.warn('slr checkconsistency fail: there are elements with ice loads where some vertices are not on the ice!')
    158156
    159         #check that  if geodetic is requested, we are a mesh3dsurface model (planet), or if we are not,
    160         #a coupler to a planet model is provided.
    161         if self.geodetic and not md.transient.iscoupler and domaintype(md.mesh) != 'mesh3dsurface':
    162             raise Exception('model is requesting geodetic computations without being a mesh3dsurface, or being coupled to one!')
     157        # Check that if geodetic is requested, we are a mesh3dsurface model (planet), or if we are not, a coupler to a planet model is provided
     158        if self.geodetic:
     159            if md.transient.iscoupler:
     160                # We are good
     161                pass
     162            else:
     163                if md.mesh.__class__.__name__ == 'mesh3dsurface':
     164                    # We are good
     165                    pass
     166                else:
     167                    raise Exception('model is requesting geodetic computations without being a mesh3dsurface, or being coupled to one!')
    163168        return md
    164169    # }}}
    … …  
    196201        WriteData(fid, prefix, 'object', self, 'fieldname', 'planetradius', 'format', 'Double')
    197202
    198     #process requested outputs
     203        # Process requested outputs
    199204        outputs = self.requested_outputs
    200205        indices = [i for i, x in enumerate(outputs) if x == 'default']

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

    r25485 r25688  
    1212
    1313class solidearth(object):
    14     '''
    15     SOLIDEARTH class definition
     14    """SOLIDEARTH class definition
    1615
    17         Usage:
    18             solidearth = solidearth()
    19     '''
     16    Usage:
     17        solidearth = solidearth()
     18    """
    2019
    2120    def __init__(self, *args):  #{{{
    … …  
    5150
    5251    def setdefaultparameters(self):  # {{{
    53         return
     52        return self
    5453    #}}}
    5554
    5655    def checkconsistency(self, md, solution, analyses):  # {{{
    57         if ('SealevelriseAnalysis' not in analyses) or (solution == 'TransientSolution' and md.transient.isslr == 0):
     56        if ('SealevelriseAnalysis' not in analyses) or (solution == 'TransientSolution' and not md.transient.isslr):
    5857            return md
    59 
    6058        md = checkfield(md, 'fieldname', 'solidearth.sealevel', 'NaN', 1, 'Inf', 1, 'size', [md.mesh.numberofvertices])
    6159        md = checkfield(md, 'fieldname', 'solidearth.requested_outputs', 'stringrow', 1)
    62 
    6360        self.settings.checkconsistency(md, solution, analyses)
    6461        self.surfaceload.checkconsistency(md, solution, analyses)
    6562        self.lovenumbers.checkconsistency(md, solution, analyses)
    6663        self.rotational.checkconsistency(md, solution, analyses)
    67 
    6864        return md
    6965    #}}}
    … …  
    9288    #}}}
    9389
    94     def extrude(self, md):  #{{{
     90    def extrude(self, md): #{{{
    9591        self.sealevel = project3d(md, 'vector', self.sealevel, 'type', 'node')
    96 
    9792        return self
    9893    #}}}

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

    r25166 r25688  
    77
    88class solidearthsettings(object):
    9     '''
    10     SOLIDEARTHSETTINGS class definition
     9    """SOLIDEARTHSETTINGS class definition
    1110
    12         Usage:
    13             solidearthsettings = solidearthsettings()
    14     '''
     11    Usage:
     12        solidearthsettings = solidearthsettings()
     13    """
    1514
    1615    def __init__(self, *args): #{{{
    … …  
    3433        s += '{}\n'.format(fielddisplay(self, 'elastic', 'elastic earth graviational potential perturbation'))
    3534        s += '{}\n'.format(fielddisplay(self, 'degacc', 'accuracy (default .01 deg) for numerical discretization of the Green\'s functions'))
    36 
    3735        return s
    3836    #}}}
    3937
    4038    def setdefaultparameters(self): # {{{
    41         #Convergence criterion: absolute, relative, and residual
     39        # Convergence criterion: absolute, relative, and residual
    4240        self.reltol = 0.01 # 1 percent
    4341        self.abstol = np.nan # default
    4442
    45         #maximum of non-linear iterations
     43        # Maximum of non-linear iterations
    4644        self.maxiter = 5
    4745
    48         #computational flags
     46        # Computational flags
    4947        self.rigid = 1
    5048        self.elastic = 1
    … …  
    5351        self.computesealevelchange = 0
    5452
    55         #numerical discetization accuracy
     53        # Numerical discretization accuracy
    5654        self.degacc = .01
    5755
    58         #how many time steps we skip before we run solidearthsettings solver during transient
     56        # How many time steps we skip before we run solidearthsettings solver during transient
    5957        self.runfrequency = 1
    6058
    61         #horizontal displacemnet? (not by default)
     59        # Horizontal displacement? (not on by default)
    6260        self.horiz = 0
     61        return self
    6362    #}}}
    6463
    6564    def checkconsistency(self, md, solution, analyses): # {{{
    66         if ('SealevelriseAnalysis' not in analyses) or (solution == 'TransientSolution' and md.transient.isslr == 0):
     65        if ('SealevelriseAnalysis' not in analyses) or (solution == 'TransientSolution' and not md.transient.isslr):
    6766            return md
    68 
    6967        md = checkfield(md, 'fieldname', 'solidearth.settings.reltol', 'size', [1])
    7068        md = checkfield(md, 'fieldname', 'solidearth.settings.abstol', 'size', [1])
    … …  
    7472        md = checkfield(md, 'fieldname', 'solidearth.settings.horiz', 'NaN', 1, 'Inf', 1, 'values', [0, 1])
    7573
    76         #a coupler to planet model is provided
     74        # A coupler to planet model is provided
    7775        if self.computesealevelchange:
    7876            if md.transient.iscoupler:
    79                 #we are good
     77                # We are good
    8078                pass
    8179            else:
    8280                if md.mesh.__class__.__name__ == 'mesh3dsurface':
    83                     #we are good
     81                    # We are good
    8482                    pass
    8583                else:
    8684                    raise Exception('model is requesting sealevel computations without being a mesh3dsurface, or being coupled to one!')
    87 
    8885        return md
    8986    #}}}

  • issm/trunk-jpl/src/m/classes/spatiallinearbasalforcings.m

    r22945 r25688  
    8181                end % }}}
    8282                function disp(self) % {{{
    83                         disp(sprintf('   basal forcings parameters:'));
     83                        disp(sprintf('   spatial linear basal forcings parameters:'));
    8484
    8585                        fielddisplay(self,'groundedice_melting_rate','basal melting rate (positive if melting) [m/yr]');
    … …  
    100100                        floatingice_melting_rate(pos)=md.basalforcings.deepwater_melting_rate(pos).*(md.geometry.base(pos)-md.basalforcings.upperwater_elevation(pos))./(md.basalforcings.deepwater_elevation(pos)-md.basalforcings.upperwater_elevation(pos));
    101101                        WriteData(fid,prefix,'name','md.basalforcings.model','data',6,'format','Integer');
    102                         WriteData(fid,prefix,'data',floatingice_melting_rate,'format','DoubleMat','name','md.basalforcings.floatingice_melting_rate','mattype',1,'scale',1./yts,'timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts)
    103                         WriteData(fid,prefix,'object',self,'fieldname','groundedice_melting_rate','format','DoubleMat','name','md.basalforcings.groundedice_melting_rate','mattype',1,'scale',1./yts,'timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts)
     102                        WriteData(fid,prefix,'data',floatingice_melting_rate,'format','DoubleMat','name','md.basalforcings.floatingice_melting_rate','mattype',1,'scale',1./yts,'timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts);
     103                        WriteData(fid,prefix,'object',self,'fieldname','groundedice_melting_rate','format','DoubleMat','name','md.basalforcings.groundedice_melting_rate','mattype',1,'scale',1./yts,'timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts);
    104104                        WriteData(fid,prefix,'object',self,'fieldname','geothermalflux','name','md.basalforcings.geothermalflux','format','DoubleMat','mattype',1,'timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts);
    105                         WriteData(fid,prefix,'object',self,'fieldname','deepwater_melting_rate','format','DoubleMat','name','md.basalforcings.deepwater_melting_rate','scale',1./yts,'mattype',1)
    106                         WriteData(fid,prefix,'object',self,'fieldname','deepwater_elevation','format','DoubleMat','name','md.basalforcings.deepwater_elevation','mattype',1)
    107                         WriteData(fid,prefix,'object',self,'fieldname','upperwater_elevation','format','DoubleMat','name','md.basalforcings.upperwater_elevation','mattype',1)
     105                        WriteData(fid,prefix,'object',self,'fieldname','deepwater_melting_rate','format','DoubleMat','name','md.basalforcings.deepwater_melting_rate','scale',1./yts,'mattype',1);
     106                        WriteData(fid,prefix,'object',self,'fieldname','deepwater_elevation','format','DoubleMat','name','md.basalforcings.deepwater_elevation','mattype',1);
     107                        WriteData(fid,prefix,'object',self,'fieldname','upperwater_elevation','format','DoubleMat','name','md.basalforcings.upperwater_elevation','mattype',1);
    108108                end % }}}
    109109        end

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

    r25023 r25688  
     1import sys
     2
    13import numpy as np
    2 import sys
     4
     5from checkfield import checkfield
     6from fielddisplay import fielddisplay
     7import MatlabFuncs as m
    38from project3d import project3d
    4 from fielddisplay import fielddisplay
    5 from checkfield import checkfield
    69from WriteData import WriteData
    7 import MatlabFuncs as m
    810
    911
    1012class stressbalance(object):
    11     """
    12     STRESSBALANCE class definition
     13    """STRESSBALANCE class definition
    1314
    14        Usage:
    15           stressbalance = stressbalance()
     15    Usage:
     16        stressbalance = stressbalance()
    1617    """
    1718
    … …  
    4142    #}}}
    4243    def __repr__(self):  # {{{
    43 
    44         string = '   StressBalance solution parameters:'
    45         string = "%s\n%s" % (string, '      Convergence criteria:')
    46         string = "%s\n%s" % (string, fielddisplay(self, 'restol', 'mechanical equilibrium residual convergence criterion'))
    47         string = "%s\n%s" % (string, fielddisplay(self, 'reltol', 'velocity relative convergence criterion, NaN: not applied'))
    48         string = "%s\n%s" % (string, fielddisplay(self, 'abstol', 'velocity absolute convergence criterion, NaN: not applied'))
    49         string = "%s\n%s" % (string, fielddisplay(self, 'isnewton', "0: Picard's fixed point, 1: Newton's method, 2: hybrid"))
    50         string = "%s\n%s" % (string, fielddisplay(self, 'maxiter', 'maximum number of nonlinear iterations'))
    51 
    52         string = "%s\n%s" % (string, '\n      boundary conditions:')
    53         string = "%s\n%s" % (string, fielddisplay(self, 'spcvx', 'x - axis velocity constraint (NaN means no constraint) [m / yr]'))
    54         string = "%s\n%s" % (string, fielddisplay(self, 'spcvy', 'y - axis velocity constraint (NaN means no constraint) [m / yr]'))
    55         string = "%s\n%s" % (string, fielddisplay(self, 'spcvz', 'z - axis velocity constraint (NaN means no constraint) [m / yr]'))
    56         string = "%s\n%s" % (string, fielddisplay(self, 'icefront', 'segments on ice front list (last column 0: Air, 1: Water, 2: Ice'))
    57 
    58         string = "%s\n%s" % (string, '\n      Rift options:')
    59         string = "%s\n%s" % (string, fielddisplay(self, 'rift_penalty_threshold', 'threshold for instability of mechanical constraints'))
    60         string = "%s\n%s" % (string, fielddisplay(self, 'rift_penalty_lock', 'number of iterations before rift penalties are locked'))
    61 
    62         string = "%s\n%s" % (string, '\n      Penalty options:')
    63         string = "%s\n%s" % (string, fielddisplay(self, 'penalty_factor', 'offset used by penalties: penalty = Kmax * 1.0**offset'))
    64         string = "%s\n%s" % (string, fielddisplay(self, 'vertex_pairing', 'pairs of vertices that are penalized'))
    65 
    66         string = "%s\n%s" % (string, '\n      Other:')
    67         string = "%s\n%s" % (string, fielddisplay(self, 'shelf_dampening', 'use dampening for floating ice ? Only for FS model'))
    68         string = "%s\n%s" % (string, fielddisplay(self, 'FSreconditioning', 'multiplier for incompressibility equation. Only for FS model'))
    69         string = "%s\n%s" % (string, fielddisplay(self, 'referential', 'local referential'))
    70         string = "%s\n%s" % (string, fielddisplay(self, 'loadingforce', 'loading force applied on each point [N / m^3]'))
    71         string = "%s\n%s" % (string, fielddisplay(self, 'requested_outputs', 'additional outputs requested'))
    72 
    73         return string
     44        s = '   StressBalance solution parameters:\n'
     45        s += '      Convergence criteria:\n'
     46        s += '{}\n'.format(fielddisplay(self, 'restol', 'mechanical equilibrium residual convergence criterion'))
     47        s += '{}\n'.format(fielddisplay(self, 'reltol', 'velocity relative convergence criterion, NaN: not applied'))
     48        s += '{}\n'.format(fielddisplay(self, 'abstol', 'velocity absolute convergence criterion, NaN: not applied'))
     49        s += '{}\n'.format(fielddisplay(self, 'isnewton', "0: Picard's fixed point, 1: Newton's method, 2: hybrid"))
     50        s += '{}\n'.format(fielddisplay(self, 'maxiter', 'maximum number of nonlinear iterations'))
     51        s += '      boundary conditions:\n'
     52        s += '{}\n'.format(fielddisplay(self, 'spcvx', 'x - axis velocity constraint (NaN means no constraint) [m / yr]'))
     53        s += '{}\n'.format(fielddisplay(self, 'spcvy', 'y - axis velocity constraint (NaN means no constraint) [m / yr]'))
     54        s += '{}\n'.format(fielddisplay(self, 'spcvz', 'z - axis velocity constraint (NaN means no constraint) [m / yr]'))
     55        s += '{}\n'.format(fielddisplay(self, 'icefront', 'segments on ice front list (last column 0: Air, 1: Water, 2: Ice'))
     56        s = "%s\n%s" % (s, '\n      Rift options:')
     57        s += '{}\n'.format(fielddisplay(self, 'rift_penalty_threshold', 'threshold for instability of mechanical constraints'))
     58        s += '{}\n'.format(fielddisplay(self, 'rift_penalty_lock', 'number of iterations before rift penalties are locked'))
     59        s += '      Penalty options:\n'
     60        s += '{}\n'.format(fielddisplay(self, 'penalty_factor', 'offset used by penalties: penalty = Kmax * 1.0**offset'))
     61        s += '{}\n'.format(fielddisplay(self, 'vertex_pairing', 'pairs of vertices that are penalized'))
     62        s += '      Other:\n'
     63        s += '{}\n'.format(fielddisplay(self, 'shelf_dampening', 'use dampening for floating ice ? Only for FS model'))
     64        s += '{}\n'.format(fielddisplay(self, 'FSreconditioning', 'multiplier for incompressibility equation. Only for FS model'))
     65        s += '{}\n'.format(fielddisplay(self, 'referential', 'local referential'))
     66        s += '{}\n'.format(fielddisplay(self, 'loadingforce', 'loading force applied on each point [N / m^3]'))
     67        s += '{}\n'.format(fielddisplay(self, 'requested_outputs', 'additional outputs requested'))
     68        return s
    7469    #}}}
    7570
    … …  
    10398        #output default:
    10499        self.requested_outputs = ['default']
    105 
    106100        return self
    107101    #}}}
    … …  
    113107            list = ['Vx', 'Vy', 'Vel', 'Pressure']
    114108        return list
    115 
    116109    #}}}
    117110
    118111    def checkconsistency(self, md, solution, analyses):  # {{{
    119         #Early return
     112        # Early return
    120113        if 'StressbalanceAnalysis' not in analyses:
     114            return md
     115        if solution == 'TransientSolution' and not md.transient.isstressbalance:
    121116            return md
    122117
    … …  
    159154            if np.any(np.logical_not(np.isnan(md.stressbalance.referential[pos, :]))):
    160155                md.checkmessage("no referential should be specified for basal vertices of grounded ice")
    161 
    162156        return md
    163157    # }}}
    … …  
    183177        WriteData(fid, prefix, 'object', self, 'class', 'stressbalance', 'fieldname', 'rift_penalty_threshold', 'format', 'Integer')
    184178        WriteData(fid, prefix, 'object', self, 'class', 'stressbalance', 'fieldname', 'referential', 'format', 'DoubleMat', 'mattype', 1)
    185 
    186179        if isinstance(self.loadingforce, (list, tuple, np.ndarray)) and np.size(self.loadingforce, 1) == 3:
    187180            WriteData(fid, prefix, 'data', self.loadingforce[:, 0], 'format', 'DoubleMat', 'mattype', 1, 'name', 'md.stressbalance.loadingforcex')
    188181            WriteData(fid, prefix, 'data', self.loadingforce[:, 1], 'format', 'DoubleMat', 'mattype', 1, 'name', 'md.stressbalance.loadingforcey')
    189182            WriteData(fid, prefix, 'data', self.loadingforce[:, 2], 'format', 'DoubleMat', 'mattype', 1, 'name', 'md.stressbalance.loadingforcez')
    190 
    191     #process requested outputs
     183        # Process requested outputs
    192184        outputs = self.requested_outputs
    193185        indices = [i for i, x in enumerate(outputs) if x == 'default']

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

    r25499 r25688  
    3131        s += '{}\n'.format(fielddisplay(self, 'waterheightchange', 'water height change: water height equivalent [mWater/yr]'))
    3232        s += '{}\n'.format(fielddisplay(self, 'other', 'other loads (sediments) [kg/m^2/yr]'))
    33 
    3433        return s
    3534    #}}}
    3635
    3736    def setdefaultparameters(self): # {{{
    38         return
     37        return self
    3938    #}}}
    4039
    4140    def checkconsistency(self, md, solution, analyses): # {{{
    42         if ('SealevelriseAnalysis' not in analyses) or (solution == 'TransientSolution' and md.transient.isslr == 0):
     41        if ('SealevelriseAnalysis' not in analyses) or (solution == 'TransientSolution' and not md.transient.isslr):
    4342            return md
    44 
    4543        if len(self.icethicknesschange):
    4644            md  = checkfield(md,'fieldname', 'solidearth.surfaceload.icethicknesschange', 'timeseries', 1, 'NaN', 1, 'Inf', 1)
    47 
    4845        if len(self.waterheightchange):
    4946            md  = checkfield(md,'fieldname', 'solidearth.surfaceload.waterheightchange', 'timeseries', 1, 'NaN', 1, 'Inf', 1)
    50 
    5147        if len(self.other):
    5248            md  = checkfield(md,'fieldname', 'solidearth.surfaceload.other', 'timeseries', 1, 'NaN', 1, 'Inf', 1)
    53 
    5449        return md
    5550    #}}}
    … …  
    5853        if len(self.icethicknesschange) == 0:
    5954            self.icethicknesschange = np.zeros((md.mesh.numberofelements + 1, ))
    60 
    6155        if len(self.waterheightchange) == 0:
    6256            self.waterheightchange = np.zeros((md.mesh.numberofelements + 1, ))
    63 
    6457        if len(self.other) == 0:
    6558            self.other = np.zeros((md.mesh.numberofelements + 1, ))
    66 
    6759        WriteData(fid, prefix, 'object', self, 'fieldname', 'icethicknesschange', 'name', 'md.solidearth.surfaceload.icethicknesschange', 'format', 'DoubleMat', 'mattype', 2, 'timeserieslength', md.mesh.numberofelements + 1, 'yts', md.constants.yts)
    6860        WriteData(fid, prefix, 'object', self, 'fieldname', 'waterheightchange', 'name', 'md.solidearth.surfaceload.waterheightchange', 'format', 'DoubleMat', 'mattype', 2, 'timeserieslength', md.mesh.numberofelements + 1, 'yts', md.constants.yts)

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

    r24213 r25688  
    11import numpy as np
    2 from project3d import project3d
    3 from WriteData import WriteData
     2
    43from checkfield import checkfield
    54from IssmConfig import IssmConfig
    65from marshallcostfunctions import marshallcostfunctions
     6from project3d import project3d
    77from supportedcontrols import *
    88from supportedcostfunctions import *
     9from WriteData import WriteData
    910
    1011
    11 class taoinversion(object):
     12class taoinversion(object): #{{{
     13    """TAOINVERSION class definition
     14
     15    Usage:
     16        inversion = taoinversion()
     17    """
     18
    1219    def __init__(self):
    1320        self.iscontrol = 0
    1421        self.incomplete_adjoint = 0
    15         self.control_parameters = float('NaN')
     22        self.control_parameters = np.nan
    1623        self.maxsteps = 0
    1724        self.maxiter = 0
    … …  
    2229        self.gttol = 0
    2330        self.algorithm = ''
    24         self.cost_functions = float('NaN')
    25         self.cost_functions_coefficients = float('NaN')
    26         self.min_parameters = float('NaN')
    27         self.max_parameters = float('NaN')
    28         self.vx_obs = float('NaN')
    29         self.vy_obs = float('NaN')
    30         self.vz_obs = float('NaN')
    31         self.vel_obs = float('NaN')
    32         self.thickness_obs = float('NaN')
    33         self.surface_obs = float('NaN')
     31        self.cost_functions = np.nan
     32        self.cost_functions_coefficients = np.nan
     33        self.min_parameters = np.nan
     34        self.max_parameters = np.nan
     35        self.vx_obs = np.nan
     36        self.vy_obs = np.nan
     37        self.vz_obs = np.nan
     38        self.vel_obs = np.nan
     39        self.thickness_obs = np.nan
     40        self.surface_obs = np.nan
     41
    3442        self.setdefaultparameters()
     43    #}}}
    3544
    3645    def __repr__(self):
    37         string = '   taoinversion parameters:'
    38         string = "%s\n%s" % (string, fieldstring(self, 'iscontrol', 'is inversion activated?'))
    39         string = "%s\n%s" % (string, fieldstring(self, 'mantle_viscosity', 'mantle viscosity constraints (NaN means no constraint) (Pa s)'))
    40         string = "%s\n%s" % (string, fieldstring(self, 'lithosphere_thickness', 'lithosphere thickness constraints (NaN means no constraint) (m)'))
    41         string = "%s\n%s" % (string, fieldstring(self, 'cross_section_shape', "1: square-edged, 2: elliptical - edged surface"))
    42         string = "%s\n%s" % (string, fieldstring(self, 'incomplete_adjoint', '1: linear viscosity, 0: non - linear viscosity'))
    43         string = "%s\n%s" % (string, fieldstring(self, 'control_parameters', 'ex: {''FrictionCoefficient''}, or {''MaterialsRheologyBbar''}'))
    44         string = "%s\n%s" % (string, fieldstring(self, 'maxsteps', 'maximum number of iterations (gradient computation)'))
    45         string = "%s\n%s" % (string, fieldstring(self, 'maxiter', 'maximum number of Function evaluation (forward run)'))
    46         string = "%s\n%s" % (string, fieldstring(self, 'fatol', 'convergence criterion: f(X) - f(X * ) (X: current iteration, X * : "true" solution, f: cost function)'))
    47         string = "%s\n%s" % (string, fieldstring(self, 'frtol', 'convergence criterion: |f(X) - f(X * )| / |f(X * )|'))
    48         string = "%s\n%s" % (string, fieldstring(self, 'gatol', 'convergence criterion: ||g(X)|| (g: gradient of the cost function)'))
    49         string = "%s\n%s" % (string, fieldstring(self, 'grtol', 'convergence criterion: ||g(X)|| / |f(X)|'))
    50         string = "%s\n%s" % (string, fieldstring(self, 'gttol', 'convergence criterion: ||g(X)|| / ||g(X0)|| (g(X0): gradient at initial guess X0)'))
    51         string = "%s\n%s" % (string, fieldstring(self, 'algorithm', 'minimization algorithm: ''tao_blmvm'', ''tao_cg'', ''tao_lmvm'''))
    52         string = "%s\n%s" % (string, fieldstring(self, 'cost_functions', 'indicate the type of response for each optimization step'))
    53         string = "%s\n%s" % (string, fieldstring(self, 'cost_functions_coefficients', 'cost_functions_coefficients applied to the misfit of each vertex and for each control_parameter'))
    54         string = "%s\n%s" % (string, fieldstring(self, 'min_parameters', 'absolute minimum acceptable value of the inversed parameter on each vertex'))
    55         string = "%s\n%s" % (string, fieldstring(self, 'max_parameters', 'absolute maximum acceptable value of the inversed parameter on each vertex'))
    56         string = "%s\n%s" % (string, fieldstring(self, 'vx_obs', 'observed velocity x component [m / yr]'))
    57         string = "%s\n%s" % (string, fieldstring(self, 'vy_obs', 'observed velocity y component [m / yr]'))
    58         string = "%s\n%s" % (string, fieldstring(self, 'vel_obs', 'observed velocity magnitude [m / yr]'))
    59         string = "%s\n%s" % (string, fieldstring(self, 'thickness_obs', 'observed thickness [m]'))
    60         string = "%s\n%s" % (string, fieldstring(self, 'surface_obs', 'observed surface elevation [m]'))
    61         string = "%s\n%s" % (string, 'Available cost functions:')
    62         string = "%s\n%s" % (string, '   101: SurfaceAbsVelMisfit')
    63         string = "%s\n%s" % (string, '   102: SurfaceRelVelMisfit')
    64         string = "%s\n%s" % (string, '   103: SurfaceLogVelMisfit')
    65         string = "%s\n%s" % (string, '   104: SurfaceLogVxVyMisfit')
    66         string = "%s\n%s" % (string, '   105: SurfaceAverageVelMisfit')
    67         string = "%s\n%s" % (string, '   201: ThicknessAbsMisfit')
    68         string = "%s\n%s" % (string, '   501: DragCoefficientAbsGradient')
    69         string = "%s\n%s" % (string, '   502: RheologyBbarAbsGradient')
    70         string = "%s\n%s" % (string, '   503: ThicknessAbsGradient')
    71         return string
     46        s = '   taoinversion parameters:\n'
     47        s += '{}'.format(fieldstring(self, 'iscontrol', 'is inversion activated?'))
     48        s += '{}'.format(fieldstring(self, 'mantle_viscosity', 'mantle viscosity constraints (NaN means no constraint) (Pa s)'))
     49        s += '{}'.format(fieldstring(self, 'lithosphere_thickness', 'lithosphere thickness constraints (NaN means no constraint) (m)'))
     50        s += '{}'.format(fieldstring(self, 'cross_section_shape', "1: square-edged, 2: elliptical - edged surface"))
     51        s += '{}'.format(fieldstring(self, 'incomplete_adjoint', '1: linear viscosity, 0: non - linear viscosity'))
     52        s += '{}'.format(fieldstring(self, 'control_parameters', 'ex: {''FrictionCoefficient''}, or {''MaterialsRheologyBbar''}'))
     53        s += '{}'.format(fieldstring(self, 'maxsteps', 'maximum number of iterations (gradient computation)'))
     54        s += '{}'.format(fieldstring(self, 'maxiter', 'maximum number of Function evaluation (forward run)'))
     55        s += '{}'.format(fieldstring(self, 'fatol', 'convergence criterion: f(X) - f(X * ) (X: current iteration, X * : "true" solution, f: cost function)'))
     56        s += '{}'.format(fieldstring(self, 'frtol', 'convergence criterion: |f(X) - f(X * )| / |f(X * )|'))
     57        s += '{}'.format(fieldstring(self, 'gatol', 'convergence criterion: ||g(X)|| (g: gradient of the cost function)'))
     58        s += '{}'.format(fieldstring(self, 'grtol', 'convergence criterion: ||g(X)|| / |f(X)|'))
     59        s += '{}'.format(fieldstring(self, 'gttol', 'convergence criterion: ||g(X)|| / ||g(X0)|| (g(X0): gradient at initial guess X0)'))
     60        s += '{}'.format(fieldstring(self, 'algorithm', 'minimization algorithm: ''tao_blmvm'', ''tao_cg'', ''tao_lmvm'''))
     61        s += '{}'.format(fieldstring(self, 'cost_functions', 'indicate the type of response for each optimization step'))
     62        s += '{}'.format(fieldstring(self, 'cost_functions_coefficients', 'cost_functions_coefficients applied to the misfit of each vertex and for each control_parameter'))
     63        s += '{}'.format(fieldstring(self, 'min_parameters', 'absolute minimum acceptable value of the inversed parameter on each vertex'))
     64        s += '{}'.format(fieldstring(self, 'max_parameters', 'absolute maximum acceptable value of the inversed parameter on each vertex'))
     65        s += '{}'.format(fieldstring(self, 'vx_obs', 'observed velocity x component [m / yr]'))
     66        s += '{}'.format(fieldstring(self, 'vy_obs', 'observed velocity y component [m / yr]'))
     67        s += '{}'.format(fieldstring(self, 'vel_obs', 'observed velocity magnitude [m / yr]'))
     68        s += '{}'.format(fieldstring(self, 'thickness_obs', 'observed thickness [m]'))
     69        s += '{}'.format(fieldstring(self, 'surface_obs', 'observed surface elevation [m]'))
     70        s += '{}'.format('Available cost functions:')
     71        s += '{}'.format('   101: SurfaceAbsVelMisfit')
     72        s += '{}'.format('   102: SurfaceRelVelMisfit')
     73        s += '{}'.format('   103: SurfaceLogVelMisfit')
     74        s += '{}'.format('   104: SurfaceLogVxVyMisfit')
     75        s += '{}'.format('   105: SurfaceAverageVelMisfit')
     76        s += '{}'.format('   201: ThicknessAbsMisfit')
     77        s += '{}'.format('   501: DragCoefficientAbsGradient')
     78        s += '{}'.format('   502: RheologyBbarAbsGradient')
     79        s += '{}'.format('   503: ThicknessAbsGradient')
     80        return s
    7281
    7382    def setdefaultparameters(self):

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

    r24567 r25688  
    11import numpy as np
     2
     3from checkfield import checkfield
     4from fielddisplay import fielddisplay
    25from project3d import project3d
    3 from fielddisplay import fielddisplay
    4 from checkfield import checkfield
    56from WriteData import WriteData
    67
    78
    89class thermal(object):
    9     """
    10     THERMAL class definition
     10    """THERMAL class definition
    1111
    1212       Usage:
    … …  
    2828        self.fe = 'P1'
    2929        self.requested_outputs = []
    30 
    31     #set defaults
    3230        self.setdefaultparameters()
    33 
    3431    #}}}
    3532
    3633    def __repr__(self):  # {{{
    37         string = '   Thermal solution parameters:'
    38         string = "%s\n%s" % (string, fielddisplay(self, 'spctemperature', 'temperature constraints (NaN means no constraint) [K]'))
    39         string = "%s\n%s" % (string, fielddisplay(self, 'stabilization', '0: no, 1: artificial_diffusivity, 2: SUPG'))
    40         string = "%s\n%s" % (string, fielddisplay(self, 'maxiter', 'maximum number of non linear iterations'))
    41         string = "%s\n%s" % (string, fielddisplay(self, 'reltol', 'relative tolerance criterion'))
    42         string = "%s\n%s" % (string, fielddisplay(self, 'penalty_lock', 'stabilize unstable thermal constraints that keep zigzagging after n iteration (default is 0, no stabilization)'))
    43         string = "%s\n%s" % (string, fielddisplay(self, 'penalty_threshold', 'threshold to declare convergence of thermal solution (default is 0)'))
    44         string = "%s\n%s" % (string, fielddisplay(self, 'isenthalpy', 'use an enthalpy formulation to include temperate ice (default is 0)'))
    45         string = "%s\n%s" % (string, fielddisplay(self, 'isdynamicbasalspc', 'enable dynamic setting of basal forcing. required for enthalpy formulation (default is 0)'))
    46         string = "%s\n%s" % (string, fielddisplay(self, 'isdrainicecolumn', 'wether waterfraction drainage is enabled for enthalpy formulation (default is 1)'))
    47         string = "%s\n%s" % (string, fielddisplay(self, 'watercolumn_upperlimit', 'upper limit of basal watercolumn for enthalpy formulation (default is 1000m)'))
    48         string = "%s\n%s" % (string, fielddisplay(self, 'requested_outputs', 'additional outputs requested'))
    49         return string
     34        s = '   Thermal solution parameters:\n'
     35        s += '{}\n'.format(fielddisplay(self, 'spctemperature', 'temperature constraints (NaN means no constraint) [K]'))
     36        s += '{}\n'.format(fielddisplay(self, 'stabilization', '0: no, 1: artificial_diffusivity, 2: SUPG'))
     37        s += '{}\n'.format(fielddisplay(self, 'maxiter', 'maximum number of non linear iterations'))
     38        s += '{}\n'.format(fielddisplay(self, 'reltol', 'relative tolerance criterion'))
     39        s += '{}\n'.format(fielddisplay(self, 'penalty_lock', 'stabilize unstable thermal constraints that keep zigzagging after n iteration (default is 0, no stabilization)'))
     40        s += '{}\n'.format(fielddisplay(self, 'penalty_threshold', 'threshold to declare convergence of thermal solution (default is 0)'))
     41        s += '{}\n'.format(fielddisplay(self, 'isenthalpy', 'use an enthalpy formulation to include temperate ice (default is 0)'))
     42        s += '{}\n'.format(fielddisplay(self, 'isdynamicbasalspc', 'enable dynamic setting of basal forcing. required for enthalpy formulation (default is 0)'))
     43        s += '{}\n'.format(fielddisplay(self, 'isdrainicecolumn', 'wether waterfraction drainage is enabled for enthalpy formulation (default is 1)'))
     44        s += '{}\n'.format(fielddisplay(self, 'watercolumn_upperlimit', 'upper limit of basal watercolumn for enthalpy formulation (default is 1000m)'))
     45        s += '{}\n'.format(fielddisplay(self, 'requested_outputs', 'additional outputs requested'))
     46        return s
    5047    #}}}
    5148
    … …  
    6461        else:
    6562            return ['Temperature', 'BasalforcingsGroundediceMeltingRate']
    66 
    6763    #}}}
    6864
    … …  
    9187        self.requested_outputs = ['default']
    9288        return self
    93 
    9489    #}}}
    9590
    9691    def checkconsistency(self, md, solution, analyses):  # {{{
    97         #Early return
     92        # Early return
    9893        if ('ThermalAnalysis' not in analyses and 'EnthalpyAnalysis' not in analyses) or (solution == 'TransientSolution' and not md.transient.isthermal):
    9994            return md
    100 
    10195        md = checkfield(md, 'fieldname', 'thermal.stabilization', 'numel', [1], 'values', [0, 1, 2])
    10296        md = checkfield(md, 'fieldname', 'thermal.spctemperature', 'Inf', 1, 'timeseries', 1)
    10397        md = checkfield(md, 'fieldname', 'thermal.requested_outputs', 'stringrow', 1)
    104 
    10598        if 'EnthalpyAnalysis' in analyses and md.thermal.isenthalpy and md.mesh.dimension() == 3:
    10699            md = checkfield(md, 'fieldname', 'thermal.isdrainicecolumn', 'numel', [1], 'values', [0, 1])
    … …  
    123116                    md.checkmessage("for a steadystate computation, thermal.reltol (relative convergence criterion) must be defined!")
    124117                md = checkfield(md, 'fieldname', 'thermal.reltol', '>', 0., 'message', "reltol must be larger than zero")
    125 
    126118        return md
    127119    # }}}

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

    r24990 r25688  
     1from checkfield import checkfield
    12from fielddisplay import fielddisplay
    2 from checkfield import checkfield
    33from WriteData import WriteData
    44
    55
    66class transient(object):
    7     """
    8     TRANSIENT class definition
     7    """TRANSIENT class definition
    98
    10        Usage:
    11           transient = transient()
     9    Usage:
     10        transient = transient()
    1211    """
    1312
    1413    def __init__(self):  # {{{
    15         self.issmb = False
    16         self.ismasstransport = False
    17         self.isstressbalance = False
    18         self.isthermal = False
    19         self.isgroundingline = False
    20         self.isgia = False
    21         self.isesa = False
    22         self.isdamageevolution = False
    23         self.ismovingfront = False
    24         self.ishydrology = False
    25         self.isslr = False
    26         self.iscoupler = False
     14        self.issmb = 0
     15        self.ismasstransport = 0
     16        self.isstressbalance = 0
     17        self.isthermal = 0
     18        self.isgroundingline = 0
     19        self.isgia = 0
     20        self.isesa = 0
     21        self.isdamageevolution = 0
     22        self.ismovingfront = 0
     23        self.ishydrology = 0
     24        self.isslr = 0
     25        self.iscoupler = 0
    2726        self.amr_frequency = 0
    28         self.isoceancoupling = False
     27        self.isoceancoupling = 0
    2928        self.requested_outputs = []
    3029
    31     #set defaults
    3230        self.setdefaultparameters()
     31    #}}}
    3332
    34     #}}}
    3533    def __repr__(self):  # {{{
    36         string = '   transient solution parameters:'
    37         string = "%s\n%s" % (string, fielddisplay(self, 'issmb', 'indicates if a surface mass balance solution is used in the transient'))
    38         string = "%s\n%s" % (string, fielddisplay(self, 'ismasstransport', 'indicates if a masstransport solution is used in the transient'))
    39         string = "%s\n%s" % (string, fielddisplay(self, 'isstressbalance', 'indicates if a stressbalance solution is used in the transient'))
    40         string = "%s\n%s" % (string, fielddisplay(self, 'isthermal', 'indicates if a thermal solution is used in the transient'))
    41         string = "%s\n%s" % (string, fielddisplay(self, 'isgroundingline', 'indicates if a groundingline migration is used in the transient'))
    42         string = "%s\n%s" % (string, fielddisplay(self, 'isgia', 'indicates if a postglacial rebound is used in the transient'))
    43         string = "%s\n%s" % (string, fielddisplay(self, 'isesa', 'indicates whether an elastic adjustment model is used in the transient'))
    44         string = "%s\n%s" % (string, fielddisplay(self, 'isdamageevolution', 'indicates whether damage evolution is used in the transient'))
    45         string = "%s\n%s" % (string, fielddisplay(self, 'ismovingfront', 'indicates whether a moving front capability is used in the transient'))
    46         string = "%s\n%s" % (string, fielddisplay(self, 'ishydrology', 'indicates whether an hydrology model is used'))
    47         string = "%s\n%s" % (string, fielddisplay(self, 'isslr', 'indicates if a sea level rise solution is used in the transient'))
    48         string = "%s\n%s" % (string, fielddisplay(self, 'isoceancoupling', 'indicates whether coupling with an ocean model is used in the transient'))
    49         string = "%s\n%s" % (string, fielddisplay(self, 'iscoupler', 'indicates whether different models are being run with need for coupling'))
    50         string = "%s\n%s" % (string, fielddisplay(self, 'amr_frequency', 'frequency at which mesh is refined in simulations with multiple time_steps'))
    51         string = "%s\n%s" % (string, fielddisplay(self, 'requested_outputs', 'list of additional outputs requested'))
    52         return string
     34        s = '   transient solution parameters:\n'
     35        s += '{}\n'.format(fielddisplay(self, 'issmb', 'indicates if a surface mass balance solution is used in the transient'))
     36        s += '{}\n'.format(fielddisplay(self, 'ismasstransport', 'indicates if a masstransport solution is used in the transient'))
     37        s += '{}\n'.format(fielddisplay(self, 'isstressbalance', 'indicates if a stressbalance solution is used in the transient'))
     38        s += '{}\n'.format(fielddisplay(self, 'isthermal', 'indicates if a thermal solution is used in the transient'))
     39        s += '{}\n'.format(fielddisplay(self, 'isgroundingline', 'indicates if a groundingline migration is used in the transient'))
     40        s += '{}\n'.format(fielddisplay(self, 'isgia', 'indicates if a postglacial rebound is used in the transient'))
     41        s += '{}\n'.format(fielddisplay(self, 'isesa', 'indicates whether an elastic adjustment model is used in the transient'))
     42        s += '{}\n'.format(fielddisplay(self, 'isdamageevolution', 'indicates whether damage evolution is used in the transient'))
     43        s += '{}\n'.format(fielddisplay(self, 'ismovingfront', 'indicates whether a moving front capability is used in the transient'))
     44        s += '{}\n'.format(fielddisplay(self, 'ishydrology', 'indicates whether an hydrology model is used'))
     45        s += '{}\n'.format(fielddisplay(self, 'isslr', 'indicates if a sea level rise solution is used in the transient'))
     46        s += '{}\n'.format(fielddisplay(self, 'isoceancoupling', 'indicates whether coupling with an ocean model is used in the transient'))
     47        s += '{}\n'.format(fielddisplay(self, 'iscoupler', 'indicates whether different models are being run with need for coupling'))
     48        s += '{}\n'.format(fielddisplay(self, 'amr_frequency', 'frequency at which mesh is refined in simulations with multiple time_steps'))
     49        s += '{}\n'.format(fielddisplay(self, 'requested_outputs', 'list of additional outputs requested'))
     50        return s
    5351    #}}}
    5452
    5553    def defaultoutputs(self, md):  # {{{
    56 
    5754        if self.issmb:
    5855            return ['SmbMassBalance']
    5956        else:
    6057            return []
    61 
    6258    #}}}
    6359
    6460    def deactivateall(self):  #{{{
    65         self.issmb = False
    66         self.ismasstransport = False
    67         self.isstressbalance = False
    68         self.isthermal = False
    69         self.isgroundingline = False
    70         self.isgia = False
    71         self.isesa = False
    72         self.isdamageevolution = False
    73         self.ismovingfront = False
    74         self.ishydrology = False
    75         self.isslr = False
    76         self.isoceancoupling = False
    77         self.iscoupler = False
     61        self.issmb = 0
     62        self.ismasstransport = 0
     63        self.isstressbalance = 0
     64        self.isthermal = 0
     65        self.isgroundingline = 0
     66        self.isgia = 0
     67        self.isesa = 0
     68        self.isdamageevolution = 0
     69        self.ismovingfront = 0
     70        self.ishydrology = 0
     71        self.isslr = 0
     72        self.isoceancoupling = 0
     73        self.iscoupler = 0
    7874        self.amr_frequency = 0
    7975
    80     #default output
     76        # Default output
    8177        self.requested_outputs = []
    8278        return self
    … …  
    8581    def setdefaultparameters(self):  # {{{
    8682        #full analysis: Stressbalance, Masstransport and Thermal but no groundingline migration for now
    87         self.issmb = True
    88         self.ismasstransport = True
    89         self.isstressbalance = True
    90         self.isthermal = True
    91         self.isgroundingline = False
    92         self.isgia = False
    93         self.isesa = False
    94         self.isdamageevolution = False
    95         self.ismovingfront = False
    96         self.ishydrology = False
    97         self.isslr = False
    98         self.isoceancoupling = False
    99         self.iscoupler = False
     83        self.issmb = 1
     84        self.ismasstransport = 1
     85        self.isstressbalance = 1
     86        self.isthermal = 1
     87        self.isgroundingline = 0
     88        self.isgia = 0
     89        self.isesa = 0
     90        self.isdamageevolution = 0
     91        self.ismovingfront = 0
     92        self.ishydrology = 0
     93        self.isslr = 0
     94        self.isoceancoupling = 0
     95        self.iscoupler = 0
    10096        self.amr_frequency = 0
    10197
    102     #default output
     98        # Default output
    10399        self.requested_outputs = ['default']
    104100        return self
    … …  
    124120        md = checkfield(md, 'fieldname', 'transient.iscoupler', 'numel', [1], 'values', [0, 1])
    125121        md = checkfield(md, 'fieldname', 'transient.amr_frequency', 'numel', [1], '>=', 0, 'NaN', 1, 'Inf', 1)
    126         md = checkfield(md, 'fieldname', 'transient.requested_outputs', 'stringrow', 1)
    127122
    128         if (solution != 'TransientSolution') and (md.transient.iscoupling):
     123        if solution != 'TransientSolution' and md.transient.iscoupling:
    129124            md.checkmessage("Coupling with ocean can only be done in transient simulations!")
    130         if (md.transient.isdamageevolution and not hasattr(md.materials, 'matdamageice')):
     125        if md.transient.isdamageevolution and not hasattr(md.materials, 'matdamageice'):
    131126            md.checkmessage("requesting damage evolution but md.materials is not of class matdamageice")
    132 
    133127        return md
    134128    # }}}
    … …  
    150144        WriteData(fid, prefix, 'object', self, 'fieldname', 'amr_frequency', 'format', 'Integer')
    151145
    152     #process requested outputs
     146        # Process requested outputs
    153147        outputs = self.requested_outputs
    154148        indices = [i for i, x in enumerate(outputs) if x == 'default']

  • issm/trunk-jpl/src/m/miscellaneous/pretty_print.m

    r25455 r25688  
    1515%   - Add an argument that allows the user to specify the number of values that
    1616%   they would like to display from the head and tail of each dimension of
    17 %   'data'.
     17%   'data' (default is 3 from each end).
     18%   - Add an argument that allows the user to designate the number of
     19%   significant figures to print for each floating point value (default is 8).
    1820
    1921if ndims(data)==1
    2022        if length(data)>6
    21                 disp(sprintf('[%d %d %d ... %d %d %d]',data(1),data(2),data(3),data(end-2),data(end-1),data(end)));
     23                output=sprintf('[%.8f %.8f %.8f ... %.8f %.8f %.8f]',data(1),data(2),data(3),data(end-2),data(end-1),data(end));
    2224        else
    23                 disp(data);
     25                output=sprintf('%.8f',data);
    2426        end
    2527elseif ndims(data)==2
    2628        shape=size(data);
    2729        if shape(1)>6
    28                 % if shape(2)==1
    29                 %       disp('NOTE: Single column printed as row!');
    30                 %       disp(sprintf('[%d %d %d ... %d %d %d]',data(1,1),data(2,1),data(3,1),data(end-2,1),data(end-1,1),data(end,1)));
    31                 % else
    3230                if shape(2)>6
    33                         disp('[');
    34                         disp(sprintf('[%d %d %d ... %d %d %d]',data(1,1),data(1,2),data(1,3),data(1,end-2),data(1,end-1),data(1,end)));
    35                         disp(sprintf('[%d %d %d ... %d %d %d]',data(2,1),data(2,2),data(2,3),data(2,end-2),data(2,end-1),data(2,end)));
    36                         disp(sprintf('[%d %d %d ... %d %d %d]',data(3,1),data(3,2),data(3,3),data(3,end-2),data(3,end-1),data(3,end)));
    37                         disp('...');
    38                         disp(sprintf('[%d %d %d ... %d %d %d]',data(end-2,1),data(end-2,2),data(end-2,3),data(end-2,end-2),data(end-2,end-1),data(end-2,end)));
    39                         disp(sprintf('[%d %d %d ... %d %d %d]',data(end-1,1),data(end-1,2),data(end-1,3),data(end-1,end-2),data(end-1,end-1),data(end-1,end)));
    40                         disp(sprintf('[%d %d %d ... %d %d %d]',data(end,1),data(end,2),data(end,3),data(end,end-2),data(end,end-1),data(end,end)));
    41                         disp(sprintf(']\n'));
     31                        output=sprintf('[[%.8f %.8f %.8f ... %.8f %.8f %.8f]\n',data(1,1),data(1,2),data(1,3),data(1,end-2),data(1,end-1),data(1,end));
     32                        output=sprintf('%s [%.8f %.8f %.8f ... %.8f %.8f %.8f]\n',data(2,1),data(2,2),data(2,3),data(2,end-2),data(2,end-1),data(2,end));
     33                        output=sprintf('%s [%.8f %.8f %.8f ... %.8f %.8f %.8f]\n',data(3,1),data(3,2),data(3,3),data(3,end-2),data(3,end-1),data(3,end));
     34                        output=sprintf('%s ...\n',output);
     35                        output=sprintf('%s [%.8f %.8f %.8f ... %.8f %.8f %.8f]\n',data(end-2,1),data(end-2,2),data(end-2,3),data(end-2,end-2),data(end-2,end-1),data(end-2,end));
     36                        output=sprintf('%s [%.8f %.8f %.8f ... %.8f %.8f %.8f]\n',data(end-1,1),data(end-1,2),data(end-1,3),data(end-1,end-2),data(end-1,end-1),data(end-1,end));
     37                        output=sprintf('%s [%.8f %.8f %.8f ... %.8f %.8f %.8f]]',data(end,1),data(end,2),data(end,3),data(end,end-2),data(end,end-1),data(end,end));
    4238                else
    43                         disp('[');
    44                         disp(data(1,:));
    45                         disp(data(2,:));
    46                         disp(data(3,:));
    47                         disp('...');
    48                         disp(data(end-2,:));
    49                         disp(data(end-1,:));
    50                         disp(data(end,:));
    51                         disp(sprintf(']\n'));
     39                        output=sprintf('[[%.8f]\n',data(1,:));
     40                        output=sprintf('%s [%.8f]\n',output,data(2,:));
     41                        output=sprintf('%s [%.8f]\n',output,data(3,:));
     42                        output=sprintf('%s ...\n',output);
     43                        output=sprintf('%s [%.8f]\n',output,data(end-2,:));
     44                        output=sprintf('%s [%.8f]\n',output,data(end-1,:));
     45                        output=sprintf('%s [%.8f]]',output,data(end,:));
    5246                end
    5347        else
    54                 % if shape(2)==1
    55                 %       data=transpose(data)
    56                 %       sprintf('% NOTE: Single column printed as row!');
    57                 %       % Get shape of transposed structure
    58                 %       shape=size(data);
    59                 %       if shape(2)>6
    60                 %               disp(sprintf('[%d %d %d ... %d %d %d]',data(1,1),data(2,1),data(3,1),data(end-2,1),data(end-1,1),data(end,1)));
    61                 %       else
    62                 %               disp(data);
    63                 %       end
    64                 % else
    6548                if shape(2)>6
    66                         disp('[');
    6749                        for i=1:shape(1)
    68                                 disp(sprintf('[%d %d %d ... %d %d %d]',data(i,1),data(i,2),data(i,3),data(i,end-2),data(i,end-1),data(i,end)));
     50                                if i==1
     51                                        output='[';
     52                                else
     53                                        output=sprintf('%s ',output);
     54                                end
     55
     56                                output=sprintf('%s[%.8f %.8f %.8f ... %.8f %.8f %.8f]',output,data(i,1),data(i,2),data(i,3),data(i,end-2),data(i,end-1),data(i,end));
     57
     58                                if i==shape(1)
     59                                        output=sprintf('%s]',output);
     60                                end
    6961                        end
    70                         disp(sprintf(']\n'));
    7162                else
    72                         disp(data);
     63                        output=sprintf('%.8f',data);
    7364                end
    7465        end
    7566else
    76         disp(data);
     67        output=sprintf('%.8f',data);
    7768end
     69
     70disp(output);

  • issm/trunk-jpl/src/m/qmu/dakota_out_parse.m

    r24759 r25688  
    9090                        [ntokens,tokens]=fltokens(fline);
    9191                        method=tokens{1}{1};
    92                         display(sprintf('Dakota method =''%s''.',method));
     92                        display(sprintf('Dakota method = ''%s''.',method));
    9393                elseif strcmp(tokens{1}{1}(7),'N')
    9494                        [fline]=findline(fidi,'methodName = ');
    9595                        [ntokens,tokens]=fltokens(fline);
    9696                        method=tokens{1}{3};
    97                         display(sprintf('Dakota methodName=''%s''.',method));
     97                        display(sprintf('Dakota methodName = ''%s''.',method));
    9898                end
    9999        end

  • issm/trunk-jpl/src/m/qmu/dakota_out_parse.py

    r24213 r25688  
    8787                [ntokens, tokens] = fltokens(fline)
    8888                method = tokens[0].strip()
    89                 print('Dakota method =\'' + method + '\'.')
     89                print('Dakota method = \'' + method + '\'.')
    9090            elif fline[6] in ['N', 'n']:
    9191                fline = findline(fidi, 'methodName = ')
    9292                [ntokens, tokens] = fltokens(fline)
    9393                method = tokens[2].strip()
    94                 print('Dakota methodName = "' + method + '".')
     94                print('Dakota methodName = \'' + method + '\'.')
    9595
    9696    #  loop through the file to find the function evaluation summary

  • issm/trunk-jpl/src/m/qmu/postqmu.py

    r25685 r25688  
     1from copy import deepcopy
    12from os import getpid, stat
    23from os.path import isfile
    … …  
    56from dakota_out_parse import *
    67from helpers import *
    7 from loadresultsfromdisk import *
     8import loadresultsfromdisk as loadresults
    89
    910
    … …  
    4849        if md.qmu.method.method == 'nond_sampling':
    4950            dakotaresults.modelresults = []
    50             md2 = copy.deepcopy(md)
     51            md2 = deepcopy(md)
    5152            md2.qmu.isdakota = 0
    5253            for i in range(md2.qmu.method.params.samples):
    53                 print('reading qmu file {}.outbin.{}'.format(md2.miscellaneous.name, i))
    54                 md2 = loadresultsfromdisk(md2, '{}.outbin.{}'.format(md2.miscellaneous.name, i))
     54                outbin_name = '{}.outbin.{}'.format(md2.miscellaneous.name, (i + 1))
     55                print('reading qmu file {}'.format(outbin_name))
     56                md2 = loadresults.loadresultsfromdisk(md2, outbin_name)
    5557                dakotaresults.modelresults.append(md2.results)
    5658
    5759    if md.qmu.statistics.method[0]['name'] != 'None':
    5860        md.qmu.isdakota = 0
    59         md = loadresultsfromdisk(md, [md.miscellaneous.name,'.stats'])
     61        md = loadresults.loadresultsfromdisk(md, '{}.stats'.format(md.miscellaneous.name))
    6062        md.qmu.isdakota = 1
    6163

  • issm/trunk-jpl/src/m/qmu/preqmu.m

    r25328 r25688  
    101101                variablepartitions{end+1}=fieldvariable.partition;
    102102                variablepartitions_npart(end+1)=qmupart2npart(fieldvariable.partition);
    103                 if isfield(struct(fieldvariable),'nsteps'),
     103                if isprop(fieldvariable,'nsteps'),
    104104                        variablepartitions_nt(end+1)=fieldvariable.nsteps;
    105105                else

  • issm/trunk-jpl/src/m/qmu/preqmu.py

    r25632 r25688  
    1111
    1212def preqmu(md, options):
    13     """QMU - apply Quantification of Margins and Uncertainties techniques
     13    """PREQMU - apply Quantification of Margins and Uncertainties techniques
    1414    to a solution sequence (like stressbalance.py, progonstic.py, etc ...),
    1515    using the Dakota software from Sandia.
    … …  
    3636    options.addfielddefault('iparams', 0)
    3737
    38     # when running in library mode, the in file needs to be called md.miscellaneous.name.qmu.in
     38    # When running in library mode, the in file needs to be called md.miscellaneous.name.qmu.in
    3939    qmufile = md.miscellaneous.name
    4040
    41     # retrieve variables and resposnes for this particular analysis.
     41    # Retrieve variables and resposnes for this particular analysis.
    4242    #print type(md.qmu.variables)
    4343    #print md.qmu.variables.__dict__
    44     #print ivar
     44    # Print ivar
    4545    variables = md.qmu.variables  #[ivar]
    4646    responses = md.qmu.responses  #[iresp]
    4747
    48     # expand variables and responses
     48    # Expand variables and responses
    4949    #print variables.__dict__
    5050    #print responses.__dict__
    … …  
    5252    responses = expandresponses(md, responses)
    5353
    54     # go through variables and responses, and check they don't have more than
    55     #   md.qmu.numberofpartitions values. Also determine numvariables and numresponses
     54    # Go through variables and responses, and check they don't have more than
     55    # md.qmu.numberofpartitions values. Also determine numvariables and
     56    # numresponses
    5657    #{{{
    5758    numvariables = 0
    … …  
    8283    #}}}
    8384
    84     # create in file for dakota
     85    # Create in file for Dakota
    8586    #dakota_in_data(md.qmu.method[imethod], variables, responses, md.qmu.params[iparams], qmufile)
    8687    dakota_in_data(md.qmu.method, variables, responses, md.qmu.params, qmufile)
    8788
    88     # build a list of variables and responses descriptors. the list is not expanded.
     89    # Build a list of variables and responses descriptors. the list is not expanded.
    8990    #{{{
    9091    variabledescriptors = []
    … …  
    115116    #}}}
    116117
    117     #build a list of variable partitions
     118    # Build a list of variable partitions
    118119    variablepartitions = []
    119120    variablepartitions_npart = []
    … …  
    125126        if type(fieldvariable) in [list, np.ndarray]:
    126127            for j in range(np.size(fieldvariable)):
    127                 if fieldvariable[j].isscaled():
     128                if fieldvariable[j].isscaled() or fieldvariable[j].isdistributed():
    128129                    variablepartitions.append(fieldvariable[j].partition)
    129130                    variablepartitions_npart.append(qmupart2npart(fieldvariable[j].partition))
    130                     variablepartitions_nt.append(fieldvariable[j].nsteps)
     131                    if hasattr(fieldvariable[j], 'nsteps'):
     132                        variablepartitions_nt.append(fieldvariable[j].nsteps)
     133                    else:
     134                        variablepartitions_nt.append(1)
    131135                else:
    132136                    variablepartitions.append([])
    … …  
    137141                variablepartitions.append(fieldvariable.partition)
    138142                variablepartitions_npart.append(qmupart2npart(fieldvariable.partition))
    139                 variablepartitions_nt.append(fieldvariable.nsteps)
     143                if hasattr(fieldvariable, 'nsteps'):
     144                    variablepartitions_nt.append(fieldvariable.nsteps)
     145                else:
     146                    variablepartitions_nt.append(1)
    140147            else:
    141148                variablepartitions.append([])
    … …  
    143150                variablepartitions_nt.append(1)
    144151
    145     #build a list of response partitions
     152    # Build a list of response partitions
    146153    responsepartitions = []
    147     responsepartitions_npart = []
     154    responsepartitions_npart = np.array([])
    148155    response_fieldnames = fieldnames(md.qmu.responses)
    149156    for i in range(len(response_fieldnames)):
    … …  
    154161                if fieldresponse[j].isscaled():
    155162                    responsepartitions.append(fieldresponse[j].partition)
    156                     responsepartitions_npart.append(qmupart2npart(fieldresponse[j].partition))
     163                    responsepartitions_npart = np.append(responsepartitions_npart, qmupart2npart(fieldresponse[j].partition))
    157164                else:
    158165                    responsepartitions.append([])
    159                     responsepartitions_npart.append(0)
     166                    responsepartitions_npart = np.append(responsepartitions_npart, 0)
    160167        else:
    161168            if fieldresponse.isscaled():
    162169                responsepartitions.append(fieldresponse.partition)
    163                 responsepartitions_npart.append(qmupart2npart(fieldresponse.partition))
     170                responsepartitions_npart = np.append(responsepartitions_npart, qmupart2npart(fieldresponse.partition))
    164171            else:
    165172                responsepartitions.append([])
    166                 responsepartitions_npart.append(0)
     173                responsepartitions_npart = np.append(responsepartitions_npart, 0)
    167174
    168     # register the fields that will be needed by the Qmu model.
     175    if responsepartitions_npart.shape[0] != 1:
     176        responsepartitions_npart = responsepartitions_npart.reshape(1, -1)
     177
     178    # Register the fields that will be needed by the Qmu model.
    169179    md.qmu.numberofresponses = numresponses
    170180    md.qmu.variabledescriptors = variabledescriptors
    … …  
    176186    md.qmu.responsepartitions_npart = responsepartitions_npart
    177187
    178     # now, we have to provide all the info necessary for the solutions to compute the
    179     # responses. For ex, if mass_flux is a response, we need a profile of points.
    180     # For a misfit, we need the observed velocity, etc ...
     188    # Now, we have to provide all the info necessary for the solutions to
     189    # compute the responses. For example, if mass_flux is a response, we need a
     190    # profile of points. For a misfit, we need the observed velocity, etc.
    181191    md = process_qmu_response_data(md)
    182192

  • issm/trunk-jpl/src/m/qmu/qmupart2npart.m

    r24989 r25688  
    11function npart=qmupart2npart(vector)
    2 
    3         %vector is full of -1 (no partition) and 0 to npart.  We need to
    4         %identify npart
    5 
     2        %vector is full of -1 (no partition) and 0 to npart. We need to identify npart=
    63        npart=max(vector)+1;
    7 
    8 

  • issm/trunk-jpl/src/m/qmu/qmupart2npart.py

    r25011 r25688  
    1 import numpy as np
    2 
    3 
    41def qmupart2npart(vector):
    52    # Vector is full of -1 (no partition) and 0 to npart. We need to identify
    63    # npart.
    7 
    8     npart = vector.max() + 1
     4    npart = int(vector.max() + 1) # cast to int as we may have a NumPy floating point type, which cannot be used as an argument to function range (see src/m/qmu/setupdesign/QmuSetupVariables.py)
    95
    106    return npart

  • issm/trunk-jpl/src/m/qmu/setupdesign/QmuSetupVariables.py

    r25031 r25688  
    11from copy import deepcopy
     2
    23from helpers import *
    34from MatlabFuncs import *
    … …  
    89
    910def QmuSetupVariables(md, variables):
    10     #get descriptor
     11    """QMUSETUPVARIABLES function
     12    """
     13
     14    # Get descriptor
    1115    descriptor = variables.descriptor
    1216
    13     #decide whether this is a distributed variable, which will drive whether we expand it into npart values,
    14     #or if we just carry it forward as is.
     17    # Decide whether this is a distributed variable, which will drive whether we expand it into npart values, or if we just carry it forward as is
    1518
     19    # Ok, key off according to type of descriptor
    1620    dvar = []
     21    if strncmp(descriptor, 'scaled_', 7):
     22        # We have a scaled variable, expand it over the partition. First recover the partition.
     23        partition = variables.partition
     24        # Figure out number of partitions
     25        npart = qmupart2npart(partition)
    1726
    18     #ok, key off according to type of descriptor:
    19     if strncmp(descriptor, 'scaled_', 7):
    20         #we have a scaled variable, expand it over the partition. First recover the partition.
    21         partition = variables.partition
    22         #figure out number of partitions
    23         npart = qmupart2npart(partition)
    24         #figure out number of time steps
     27        # Figure out number of time steps
    2528        nt = variables.nsteps
    2629
    … …  
    4447                raise RuntimeError('QmuSetupVariables error message: stddev and mean fields should have the same number of cols as the number of time steps')
    4548
    46         #ok, dealing with semi-discrete distributed variable. Distribute according to how many
    47         #partitions we want, and number of time steps
     49        # Ok, dealing with semi-discrete distributed variable. Distribute according to how many partitions we want, and number of time steps
    4850        if nt == 1:
    4951            for j in range(npart):
    5052                dvar.append(deepcopy(variables))
    5153
    52                 # text parsing in dakota requires literal "'identifier'" not just "identifier"
     54                # Text parsing in Dakota requires literal "'identifier'" not just "identifier"
    5355                dvar[-1].descriptor = "'" + str(variables.descriptor) + '_' + str(j + 1) + "'"
    5456
    … …  
    6466                    dvar.append(deepcopy(variables))
    6567
    66                     # text parsing in dakota requires literal "'identifier'" not just "identifier"
     68                    # Text parsing in Dakota requires literal "'identifier'" not just "identifier"
    6769                    dvar[-1].descriptor = "'" + str(variables.descriptor) + '_' + str(j + 1) + '_' + str(k + 1) + "'"
    6870
    … …  
    7678        dvar.append(deepcopy(variables))
    7779
    78         # text parsing in dakota requires literal "'identifier'" not just "identifier"
     80        # text parsing in Dakota requires literal "'identifier'" not just "identifier"
    7981        dvar[-1].descriptor = "'" + str(variables.descriptor) + "'"
    8082

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

    r25515 r25688  
    164164        for i in range(s[0]):
    165165            if np.ndim(data) == 1:
    166                 fid.write(pack('d', float(data[i])))  #get to the "c" convention, hence the transpose
     166                fid.write(pack('d', float(data[i])))
    167167            else:
    168168                for j in range(s[1]):
    169                     fid.write(pack('d', float(data[i][j])))  #get to the "c" convention, hence the transpose
     169                    fid.write(pack('d', float(data[i][j])))
    170170    # }}}
    171171
    … …  
    352352
    353353def FormatToCode(datatype):  # {{{
    354     """
    355     This routine takes the datatype string, and hardcodes it into an integer, which
    356     is passed along the record, in order to identify the nature of the dataset being
    357     sent.
     354    """ FORMATTOCODE - This routine takes the datatype string, and hardcodes it
     355    into an integer, which is passed along the record, in order to identify the
     356    nature of the dataset being sent.
    358357    """
    359358    if datatype == 'Boolean':

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

    r25685 r25688  
    4343        if md.qmu.output and md.qmu.statistics.method[0]['name'] == 'None':
    4444            if md.qmu.method.method == 'nond_sampling':
    45                 for i in range(len(md.qmu.method.params.samples)):
    46                     filelist.append(md.miscellaneous.name + '.outbin' + str(i))
     45                for i in range(md.qmu.method.params.samples):
     46                    filelist.append(md.miscellaneous.name + '.outbin.' + str(i + 1))
    4747        if md.qmu.statistics.method[0]['name'] != 'None':
    4848            filelist.append(md.miscellaneous.name + '.stats')

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

    r25627 r25688  
    6262
    6363        if ~isempty(md.results.(structure(1).SolutionType)(1).errlog),
    64                 disp(['loadresultsfromcluster info message: error during solution. Check your errlog and outlog model fields']);
     64                disp(['loadresultsfromdisk info message: error during solution. Check your errlog and outlog model fields']);
    6565        end
    6666

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

    r25685 r25688  
    11import os
     2
     3import numpy as np # Remove: used temporarily for debugging
    24
    35from helpers import fieldnames
    … …  
    4042        if not structure:
    4143            raise RuntimeError("No result found in binary file '{}'. Check for solution crash.".format(filename))
    42         if not structure[0].SolutionType:
    43             if structure[-1].SolutionType:
     44        if not hasattr(structure[0], 'SolutionType'):
     45            if hasattr(structure[-1], 'SolutionType'):
    4446                structure[0].SolutionType = structure[-1].SolutionType
    4547            else:
    … …  
    6567
    6668        if getattr(md.results, structure[0].SolutionType)[0].errlog:
    67             print("loadresultsfromcluster info message: error during solution. Check your errlog and outlog model fields.")
     69            print("loadresultsfromdisk info message: error during solution. Check your errlog and outlog model fields.")
    6870
    6971        # If only one solution, extract it from list for user friendliness

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

    r25499 r25688  
    1818end
    1919
    20 %Go through all model fields: check that it is a class and call checkconsistency
    21 fields=properties('model');
     20% Go through all model fields: check that it is a class and call checkconsistency
     21fields=sort(properties('model')); %sort fields so that comparison of binary files is easier
    2222for i=1:length(fields),
    2323        field=fields{i};
    … …  
    4343%close file
    4444st=fclose(fid);
    45 % % Uncomment the following to make a copy of the binary input file for debugging
    46 % % purposes (can be fed into scripts/ReadBin.py).
    47 copyfile([md.miscellaneous.name '.bin'], [md.miscellaneous.name '.m.bin'])
     45
     46% Uncomment the following to make a copy of the binary input file for debugging
     47% purposes (can be fed into scripts/ReadBin.py).
     48% copyfile([md.miscellaneous.name '.bin'], [md.miscellaneous.name '.m.bin'])
     49
    4850if st==-1,
    4951        error(['marshall error message: could not close file ' [md.miscellaneous.name '.bin']]);

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

    r25499 r25688  
    55
    66def marshall(md):
    7     '''
    8     MARSHALL - outputs a compatible binary file from @model md, for certain solution type.
     7    """MARSHALL - outputs a compatible binary file from @model md, for certain solution type.
    98
    10         The routine creates a compatible binary file from @model md
    11         This binary file will be used for parallel runs in JPL-package
     9    The routine creates a compatible binary file from @model md
     10    his binary file will be used for parallel runs in JPL-package
    1211
    13         Usage:
    14             marshall(md)
    15     '''
     12    Usage:
     13        marshall(md)
     14    """
     15
    1616    if md.verbose.solution:
    1717        print("marshalling file {}.bin".format(md.miscellaneous.name))
    … …  
    2323        raise IOError("marshall error message: could not open '%s.bin' file for binary writing. Due to: ".format(md.miscellaneous.name), e)
    2424
    25     for field in md.properties():
     25    fields = md.properties()
     26    fields.sort() # sort fields so that comparison of binary files is easier
     27    for field in fields:
    2628        #Some properties do not need to be marshalled
    2729        if field in ['results', 'radaroverlay', 'toolkits', 'cluster', 'private']:
    … …  
    4244    try:
    4345        fid.close()
    44         # # Uncomment the following to make a copy of the binary input file for
    45         # # debugging purposes (can be fed into scripts/ReadBin.py).
    46         copy_cmd = "cp {}.bin {}.py.bin".format(md.miscellaneous.name, md.miscellaneous.name)
    47         subprocess.call(copy_cmd, shell=True)
     46
     47        # Uncomment the following to make a copy of the binary input file for
     48        # debugging purposes (can be fed into scripts/ReadBin.py).
     49        # copy_cmd = "cp {}.bin {}.py.bin".format(md.miscellaneous.name, md.miscellaneous.name)
     50        # subprocess.call(copy_cmd, shell=True)
     51
    4852    except IOError as e:
    4953        print("marshall error message: could not close file '{}.bin' due to:".format(md.miscellaneous.name), e)

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

    r25627 r25688  
    1 function results=parseresultsfromdisk(md,filename,iosplit)
     1function results=parseresultsfromdisk(md,filename,iosplit) % {{{
    22
    33if iosplit,
    … …  
    77        %results=parseresultsfromdiskioserialsequential(md,filename);
    88end
    9 
    10 function results=parseresultsfromdiskioserialsequential(md,filename) % {{{
     9% }}}
     10function results=parseresultsfromdiskiosplit(md,filename) % {{{
    1111
    1212%Open file
    … …  
    1515        error(['loadresultsfromdisk error message: could not open ',filename,' for binary reading']);
    1616end
    17 
    18 %first pass to figure out the steps we have:
    19 steps=[];
    20 while  1,
    21         result  = ReadDataDimensions(fid);
    22         if isempty(result),
    23                 break;
    24         end
    25         if result.step~=-9999,
    26                 steps=[steps result.step];
    27         end
    28 end
    29 
    30 steps=unique(steps);
    31 
    32 %create structure:
    33 results=struct('step',num2cell(steps));
    34 
    35 %second pass to fill the steps we have:
     17results=struct();
     18
     19%if we have done split I/O, ie, we have results that are fragmented across patches,
     20%do a first pass, and figure out the structure of results
     21result=ReadDataDimensions(fid);
     22while ~isempty(result),
     23
     24        %Get time and step
     25        results(result.step).step=result.step;
     26        if result.time~=-9999,
     27                results(result.step).time=result.time;
     28        end
     29
     30        %Add result
     31        results(result.step).(result.fieldname)=NaN;
     32
     33        %read next result
     34        result=ReadDataDimensions(fid);
     35end
     36
     37%do a second pass, and figure out the size of the patches
    3638fseek(fid,0,-1); %rewind
    37 while  1,
    38         result  = ReadData(fid,md);
    39         if isempty(result),
    40                 break;
    41         end
    42         if result.step==-9999,
    43                 result.step=1;
    44                 result.time=0;
    45         end
    46         %find where we are putting this step:
    47         ind=find(steps==result.step);
    48         if isempty(ind),
    49                 error('could not find a step in our pre-structure! Something is very off!');
    50         end
    51 
    52         %plug:
    53         results(ind).time=result.time;
    54         results(ind).(result.fieldname)=result.field;
    55 end
     39result=ReadDataDimensions(fid);
     40while ~isempty(result),
     41        %read next result
     42        result=ReadDataDimensions(fid);
     43end
     44
     45%third pass, this time to read the real information
     46fseek(fid,0,-1); %rewind
     47result=ReadData(fid,md);
     48while ~isempty(result),
     49
     50        %Get time and step
     51        results(result.step).step=result.step;
     52        if result.time~=-9999,
     53                results(result.step).time=result.time;
     54        end
     55
     56        %Add result
     57        results(result.step).(result.fieldname)=result.field;
     58
     59        %read next result
     60        try,
     61                result=ReadData(fid,md);
     62        catch me,
     63                disp('WARNING: file corrupted, results partial recovery');
     64                result=[];
     65        end
     66
     67end
     68
     69%close file
    5670fclose(fid);
    5771% }}}
    … …  
    7185        %read next result
    7286        try,
    73                 result  = ReadData(fid,md);
     87                result = ReadData(fid,md);
    7488        catch me,
    7589                disp('WARNING: file corrupted, trying partial recovery');
    … …  
    92106fclose(fid);
    93107
    94 %Now, process all results and find how many steps we have
     108%Now, process all results and find out how many steps we have
    95109numresults = numel(allresults);
    96110allsteps   = zeros(numresults,1);
    … …  
    105119for i=1:numresults
    106120        result = allresults{i};
    107 
    108121        index = 1;
    109         if result.step ~= -9999
     122        if(result.step ~= -9999)
    110123                index = find(result.step == allsteps);
    111         end
    112 
    113         if(result.step~=-9999) results(index).step=result.step; end
    114         if(result.time~=-9999) results(index).time=result.time; end
     124                results(index).step=result.step;
     125        end
     126        if(result.time~=-9999)
     127                results(index).time=result.time;
     128        end
    115129        results(index).(result.fieldname)=result.field;
    116130end
    117131% }}}
    118 function results=parseresultsfromdiskiosplit(md,filename) % {{{
     132function results=parseresultsfromdiskioserialsequential(md,filename) % {{{
    119133
    120134%Open file
    … …  
    123137        error(['loadresultsfromdisk error message: could not open ',filename,' for binary reading']);
    124138end
    125 results=struct();
    126 
    127 %if we have done split I/O, ie, we have results that are fragmented across patches,
    128 %do a first pass, and figure out the structure of results
    129 result=ReadDataDimensions(fid);
    130 while ~isempty(result),
    131 
    132         %Get time and step
    133         results(result.step).step=result.step;
    134         if result.time~=-9999,
    135                 results(result.step).time=result.time;
    136         end
    137 
    138         %Add result
    139         results(result.step).(result.fieldname)=NaN;
    140 
    141         %read next result
    142         result=ReadDataDimensions(fid);
    143 end
    144 
    145 %do a second pass, and figure out the size of the patches
     139
     140%first pass to figure out the steps we have:
     141steps=[];
     142while  1,
     143        result  = ReadDataDimensions(fid);
     144        if isempty(result),
     145                break;
     146        end
     147        if result.step~=-9999,
     148                steps=[steps result.step];
     149        end
     150end
     151
     152steps=unique(steps);
     153
     154%create structure:
     155results=struct('step',num2cell(steps));
     156
     157%second pass to fill the steps we have:
    146158fseek(fid,0,-1); %rewind
    147 result=ReadDataDimensions(fid);
    148 while ~isempty(result),
    149         %read next result
    150         result=ReadDataDimensions(fid);
    151 end
    152 
    153 %third pass, this time to read the real information
    154 fseek(fid,0,-1); %rewind
    155 result=ReadData(fid,md);
    156 while ~isempty(result),
    157 
    158         %Get time and step
    159         results(result.step).step=result.step;
    160         if result.time~=-9999,
    161                 results(result.step).time=result.time;
    162         end
    163 
    164         %Add result
    165         results(result.step).(result.fieldname)=result.field;
    166 
    167         %read next result
    168         try,
    169                 result=ReadData(fid,md);
    170         catch me,
    171                 disp('WARNING: file corrupted, results partial recovery');
    172                 result=[];
    173         end
    174 
    175 end
    176 
    177 %close file
     159while  1,
     160        result  = ReadData(fid,md);
     161        if isempty(result),
     162                break;
     163        end
     164        if result.step==-9999,
     165                result.step=1;
     166                result.time=0;
     167        end
     168        %find where we are putting this step:
     169        ind=find(steps==result.step);
     170        if isempty(ind),
     171                error('could not find a step in our pre-structure! Something is very off!');
     172        end
     173
     174        %plug:
     175        results(ind).time=result.time;
     176        results(ind).(result.fieldname)=result.field;
     177end
    178178fclose(fid);
    179         % }}}
     179%}}}
    180180function result=ReadData(fid,md) % {{{
    181181
    … …  
    312312        end
    313313
     314        if time~=-9999,
     315                time=time/yts;
     316        end
     317
    314318        result.fieldname=fieldname;
    315319        result.time=time;
    316         if result.time~=-9999,
    317                 result.time=time/yts;
    318         end
    319320        result.step=step;
    320321        result.field=field;

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

    r25303 r25688  
     1from collections import OrderedDict
    12import struct
     3
    24import numpy as np
    3 from collections import OrderedDict
     5
    46import results as resultsclass
    57
    68
    7 def parseresultsfromdisk(md, filename, iosplit):
     9def parseresultsfromdisk(md, filename, iosplit): #{{{
    810    if iosplit:
    911        saveres = parseresultsfromdiskiosplit(md, filename)
    1012    else:
    1113        saveres = parseresultsfromdiskioserial(md, filename)
    12     return saveres
    13 
     14        #saveres = parseresultsfromdiskioserialsequential(md, filename)
     15    return saveres
     16#}}}
     17
     18def parseresultsfromdiskiosplit(md, filename):  # {{{
     19    #Open file
     20    try:
     21        fid = open(filename, 'rb')
     22    except IOError as e:
     23        raise IOError("parseresultsfromdisk error message: could not open '{}' for binary reading.".format(filename))
     24
     25    saveres = []
     26
     27    #if we have done split I / O, ie, we have results that are fragmented across patches,
     28    #do a first pass, and figure out the structure of results
     29    loadres = ReadDataDimensions(fid)
     30    while loadres:
     31
     32        #Get time and step
     33        if loadres['step'] > len(saveres):
     34            for i in range(len(saveres), loadres['step'] - 1):
     35                saveres.append(None)
     36            saveres.append(resultsclass.results())
     37        setattr(saveres[loadres['step'] - 1], 'step', loadres['step'])
     38        setattr(saveres[loadres['step'] - 1], 'time', loadres['time'])
     39
     40    #Add result
     41        setattr(saveres[loadres['step'] - 1], loadres['fieldname'], float('NaN'))
     42
     43    #read next result
     44        loadres = ReadDataDimensions(fid)
     45
     46    #do a second pass, and figure out the size of the patches
     47    fid.seek(0)  #rewind
     48    loadres = ReadDataDimensions(fid)
     49    while loadres:
     50
     51        #read next result
     52        loadres = ReadDataDimensions(fid)
     53
     54    #third pass, this time to read the real information
     55    fid.seek(0)  #rewind
     56    loadres = ReadData(fid, md)
     57    while loadres:
     58
     59        #Get time and step
     60        if loadres['step'] > len(saveres):
     61            for i in range(len(saveres), loadres['step'] - 1):
     62                saveres.append(None)
     63            saveres.append(saveresclass.saveres())
     64        setattr(saveres[loadres['step'] - 1], 'step', loadres['step'])
     65        setattr(saveres[loadres['step'] - 1], 'time', loadres['time'])
     66
     67    #Add result
     68        setattr(saveres[loadres['step'] - 1], loadres['fieldname'], loadres['field'])
     69
     70    #read next result
     71        loadres = ReadData(fid, md)
     72
     73    #close file
     74    fid.close()
     75
     76    return saveres
     77# }}}
    1478
    1579def parseresultsfromdiskioserial(md, filename):  # {{{
    … …  
    73137
    74138    return saveres
    75     # }}}
    76 
    77 
    78 def parseresultsfromdiskiosplit(md, filename):  # {{{
    79     #Open file
    80     try:
    81         fid = open(filename, 'rb')
    82     except IOError as e:
    83         raise IOError("loadresultsfromdisk error message: could not open '{}' for binary reading.".format(filename))
    84 
    85     saveres = []
    86 
    87     #if we have done split I / O, ie, we have results that are fragmented across patches,
    88     #do a first pass, and figure out the structure of results
    89     loadres = ReadDataDimensions(fid)
    90     while loadres:
    91 
    92         #Get time and step
    93         if loadres['step'] > len(saveres):
    94             for i in range(len(saveres), loadres['step'] - 1):
    95                 saveres.append(None)
    96             saveres.append(resultsclass.results())
    97         setattr(saveres[loadres['step'] - 1], 'step', loadres['step'])
    98         setattr(saveres[loadres['step'] - 1], 'time', loadres['time'])
    99 
    100     #Add result
    101         setattr(saveres[loadres['step'] - 1], loadres['fieldname'], float('NaN'))
    102 
    103     #read next result
    104         loadres = ReadDataDimensions(fid)
    105 
    106     #do a second pass, and figure out the size of the patches
    107     fid.seek(0)  #rewind
    108     loadres = ReadDataDimensions(fid)
    109     while loadres:
    110 
    111         #read next result
    112         loadres = ReadDataDimensions(fid)
    113 
    114     #third pass, this time to read the real information
    115     fid.seek(0)  #rewind
    116     loadres = ReadData(fid, md)
    117     while loadres:
    118 
    119         #Get time and step
    120         if loadres['step'] > len(saveres):
    121             for i in range(len(saveres), loadres['step'] - 1):
    122                 saveres.append(None)
    123             saveres.append(saveresclass.saveres())
    124         setattr(saveres[loadres['step'] - 1], 'step', loadres['step'])
    125         setattr(saveres[loadres['step'] - 1], 'time', loadres['time'])
    126 
    127     #Add result
    128         setattr(saveres[loadres['step'] - 1], loadres['fieldname'], loadres['field'])
    129 
    130     #read next result
    131         loadres = ReadData(fid, md)
    132 
    133     #close file
    134     fid.close()
    135 
    136     return saveres
    137     # }}}
    138 
     139# }}}
    139140
    140141def ReadData(fid, md):  # {{{
    141     '''
    142     READDATA -
    143 
    144         Usage:
    145            field = ReadData(fid, md)
    146     '''
     142    """READDATA
     143
     144    Usage:
     145        field = ReadData(fid, md)
     146    """
    147147
    148148    #read field
    … …  
    291291
    292292    return saveres
    293     # }}}
     293# }}}
    294294
    295295
    296296def ReadDataDimensions(fid):  # {{{
     297    """READDATADIMENSIONS - read data dimensions, step and time, but not the data itself.
     298
     299    Usage:
     300        field = ReadDataDimensions(fid)
    297301    """
    298     READDATADIMENSIONS - read data dimensions, step and time, but not the data itself.
    299 
    300         Usage:
    301            field = ReadDataDimensions(fid)
    302     """
    303 
    304     #read field
     302
     303    # Read field
    305304    try:
    306305        length = struct.unpack('i', fid.read(struct.calcsize('i')))[0]
    … …  
    310309        datatype = struct.unpack('i', fid.read(struct.calcsize('i')))[0]
    311310        M = struct.unpack('i', fid.read(struct.calcsize('i')))[0]
    312         N = 1  #default
     311        N = 1 # default
    313312        if datatype == 1:
    314313            fid.seek(M * 8, 1)
    … …  
    333332
    334333    return saveres
    335     # }}}
     334# }}}

  • issm/trunk-jpl/src/wrappers/CoordTransform

    • Property svn:ignore set to
      .deps

  • issm/trunk-jpl/test

  • issm/trunk-jpl/test/NightlyRun

    • Property svn:ignore

      •  

        old new  
        1818run.old
        1919run_matlab
         20test218.qmu.in
         21test218.qmu.out
         22test218.qmu.err

  • issm/trunk-jpl/test/NightlyRun/runme.m

    r25670 r25688  
    107107%Process Ids according to benchmarks{{{
    108108if strcmpi(benchmark,'nightly'),
    109         test_ids=intersect(test_ids,[1:999]);
     109        test_ids=intersect(test_ids,union([1:999],[2001:2500]));
    110110elseif strcmpi(benchmark,'validation'),
    111111        test_ids=intersect(test_ids,[1001:1999]);
    … …  
    219219                                        archive_cell=archread(['../Archives/' archive_name '.arch'],[archive_name '_field' num2str(k)]);
    220220                                        archive=archive_cell{1};
    221                                         error_diff=full(max(abs(archive(:)-field(:)))/(max(abs(archive(:)))+eps));
    222 
    223                                         %disp test result
     221                                        error_diff=full(max(abs(archive(:)-field(:)))/(max(abs(archive(:)))+eps));                                      %disp test result
    224222                                        if (error_diff>tolerance | isnan(error_diff));
    225223                                                disp(sprintf(['ERROR   difference: %-7.2g > %7.2g test id: %i test name: %s field: %s'],...

  • issm/trunk-jpl/test/NightlyRun/runme.py

    r25670 r25688  
    114114    #Process Ids according to benchmarks {{{
    115115    if benchmark == 'nightly':
    116         test_ids = test_ids.intersection(set(range(1, 1000)))
     116        test_ids = test_ids.intersection(set(range(1, 1000)).union(set(range(2001, 2500))))
    117117    elif benchmark == 'validation':
    118118        test_ids = test_ids.intersection(set(range(1001, 2000)))

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

    r25181 r25688  
    5252
    5353#make sure wherever there is an ice load, that the mask is set to ice:
    54 #pos = np.nonzero(md.solidearth.surfaceload.icethicknesschange)[0] # Do we need to do this twice?
     54#pos = np.nonzero(md.solidearth.surfaceload.icethicknesschange)[0] # TODO: Do we need to do this twice?
    5555md.mask.ice_levelset[md.mesh.elements[pos, :] - 1] = -1
    5656# }}}
    … …  
    7272md.solidearth.settings.computesealevelchange = 1
    7373
    74 #max number of iteration reverted back to 10 (i.e., the original default value)
     74#max number of iterations reverted back to 10 (i.e., the original default value)
    7575md.solidearth.settings.maxiter = 10
    7676

  • issm/trunk-jpl/test/NightlyRun/test2005.m

    r25627 r25688  
    4040
    4141%make sure wherever there is an ice load, that the mask is set to ice:
    42 pos=find(md.solidearth.surfaceload.icethicknesschange);
     42%pos=find(md.solidearth.surfaceload.icethicknesschange); % TODO: Do we need to do this twice?
    4343md.mask.ice_levelset(md.mesh.elements(pos,:))=-1;
    4444% }}}
    … …  
    5151%Materials:
    5252md.materials=materials('hydro');
    53 
    5453
    5554%Miscellaneous
    … …  
    8685deltathickness(end,:)=0:1:9;
    8786md.solidearth.surfaceload.icethicknesschange=deltathickness;
    88 
    8987%hack:
    9088md.geometry.surface=zeros(md.mesh.numberofvertices,1);
    … …  
    103101%Fields and tolerances to track changes
    104102field_names={'Sealevel1','Sealevel5','Sealevel10','Seustatic10'};
    105 field_tolerances={1e-13,1e-13,1e-13,1e-13,1e-13};
     103field_tolerances={1e-13,1e-13,1e-13,1e-13};
    106104field_values={S1,S5,S10,Seus10};
    107105

  • issm/trunk-jpl/test/NightlyRun/test2006.m

    r25627 r25688  
    5252%Materials:
    5353md.materials=materials('hydro');
    54 
    5554
    5655%Miscellaneous
    … …  
    138137end
    139138% }}}
    140         %algorithm:  % {{{
    141         md.qmu.method     =dakota_method('nond_samp');
    142         md.qmu.method(end)=dmeth_params_set(md.qmu.method(end),...
    143         'seed',1234,...
    144         'samples',10,...
    145         'sample_type','random');
    146         md.qmu.output=1;
    147         %}}}
    148         %parameters % {{{
    149         md.qmu.params.direct=true;
    150         md.qmu.params.interval_type='forward';
    151         md.qmu.params.analysis_driver='matlab';
    152         md.qmu.params.evaluation_scheduling='master';
    153         md.qmu.params.processors_per_evaluation=2;
    154         md.qmu.params.tabular_graphics_data=true;
    155         md.qmu.isdakota=1;
    156         md.verbose=verbose(0); md.verbose.qmu=1;
    157         % }}}
    158         %qmu statistics %{{{
    159         md.qmu.statistics.nfiles_per_directory=2;
    160         md.qmu.statistics.ndirectories=5;
    161        
    162         md.qmu.statistics.method(1).name='Histogram';
    163         md.qmu.statistics.method(1).fields={'Sealevel','BslrIce'};
    164         md.qmu.statistics.method(1).steps=1:10;
    165         md.qmu.statistics.method(1).nbins=20;
    166 
    167         md.qmu.statistics.method(2).name='MeanVariance';
    168         md.qmu.statistics.method(2).fields={'Sealevel','BslrIce'};
    169         md.qmu.statistics.method(2).steps=[1:10];
    170 
    171         md.qmu.statistics.method(3).name='SampleSeries';
    172         md.qmu.statistics.method(3).fields={'Sealevel','BslrIce'};
    173         md.qmu.statistics.method(3).steps=[1:10];
    174         md.qmu.statistics.method(3).indices=locations;
    175         %}}}
     139
     140%algorithm:  % {{{
     141md.qmu.method     =dakota_method('nond_samp');
     142md.qmu.method(end)=dmeth_params_set(md.qmu.method(end),...
     143'seed',1234,...
     144'samples',10,...
     145'sample_type','random');
     146md.qmu.output=1;
     147%}}}
     148%parameters % {{{
     149md.qmu.params.direct=true;
     150md.qmu.params.interval_type='forward';
     151md.qmu.params.analysis_driver='matlab';
     152md.qmu.params.evaluation_scheduling='master';
     153md.qmu.params.processors_per_evaluation=2;
     154md.qmu.params.tabular_graphics_data=true;
     155md.qmu.isdakota=1;
     156md.verbose=verbose(0); md.verbose.qmu=1;
     157% }}}
     158%qmu statistics %{{{
     159md.qmu.statistics.nfiles_per_directory=2;
     160md.qmu.statistics.ndirectories=5;
     161
     162md.qmu.statistics.method(1).name='Histogram';
     163md.qmu.statistics.method(1).fields={'Sealevel','BslrIce'};
     164md.qmu.statistics.method(1).steps=[1:10];
     165md.qmu.statistics.method(1).nbins=20;
     166
     167md.qmu.statistics.method(2).name='MeanVariance';
     168md.qmu.statistics.method(2).fields={'Sealevel','BslrIce'};
     169md.qmu.statistics.method(2).steps=[1:10];
     170
     171md.qmu.statistics.method(3).name='SampleSeries';
     172md.qmu.statistics.method(3).fields={'Sealevel','BslrIce'};
     173md.qmu.statistics.method(3).steps=[1:10];
     174md.qmu.statistics.method(3).indices=locations;
     175%}}}
    176176
    177177%run transient dakota solution:
    … …  
    182182md=solve(md,'Transient');
    183183
    184 %compare statistics with our own here:
     184disp(mds.results.StatisticsSolution)
     185
     186%compare statistics with our own here:
    185187svalues=mds.results.StatisticsSolution(end).SealevelSamples; %all values at locations.
    186188
    … …  
    190192end
    191193
     194disp(dvalues)
     195disp(size(dvalues))
     196
    192197samplesnorm=norm(dvalues-svalues,'fro');
    193198
Note: See TracChangeset for help on using the changeset viewer.