Changeset 25125


Ignore:
Timestamp:
06/23/20 15:42:10 (5 years ago)
Author:
jdquinn
Message:

CHG: Saving progress on MATLAB -> Python translations; cleanup

Location:
issm/trunk-jpl
Files:
6 added
44 edited

Legend:

Unmodified
Added
Removed

  • issm/trunk-jpl/etc/environment.sh

    r24877 r25125  
    372372fi
    373373
     374NETCDF_PYTHON_DIR="${ISSM_DIR}/externalpackages/netcdf-python/install"
     375if [ -d "${NETCDF_PYTHON_DIR}" ]; then
     376        if [ -d "${NETCDF_PYTHON_DIR}/lib/python2.7/site-packages" ]; then
     377                ld_library_path_append "${NETCDF_PYTHON_DIR}/lib/python2.7/site-packages"
     378        fi
     379fi
     380
    374381HDF5_DIR="${ISSM_DIR}/externalpackages/hdf5/install"
    375382if [ -d "${HDF5_DIR}" ]; then

  • issm/trunk-jpl/externalpackages/netcdf-python/install.sh

    r23435 r25125  
    22set -eu
    33
    4 #Some cleanup
    5 rm -rf install netCDF4-1.0
    6 mkdir install
    74
    8 #Download from ISSM server
    9 $ISSM_DIR/scripts/DownloadExternalPackage.sh "https://issm.ess.uci.edu/files/externalpackages/netCDF4-1.0.tar.gz" "netCDF4-1.0.tar.gz"
     5# Constants
     6#
     7VER="1.5.3"
    108
    11 #Untar
    12 tar -zxvf  netCDF4-1.0.tar.gz
     9# Environment
     10#
     11export HDF5_DIR="${ISSM_DIR}/externalpackages/petsc/install"
     12export MPIINC_DIR="${ISSM_DIR}/externalpackages/petsc/install/include"
     13export NETCDF_DIR="${ISSM_DIR}/externalpackages/netcdf/install"
     14export PYTHONUSERBASE="${ISSM_DIR}/externalpackages/netcdf-python/install" # This variable and '--user' option supplied to 'setup.py install' are required to install to custom location (source: https://setuptools.readthedocs.io/en/latest/easy_install.html#custom-installation-locations)
    1315
    14 #for later:
    15 rm -rf ISSMDIR
    16 echo $ISSM_DIR | sed 's/\//\\\//g' > ISSMDIR
    17 ISSMDIR=`cat ISSMDIR` && rm -rf ISSMDIR
     16# Download source
     17$ISSM_DIR/scripts/DownloadExternalPackage.sh "https://issm.ess.uci.edu/files/externalpackages/netcdf4-python-${VER}rel.tar.gz" "netcdf4-python-${VER}rel.tar.gz"
    1818
    19 #Move netCDF4 to install directory
    20 rm -rf install/*
    21 mv netCDF4-1.0/* install/
    22 rm -rf netCDF4-1.0
     19# Unpack source
     20tar -zxvf netcdf4-python-${VER}rel.tar.gz
    2321
    24 #Configur and compile
    25 cd install
    26 #edit setup.cfg to point to hdf5 and netcdf
    27 cat setup.cfg.template  | sed "s/\#netCDF4_dir = \/usr\/local/netCDF4_dir = $ISSMDIR\/externalpackages\/netcdf\/install/g"  | sed "s/\#HDF5_dir = \/usr\/local/HDF5_DIR = $ISSMDIR\/externalpackages\/hdf5\/install/g" > setup.cfg
     22# Cleanup
     23rm -rf install src
     24mkdir install src
    2825
     26# Move source to 'src' directory
     27mv netcdf4-python-${VER}rel/* src/
     28rm -rf netcdf4-python-${VER}rel
     29
     30# Compile and install
     31cd src
    2932python setup.py build
    30 python setup.py install
     33python setup.py install --user
     34
     35# Unzip eggs
     36for egg in $(find "${PYTHONUSERBASE}/lib" -name *.egg); do
     37        parent_dir=$(dirname ${egg})
     38        filename=$(basename -- ${egg})
     39        extension="${filename##*.}"
     40        filename="${filename%.*}"
     41        src_dir="${parent_dir}/${filename}"
     42        unzip ${egg} ${src_dir}
     43done

  • issm/trunk-jpl/externalpackages/netcdf/install-4.7-parallel-static-with_tests.sh

    r25073 r25125  
    5353        --enable-fast-install \
    5454        --disable-doxygen \
     55        --enable-netcdf4 \
    5556        --disable-filter-testing \
    5657        --disable-dap-remote-tests \

  • issm/trunk-jpl/externalpackages/netcdf/install-4.7-parallel-static.sh

    r25073 r25125  
    5353        --enable-fast-install \
    5454        --disable-doxygen \
     55        --enable-netcdf4 \
    5556        --disable-testsets \
    5657        --disable-examples

  • issm/trunk-jpl/externalpackages/netcdf/install-4.7-parallel-with_tests.sh

    r25073 r25125  
    5050        --enable-fast-install \
    5151        --disable-doxygen \
     52        --enable-netcdf4 \
    5253        --disable-filter-testing \
    5354        --disable-dap-remote-tests \

  • issm/trunk-jpl/externalpackages/netcdf/install-4.7-parallel.sh

    r25073 r25125  
    5151        --disable-doxygen \
    5252        --disable-testsets \
    53         --disable-examples
     53        --disable-examples \
     54        --enable-netcdf4
    5455
    5556# Compile and install

  • issm/trunk-jpl/externalpackages/netcdf/install-4.7.sh

    r25073 r25125  
    5050        --disable-doxygen \
    5151        --disable-testsets \
    52         --disable-examples
     52        --disable-examples \
     53        --enable-netcdf4
    5354
    5455# Compile and install

  • issm/trunk-jpl/jenkins/ross-debian_linux-solid_earth

    r25095 r25125  
    3131
    3232EXTERNALPACKAGES="
    33         autotools       install-debian-linux.sh
    34         cmake           install.sh
    35         petsc           install-3.12-linux-solid_earth.sh
    36         triangle        install-linux.sh
    37         chaco           install.sh
    38         m1qn3           install.sh
    39         semic           install.sh
    40         boost           install-1.7-linux.sh
    41         curl            install-7.67.sh
    42         netcdf          install-4.7.sh
    43         proj            install-6.2.sh
    44         gdal            install-3.0-python-netcdf.sh
    45         gshhg           install.sh
    46         gmt                     install-6.0-linux.sh
    47         gmsh            install-4.sh
    48         shell2junit     install.sh
     33        autotools               install-debian-linux.sh
     34        cmake                   install.sh
     35        petsc                   install-3.12-linux-solid_earth.sh
     36        triangle                install-linux.sh
     37        chaco                   install.sh
     38        m1qn3                   install.sh
     39        semic                   install.sh
     40        boost                   install-1.7-linux.sh
     41        curl                    install-7.67.sh
     42        netcdf                  install-4.7-parallel.sh
     43        netcdf-python   install.sh
     44        proj                    install-6.2.sh
     45        gdal                    install-3.0-python-netcdf.sh
     46        gshhg                   install.sh
     47        gmt                             install-6.0-linux.sh
     48        gmsh                    install-4.sh
     49        shell2junit             install.sh
    4950"
    5051
    … …  
    5455
    5556# Test suites
    56 MATLAB_TEST=1
     57MATLAB_TEST=0
    5758PYTHON_TEST=1
    5859JAVASCRIPT_TEST=0

  • issm/trunk-jpl/scripts/devpath.py

    r24593 r25125  
    22
    33# NOTE: This script is a stripped-down version of
    4 #       $ISSM_DIR/src/m/dev/devpath.py and is intended only for loading ISSM in
    5 #       order to test our distributable packages. It assumes the following is
    6 #       set before $ISSM_DIR/test/NightlyRun/runme.py is called,
     4#       $ISSM_DIR/src/m/dev/devpath.py and is intended only for loading ISSM in 
     5#       order to test our distributable packages. It assumes the following is 
     6#       set before $ISSM_DIR/test/NightlyRun/runme.py is called, 
    77#
    88#           export ISSM_DIR=</path/to/ISSM>
    9 #                       export PATH="${PATH}:${ISSM_DIR}/bin"
     9#           export PATH="${PATH}:${ISSM_DIR}/bin"
    1010#           export PYTHONPATH="${ISSM_DIR}/scripts"
    1111#           export PYTHONSTARTUP="${PYTHONPATH}/devpath.py"

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

    r25065 r25125  
    11import math
     2import numpy as np
    23
    34from boundary import boundary
    45from epsg2proj import epsg2proj
     6from gdaltransform import gdaltransform
    57from helpers import fileparts
    68from laea import laea
    … …  
    7274
    7375    def isnameany(self, *args): #{{{
     76        boolean = 0
    7477        for arg in args:
    75             if arg == self.name:
    76                 return 1
    77         return 0
     78            if type(arg) in [np.ndarray, list]:
     79                for name in arg:
     80                    if name == self.name:
     81                        boolean = 1
     82                        break
     83            elif arg == self.name:
     84                boolean = 1
     85                break
     86        return boolean
    7887    #}}}
    7988
    8089    def iscontinentany(self, *args): #{{{
     90        boolean = 0
    8191        for arg in args:
    82             if arg == self.continent:
    83                 return 1
    84         return 0
     92            if type(arg) in [np.ndarray, list]:
     93                for continent in arg:
     94                    if continent == self.continent:
     95                        boolean = 1
     96                        break
     97            elif arg == self.continent:
     98                boolean = 1
     99                break
     100        return boolean
    85101    #}}}
    86102
    … …  
    130146            y.append(y[0])
    131147
    132         setattr(out, 'x', x)
    133         setattr(out, 'y', y)
    134         setattr(out, 'nods', len(x))
     148        out = {}
     149        out['x'] = x
     150        out['y'] = y
     151        out['nods'] = len(x)
    135152
    136153        return out
    … …  
    198215        ba = self.contour()
    199216        flags = np.zeros(len(contours))
    200         for i in range(len(contours))
     217        for i in range(len(contours)):
    201218            h = contours[i]
    202219            isin = inpolygon(h.x, h.y, ba.x, ba.y)

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

    r25068 r25125  
    7474                        end
    7575                        output=shpread([self.shppath '/' name '.' ext]);
    76 
     76                       
    7777                        %do we reverse?
    7878                        if strcmpi(self.orientation,'reverse'),

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

    r25065 r25125  
    44
    55from epsg2proj import epsg2proj
     6from helpers import *
    67from pairoptions import pairoptions
    78from shpread import shpread
    … …  
    7576    def edges(self): #{{{
    7677        #read domain
    77         path, name, ext = fileparts(shp.shpfilename)
     78        path, name, ext = fileparts(self.shpfilename)
    7879        if ext != 'shp':
    7980            ext = 'shp'
    … …  
    8283        #do we reverse?
    8384        if self.orientation == 'reverse':
    84             output.x = np.flipud(output.x)
    85             output.y = np.flipud(output.y)
     85            for i in range(len(output)):
     86                output[i].x = np.flipud(output[i].x)
     87                output[i].y = np.flipud(output[i].y)
    8688    #}}}
    8789
    … …  
    104106
    105107        #read domain
    106         path, name, ext = fileparts(shp.shpfilename)
     108        path, name, ext = fileparts(self.shpfilename)
    107109        if ext != 'shp':
    108110            ext = 'shp'
    … …  
    167169
    168170        #read domain
    169         path, name, ext = fileparts(shp.shpfilename)
     171        path, name, ext = fileparts(self.shpfilename)
    170172        if ext != 'shp':
    171173            ext = 'shp'

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

    r24261 r25125  
    3434        string = "%s\n%s" % (string, fielddisplay(self, 'sh_nmin', 'minimum spherical harmonic degree (default 1)'))
    3535        string = "%s\n%s" % (string, fielddisplay(self, 'g0', 'adimensioning constant for gravity (default 10) [m / s^2]'))
    36         string = "%s\n%s" % (string, fielddisplay(self, 'r0', 'adimensioning constant for radius (default 6378 * 1.0e3) [m]'))
     36        string = "%s\n%s" % (string, fielddisplay(self, 'r0', 'adimensioning constant for radius (default 6378e3) [m]'))
    3737        string = "%s\n%s" % (string, fielddisplay(self, 'mu0', 'adimensioning constant for stress (default 1.0e11) [Pa]'))
    3838        string = "%s\n%s" % (string, fielddisplay(self, 'allow_layer_deletion', 'allow for migration of the integration boundary with increasing spherical harmonics degree (default 1)'))
    … …  
    6666        self.sh_nmin = 1
    6767        self.g0 = 10  # m / s^2
    68         self.r0 = 6378 * 1e3  #m
     68        self.r0 = 6378e3  #m
    6969        self.mu0 = 1e11  # Pa
    7070        self.allow_layer_deletion = 1

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

    r24546 r25125  
     1import numpy as np
     2
     3from checkfield import checkfield
    14from fielddisplay import fielddisplay
    2 from checkfield import checkfield
    35from WriteData import WriteData
    4 import numpy as np
    56
    67

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

    r25118 r25125  
    3131                        %tidal love numbers:
    3232                        self.th=0.6149; %degree 2
    33                         self.tk=0.3055; % degree 2
     33                        self.tk=0.3055; %degree 2
    3434
    3535                        %secular fluid love number:
    … …  
    7575                end % }}}
    7676                function marshall(self,prefix,md,fid) % {{{
    77                        
     77
    7878                        WriteData(fid,prefix,'name','md.solidearth.love.model','data',1,'format','Integer');
    7979
    … …  
    8585                        WriteData(fid,prefix,'object',self,'fieldname','tk','name','md.solidearth.love.tk','format','Double');
    8686                        WriteData(fid,prefix,'object',self,'fieldname','tk2secular','name','md.solidearth.love.tk2secular','format','Double');
    87 
    8887
    8988                end % }}}

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

    r24213 r25125  
    8686        self.heatcapacity = 2093.
    8787        #ice latent heat of fusion L (J / kg)
    88         self.latentheat = 3.34 * 1.0e5
     88        self.latentheat = 3.34e5
    8989        #ice thermal conductivity (W / m / K)
    9090        self.thermalconductivity = 2.4
    … …  
    9696        self.meltingpoint = 273.15
    9797        #rate of change of melting point with pressure (K / Pa)
    98         self.beta = 9.8 * 1.0e-8
     98        self.beta = 9.8e-8
    9999        #mixed layer (ice-water interface) heat capacity (J / kg / K)
    100100        self.mixed_layer_capacity = 3974.
    101101        #thermal exchange velocity (ice-water interface) (m / s)
    102         self.thermal_exchange_velocity = 1.00 * 1.0e-4
     102        self.thermal_exchange_velocity = 1.00e-4
    103103        #Rheology law: what is the temperature dependence of B with T
    104104        #available: none, paterson and arrhenius
    … …  
    106106
    107107        # GIA:
    108         self.lithosphere_shear_modulus = 6.7 * 1.0e10  # (Pa)
     108        self.lithosphere_shear_modulus = 6.7e10  # (Pa)
    109109        self.lithosphere_density = 3.32  # (g / cm^ - 3)
    110         self.mantle_shear_modulus = 1.45 * 1.0e11  # (Pa)
     110        self.mantle_shear_modulus = 1.45e11  # (Pa)
    111111        self.mantle_density = 3.34  # (g / cm^ - 3)
    112112

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

    r24756 r25125  
    8080                self.heatcapacity = 2093.
    8181                #ice latent heat of fusion L (J / kg)
    82                 self.latentheat = 3.34 * 1.0e5
     82                self.latentheat = 3.34e5
    8383                #ice thermal conductivity (W / m / K)
    8484                self.thermalconductivity = 2.4
    … …  
    8888                self.meltingpoint = 273.15
    8989                #rate of change of melting point with pressure (K / Pa)
    90                 self.beta = 9.8 * 1.0e-8
     90                self.beta = 9.8e-8
    9191                #mixed layer (ice-water interface) heat capacity (J / kg / K)
    9292                self.mixed_layer_capacity = 3974.
    9393                #thermal exchange velocity (ice-water interface) (m / s)
    94                 self.thermal_exchange_velocity = 1.00 * 1.0e-4
     94                self.thermal_exchange_velocity = 1.00e-4
    9595                #Rheology law: what is the temperature dependence of B with T
    9696                #available: none, paterson and arrhenius
    … …  
    101101                #center of the earth (approximation, must not be 0), then the lab (lithosphere / asthenosphere boundary) then the surface
    102102                #(with 1d3 to avoid numerical singularities)
    103                 self.radius = [1e3, 6278 * 1e3, 6378 * 1e3]
     103                self.radius = [1e3, 6278e3, 6378e3]
    104104                self.viscosity = [1e21, 1e40]  #mantle and lithosphere viscosity (respectively) [Pa.s]
    105                 self.lame_mu = [1.45 * 1e11, 6.7 * 1e10]  # (Pa)  #lithosphere and mantle shear modulus (respectively) [Pa]
     105                self.lame_mu = [1.45e11, 6.7e10]  # (Pa)  #lithosphere and mantle shear modulus (respectively) [Pa]
    106106                self.lame_lambda = self.lame_mu  # (Pa)  #mantle and lithosphere lamba parameter (respectively) [Pa]
    107107                self.burgers_viscosity = [np.nan, np.nan]
    108108                self.burgers_mu = [np.nan, np.nan]
    109109                self.isburgers = [0, 0]
    110                 self.density = [5.51 * 1e3, 5.50 * 1e3]  # (Pa)  #mantle and lithosphere density [kg / m^3]
     110                self.density = [5.51e3, 5.50e3]  # (Pa)  #mantle and lithosphere density [kg / m^3]
    111111                self.issolid = [1, 1]  # is layer solid or liquid.
    112112            elif nat == 'hydro':

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

    r24213 r25125  
    9191        self.heatcapacity = 2093.
    9292        #ice latent heat of fusion L (J / kg)
    93         self.latentheat = 3.34 * 1.0e5
     93        self.latentheat = 3.34e5
    9494        #ice thermal conductivity (W / m / K)
    9595        self.thermalconductivity = 2.4
    … …  
    101101        self.meltingpoint = 273.15
    102102        #rate of change of melting point with pressure (K / Pa)
    103         self.beta = 9.8 * 1.0e-8
     103        self.beta = 9.8e-8
    104104        #mixed layer (ice-water interface) heat capacity (J / kg / K)
    105105        self.mixed_layer_capacity = 3974.
    106106        #thermal exchange velocity (ice-water interface) (m / s)
    107         self.thermal_exchange_velocity = 1.00 * 1.0e-4
     107        self.thermal_exchange_velocity = 1.00e-4
    108108        #Rheology law: what is the temperature dependence of B with T
    109109        #available: none, paterson and arrhenius
    110110        self.rheology_law = 'Paterson'
    111111    # GIA:
    112         self.lithosphere_shear_modulus = 6.7 * 1.0e10  # (Pa)
     112        self.lithosphere_shear_modulus = 6.7e10  # (Pa)
    113113        self.lithosphere_density = 3.32  # (g / cm^ - 3)
    114         self.mantle_shear_modulus = 1.45 * 1.0e11  # (Pa)
     114        self.mantle_shear_modulus = 1.45e11  # (Pa)
    115115        self.mantle_density = 3.34  # (g / cm^ - 3)
    116116    #SLR

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

    r24213 r25125  
    8787        self.heatcapacity = 2093.
    8888        #ice latent heat of fusion L (J / kg)
    89         self.latentheat = 3.34 * 1.0e5
     89        self.latentheat = 3.34e5
    9090        #ice thermal conductivity (W / m / K)
    9191        self.thermalconductivity = 2.4
    … …  
    9797        self.meltingpoint = 273.15
    9898        #rate of change of melting point with pressure (K / Pa)
    99         self.beta = 9.8 * 1.0e-8
     99        self.beta = 9.8e-8
    100100        #mixed layer (ice-water interface) heat capacity (J / kg / K)
    101101        self.mixed_layer_capacity = 3974.
    102102        #thermal exchange velocity (ice-water interface) (m / s)
    103         self.thermal_exchange_velocity = 1.00 * 1.0e-4
     103        self.thermal_exchange_velocity = 1.00e-4
    104104        #Rheology law: what is the temperature dependence of B with T
    105105        #available: none, paterson and arrhenius
    … …  
    107107
    108108        # GIA:
    109         self.lithosphere_shear_modulus = 6.7 * 1.0e10  # (Pa)
     109        self.lithosphere_shear_modulus = 6.7e10  # (Pa)
    110110        self.lithosphere_density = 3.32  # (g / cm^ - 3)
    111         self.mantle_shear_modulus = 1.45 * 1.0e11  # (Pa)
     111        self.mantle_shear_modulus = 1.45e11  # (Pa)
    112112        self.mantle_density = 3.34  # (g / cm^ - 3)
    113113

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

    r25094 r25125  
    2929from initialization import initialization
    3030from rifts import rifts
    31 from slr import slr
     31from solidearth import solidearth
    3232from dsl import dsl
    3333from debug import debug
    … …  
    9898        self.rifts = rifts()
    9999        self.dsl = dsl()
    100         self.slr = slr()
     100        self.solidearth = solidearth()
    101101
    102102        self.debug = debug()
    … …  
    147147                'initialization',
    148148                'rifts',
    149                 'slr',
     149                'solidearth',
    150150                'dsl',
    151151                'debug',
    … …  
    194194        string = "%s\n%s" % (string, "%19s: % - 22s - -  %s" % ("initialization", "[%s, %s]" % ("1x1", obj.initialization.__class__.__name__), "initial guess / state"))
    195195        string = "%s\n%s" % (string, "%19s: % - 22s - -  %s" % ("rifts", "[%s, %s]" % ("1x1", obj.rifts.__class__.__name__), "rifts properties"))
    196         string = "%s\n%s" % (string, "%19s: % - 22s - -  %s" % ("slr", "[%s, %s]" % ("1x1", obj.slr.__class__.__name__), "slr forcings"))
     196        string = "%s\n%s" % (string, "%19s: % - 22s - -  %s" % ("solidearth", "[%s, %s]" % ("1x1", obj.solidearth.__class__.__name__), "solid earth inputs and settings"))
    197197        string = "%s\n%s" % (string, "%19s: % - 22s - -  %s" % ("dsl", "[%s, %s]" % ("1x1", obj.dsl.__class__.__name__), "dynamic sea level"))
    198198        string = "%s\n%s" % (string, "%19s: % - 22s - -  %s" % ("debug", "[%s, %s]" % ("1x1", obj.debug.__class__.__name__), "debugging tools (valgrind, gprof)"))

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

    r25012 r25125  
    4343        #check length of input
    4444        if len(arg) % 2:
    45             raise TypeError('Invalid parameter / value pair arguments')
     45            raise TypeError('Invalid parameter/value pair arguments')
    4646        numoptions = int(len(arg) / 2)
    4747

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

    r25118 r25125  
    2222               
    2323                %moment of inertia:
    24                 self.equatorialmoi=8.0077*10^37; % [kg m^2]
    25                 self.polarmoi            =8.0345*10^37; % [kg m^2]
     24                self.equatorialmoi      =8.0077*10^37; % [kg m^2]
     25                self.polarmoi           =8.0345*10^37; % [kg m^2]
    2626
    2727                % mean rotational velocity of earth

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

    r25065 r25125  
    220220                        if iscell(continent),
    221221                                if length(continent)==1,
    222                                          if strcmpi(continent{1},'all'),
    223                                                  %need to transform this into a list of continents:
    224                                                  continent={};
    225                                                  for i=1:length(self.basins),
    226                                                          continent{end+1}=self.basins{i}.continent;
    227                                                  end
    228                                                  continent=unique(continent);
    229                                          end
     222                                        if strcmpi(continent{1},'all'),
     223                                                %need to transform this into a list of continents:
     224                                                continent={};
     225                                                for i=1:length(self.basins),
     226                                                        continent{end+1}=self.basins{i}.continent;
     227                                                end
     228                                                continent=unique(continent);
     229                                        end
    230230                                else
    231231                                        %nothing to do, we have a list of continents
    … …  
    236236                                        continent={};
    237237                                        for i=1:length(self.basins),
    238                                                  continent{end+1}=self.basins{i}.continent;
     238                                                continent{end+1}=self.basins{i}.continent;
    239239                                        end
    240240                                        continent=unique(continent);
    … …  
    247247                        if iscell(bas),
    248248                                if length(bas)==1,
    249                                          if strcmpi(bas{1},'all'),
    250                                                  %need to transform this into a list of basins:
    251                                                  baslist=[];
    252                                                  for i=1:length(self.basins),
    253                                                          if self.basins{i}.iscontinentany(continent{:}),
    254                                                                  baslist(end+1)=i;
    255                                                          end
    256                                                  end
    257                                                  baslist=unique(baslist);
    258                                          else
    259                                                  bas=bas{1};
    260                                                  baslist=[];
    261                                                  for i=1:length(self.basins),
    262                                                          if self.basins{i}.iscontinentany(continent{:}),
    263                                                                  if self.basins{i}.isnameany(bas),
    264                                                                          baslist(end+1)=i;
    265                                                                  end
    266                                                          end
    267                                                  end
    268 
    269                                          end
     249                                        if strcmpi(bas{1},'all'),
     250                                                %need to transform this into a list of basins:
     251                                                baslist=[];
     252                                                for i=1:length(self.basins),
     253                                                        if self.basins{i}.iscontinentany(continent{:}),
     254                                                                baslist(end+1)=i;
     255                                                        end
     256                                                end
     257                                                baslist=unique(baslist);
     258                                        else
     259                                        bas=bas{1};
     260                                        baslist=[];
     261                                        for i=1:length(self.basins),
     262                                                if self.basins{i}.iscontinentany(continent{:}),
     263                                                        if self.basins{i}.isnameany(bas),
     264                                                                baslist(end+1)=i;
     265                                                        end
     266                                                end
     267                                        end
     268
     269                                        end
    270270                                else
    271271                                        %we have a list of basin names:

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

    r25065 r25125  
    4343
    4444        # Create a default object
    45         self.setdefaultparameters()
     45        self = self.setdefaultparameters()
    4646
    4747        if len(args):
    … …  
    7777        self.eltransitions  = []
    7878        self.planet         = 'earth'
    79 
    80         return self
    8179    #}}}
    8280
    … …  
    191189            self.eltransitions.append(meshintersect3d(xe, ye, ze, xei, yei, zei, 'force', force))
    192190    #}}}
     191
     192    def checkintersections(self): #{{{
     193        flags = np.zeros(self.earth.mesh.numberofvertices, 1)
     194        for i in range(len(self.basins)):
     195            flags[self.transitions[i]] = i
     196        plotmodel(self.earth, 'data', flags, 'coastline', 'on')
     197    #}}}
     198
     199    def checkbasinconsistency(self): #{{{
     200        for i in range(len(self.basins)):
     201            self.basins[i].checkconsistency()
     202    #}}}
     203
     204    def basinindx(self, *args): #{{{
     205        options = pairoptions(*args)
     206        continent = options.getfieldvalue('continent', 'all')
     207        bas = options.getfieldvalue('basin', 'all')
     208
     209        #expand continent list #{{{
     210        if type(continent) == np.ndarray:
     211            if continent.shape[1] == 1:
     212                if continent[0] == 'all':
     213                    #need to transform this into a list of continents
     214                    continent = []
     215                    for i in range(len(self.basins)):
     216                        continent.append(self.basins[i].continent)
     217                    continent = np.unique(continent)
     218            else:
     219                pass #nothing to do, we have a list of continents
     220        else:
     221            if continent == 'all':
     222                #need to transform this into a list of continents
     223                continent = []
     224                for i in range(len(self.basins)):
     225                    continent.append(self.basins[i].continent)
     226                continent = np.unique(continent)
     227            else:
     228                continent = [continent]
     229        #}}}
     230        #expand basins list using the continent list above and the extra bas discriminator #{{{
     231        if type(bas) == np.ndarray:
     232            if bas.shape[1] == 1:
     233                if bas[0] == 'all':
     234                    #need to transform this into a list of basins
     235                    baslist = []
     236                    for i in range(len(self.basins)):
     237                        if self.basins[i].iscontinentany(continent):
     238                            baslist.append(i)
     239                    baslist = np.unique(baslist)
     240                else:
     241                    bas = bas[0]
     242                    baslist = []
     243                    for i in range(len(self.basins)):
     244                        if self.basins[i].iscontinentany(continent):
     245                            if self.basins[i].isnameany(bas):
     246                                baslist.append(i)
     247            else:
     248                #we have a list of basin names
     249                baslist = []
     250                for i in range(len(bas)):
     251                    basname = bas[i]
     252                    for j in range(len(self.basins)):
     253                        if self.basins[j].iscontinentany(continent):
     254                            if self.basins[j].isnameany(basname):
     255                                baslist.append(j)
     256                    baslist = np.unique(baslist)
     257        else:
     258            if bas == 'all':
     259                baslist = []
     260                for i in range(len(self.basins)):
     261                    if self.basins[i].iscontinentany(continent):
     262                        baslist.append(i)
     263                baslist = np.unique(baslist)
     264            else:
     265                baslist = []
     266                for i in range(len(self.basins)):
     267                    if self.basins[i].iscontinentany(continent):
     268                        if self.basins[i].isnameany(bas):
     269                            baslist.append(i)
     270                baslist = np.unique(baslist)
     271
     272        return baslist
     273        #}}}
     274    #}}}

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

    r25014 r25125  
    104104        self.fluid_love = 0.942
    105105        #moment of inertia:
    106         self.equatorial_moi = 8.0077 * 1.0e37  # [kg m^2]
    107         self.polar_moi = 8.0345 * 1.0e37  # [kg m^2]
     106        self.equatorial_moi = 8.0077e37  # [kg m^2]
     107        self.polar_moi = 8.0345e37  # [kg m^2]
    108108        #mean rotational velocity of earth
    109         self.angular_velocity = 7.2921 * 1.0e-5  # [s^ - 1]
     109        self.angular_velocity = 7.2921e-5  # [s^ - 1]
    110110        #numerical discretization accuracy
    111111        self.degacc = 0.01

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

    r25118 r25125  
    4646                self.degacc=.01;
    4747
    48                 %how many time steps we skip before we run SOLIDEARTHSETTINGS solver during transient
     48                %how many time steps we skip before we run solidearthsettings solver during transient
    4949                self.runfrequency=1;
    5050       
    … …  
    8686                        fielddisplay(self,'ocean_area_scaling','correction for model representation of ocean area [default: No correction]');
    8787                        fielddisplay(self,'computesealevelchange','compute sealevel change from GRD in addition to steric?) default 0');
    88                         fielddisplay(self,'runfrequency','how many time steps we skip before we run SOLIDEARTHSETTINGS solver during transient (default: 1)');
     88                        fielddisplay(self,'runfrequency','how many time steps we skip before we run solidearthsettings solver during transient (default: 1)');
    8989                        fielddisplay(self,'rigid','rigid earth graviational potential perturbation');
    9090                        fielddisplay(self,'elastic','elastic earth graviational potential perturbation');

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

    r25118 r25125  
    6060                function marshall(self,prefix,md,fid) % {{{
    6161
    62                         if isempty(self.icethicknesschange), self.icethicknesschange=zeros(md.mesh.numberofelements+1,1); end
    63                         if isempty(self.waterheightchange), self.waterheightchange=zeros(md.mesh.numberofelements+1,1); end
    64                         if isempty(self.other), self.other=zeros(md.mesh.numberofelements+1,1); end
     62                        if isempty(self.icethicknesschange),
     63                                self.icethicknesschange=zeros(md.mesh.numberofelements+1,1);
     64                        end
     65                        if isempty(self.waterheightchange),
     66                                self.waterheightchange=zeros(md.mesh.numberofelements+1,1);
     67                        end
     68                        if isempty(self.other),
     69                                self.other=zeros(md.mesh.numberofelements+1,1);
     70                        end
    6571                        WriteData(fid,prefix,'object',self,'fieldname','icethicknesschange','name','md.solidearth.surfaceload.icethicknesschange',...
    6672                                'format','DoubleMat','mattype',2,'timeserieslength',md.mesh.numberofelements+1,'yts',md.constants.yts);

  • issm/trunk-jpl/src/m/coordsystems/gdaltransform.py

    r25065 r25125  
    4343    filename_out = file_out.name
    4444
    45     file_in.write(b'%8g %8g\n' % (x.flatten(1) y.flatten(1)))
     45    file_in.write(b'%8g %8g\n' % (x.flatten(1), y.flatten(1)))
    4646    file_in.close()
    4747

  • issm/trunk-jpl/src/m/dev/devpath.py

    r24254 r25125  
    1111    raise NameError('"ISSM_DIR" environment variable is empty! You should define ISSM_DIR in your .cshrc or .bashrc!')
    1212
    13     #Go through src/m and append any directory that contains a *.py file to PATH
     13#Get paths to locally installed Python libraries
     14tpls = [
     15    'cftime', # Needed for netCDF4
     16    'netCDF4'
     17]
     18for tpl in os.listdir(ISSM_DIR + '/externalpackages'):
     19    tpl_path = ISSM_DIR + '/externalpackages/' + tpl
     20    for dirpath, dirnames, filenames in os.walk(tpl_path):
     21        install_lib_path = dirpath + '/install/lib'
     22        for dirpath, dirnames, filenames in os.walk(install_lib_path):
     23            for dirname in dirnames:
     24                if dirname in tpls:
     25                    dir = dirpath + '/' + dirname
     26                    for file in os.listdir(dir):
     27                        if file.find(".py") != -1:
     28                            if file.find(".pyc") == -1:
     29                                if dirpath not in sys.path:
     30                                    sys.path.append(dirpath)
     31
     32#Go through src/m and append any directory that contains a *.py file to PATH
    1433for root, dirs, files in os.walk(ISSM_DIR + '/src/m'):
    1534    if '.svn' in dirs:
    … …  
    2140                    sys.path.append(root)
    2241
    23     #Also add the Nightly run directory
     42#Also add the Nightly run directory
    2443if ISSM_DIR + '/test/NightlyRun' not in sys.path:
    2544    sys.path.append(ISSM_DIR + '/test/NightlyRun')
    … …  
    2847if ISSM_DIR + '/src/wrappers/python/.libs' not in sys.path:
    2948    sys.path.append(ISSM_DIR + '/src/wrappers/python/.libs')
    30     # If using clusters, we need to have the path to the cluster settings directory
     49
     50# If using clusters, we need to have the path to the cluster settings directory
    3151if JPL_SVN is not None:
    3252    jpl_path = JPL_SVN + '/usr/' + USERNAME

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

    r25065 r25125  
    44
    55
    6 def polyarea: #{{{
     6def polyarea(x, y): #{{{
    77    '''
    88    POLYAREA - returns the area of the 2-D polygon defined by the vertices in

  • issm/trunk-jpl/src/m/io/loadvars.py

    r25092 r25125  
     1from collections import OrderedDict
     2from netCDF4 import Dataset
     3import numpy as np
     4from re import findall
    15import shelve
    2 import numpy as np
    3 from netCDF4 import Dataset
    4 from re import findall
    5 from collections import OrderedDict
     6
    67from model import *
    78#hack to keep python 2 compatibility
    … …  
    1516
    1617def loadvars(*args, OL):
    17     """
    18     LOADVARS - function to load variables to a file.
    19 
    20     This function loads one or more variables from a file.  The names of the variables
    21     must be supplied.  If more than one variable is specified, it may be done with
    22     a list of names or a dictionary of name as keys.  The output type will correspond
    23     to the input type.  All the variables in the file may be loaded by specifying only
    24     the file name.
     18    '''
     19    LOADVARS - function to load variables from a file.
     20
     21    This function loads one or more variables from a file. The names of the
     22    variables must be supplied. If more than one variable is specified, it may
     23    be done with a list of names or a dictionary of name as keys. The output
     24    type will correspond to the input type. All the variables in the file may
     25    be loaded by specifying only the file name.
    2526
    2627    Usage:
    … …  
    2930        nvdict = loadvars('shelve.dat', {'a':None, 'b':None})
    3031        nvdict = loadvars('shelve.dat')
    31 
    32     """
     32    '''
    3333
    3434    filename = ''
    … …  
    4040        if not filename:
    4141            filename = '/tmp/shelve.dat'
    42 
    4342    else:
    4443        raise TypeError("Missing file name.")
    … …  
    4746        for name in args[1:]:
    4847            nvdict[name] = None
    49 
    5048    elif len(args) == 2 and isinstance(args[1], list):  # (filename, [names])
    5149        for name in args[1]:
    5250            nvdict[name] = None
    53 
    5451    elif len(args) == 2 and isinstance(args[1], dict):  # (filename, {names:values})
    5552        nvdict = args[1]
    56 
    5753    elif len(args) == 1:  #  (filename)
    5854        pass
    59 
    6055    else:
    6156        raise TypeError("Unrecognized input arguments.")
    … …  
    7873                print(("Variable '%s' loaded." % name))
    7974        my_shelf.close()
    80 
    8175    else:  #We used netcdf for the save
    8276        try:

  • issm/trunk-jpl/src/m/materials/nye.py

    r24794 r25125  
    4040        n = 8.                  # Glen's exponent
    4141    elif ice_type == 2:         # H2O ice
    42         A_const = 9 * 1.0e4       # s^-1 MPa
     42        A_const = 9e4       # s^-1 MPa
    4343        Q = 60000.              #  J mol^-1
    4444        n = 3.                  # Glen's exponent

  • issm/trunk-jpl/src/m/plot/applyoptions.py

    r25065 r25125  
    55    from matplotlib.ticker import MaxNLocator
    66except ImportError:
    7     print("could not import pylab, matplotlib has not been installed, no plotting capabilities enabled")
     7    print("could not import pyplot, matplotlib has not been installed, no plotting capabilities enabled")
    88import numpy as np
    99

  • issm/trunk-jpl/src/m/plot/plot_unit.py

    r25065 r25125  
    77    from mpl_toolkits.mplot3d.art3d import Poly3DCollection
    88except ImportError:
    9     print("could not import pylab, matplotlib has not been installed, no plotting capabilities enabled")
     9    print("could not import pyplot, matplotlib has not been installed, no plotting capabilities enabled")
    1010import numpy as np
    1111

  • issm/trunk-jpl/src/m/plot/plotmodel.py

    r25065 r25125  
    55    from mpl_toolkits.axes_grid1 import ImageGrid
    66except ImportError:
    7     print("could not import pylab, matplotlib has not been installed, no plotting capabilities enabled")
     7    print("could not import pyplot, matplotlib has not been installed, no plotting capabilities enabled")
    88import numpy as np
    99

  • issm/trunk-jpl/src/m/shp/shpread.m

    r25065 r25125  
    2121%recover options
    2222options=pairoptions(varargin{:});
     23invert=getfieldvalue(options,'invert',0);
    2324
    2425%some checks
    … …  
    9798end
    9899
    99 invert=getfieldvalue(options,'invert',0);
    100100if invert,
    101101        for i=1:length(Struct),

  • issm/trunk-jpl/src/m/shp/shpread.py

    r25065 r25125  
    55    print("could not import shapefile, PyShp has not been installed, no shapefile reading capabilities enabled")
    66
     7from helpers import OrderedStruct
    78from pairoptions import pairoptions
    89
    … …  
    4748    Sources:
    4849    - https://github.com/GeospatialPython/pyshp
     50
     51    TODO:
     52    - Create class that can be used to store and pretty print shape structs
     53      (ala OrderedStruct from src/m/qmu/helpers.py).
    4954    '''
    5055
    … …  
    5964    sf = shapefile.Reader(filename)
    6065
    61     Dicts = []
    62 
    63     #retrieve ID (if it exists)
    64     fields = sf.fields
    65     for i in range(1, len(fields)): # skip over first field, which is "DeletionFlag"
    66         if fields[i][0] == 'id':
    67             name = str(sf.record(i - 1)[0]) # record index is offset by one, again, because of "DeletionFlag"
    68             break
    69 
     66    Structs = []
    7067    shapes = sf.shapes()
    7168    for i in range(len(shapes)):
    72         Dict = {}
     69        Struct = OrderedStruct()
    7370        shape = shapes[i]
    7471        if shape.shapeType == shapefile.POINT:
    75             Dict['x'] = shape.points[0][0]
    76             Dict['y'] = shape.points[0][1]
    77             Dict['density'] = 1
    78             Dict['Geometry'] = 'Point'
     72            Struct.x = shape.points[0][0]
     73            Struct.y = shape.points[0][1]
     74            Struct.density = 1
     75            Struct.Geometry = 'Point'
    7976        elif shape.shapeType == shapefile.POLYLINE:
    8077            num_points = len(shape.points)
    … …  
    8582                x.append(point[0])
    8683                y.append(point[1])
    87             Dict['x'] = x
    88             Dict['y'] = y
    89             Dict['nods'] = num_points
    90             Dict['density'] = 1
    91             Dict['closed'] = 1
    92             Dict['BoundingBox'] = shape.bbox
    93             Dict['Geometry'] = 'Line'
     84            Struct.x = x
     85            Struct.y = y
     86            Struct.nods = num_points
     87            Struct.density = 1
     88            Struct.closed = 1
     89            Struct.BoundingBox = shape.bbox
     90            Struct.Geometry = 'Line'
    9491        elif shape.shapeType == shapefile.POLYGON:
    9592            num_points = len(shape.points)
    … …  
    10097                x.append(point[0])
    10198                y.append(point[1])
    102             Dict['x'] = x
    103             Dict['y'] = y
    104             Dict['nods'] = num_points
    105             Dict['density'] = 1
    106             Dict['closed'] = 1
    107             Dict['BoundingBox'] = shape.bbox
    108             Dict['Geometry'] = 'Polygon'
     99            Struct.x = x
     100            Struct.y = y
     101            Struct.nods = num_points
     102            Struct.density = 1
     103            Struct.closed = 1
     104            Struct.BoundingBox = shape.bbox
     105            Struct.Geometry = 'Polygon'
    109106        else:
    110107            # NOTE: We could do this once before looping over shapes as all
    … …  
    115112            #
    116113            raise Exception('shpread error: geometry {} is not currently supported'.format(shape.shapeTypeName))
     114
    117115        name = ''
    118116        fields = sf.fields
    … …  
    121119            # 'id' field gets special treatment
    122120            if fieldname == 'id':
    123                 name = str(sf.record(j - 1)[0]) # record index is offset by one, again, because of "DeletionFlag"
     121                name = str(sf.record(i)[j - 1]) # record index is offset by one, again, because of "DeletionFlag"
    124122            else:
    125                 Dict[fieldname] = sf.record(j - 1)[0]
    126         Dict['name'] = name
    127         Dicts.append(Dict)
     123                setattr(Struct, str(fieldname), sf.record(i)[j - 1]) # cast to string removes "u" from "u'fieldname'"
     124        Struct.name = name
     125        Structs.append(Struct)
    128126
    129127    invert = options.getfieldvalue('invert', 0)
    130128    if invert:
    131         for i in range(len(Dicts)):
    132             Dicts[i].x = np.flipud(Dicts[i].x)
    133             Dicts[i].y = np.flipud(Dicts[i].y)
     129        for i in range(len(Structs)):
     130            Structs[i].x = np.flipud(Structs[i].x)
     131            Structs[i].y = np.flipud(Structs[i].y)
    134132
    135     return Dicts
     133    return Structs
    136134#}}}

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

    r24214 r25125  
    11#Test Name: SquareSheetConstrainedGia2d
     2from socket import gethostname
     3
    24import numpy as np
     5
    36from model import *
    4 from socket import gethostname
    5 from triangle import *
    6 from setmask import *
    77from parameterize import *
    88from setflowequation import *
     9from setmask import *
    910from solve import *
     11from triangle import *
    1012
    1113
    12 #Define a model
    13 md = model()
    14 md = triangle(md, '../Exp/Square.exp', 100000.)
     14md = triangle(model(), '../Exp/Square.exp', 100000.)
    1515md = setmask(md, '', '')
    1616md = parameterize(md, '../Par/SquareSheetConstrained.py')
    1717
     18#GIA
     19md.gia = giaivins()
     20md.gia.lithosphere_thickness = 100. * np.ones(md.mesh.numberofvertices) # in km
     21md.gia.mantle_viscosity = 1.0e21 * np.ones(md.mesh.numberofvertices) # in Pa.s
     22md.materials.lithosphere_shear_modulus = 6.7e10 # in Pa
     23md.materials.lithosphere_density = 3.32 # in g/cm^3
     24md.materials.mantle_shear_modulus = 1.45e11 # in Pa
     25md.materials.mantle_density = 3.34 # in g/cm^3
     26
    1827#Indicate what you want to compute
    19 md.gia.cross_section_shape = 1  # for square-edged x - section
     28md.gia.cross_section_shape = 1 # for square-edged x-section
    2029
    2130#Define loading history (see test2001.m for the description)
    … …  
    3544field_names = ['GiaW', 'GiadWdt']
    3645field_tolerances = [1e-13, 1e-13]
    37 field_values = [md.results.GiaSolution.GiaW,
    38                 md.results.GiaSolution.GiadWdt]
     46field_values = [md.results.GiaSolution.UGia,
     47                md.results.GiaSolution.UGiaRate]

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

    r24914 r25125  
    1 #Test Name: EarthSlr
     1#Test Name: EarthSolidearth
    22import numpy as np
    33
    … …  
    1515md.mesh = gmshplanet('radius', 6.371012 * 10**3, 'resolution', 700.)  #500 km resolution mesh
    1616
    17 #parameterize slr solution:
    18 #slr loading:
    19 md.slr.deltathickness = np.zeros((md.mesh.numberofelements))
    20 md.slr.sealevel = np.zeros((md.mesh.numberofvertices))
     17#parameterize solidearth solution:
     18#solidearth loading:
     19md.solidearth.deltathickness = np.zeros((md.mesh.numberofelements))
     20md.solidearth.sealevel = np.zeros((md.mesh.numberofvertices))
    2121md.dsl.global_average_thermosteric_sea_level_change=np.zeros((2, ))
    2222md.dsl.sea_surface_height_change_above_geoid=np.zeros((md.mesh.numberofvertices+1, ))
    … …  
    2828longe = np.sum(md.mesh.long[md.mesh.elements - 1], axis=1) / 3
    2929pos = np.where(late < -80)
    30 md.slr.deltathickness[pos] = -100
     30md.solidearth.deltathickness[pos] = -100
    3131#greenland
    3232pos = np.where(np.logical_and.reduce((late > 70, late < 80, longe > -60, longe < -30)))
    33 md.slr.deltathickness[pos] = -100
     33md.solidearth.deltathickness[pos] = -100
    3434
    3535#elastic loading from love numbers:
    3636nlov = 101
    37 md.slr.love_h = love_numbers('h')[:nlov]
    38 md.slr.love_k = love_numbers('k')[:nlov]
    39 md.slr.love_l = love_numbers('l')[:nlov]
     37md.solidearth.love_h = love_numbers('h')[:nlov]
     38md.solidearth.love_k = love_numbers('k')[:nlov]
     39md.solidearth.love_l = love_numbers('l')[:nlov]
    4040
    4141#mask:
    … …  
    5454
    5555#make sure that the elements that have loads are fully grounded:
    56 pos = np.nonzero(md.slr.deltathickness)[0]
     56pos = np.nonzero(md.solidearth.deltathickness)[0]
    5757md.mask.ocean_levelset[md.mesh.elements[pos, :] - 1] = 1
    5858
    … …  
    6161md.mask.ice_levelset = icemask
    6262
    63 md.slr.ocean_area_scaling = 0
     63md.solidearth.ocean_area_scaling = 0
    6464
    6565#geometry for the bed, arbitrary
    … …  
    7373
    7474#New stuff
    75 md.slr.spcthickness = np.nan * np.ones((md.mesh.numberofvertices, ))
    76 md.slr.Ngia = np.zeros((md.mesh.numberofvertices, ))
    77 md.slr.Ugia = np.zeros((md.mesh.numberofvertices, ))
    78 md.slr.hydro_rate = np.zeros((md.mesh.numberofvertices, ))
     75md.solidearth.spcthickness = np.nan * np.ones((md.mesh.numberofvertices, ))
     76md.solidearth.Ngia = np.zeros((md.mesh.numberofvertices, ))
     77md.solidearth.Ugia = np.zeros((md.mesh.numberofvertices, ))
     78md.solidearth.hydro_rate = np.zeros((md.mesh.numberofvertices, ))
    7979
    8080#Miscellaneous
    … …  
    8282
    8383#Solution parameters
    84 md.slr.reltol = np.nan
    85 md.slr.abstol = 1e-3
    86 md.slr.geodetic = 1
     84md.solidearth.reltol = np.nan
     85md.solidearth.abstol = 1e-3
     86md.solidearth.geodetic = 1
    8787
    8888#max number of iteration reverted back to 10 (i.e., the original default value)
    89 md.slr.maxiter = 10
     89md.solidearth.maxiter = 10
    9090
    9191#eustatic run:
    92 md.slr.rigid = 0
    93 md.slr.elastic = 0
    94 md.slr.rotation = 0
     92md.solidearth.rigid = 0
     93md.solidearth.elastic = 0
     94md.solidearth.rotation = 0
    9595md = solve(md, 'Sealevelrise')
    9696Seustatic = md.results.SealevelriseSolution.Sealevel
    9797
    9898#eustatic + rigid run:
    99 md.slr.rigid = 1
    100 md.slr.elastic = 0
    101 md.slr.rotation = 0
     99md.solidearth.rigid = 1
     100md.solidearth.elastic = 0
     101md.solidearth.rotation = 0
    102102md = solve(md, 'Sealevelrise')
    103103Srigid = md.results.SealevelriseSolution.Sealevel
    104104
    105105#eustatic + rigid + elastic run:
    106 md.slr.rigid = 1
    107 md.slr.elastic = 1
    108 md.slr.rotation = 0
     106md.solidearth.rigid = 1
     107md.solidearth.elastic = 1
     108md.solidearth.rotation = 0
    109109md = solve(md, 'Sealevelrise')
    110110Selastic = md.results.SealevelriseSolution.Sealevel
    111111
    112112#eustatic + rigid + elastic + rotation run:
    113 md.slr.rigid = 1
    114 md.slr.elastic = 1
    115 md.slr.rotation = 1
     113md.solidearth.rigid = 1
     114md.solidearth.elastic = 1
     115md.solidearth.rotation = 1
    116116md = solve(md, 'Sealevelrise')
    117117Srotation = md.results.SealevelriseSolution.Sealevel

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

    r24919 r25125  
    1 #Test Name: EarthSlr_rotationalFeedback
     1#Test Name: EarthSolidearth_rotationalFeedback
    22import numpy as np
    33from model import *
    … …  
    1313md.mesh = gmshplanet('radius', 6.371012 * 1e3, 'resolution', 1000.)  #500 km resolution mesh
    1414
    15 #parameterize slr solution:
    16 #slr loading:  {{{
    17 md.slr.deltathickness = np.zeros((md.mesh.numberofelements, ))
    18 md.slr.sealevel = np.zeros((md.mesh.numberofvertices, ))
     15#parameterize solidearth solution:
     16#solidearth loading:  {{{
     17md.solidearth.deltathickness = np.zeros((md.mesh.numberofelements, ))
     18md.solidearth.sealevel = np.zeros((md.mesh.numberofvertices, ))
    1919md.dsl.global_average_thermosteric_sea_level_change=np.zeros((2, ))
    2020md.dsl.sea_surface_height_change_above_geoid=np.zeros((md.mesh.numberofvertices+1, ))
    … …  
    2828longe = sum(md.mesh.long[md.mesh.elements - 1], 1) / 3
    2929pos = np.intersect1d(np.array(np.where(late < -75)), np.array(np.where(longe < 0)))
    30 md.slr.deltathickness[pos] = -1
     30md.solidearth.deltathickness[pos] = -1
    3131
    3232#elastic loading from love numbers:
    3333nlov = 1000
    34 md.slr.love_h = np.array(love_numbers('h'))
    35 md.slr.love_h = np.resize(md.slr.love_h, nlov + 1)
    36 md.slr.love_k = np.array(love_numbers('k'))
    37 md.slr.love_k = np.resize(md.slr.love_k, nlov + 1)
    38 md.slr.love_l = np.array(love_numbers('l'))
    39 md.slr.love_l = np.resize(md.slr.love_l, nlov + 1)
     34md.solidearth.love_h = np.array(love_numbers('h'))
     35md.solidearth.love_h = np.resize(md.solidearth.love_h, nlov + 1)
     36md.solidearth.love_k = np.array(love_numbers('k'))
     37md.solidearth.love_k = np.resize(md.solidearth.love_k, nlov + 1)
     38md.solidearth.love_l = np.array(love_numbers('l'))
     39md.solidearth.love_l = np.resize(md.solidearth.love_l, nlov + 1)
    4040#}}}
    4141
    … …  
    5656
    5757#make sure that the elements that have loads are fully grounded:
    58 pos = np.nonzero(md.slr.deltathickness)[0]
     58pos = np.nonzero(md.solidearth.deltathickness)[0]
    5959md.mask.ocean_levelset[md.mesh.elements[pos, :] - 1] = 1
    6060
    … …  
    6565
    6666# use model representation of ocea area (not the ture area)
    67 md.slr.ocean_area_scaling = 0
     67md.solidearth.ocean_area_scaling = 0
    6868
    6969#geometry
    … …  
    8383
    8484#New stuff
    85 md.slr.spcthickness = np.nan * np.ones((md.mesh.numberofvertices, ))
    86 md.slr.Ngia = np.zeros((md.mesh.numberofvertices, ))
    87 md.slr.Ugia = np.zeros((md.mesh.numberofvertices, ))
    88 md.slr.hydro_rate = np.zeros((md.mesh.numberofvertices, ))
     85md.solidearth.spcthickness = np.nan * np.ones((md.mesh.numberofvertices, ))
     86md.solidearth.Ngia = np.zeros((md.mesh.numberofvertices, ))
     87md.solidearth.Ugia = np.zeros((md.mesh.numberofvertices, ))
     88md.solidearth.hydro_rate = np.zeros((md.mesh.numberofvertices, ))
    8989
    9090#Solution parameters
    91 md.slr.reltol = float('NaN')
    92 md.slr.abstol = 1e-3
    93 md.slr.geodetic = 1
     91md.solidearth.reltol = float('NaN')
     92md.solidearth.abstol = 1e-3
     93md.solidearth.geodetic = 1
    9494
    9595#eustatic + rigid + elastic run:
    96 md.slr.rigid = 1
    97 md.slr.elastic = 1
    98 md.slr.rotation = 0
     96md.solidearth.rigid = 1
     97md.solidearth.elastic = 1
     98md.solidearth.rotation = 0
    9999md.cluster = generic('name', gethostname(), 'np', 3)
    100100#md.verbose = verbose('111111111')
    … …  
    108108
    109109#eustatic + rigid + elastic + rotation run:
    110 md.slr.rigid = 1
    111 md.slr.elastic = 1
    112 md.slr.rotation = 1
     110md.solidearth.rigid = 1
     111md.solidearth.elastic = 1
     112md.solidearth.rotation = 1
    113113md.cluster = generic('name', gethostname(), 'np', 3)
    114114#md.verbose = verbose('111111111')

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

    r25124 r25125  
    11%Test Name: Sea-Level-Partitions
     2testagainst2002=0;
     3
    24%Data paths: {{{
    3 modeldatapath='/Users/larour/ModelData';
    45shppath='../Data/shp/';
    56gshhsshapefile=[shppath 'GSHHS_c_L1-NightlyRun.shp'];
    67%}}}
    7                
    8 testagainst2002=0;
    9 
    10 %Data paths: {{{
    11 modeldatapath='/Users/larour/ModelData';
    12 shppath='../Data/shp/';
    13 gshhsshapefile=[shppath 'GSHHS_c_L1-NightlyRun.shp'];
    14 %}}}
    158
    169%create sealevel model to hold our information:
    17         sl=sealevelmodel();
     10sl=sealevelmodel();
    1811
    1912%Create basins using boundaries from shapefile:
    … …  
    3326        boundary('shppath',shppath,'shpfilename','SouthAntarctica','proj',proj3031)...
    3427        }));
    35         %}}}
    36         %Ross: {{{
    37         sl.addbasin(basin('continent','antarctica','name','ross','proj',proj3031,'boundaries',{...
    38                 boundary('shppath',shppath,'shpfilename','SouthAntarctica','proj',proj3031),...
    39                 boundary('shppath',shppath,'shpfilename','RossIceShelf','proj',proj3031),...
    40                 boundary('shppath',shppath,'shpfilename','RossWestFront','proj',proj3031),...
    41                 boundary('shppath',shppath,'shpfilename','RossFront','proj',proj3031,'orientation','reverse')...
    42                 }));
    43 %}}}
    44         %HemisphereEast: {{{
    45         sl.addbasin(basin('continent','hemisphereeast','name','hemisphereeast','proj',laea(0,+90),'boundaries',{... %Australian projection lat,long
    46                 boundary('shppath',shppath,'shpfilename','HemisphereSplit','proj',proj4326),...
    47                 boundary('shppath',shppath,'shpfilename','SouthAntarctica','proj',proj3031),...
    48                 boundary('shppath',shppath,'shpfilename','RossFront','proj',proj3031),...
    49                 boundary('shppath',shppath,'shpfilename','RossWestFront','proj',proj3031),...
    50                 boundary('shppath',shppath,'shpfilename','EastAntarctica2','proj',proj3031,'orientation','reverse'),...
    51                 boundary('shppath',shppath,'shpfilename','NorthAntarctica','proj',proj3031)...
    52                 }));
     28%}}}
     29%Ross: {{{
     30sl.addbasin(basin('continent','antarctica','name','ross','proj',proj3031,'boundaries',{...
     31        boundary('shppath',shppath,'shpfilename','SouthAntarctica','proj',proj3031),...
     32        boundary('shppath',shppath,'shpfilename','RossIceShelf','proj',proj3031),...
     33        boundary('shppath',shppath,'shpfilename','RossWestFront','proj',proj3031),...
     34        boundary('shppath',shppath,'shpfilename','RossFront','proj',proj3031,'orientation','reverse')...
     35        }));
     36%}}}
     37%HemisphereEast: {{{
     38sl.addbasin(basin('continent','hemisphereeast','name','hemisphereeast','proj',laea(0,+90),'boundaries',{... %Australian projection lat,long
     39        boundary('shppath',shppath,'shpfilename','HemisphereSplit','proj',proj4326),...
     40        boundary('shppath',shppath,'shpfilename','SouthAntarctica','proj',proj3031),...
     41        boundary('shppath',shppath,'shpfilename','RossFront','proj',proj3031),...
     42        boundary('shppath',shppath,'shpfilename','RossWestFront','proj',proj3031),...
     43        boundary('shppath',shppath,'shpfilename','EastAntarctica2','proj',proj3031,'orientation','reverse'),...
     44        boundary('shppath',shppath,'shpfilename','NorthAntarctica','proj',proj3031)...
     45        }));
    5346%}}}
    5447%Antarctica excluding Ronne: {{{
    … …  
    5952        boundary('shppath',shppath,'shpfilename','RossIceShelf','proj',proj3031,'orientation','reverse'),...
    6053        boundary('shppath',shppath,'shpfilename','SouthAntarctica','proj',proj3031),...
    61         boundary('shppath',shppath,'shpfilename','WestAntarctica2','proj',proj3031)...
    62         boundary('shppath',shppath,'shpfilename','RonneWestSummit','proj',proj3031)...
    63         boundary('shppath',shppath,'shpfilename','RonneIceShelf','proj',proj3031)...
    64         boundary('shppath',shppath,'shpfilename','RonneEastSummit','proj',proj3031)...
     54        boundary('shppath',shppath,'shpfilename','WestAntarctica2','proj',proj3031),...
     55        boundary('shppath',shppath,'shpfilename','RonneWestSummit','proj',proj3031),...
     56        boundary('shppath',shppath,'shpfilename','RonneIceShelf','proj',proj3031),...
     57        boundary('shppath',shppath,'shpfilename','RonneEastSummit','proj',proj3031),...
    6558        boundary('shppath',shppath,'shpfilename','RonneBrunt','proj',proj3031)...
    6659        }));
    67         %}}}
    68         %Ronne: {{{
    69         sl.addbasin(basin('continent','antarctica','name','ronne','proj',proj3031,'boundaries',{...
    70                 boundary('shppath',shppath,'shpfilename','RonneWestSummit','proj',proj3031),...
    71                 boundary('shppath',shppath,'shpfilename','RonneIceShelf','proj',proj3031),...
    72                 boundary('shppath',shppath,'shpfilename','RonneEastSummit','proj',proj3031),...
    73                 boundary('shppath',shppath,'shpfilename','RonneFront','proj',proj3031,'orientation','reverse')...
    74                 }));
     60%}}}
     61%Ronne: {{{
     62sl.addbasin(basin('continent','antarctica','name','ronne','proj',proj3031,'boundaries',{...
     63        boundary('shppath',shppath,'shpfilename','RonneWestSummit','proj',proj3031),...
     64        boundary('shppath',shppath,'shpfilename','RonneIceShelf','proj',proj3031),...
     65        boundary('shppath',shppath,'shpfilename','RonneEastSummit','proj',proj3031),...
     66        boundary('shppath',shppath,'shpfilename','RonneFront','proj',proj3031,'orientation','reverse')...
     67        }));
    7568%}}}
    7669
    … …  
    372365md.timestepping.time_step=1;
    373366
    374 % max number of iteration reverted back to 10 (i.e., the original default value)
     367%max number of iterations reverted back to 10 (i.e. the original default value)
    375368md.solidearth.settings.maxiter=10;
    376369

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

    r24925 r25125  
    1414md.mesh = gmshplanet('radius', rad_e, 'resolution', 1000.0)  # km resolution
    1515
    16 #parameterize slr solution:
    17 #slr loading:  {{{
     16#parameterize solidearth solution:
     17#solidearth loading:  {{{
    1818late = sum(md.mesh.lat[md.mesh.elements - 1], 1) / 3
    1919longe = sum(md.mesh.long[md.mesh.elements - 1], 1) / 3
    2020
    21 md.slr.deltathickness = np.zeros((md.mesh.numberofelements, ))
     21md.solidearth.deltathickness = np.zeros((md.mesh.numberofelements, ))
    2222pos = np.intersect1d(np.array(np.where(late < -75)), np.array(np.where(longe > 0)))
    2323#python does not include last element in array slices, (6:7) -> [5:7]
    24 md.slr.deltathickness[pos[5:7]] = -1
     24md.solidearth.deltathickness[pos[5:7]] = -1
    2525
    26 md.slr.sealevel = np.zeros((md.mesh.numberofvertices, ))
     26md.solidearth.sealevel = np.zeros((md.mesh.numberofvertices, ))
    2727md.dsl.global_average_thermosteric_sea_level_change=np.zeros((2, ))
    2828md.dsl.sea_surface_height_change_above_geoid=np.zeros((md.mesh.numberofvertices+1, ))
    2929md.dsl.sea_water_pressure_change_at_sea_floor=np.zeros((md.mesh.numberofvertices+1, ))
    3030
    31 md.slr.ocean_area_scaling = 1
     31md.solidearth.ocean_area_scaling = 1
    3232
    3333#elastic loading from love numbers:
    3434nlov = 1000
    35 md.slr.love_h = np.array(love_numbers('h'))
    36 md.slr.love_h = np.resize(md.slr.love_h, nlov + 1)
    37 md.slr.love_k = np.array(love_numbers('k'))
    38 md.slr.love_k = np.resize(md.slr.love_k, nlov + 1)
    39 md.slr.love_l = np.array(love_numbers('l'))
    40 md.slr.love_l = np.resize(md.slr.love_l, nlov + 1)
     35md.solidearth.love_h = np.array(love_numbers('h'))
     36md.solidearth.love_h = np.resize(md.solidearth.love_h, nlov + 1)
     37md.solidearth.love_k = np.array(love_numbers('k'))
     38md.solidearth.love_k = np.resize(md.solidearth.love_k, nlov + 1)
     39md.solidearth.love_l = np.array(love_numbers('l'))
     40md.solidearth.love_l = np.resize(md.solidearth.love_l, nlov + 1)
    4141
    4242#}}}
    … …  
    5656
    5757#make sure that the elements that have loads are fully grounded:
    58 pos = np.nonzero(md.slr.deltathickness)[0]
     58pos = np.nonzero(md.solidearth.deltathickness)[0]
    5959md.mask.ocean_levelset[md.mesh.elements[pos, :] - 1] = 1
    6060
    … …  
    7979# }}}
    8080#Solution parameters {{{
    81 md.slr.reltol = float('NaN')
    82 md.slr.abstol = 1e-3
    83 md.slr.geodetic = 1
     81md.solidearth.reltol = float('NaN')
     82md.solidearth.abstol = 1e-3
     83md.solidearth.geodetic = 1
    8484# }}}
    8585
    8686#New stuff
    87 md.slr.spcthickness = np.nan * np.ones((md.mesh.numberofvertices, ))
    88 md.slr.Ngia = np.zeros((md.mesh.numberofvertices, ))
    89 md.slr.Ugia = np.zeros((md.mesh.numberofvertices, ))
    90 md.slr.hydro_rate = np.zeros((md.mesh.numberofvertices, ))
     87md.solidearth.spcthickness = np.nan * np.ones((md.mesh.numberofvertices, ))
     88md.solidearth.Ngia = np.zeros((md.mesh.numberofvertices, ))
     89md.solidearth.Ugia = np.zeros((md.mesh.numberofvertices, ))
     90md.solidearth.hydro_rate = np.zeros((md.mesh.numberofvertices, ))
    9191
    9292#eustatic + rigid + elastic run:
    93 md.slr.rigid = 1
    94 md.slr.elastic = 1
    95 md.slr.rotation = 1
     93md.solidearth.rigid = 1
     94md.solidearth.elastic = 1
     95md.solidearth.rotation = 1
    9696md.cluster = generic('name', gethostname(), 'np', 3)
    9797
    … …  
    9999md = solve(md, 'Sealevelrise')
    100100eus = md.results.SealevelriseSolution.SealevelRSLEustatic
    101 slr = md.results.SealevelriseSolution.Sealevel
     101solidearth = md.results.SealevelriseSolution.Sealevel
    102102moixz = md.results.SealevelriseSolution.SealevelInertiaTensorXZ
    103103moiyz = md.results.SealevelriseSolution.SealevelInertiaTensorYZ
    … …  
    105105
    106106# analytical moi = > just checking FOR ICE only!!! {{{
    107 # ...have to mute**slr induced MOI in Tria.cpp**prior to the comparison
     107# ...have to mute**solidearth induced MOI in Tria.cpp**prior to the comparison
    108108#rad_e = rad_e * 1e3  # now in meters
    109109#areas = GetAreasSphericalTria(md.mesh.elements, md.mesh.lat, md.mesh.long, rad_e)
    110110#lat = late * pi / 180 lon = longe * pi / 180
    111 #moi_xz = sum(-md.materials.rho_freshwater. * md.slr.deltathickness. * areas. * rad_e^2. * sin(lat). * cos(lat). * cos(lon))
    112 #moi_yz = sum(-md.materials.rho_freshwater. * md.slr.deltathickness. * areas. * rad_e^2. * sin(lat). * cos(lat). * sin(lon))
     111#moi_xz = sum(-md.materials.rho_freshwater. * md.solidearth.deltathickness. * areas. * rad_e^2. * sin(lat). * cos(lat). * cos(lon))
     112#moi_yz = sum(-md.materials.rho_freshwater. * md.solidearth.deltathickness. * areas. * rad_e^2. * sin(lat). * cos(lat). * sin(lon))
    113113# }}}
    114114
    115115#Fields and tolerances to track changes
    116 field_names = ['eus', 'slr', 'moixz', 'moiyz', 'moizz']
     116field_names = ['eus', 'solidearth', 'moixz', 'moiyz', 'moizz']
    117117field_tolerances = [1e-13, 1e-13, 1e-13, 1e-13, 1e-13]
    118 field_values = [eus, slr, moixz, moiyz, moizz]
     118field_values = [eus, solidearth, moixz, moiyz, moizz]

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

    r24862 r25125  
    3030
    3131md.timestepping.time_step = .1
    32 md.slr.sealevel = newforcing(md.timestepping.start_time, md.timestepping.final_time,
     32md.solidearth.sealevel = newforcing(md.timestepping.start_time, md.timestepping.final_time,
    3333                             md.timestepping.time_step, -200., 200., md.mesh.numberofvertices)
    3434

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

    r24261 r25125  
    3939md.geometry.bed = md.geometry.bed + 1000
    4040md.geometry.surface = md.geometry.surface + 1000
    41 md.slr.sealevel = 1000 * np.ones((md.mesh.numberofvertices, ))
     41md.solidearth.sealevel = 1000 * np.ones((md.mesh.numberofvertices, ))
    4242
    4343md = solve(md, 'Transient', 'checkconsistency', 'no')
Note: See TracChangeset for help on using the changeset viewer.