Index: /issm/trunk-jpl/etc/environment.sh =================================================================== --- /issm/trunk-jpl/etc/environment.sh (revision 25124) +++ /issm/trunk-jpl/etc/environment.sh (revision 25125) @@ -372,4 +372,11 @@ fi +NETCDF_PYTHON_DIR="${ISSM_DIR}/externalpackages/netcdf-python/install" +if [ -d "${NETCDF_PYTHON_DIR}" ]; then + if [ -d "${NETCDF_PYTHON_DIR}/lib/python2.7/site-packages" ]; then + ld_library_path_append "${NETCDF_PYTHON_DIR}/lib/python2.7/site-packages" + fi +fi + HDF5_DIR="${ISSM_DIR}/externalpackages/hdf5/install" if [ -d "${HDF5_DIR}" ]; then Index: /issm/trunk-jpl/externalpackages/netcdf-python/install.sh =================================================================== --- /issm/trunk-jpl/externalpackages/netcdf-python/install.sh (revision 25124) +++ /issm/trunk-jpl/externalpackages/netcdf-python/install.sh (revision 25125) @@ -2,29 +2,42 @@ set -eu -#Some cleanup -rm -rf install netCDF4-1.0 -mkdir install -#Download from ISSM server -$ISSM_DIR/scripts/DownloadExternalPackage.sh "https://issm.ess.uci.edu/files/externalpackages/netCDF4-1.0.tar.gz" "netCDF4-1.0.tar.gz" +# Constants +# +VER="1.5.3" -#Untar -tar -zxvf netCDF4-1.0.tar.gz +# Environment +# +export HDF5_DIR="${ISSM_DIR}/externalpackages/petsc/install" +export MPIINC_DIR="${ISSM_DIR}/externalpackages/petsc/install/include" +export NETCDF_DIR="${ISSM_DIR}/externalpackages/netcdf/install" +export 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) -#for later: -rm -rf ISSMDIR -echo $ISSM_DIR | sed 's/\//\\\//g' > ISSMDIR -ISSMDIR=`cat ISSMDIR` && rm -rf ISSMDIR +# Download source +$ISSM_DIR/scripts/DownloadExternalPackage.sh "https://issm.ess.uci.edu/files/externalpackages/netcdf4-python-${VER}rel.tar.gz" "netcdf4-python-${VER}rel.tar.gz" -#Move netCDF4 to install directory -rm -rf install/* -mv netCDF4-1.0/* install/ -rm -rf netCDF4-1.0 +# Unpack source +tar -zxvf netcdf4-python-${VER}rel.tar.gz -#Configur and compile -cd install -#edit setup.cfg to point to hdf5 and netcdf -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 +# Cleanup +rm -rf install src +mkdir install src +# Move source to 'src' directory +mv netcdf4-python-${VER}rel/* src/ +rm -rf netcdf4-python-${VER}rel + +# Compile and install +cd src python setup.py build -python setup.py install +python setup.py install --user + +# Unzip eggs +for egg in $(find "${PYTHONUSERBASE}/lib" -name *.egg); do + parent_dir=$(dirname ${egg}) + filename=$(basename -- ${egg}) + extension="${filename##*.}" + filename="${filename%.*}" + src_dir="${parent_dir}/${filename}" + unzip ${egg} ${src_dir} +done Index: /issm/trunk-jpl/externalpackages/netcdf/install-4.7-parallel-static-with_tests.sh =================================================================== --- /issm/trunk-jpl/externalpackages/netcdf/install-4.7-parallel-static-with_tests.sh (revision 25124) +++ /issm/trunk-jpl/externalpackages/netcdf/install-4.7-parallel-static-with_tests.sh (revision 25125) @@ -53,4 +53,5 @@ --enable-fast-install \ --disable-doxygen \ + --enable-netcdf4 \ --disable-filter-testing \ --disable-dap-remote-tests \ Index: /issm/trunk-jpl/externalpackages/netcdf/install-4.7-parallel-static.sh =================================================================== --- /issm/trunk-jpl/externalpackages/netcdf/install-4.7-parallel-static.sh (revision 25124) +++ /issm/trunk-jpl/externalpackages/netcdf/install-4.7-parallel-static.sh (revision 25125) @@ -53,4 +53,5 @@ --enable-fast-install \ --disable-doxygen \ + --enable-netcdf4 \ --disable-testsets \ --disable-examples Index: /issm/trunk-jpl/externalpackages/netcdf/install-4.7-parallel-with_tests.sh =================================================================== --- /issm/trunk-jpl/externalpackages/netcdf/install-4.7-parallel-with_tests.sh (revision 25124) +++ /issm/trunk-jpl/externalpackages/netcdf/install-4.7-parallel-with_tests.sh (revision 25125) @@ -50,4 +50,5 @@ --enable-fast-install \ --disable-doxygen \ + --enable-netcdf4 \ --disable-filter-testing \ --disable-dap-remote-tests \ Index: /issm/trunk-jpl/externalpackages/netcdf/install-4.7-parallel.sh =================================================================== --- /issm/trunk-jpl/externalpackages/netcdf/install-4.7-parallel.sh (revision 25124) +++ /issm/trunk-jpl/externalpackages/netcdf/install-4.7-parallel.sh (revision 25125) @@ -51,5 +51,6 @@ --disable-doxygen \ --disable-testsets \ - --disable-examples + --disable-examples \ + --enable-netcdf4 # Compile and install Index: /issm/trunk-jpl/externalpackages/netcdf/install-4.7.sh =================================================================== --- /issm/trunk-jpl/externalpackages/netcdf/install-4.7.sh (revision 25124) +++ /issm/trunk-jpl/externalpackages/netcdf/install-4.7.sh (revision 25125) @@ -50,5 +50,6 @@ --disable-doxygen \ --disable-testsets \ - --disable-examples + --disable-examples \ + --enable-netcdf4 # Compile and install Index: /issm/trunk-jpl/jenkins/ross-debian_linux-solid_earth =================================================================== --- /issm/trunk-jpl/jenkins/ross-debian_linux-solid_earth (revision 25124) +++ /issm/trunk-jpl/jenkins/ross-debian_linux-solid_earth (revision 25125) @@ -31,20 +31,21 @@ EXTERNALPACKAGES=" - autotools install-debian-linux.sh - cmake install.sh - petsc install-3.12-linux-solid_earth.sh - triangle install-linux.sh - chaco install.sh - m1qn3 install.sh - semic install.sh - boost install-1.7-linux.sh - curl install-7.67.sh - netcdf install-4.7.sh - proj install-6.2.sh - gdal install-3.0-python-netcdf.sh - gshhg install.sh - gmt install-6.0-linux.sh - gmsh install-4.sh - shell2junit install.sh + autotools install-debian-linux.sh + cmake install.sh + petsc install-3.12-linux-solid_earth.sh + triangle install-linux.sh + chaco install.sh + m1qn3 install.sh + semic install.sh + boost install-1.7-linux.sh + curl install-7.67.sh + netcdf install-4.7-parallel.sh + netcdf-python install.sh + proj install-6.2.sh + gdal install-3.0-python-netcdf.sh + gshhg install.sh + gmt install-6.0-linux.sh + gmsh install-4.sh + shell2junit install.sh " @@ -54,5 +55,5 @@ # Test suites -MATLAB_TEST=1 +MATLAB_TEST=0 PYTHON_TEST=1 JAVASCRIPT_TEST=0 Index: /issm/trunk-jpl/scripts/devpath.py =================================================================== --- /issm/trunk-jpl/scripts/devpath.py (revision 25124) +++ /issm/trunk-jpl/scripts/devpath.py (revision 25125) @@ -2,10 +2,10 @@ # NOTE: This script is a stripped-down version of -# $ISSM_DIR/src/m/dev/devpath.py and is intended only for loading ISSM in -# order to test our distributable packages. It assumes the following is -# set before $ISSM_DIR/test/NightlyRun/runme.py is called, +# $ISSM_DIR/src/m/dev/devpath.py and is intended only for loading ISSM in +# order to test our distributable packages. It assumes the following is +# set before $ISSM_DIR/test/NightlyRun/runme.py is called, # # export ISSM_DIR= -# export PATH="${PATH}:${ISSM_DIR}/bin" +# export PATH="${PATH}:${ISSM_DIR}/bin" # export PYTHONPATH="${ISSM_DIR}/scripts" # export PYTHONSTARTUP="${PYTHONPATH}/devpath.py" Index: /issm/trunk-jpl/src/m/classes/basin.py =================================================================== --- /issm/trunk-jpl/src/m/classes/basin.py (revision 25124) +++ /issm/trunk-jpl/src/m/classes/basin.py (revision 25125) @@ -1,6 +1,8 @@ import math +import numpy as np from boundary import boundary from epsg2proj import epsg2proj +from gdaltransform import gdaltransform from helpers import fileparts from laea import laea @@ -72,15 +74,29 @@ def isnameany(self, *args): #{{{ + boolean = 0 for arg in args: - if arg == self.name: - return 1 - return 0 + if type(arg) in [np.ndarray, list]: + for name in arg: + if name == self.name: + boolean = 1 + break + elif arg == self.name: + boolean = 1 + break + return boolean #}}} def iscontinentany(self, *args): #{{{ + boolean = 0 for arg in args: - if arg == self.continent: - return 1 - return 0 + if type(arg) in [np.ndarray, list]: + for continent in arg: + if continent == self.continent: + boolean = 1 + break + elif arg == self.continent: + boolean = 1 + break + return boolean #}}} @@ -130,7 +146,8 @@ y.append(y[0]) - setattr(out, 'x', x) - setattr(out, 'y', y) - setattr(out, 'nods', len(x)) + out = {} + out['x'] = x + out['y'] = y + out['nods'] = len(x) return out @@ -198,5 +215,5 @@ ba = self.contour() flags = np.zeros(len(contours)) - for i in range(len(contours)) + for i in range(len(contours)): h = contours[i] isin = inpolygon(h.x, h.y, ba.x, ba.y) Index: /issm/trunk-jpl/src/m/classes/boundary.m =================================================================== --- /issm/trunk-jpl/src/m/classes/boundary.m (revision 25124) +++ /issm/trunk-jpl/src/m/classes/boundary.m (revision 25125) @@ -74,5 +74,5 @@ end output=shpread([self.shppath '/' name '.' ext]); - + %do we reverse? if strcmpi(self.orientation,'reverse'), Index: /issm/trunk-jpl/src/m/classes/boundary.py =================================================================== --- /issm/trunk-jpl/src/m/classes/boundary.py (revision 25124) +++ /issm/trunk-jpl/src/m/classes/boundary.py (revision 25125) @@ -4,4 +4,5 @@ from epsg2proj import epsg2proj +from helpers import * from pairoptions import pairoptions from shpread import shpread @@ -75,5 +76,5 @@ def edges(self): #{{{ #read domain - path, name, ext = fileparts(shp.shpfilename) + path, name, ext = fileparts(self.shpfilename) if ext != 'shp': ext = 'shp' @@ -82,6 +83,7 @@ #do we reverse? if self.orientation == 'reverse': - output.x = np.flipud(output.x) - output.y = np.flipud(output.y) + for i in range(len(output)): + output[i].x = np.flipud(output[i].x) + output[i].y = np.flipud(output[i].y) #}}} @@ -104,5 +106,5 @@ #read domain - path, name, ext = fileparts(shp.shpfilename) + path, name, ext = fileparts(self.shpfilename) if ext != 'shp': ext = 'shp' @@ -167,5 +169,5 @@ #read domain - path, name, ext = fileparts(shp.shpfilename) + path, name, ext = fileparts(self.shpfilename) if ext != 'shp': ext = 'shp' Index: /issm/trunk-jpl/src/m/classes/fourierlove.py =================================================================== --- /issm/trunk-jpl/src/m/classes/fourierlove.py (revision 25124) +++ /issm/trunk-jpl/src/m/classes/fourierlove.py (revision 25125) @@ -34,5 +34,5 @@ string = "%s\n%s" % (string, fielddisplay(self, 'sh_nmin', 'minimum spherical harmonic degree (default 1)')) string = "%s\n%s" % (string, fielddisplay(self, 'g0', 'adimensioning constant for gravity (default 10) [m / s^2]')) - string = "%s\n%s" % (string, fielddisplay(self, 'r0', 'adimensioning constant for radius (default 6378 * 1.0e3) [m]')) + string = "%s\n%s" % (string, fielddisplay(self, 'r0', 'adimensioning constant for radius (default 6378e3) [m]')) string = "%s\n%s" % (string, fielddisplay(self, 'mu0', 'adimensioning constant for stress (default 1.0e11) [Pa]')) string = "%s\n%s" % (string, fielddisplay(self, 'allow_layer_deletion', 'allow for migration of the integration boundary with increasing spherical harmonics degree (default 1)')) @@ -66,5 +66,5 @@ self.sh_nmin = 1 self.g0 = 10 # m / s^2 - self.r0 = 6378 * 1e3 #m + self.r0 = 6378e3 #m self.mu0 = 1e11 # Pa self.allow_layer_deletion = 1 Index: /issm/trunk-jpl/src/m/classes/linearbasalforcings.py =================================================================== --- /issm/trunk-jpl/src/m/classes/linearbasalforcings.py (revision 25124) +++ /issm/trunk-jpl/src/m/classes/linearbasalforcings.py (revision 25125) @@ -1,6 +1,7 @@ +import numpy as np + +from checkfield import checkfield from fielddisplay import fielddisplay -from checkfield import checkfield from WriteData import WriteData -import numpy as np Index: /issm/trunk-jpl/src/m/classes/loadinglove.m =================================================================== --- /issm/trunk-jpl/src/m/classes/loadinglove.m (revision 25124) +++ /issm/trunk-jpl/src/m/classes/loadinglove.m (revision 25125) @@ -31,5 +31,5 @@ %tidal love numbers: self.th=0.6149; %degree 2 - self.tk=0.3055; % degree 2 + self.tk=0.3055; %degree 2 %secular fluid love number: @@ -75,5 +75,5 @@ end % }}} function marshall(self,prefix,md,fid) % {{{ - + WriteData(fid,prefix,'name','md.solidearth.love.model','data',1,'format','Integer'); @@ -85,5 +85,4 @@ WriteData(fid,prefix,'object',self,'fieldname','tk','name','md.solidearth.love.tk','format','Double'); WriteData(fid,prefix,'object',self,'fieldname','tk2secular','name','md.solidearth.love.tk2secular','format','Double'); - end % }}} Index: /issm/trunk-jpl/src/m/classes/loadinglove.py =================================================================== --- /issm/trunk-jpl/src/m/classes/loadinglove.py (revision 25125) +++ /issm/trunk-jpl/src/m/classes/loadinglove.py (revision 25125) @@ -0,0 +1,95 @@ +import numpy as np + +from checkfield import checkfield +from fielddisplay import fielddisplay +from WriteData import WriteData + + +class loadinglove(object): + ''' + LOADINGLOVE class definition + + Usage: + love = loadinglove() + ''' + + def __init__(self, *args): #{{{ + #regular love numbers + self.h = 0 #provided by PREM model + self.k = 0 #idem + self.l = 0 #idem + + #tidal love numbers for computing rotational feedback + self.th = 0 + self.tk = 0 + self.tk2secular = 0 + + nargin = len(args) + + if nargin == 0: + self.setdefaultparameters() + else: + raise Exception('constructor not supported') + #}}} + + def __repr__(self): # {{{ + s = ' loadinglove parameters:\n' + s += '{}\n'.format(fielddisplay(self, 'h', 'load Love number for radial displacement')) + s += '{}\n'.format(fielddisplay(self, 'k', 'load Love number for gravitational potential perturbation')) + s += '{}\n'.format(fielddisplay(self, 'l', 'load Love number for horizontal displacements')) + s += '{}\n'.format(fielddisplay(self, 'th', 'tidal load Love number (deg 2)')) + s += '{}\n'.format(fielddisplay(self, 'tk', 'tidal load Love number (deg 2)')) + s += '{}\n'.format(fielddisplay(self, 'tk2secular', 'secular fluid Love number')) + + return s + #}}} + + def setdefaultparameters(self): # {{{ + #tidal love numbers + self.th = 0.6149 #degree 2 + self.tk = 0.3055 #degree 2 + + #secular fluid love number + self.tk2secular = 0.942 + + return self + #}}} + + def checkconsistency(self, md, solution, analyses): # {{{ + if not 'SealevelriseAnalysis' in analyses or solution == 'TransientSolution' and md.transient.islove == 0: + return md + + md = checkfield(md, 'fieldname', 'solidearth.love.h', 'NaN', 1, 'Inf', 1) + md = checkfield(md, 'fieldname', 'solidearth.love.k', 'NaN', 1, 'Inf', 1) + md = checkfield(md, 'fieldname', 'solidearth.love.l', 'NaN', 1, 'Inf', 1) + + md = checkfield(md, 'fieldname', 'solidearth.love.th', 'NaN', 1, 'Inf', 1) + md = checkfield(md, 'fieldname', 'solidearth.love.tk', 'NaN', 1, 'Inf', 1) + md = checkfield(md, 'fieldname', 'solidearth.love.tk2secular', 'NaN', 1, 'Inf', 1) + + #check that love nunbers are provided at the same level of accuracy + if np.shape(self.h)[0] != np.shape(self.k)[0] or np.shape(self.h)[0] != np.shape(self.kl[0]): + raise Exception('love error message: love numbers should be provided at the same level of accuracy') + + return md + #}}} + + def defaultoutputs(self, md): #{{{ + return [] + #}}} + + def marshall(self, prefix, md, fid): #{{{ + WriteData(fid, prefix, 'name', 'md.solidearth.love.model', 'data', 1, 'format', 'Integer') + + WriteData(fid, prefix, 'object', self, 'fieldname', 'h', 'name', 'md.solidearth.love.h', 'format', 'DoubleMat', 'mattype', 1) + WriteData(fid, prefix, 'object', self, 'fieldname', 'k', 'name', 'md.solidearth.love.k', 'format', 'DoubleMat', 'mattype', 1) + WriteData(fid, prefix, 'object', self, 'fieldname', 'l', 'name', 'md.solidearth.love.l', 'format', 'DoubleMat', 'mattype', 1) + + WriteData(fid, prefix, 'object', self, 'fieldname', 'th', 'name', 'md.solidearth.love.th', 'format', 'Double') + WriteData(fid, prefix, 'object', self, 'fieldname', 'tk', 'name', 'md.solidearth.love.tk', 'format', 'Double') + WriteData(fid, prefix, 'object', self, 'fieldname', 'tk2secular', 'name', 'md.solidearth.love.tk2secular', 'format', 'Double') + #}}} + + def extrude(self, md): #{{{ + return self + #}}} Index: /issm/trunk-jpl/src/m/classes/matdamageice.py =================================================================== --- /issm/trunk-jpl/src/m/classes/matdamageice.py (revision 25124) +++ /issm/trunk-jpl/src/m/classes/matdamageice.py (revision 25125) @@ -86,5 +86,5 @@ self.heatcapacity = 2093. #ice latent heat of fusion L (J / kg) - self.latentheat = 3.34 * 1.0e5 + self.latentheat = 3.34e5 #ice thermal conductivity (W / m / K) self.thermalconductivity = 2.4 @@ -96,9 +96,9 @@ self.meltingpoint = 273.15 #rate of change of melting point with pressure (K / Pa) - self.beta = 9.8 * 1.0e-8 + self.beta = 9.8e-8 #mixed layer (ice-water interface) heat capacity (J / kg / K) self.mixed_layer_capacity = 3974. #thermal exchange velocity (ice-water interface) (m / s) - self.thermal_exchange_velocity = 1.00 * 1.0e-4 + self.thermal_exchange_velocity = 1.00e-4 #Rheology law: what is the temperature dependence of B with T #available: none, paterson and arrhenius @@ -106,7 +106,7 @@ # GIA: - self.lithosphere_shear_modulus = 6.7 * 1.0e10 # (Pa) + self.lithosphere_shear_modulus = 6.7e10 # (Pa) self.lithosphere_density = 3.32 # (g / cm^ - 3) - self.mantle_shear_modulus = 1.45 * 1.0e11 # (Pa) + self.mantle_shear_modulus = 1.45e11 # (Pa) self.mantle_density = 3.34 # (g / cm^ - 3) Index: /issm/trunk-jpl/src/m/classes/materials.py =================================================================== --- /issm/trunk-jpl/src/m/classes/materials.py (revision 25124) +++ /issm/trunk-jpl/src/m/classes/materials.py (revision 25125) @@ -80,5 +80,5 @@ self.heatcapacity = 2093. #ice latent heat of fusion L (J / kg) - self.latentheat = 3.34 * 1.0e5 + self.latentheat = 3.34e5 #ice thermal conductivity (W / m / K) self.thermalconductivity = 2.4 @@ -88,9 +88,9 @@ self.meltingpoint = 273.15 #rate of change of melting point with pressure (K / Pa) - self.beta = 9.8 * 1.0e-8 + self.beta = 9.8e-8 #mixed layer (ice-water interface) heat capacity (J / kg / K) self.mixed_layer_capacity = 3974. #thermal exchange velocity (ice-water interface) (m / s) - self.thermal_exchange_velocity = 1.00 * 1.0e-4 + self.thermal_exchange_velocity = 1.00e-4 #Rheology law: what is the temperature dependence of B with T #available: none, paterson and arrhenius @@ -101,12 +101,12 @@ #center of the earth (approximation, must not be 0), then the lab (lithosphere / asthenosphere boundary) then the surface #(with 1d3 to avoid numerical singularities) - self.radius = [1e3, 6278 * 1e3, 6378 * 1e3] + self.radius = [1e3, 6278e3, 6378e3] self.viscosity = [1e21, 1e40] #mantle and lithosphere viscosity (respectively) [Pa.s] - self.lame_mu = [1.45 * 1e11, 6.7 * 1e10] # (Pa) #lithosphere and mantle shear modulus (respectively) [Pa] + self.lame_mu = [1.45e11, 6.7e10] # (Pa) #lithosphere and mantle shear modulus (respectively) [Pa] self.lame_lambda = self.lame_mu # (Pa) #mantle and lithosphere lamba parameter (respectively) [Pa] self.burgers_viscosity = [np.nan, np.nan] self.burgers_mu = [np.nan, np.nan] self.isburgers = [0, 0] - self.density = [5.51 * 1e3, 5.50 * 1e3] # (Pa) #mantle and lithosphere density [kg / m^3] + self.density = [5.51e3, 5.50e3] # (Pa) #mantle and lithosphere density [kg / m^3] self.issolid = [1, 1] # is layer solid or liquid. elif nat == 'hydro': Index: /issm/trunk-jpl/src/m/classes/matestar.py =================================================================== --- /issm/trunk-jpl/src/m/classes/matestar.py (revision 25124) +++ /issm/trunk-jpl/src/m/classes/matestar.py (revision 25125) @@ -91,5 +91,5 @@ self.heatcapacity = 2093. #ice latent heat of fusion L (J / kg) - self.latentheat = 3.34 * 1.0e5 + self.latentheat = 3.34e5 #ice thermal conductivity (W / m / K) self.thermalconductivity = 2.4 @@ -101,16 +101,16 @@ self.meltingpoint = 273.15 #rate of change of melting point with pressure (K / Pa) - self.beta = 9.8 * 1.0e-8 + self.beta = 9.8e-8 #mixed layer (ice-water interface) heat capacity (J / kg / K) self.mixed_layer_capacity = 3974. #thermal exchange velocity (ice-water interface) (m / s) - self.thermal_exchange_velocity = 1.00 * 1.0e-4 + self.thermal_exchange_velocity = 1.00e-4 #Rheology law: what is the temperature dependence of B with T #available: none, paterson and arrhenius self.rheology_law = 'Paterson' # GIA: - self.lithosphere_shear_modulus = 6.7 * 1.0e10 # (Pa) + self.lithosphere_shear_modulus = 6.7e10 # (Pa) self.lithosphere_density = 3.32 # (g / cm^ - 3) - self.mantle_shear_modulus = 1.45 * 1.0e11 # (Pa) + self.mantle_shear_modulus = 1.45e11 # (Pa) self.mantle_density = 3.34 # (g / cm^ - 3) #SLR Index: /issm/trunk-jpl/src/m/classes/matice.py =================================================================== --- /issm/trunk-jpl/src/m/classes/matice.py (revision 25124) +++ /issm/trunk-jpl/src/m/classes/matice.py (revision 25125) @@ -87,5 +87,5 @@ self.heatcapacity = 2093. #ice latent heat of fusion L (J / kg) - self.latentheat = 3.34 * 1.0e5 + self.latentheat = 3.34e5 #ice thermal conductivity (W / m / K) self.thermalconductivity = 2.4 @@ -97,9 +97,9 @@ self.meltingpoint = 273.15 #rate of change of melting point with pressure (K / Pa) - self.beta = 9.8 * 1.0e-8 + self.beta = 9.8e-8 #mixed layer (ice-water interface) heat capacity (J / kg / K) self.mixed_layer_capacity = 3974. #thermal exchange velocity (ice-water interface) (m / s) - self.thermal_exchange_velocity = 1.00 * 1.0e-4 + self.thermal_exchange_velocity = 1.00e-4 #Rheology law: what is the temperature dependence of B with T #available: none, paterson and arrhenius @@ -107,7 +107,7 @@ # GIA: - self.lithosphere_shear_modulus = 6.7 * 1.0e10 # (Pa) + self.lithosphere_shear_modulus = 6.7e10 # (Pa) self.lithosphere_density = 3.32 # (g / cm^ - 3) - self.mantle_shear_modulus = 1.45 * 1.0e11 # (Pa) + self.mantle_shear_modulus = 1.45e11 # (Pa) self.mantle_density = 3.34 # (g / cm^ - 3) Index: /issm/trunk-jpl/src/m/classes/model.py =================================================================== --- /issm/trunk-jpl/src/m/classes/model.py (revision 25124) +++ /issm/trunk-jpl/src/m/classes/model.py (revision 25125) @@ -29,5 +29,5 @@ from initialization import initialization from rifts import rifts -from slr import slr +from solidearth import solidearth from dsl import dsl from debug import debug @@ -98,5 +98,5 @@ self.rifts = rifts() self.dsl = dsl() - self.slr = slr() + self.solidearth = solidearth() self.debug = debug() @@ -147,5 +147,5 @@ 'initialization', 'rifts', - 'slr', + 'solidearth', 'dsl', 'debug', @@ -194,5 +194,5 @@ string = "%s\n%s" % (string, "%19s: % - 22s - - %s" % ("initialization", "[%s, %s]" % ("1x1", obj.initialization.__class__.__name__), "initial guess / state")) string = "%s\n%s" % (string, "%19s: % - 22s - - %s" % ("rifts", "[%s, %s]" % ("1x1", obj.rifts.__class__.__name__), "rifts properties")) - string = "%s\n%s" % (string, "%19s: % - 22s - - %s" % ("slr", "[%s, %s]" % ("1x1", obj.slr.__class__.__name__), "slr forcings")) + string = "%s\n%s" % (string, "%19s: % - 22s - - %s" % ("solidearth", "[%s, %s]" % ("1x1", obj.solidearth.__class__.__name__), "solid earth inputs and settings")) string = "%s\n%s" % (string, "%19s: % - 22s - - %s" % ("dsl", "[%s, %s]" % ("1x1", obj.dsl.__class__.__name__), "dynamic sea level")) string = "%s\n%s" % (string, "%19s: % - 22s - - %s" % ("debug", "[%s, %s]" % ("1x1", obj.debug.__class__.__name__), "debugging tools (valgrind, gprof)")) Index: /issm/trunk-jpl/src/m/classes/pairoptions.py =================================================================== --- /issm/trunk-jpl/src/m/classes/pairoptions.py (revision 25124) +++ /issm/trunk-jpl/src/m/classes/pairoptions.py (revision 25125) @@ -43,5 +43,5 @@ #check length of input if len(arg) % 2: - raise TypeError('Invalid parameter / value pair arguments') + raise TypeError('Invalid parameter/value pair arguments') numoptions = int(len(arg) / 2) Index: /issm/trunk-jpl/src/m/classes/rotational.m =================================================================== --- /issm/trunk-jpl/src/m/classes/rotational.m (revision 25124) +++ /issm/trunk-jpl/src/m/classes/rotational.m (revision 25125) @@ -22,6 +22,6 @@ %moment of inertia: - self.equatorialmoi=8.0077*10^37; % [kg m^2] - self.polarmoi =8.0345*10^37; % [kg m^2] + self.equatorialmoi =8.0077*10^37; % [kg m^2] + self.polarmoi =8.0345*10^37; % [kg m^2] % mean rotational velocity of earth Index: /issm/trunk-jpl/src/m/classes/rotational.py =================================================================== --- /issm/trunk-jpl/src/m/classes/rotational.py (revision 25125) +++ /issm/trunk-jpl/src/m/classes/rotational.py (revision 25125) @@ -0,0 +1,72 @@ +import numpy as np + +from checkfield import checkfield +from fielddisplay import fielddisplay +from WriteData import WriteData + + +class rotational(object): + ''' + ROTATIONAL class definition + + Usage: + rotational = rotational() + ''' + + def __init__(self, *args): #{{{ + self.equatorialmoi = 0 + self.polarmoi = 0 + self.langularvelocity = 0 + + nargin = len(args) + + if nargin == 0: + self.setdefaultparameters() + else: + raise Exception('constructor not supported') + #}}} + + def __repr__(self): # {{{ + s = ' rotational parameters:\n' + s += '{}\n'.format(fielddisplay(self, 'equatorialmoi', 'mean equatorial moment of inertia [kg m^2]')) + s += '{}\n'.format(fielddisplay(self, 'polarmoi', 'polar moment of inertia [kg m^2]')) + s += '{}\n'.format(fielddisplay(self, 'angularvelocity', 'mean rotational velocity of earth [per second]')) + + return s + #}}} + + def setdefaultparameters(self): # {{{ + #moment of inertia + self.equatorialmoi = 8.0077e37 # [kg m^2] + self.polarmoi = 8.0345e37 # [kg m^2] + + #mean rotational velocity of earth + self.angularvelocity = 7.2921e-5 # [s^-1] + + return self + #}}} + + def checkconsistency(self, md, solution, analyses): # {{{ + if not 'SealevelriseAnalysis' in analyses or solution == 'TransientSolution' and md.transient.isrotational == 0: + return md + + md = checkfield(md, 'fieldname', 'solidearth.rotational.equatorialmoi', 'NaN', 1, 'Inf', 1) + md = checkfield(md, 'fieldname', 'solidearth.rotational.polarmoi', 'NaN', 1, 'Inf', 1) + md = checkfield(md, 'fieldname', 'solidearth.rotational.angularvelocity', 'NaN', 1, 'Inf', 1) + + return md + #}}} + + def defaultoutputs(self, md): #{{{ + return [] + #}}} + + def marshall(self, prefix, md, fid): #{{{ + WriteData(fid, prefix, 'object', self, 'fieldname', 'equatorialmoi', 'name', 'md.solidearth.rotational.equatorialmoi', 'format', 'Double') + WriteData(fid, prefix, 'object', self, 'fieldname', 'polarmoi', 'name', 'md.solidearth.rotational.polarmoi', 'format', 'Double') + WriteData(fid, prefix, 'object', self, 'fieldname', 'angularvelocity', 'name', 'md.solidearth.rotational.angularvelocity', 'format', 'Double') + #}}} + + def extrude(self, md): #{{{ + return self + #}}} Index: /issm/trunk-jpl/src/m/classes/sealevelmodel.m =================================================================== --- /issm/trunk-jpl/src/m/classes/sealevelmodel.m (revision 25124) +++ /issm/trunk-jpl/src/m/classes/sealevelmodel.m (revision 25125) @@ -220,12 +220,12 @@ if iscell(continent), if length(continent)==1, - if strcmpi(continent{1},'all'), - %need to transform this into a list of continents: - continent={}; - for i=1:length(self.basins), - continent{end+1}=self.basins{i}.continent; - end - continent=unique(continent); - end + if strcmpi(continent{1},'all'), + %need to transform this into a list of continents: + continent={}; + for i=1:length(self.basins), + continent{end+1}=self.basins{i}.continent; + end + continent=unique(continent); + end else %nothing to do, we have a list of continents @@ -236,5 +236,5 @@ continent={}; for i=1:length(self.basins), - continent{end+1}=self.basins{i}.continent; + continent{end+1}=self.basins{i}.continent; end continent=unique(continent); @@ -247,25 +247,25 @@ if iscell(bas), if length(bas)==1, - if strcmpi(bas{1},'all'), - %need to transform this into a list of basins: - baslist=[]; - for i=1:length(self.basins), - if self.basins{i}.iscontinentany(continent{:}), - baslist(end+1)=i; - end - end - baslist=unique(baslist); - else - bas=bas{1}; - baslist=[]; - for i=1:length(self.basins), - if self.basins{i}.iscontinentany(continent{:}), - if self.basins{i}.isnameany(bas), - baslist(end+1)=i; - end - end - end - - end + if strcmpi(bas{1},'all'), + %need to transform this into a list of basins: + baslist=[]; + for i=1:length(self.basins), + if self.basins{i}.iscontinentany(continent{:}), + baslist(end+1)=i; + end + end + baslist=unique(baslist); + else + bas=bas{1}; + baslist=[]; + for i=1:length(self.basins), + if self.basins{i}.iscontinentany(continent{:}), + if self.basins{i}.isnameany(bas), + baslist(end+1)=i; + end + end + end + + end else %we have a list of basin names: Index: /issm/trunk-jpl/src/m/classes/sealevelmodel.py =================================================================== --- /issm/trunk-jpl/src/m/classes/sealevelmodel.py (revision 25124) +++ /issm/trunk-jpl/src/m/classes/sealevelmodel.py (revision 25125) @@ -43,5 +43,5 @@ # Create a default object - self.setdefaultparameters() + self = self.setdefaultparameters() if len(args): @@ -77,6 +77,4 @@ self.eltransitions = [] self.planet = 'earth' - - return self #}}} @@ -191,2 +189,86 @@ self.eltransitions.append(meshintersect3d(xe, ye, ze, xei, yei, zei, 'force', force)) #}}} + + def checkintersections(self): #{{{ + flags = np.zeros(self.earth.mesh.numberofvertices, 1) + for i in range(len(self.basins)): + flags[self.transitions[i]] = i + plotmodel(self.earth, 'data', flags, 'coastline', 'on') + #}}} + + def checkbasinconsistency(self): #{{{ + for i in range(len(self.basins)): + self.basins[i].checkconsistency() + #}}} + + def basinindx(self, *args): #{{{ + options = pairoptions(*args) + continent = options.getfieldvalue('continent', 'all') + bas = options.getfieldvalue('basin', 'all') + + #expand continent list #{{{ + if type(continent) == np.ndarray: + if continent.shape[1] == 1: + if continent[0] == 'all': + #need to transform this into a list of continents + continent = [] + for i in range(len(self.basins)): + continent.append(self.basins[i].continent) + continent = np.unique(continent) + else: + pass #nothing to do, we have a list of continents + else: + if continent == 'all': + #need to transform this into a list of continents + continent = [] + for i in range(len(self.basins)): + continent.append(self.basins[i].continent) + continent = np.unique(continent) + else: + continent = [continent] + #}}} + #expand basins list using the continent list above and the extra bas discriminator #{{{ + if type(bas) == np.ndarray: + if bas.shape[1] == 1: + if bas[0] == 'all': + #need to transform this into a list of basins + baslist = [] + for i in range(len(self.basins)): + if self.basins[i].iscontinentany(continent): + baslist.append(i) + baslist = np.unique(baslist) + else: + bas = bas[0] + baslist = [] + for i in range(len(self.basins)): + if self.basins[i].iscontinentany(continent): + if self.basins[i].isnameany(bas): + baslist.append(i) + else: + #we have a list of basin names + baslist = [] + for i in range(len(bas)): + basname = bas[i] + for j in range(len(self.basins)): + if self.basins[j].iscontinentany(continent): + if self.basins[j].isnameany(basname): + baslist.append(j) + baslist = np.unique(baslist) + else: + if bas == 'all': + baslist = [] + for i in range(len(self.basins)): + if self.basins[i].iscontinentany(continent): + baslist.append(i) + baslist = np.unique(baslist) + else: + baslist = [] + for i in range(len(self.basins)): + if self.basins[i].iscontinentany(continent): + if self.basins[i].isnameany(bas): + baslist.append(i) + baslist = np.unique(baslist) + + return baslist + #}}} + #}}} Index: /issm/trunk-jpl/src/m/classes/slr.py =================================================================== --- /issm/trunk-jpl/src/m/classes/slr.py (revision 25124) +++ /issm/trunk-jpl/src/m/classes/slr.py (revision 25125) @@ -104,8 +104,8 @@ self.fluid_love = 0.942 #moment of inertia: - self.equatorial_moi = 8.0077 * 1.0e37 # [kg m^2] - self.polar_moi = 8.0345 * 1.0e37 # [kg m^2] + self.equatorial_moi = 8.0077e37 # [kg m^2] + self.polar_moi = 8.0345e37 # [kg m^2] #mean rotational velocity of earth - self.angular_velocity = 7.2921 * 1.0e-5 # [s^ - 1] + self.angular_velocity = 7.2921e-5 # [s^ - 1] #numerical discretization accuracy self.degacc = 0.01 Index: /issm/trunk-jpl/src/m/classes/solidearth.py =================================================================== --- /issm/trunk-jpl/src/m/classes/solidearth.py (revision 25125) +++ /issm/trunk-jpl/src/m/classes/solidearth.py (revision 25125) @@ -0,0 +1,102 @@ +import numpy as np + +from checkfield import checkfield +from fielddisplay import fielddisplay +from loadinglove import loadinglove +from planetradius import planetradius +from project3d import project3d +from rotational import rotational +from solidearthsettings import solidearthsettings +from surfaceload import surfaceload +from WriteData import WriteData + + +class solidearth(object): + ''' + SOLIDEARTH class definition + + Usage: + solidearth = solidearth() + ''' + + def __init__(self, *args): #{{{ + self.sealevel = np.nan + self.settings = solidearthsettings() + self.surfaceload = surfaceload() + self.love = loadinglove() + self.rotational = rotational() + self.planetradius = planetradius('earth') + self.requested_outputs = ['default'] + self.transitions = [] + + nargin = len(args) + + if nargin == 0: + self.setdefaultparameters() + else: + raise Exception('constructor not supported') + #}}} + + def __repr__(self): # {{{ + s = ' solidearthinputs, forcings and settings:\n' + s += '{}\n'.format(fielddisplay(self, 'sealevel', 'current sea level (prior to computation) [m]')) + s += '{}\n'.format(fielddisplay(self, 'planetradius', 'planet radius [m]')) + s += '{}\n'.format(fielddisplay(self, 'transitions', 'indices into parts of the mesh that will be icecaps')) + s += '{}\n'.format(fielddisplay(self, 'requested_outputs', 'additional outputs requested')) + s += '{}\n'.format(fielddisplay(self, 'miscellaneous', 'miscellaneous fields')) + s += '{}\n'.format(self.settings) + s += '{}\n'.format(self.surfaceload) + s += '{}\n'.format(self.love) + s += '{}\n'.format(self.rotational) + return s + #}}} + + def setdefaultparameters(self): # {{{ + return self + #}}} + + def checkconsistency(self, md, solution, analyses): # {{{ + if not 'SealevelriseAnalysis' in analyses or solution == 'TransientSolution' and md.transient.isslr == 0: + return md + + md = checkfield(md, 'fieldname', 'solidearth.sealevel', 'NaN', 1, 'Inf', 1, 'size', [md.mesh.numberofvertices]) + md = checkfield(md, 'fieldname', 'solidearth.requested_outputs', 'stringrow', 1) + + self.settings.checkconsistency(md, solution, analyses) + self.surfaceload.checkconsistency(md, solution, analyses) + self.love.checkconsistency(md, solution, analyses) + self.rotational.checkconsistency(md, solution, analyses) + + return md + #}}} + + def defaultoutputs(self, md): #{{{ + return ['Sealevel'] + #}}} + + def marshall(self, prefix, md, fid): #{{{ + WriteData(fid, prefix, 'object', self, 'fieldname', 'sealevel', 'mattype', 1, 'format', 'DoubleMat', 'timeserieslength', md.mesh.numberofvertices + 1, 'yts', md.constants.yts) + WriteData(fid, prefix, 'object', self, 'fieldname', 'planetradius', 'format', 'Double') + WriteData(fid, prefix, 'object', self, 'fieldname', 'transitions', 'format', 'MatArray') + + self.settings.marshall(prefix, md, fid) + self.surfaceload.marshall(prefix, md, fid) + self.love.marshall(prefix, md, fid) + self.rotational.marshall(prefix, md, fid) + + #process requested outputs + outputs = self.requested_outputs + pos = np.where(np.asarray(outputs) == 'default')[0] + for i in pos: + outputs[i] = '' #remove 'default' from outputs + outputs.extend(self.defaultoutputs(md)) #add defaults + + WriteData(fid, prefix, 'data', outputs, 'name', 'md.solidearth.requested_outputs', 'format', 'StringArray') + #}}} + + def extrude(self, md): #{{{ + self.sealevel = project3d(md, 'vector', self.sealevel, 'type', 'node') + + return self + #}}} + Index: /issm/trunk-jpl/src/m/classes/solidearthsettings.m =================================================================== --- /issm/trunk-jpl/src/m/classes/solidearthsettings.m (revision 25124) +++ /issm/trunk-jpl/src/m/classes/solidearthsettings.m (revision 25125) @@ -46,5 +46,5 @@ self.degacc=.01; - %how many time steps we skip before we run SOLIDEARTHSETTINGS solver during transient + %how many time steps we skip before we run solidearthsettings solver during transient self.runfrequency=1; @@ -86,5 +86,5 @@ fielddisplay(self,'ocean_area_scaling','correction for model representation of ocean area [default: No correction]'); fielddisplay(self,'computesealevelchange','compute sealevel change from GRD in addition to steric?) default 0'); - fielddisplay(self,'runfrequency','how many time steps we skip before we run SOLIDEARTHSETTINGS solver during transient (default: 1)'); + fielddisplay(self,'runfrequency','how many time steps we skip before we run solidearthsettings solver during transient (default: 1)'); fielddisplay(self,'rigid','rigid earth graviational potential perturbation'); fielddisplay(self,'elastic','elastic earth graviational potential perturbation'); Index: /issm/trunk-jpl/src/m/classes/solidearthsettings.py =================================================================== --- /issm/trunk-jpl/src/m/classes/solidearthsettings.py (revision 25125) +++ /issm/trunk-jpl/src/m/classes/solidearthsettings.py (revision 25125) @@ -0,0 +1,109 @@ +import numpy as np + +from checkfield import checkfield +from fielddisplay import fielddisplay +from WriteData import WriteData + + +class solidearthsettings(object): + ''' + SOLIDEARTHSETTINGS class definition + + Usage: + solidearthsettings = solidearthsettings() + ''' + + def __init__(self, *args): #{{{ + nargin = len(args) + + if nargin == 0: + self.setdefaultparameters() + else: + raise Exception('constructor not supported') + #}}} + + def __repr__(self): # {{{ + s = ' solidearth settings:\n' + s += '{}\n'.format(fielddisplay(self, 'reltol', 'sea level rise relative convergence criterion, (default, NaN: not applied)')) + s += '{}\n'.format(fielddisplay(self, 'abstol', 'sea level rise absolute convergence criterion, NaN: not applied')) + s += '{}\n'.format(fielddisplay(self, 'maxiter', 'maximum number of nonlinear iterations')) + s += '{}\n'.format(fielddisplay(self, 'ocean_area_scaling', 'correction for model representation of ocean area [default: No correction]')) + s += '{}\n'.format(fielddisplay(self, 'computesealevelchange', 'compute sealevel change from GRD in addition to steric?) default 0')) + s += '{}\n'.format(fielddisplay(self, 'runfrequency', 'how many time steps we skip before we run solidearthsettings solver during transient (default: 1)')) + s += '{}\n'.format(fielddisplay(self, 'rigid', 'rigid earth graviational potential perturbation')) + s += '{}\n'.format(fielddisplay(self, 'elastic', 'elastic earth graviational potential perturbation')) + s += '{}\n'.format(fielddisplay(self, 'degacc', 'accuracy (default .01 deg) for numerical discretization of the Green\'s functions')) + + return s + #}}} + + def setdefaultparameters(self): # {{{ + #Convergence criterion: absolute, relative, and residual + self.reltol = 0.01 # 1 percent + self.abstol = np.nan # default + + #maximum of non-linear iterations + self.maxiter = 5 + + #computational flags + self.rigid = 1 + self.elastic = 1 + self.rotation = 1 + self.ocean_area_scaling = 0 + self.computesealevelchange = 0 + + #numerical discetization accuracy + self.degacc = .01 + + #how many time steps we skip before we run solidearthsettings solver during transient + self.runfrequency = 1 + + #horizontal displacemnet? (not by default) + self.horiz = 0 + + return self + #}}} + + def checkconsistency(self, md, solution, analyses): # {{{ + if not 'SealevelriseAnalysis' in analyses or solution == 'TransientSolution' and md.transient.issolidearthsettings == 0: + return md + + md = checkfield(md, 'fieldname', 'solidearth.settings.reltol', 'size', [1, 1]) + md = checkfield(md, 'fieldname', 'solidearth.settings.abstol', 'size', [1, 1]) + md = checkfield(md, 'fieldname', 'solidearth.settings.maxiter', 'size', [1, 1], '>=', 1) + md = checkfield(md, 'fieldname', 'solidearth.settings.runfrequency', 'size', [1, 1], '>=', 1) + md = checkfield(md, 'fieldname', 'solidearth.settings.degacc', 'size', [1, 1], '>=', 1e-10) + md = checkfield(md, 'fieldname', 'solidearth.settings.horiz', 'NaN', 1, 'Inf', 1, 'values', [0, 1]) + + #a coupler to planet model is provided + if self.computesealevelchange: + if md.transient.iscoupler: + #we are good + pass + else: + if md.mesh.__class__.__name__ == 'mesh3dsurface': + #we are good + pass + else: + raise Exception('model is requesting sealevel computations without being a mesh3dsurface, or being coupled to one!') + + return md + #}}} + + def marshall(self, prefix, md, fid): #{{{ + WriteData(fid, prefix, 'object', self, 'fieldname', 'reltol', 'name', 'md.solidearth.settings.reltol', 'format', 'Double') + WriteData(fid,prefix,'object', self,'fieldname', 'abstol', 'name', 'md.solidearth.settings.abstol', 'format', 'Double') + WriteData(fid,prefix,'object', self,'fieldname', 'maxiter', 'name', 'md.solidearth.settings.maxiter', 'format', 'Integer') + WriteData(fid,prefix,'object', self,'fieldname', 'rigid', 'name', 'md.solidearth.settings.rigid', 'format', 'Boolean') + WriteData(fid,prefix,'object', self,'fieldname', 'elastic', 'name', 'md.solidearth.settings.elastic', 'format', 'Boolean') + WriteData(fid,prefix,'object', self,'fieldname', 'rotation', 'name', 'md.solidearth.settings.rotation', 'format', 'Boolean') + WriteData(fid,prefix,'object', self,'fieldname', 'ocean_area_scaling', 'name', 'md.solidearth.settings.ocean_area_scaling', 'format', 'Boolean') + WriteData(fid,prefix,'object', self,'fieldname', 'runfrequency', 'name', 'md.solidearth.settings.runfrequency', 'format', 'Integer') + WriteData(fid,prefix,'object', self,'fieldname', 'degacc', 'name', 'md.solidearth.settings.degacc', 'format', 'Double') + WriteData(fid,prefix,'object', self,'fieldname', 'horiz', 'name', 'md.solidearth.settings.horiz', 'format', 'Integer') + WriteData(fid,prefix,'object', self,'fieldname', 'computesealevelchange', 'name', 'md.solidearth.settings.computesealevelchange', 'format', 'Integer') + #}}} + + def extrude(self, md): #{{{ + return self + #}}} Index: /issm/trunk-jpl/src/m/classes/surfaceload.m =================================================================== --- /issm/trunk-jpl/src/m/classes/surfaceload.m (revision 25124) +++ /issm/trunk-jpl/src/m/classes/surfaceload.m (revision 25125) @@ -60,7 +60,13 @@ function marshall(self,prefix,md,fid) % {{{ - if isempty(self.icethicknesschange), self.icethicknesschange=zeros(md.mesh.numberofelements+1,1); end - if isempty(self.waterheightchange), self.waterheightchange=zeros(md.mesh.numberofelements+1,1); end - if isempty(self.other), self.other=zeros(md.mesh.numberofelements+1,1); end + if isempty(self.icethicknesschange), + self.icethicknesschange=zeros(md.mesh.numberofelements+1,1); + end + if isempty(self.waterheightchange), + self.waterheightchange=zeros(md.mesh.numberofelements+1,1); + end + if isempty(self.other), + self.other=zeros(md.mesh.numberofelements+1,1); + end 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); Index: /issm/trunk-jpl/src/m/classes/surfaceload.py =================================================================== --- /issm/trunk-jpl/src/m/classes/surfaceload.py (revision 25125) +++ /issm/trunk-jpl/src/m/classes/surfaceload.py (revision 25125) @@ -0,0 +1,75 @@ +import numpy as np + +from checkfield import checkfield +from fielddisplay import fielddisplay +from WriteData import WriteData + + +class surfaceload(object): + ''' + SURFACELOAD class definition + + Usage: + surfaceload = surfaceload() + ''' + + def __init__(self, *args): #{{{ + self.icethicknesschange = [] + self.waterheightchange = [] + self.other = [] + + nargin = len(args) + + if nargin == 0: + self.setdefaultparameters() + else: + raise Exception('constructor not supported') + #}}} + + def __repr__(self): # {{{ + s = ' surfaceload:\n' + s += '{}\n'.format(fielddisplay(self, 'icethicknesschange', 'thickness change: ice height equivalent [mIce/yr]')) + s += '{}\n'.format(fielddisplay(self, 'waterheightchange', 'water height change: water height equivalent [mWater/yr]')) + s += '{}\n'.format(fielddisplay(self, 'other', 'other loads (sediments) [kg/m^2/yr]')) + + return s + #}}} + + def setdefaultparameters(self): # {{{ + return self + #}}} + + def checkconsistency(self, md, solution, analyses): # {{{ + if not 'SealevelriseAnalysis' in analyses or solution == 'TransientSolution' and md.transient.issurfaceload == 0: + return md + + if len(self.icethicknesschange) > 0: + md = checkfield(md,'fieldname', 'solidearth.surfaceload.icethicknesschange', 'timeseries', 1, 'NaN', 1, 'Inf', 1) + + if len(self.waterheightchange) > 0: + md = checkfield(md,'fieldname', 'solidearth.surfaceload.waterheightchange', 'timeseries', 1, 'NaN', 1, 'Inf', 1) + + if len(self.other) > 0: + md = checkfield(md,'fieldname', 'solidearth.surfaceload.other', 'timeseries', 1, 'NaN', 1, 'Inf', 1) + + return md + #}}} + + def marshall(self, prefix, md, fid): #{{{ + if len(self.icethicknesschange) == 0: + self.icethicknesschange = np.zeros(md.mesh.numberofelements + 1) + + if len(self.waterheightchange) == 0: + self.waterheightchange = np.zeros(md.mesh.numberofelements + 1) + + if len(self.other) == 0: + self.other = np.zeros(md.mesh.numberofelements + 1) + + 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) + 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) + WriteData(fid, prefix, 'object', self, 'fieldname', 'other', 'name', 'md.solidearth.surfaceload.other', 'format', 'DoubleMat', 'mattype', 2, 'timeserieslength', md.mesh.numberofelements + 1, 'yts', md.constants.yts) + #}}} + + def extrude(self, md): #{{{ + return self + #}}} Index: /issm/trunk-jpl/src/m/classes/tidallove.py =================================================================== --- /issm/trunk-jpl/src/m/classes/tidallove.py (revision 25125) +++ /issm/trunk-jpl/src/m/classes/tidallove.py (revision 25125) @@ -0,0 +1,74 @@ +import numpy as np + +from checkfield import checkfield +from fielddisplay import fielddisplay +from WriteData import WriteData + + +class tidallove(object): + ''' + TIDALLOVE class definition + + Usage: + tidallove = tidallove() + ''' + + def __init__(self, *args): #{{{ + self.h = 0 + self.k = 0 + self.secular = 0 + + nargin = len(args) + + if nargin == 0: + self.setdefaultparameters() + else: + raise Exception('constructor not supported') + #}}} + + def __repr__(self): # {{{ + s = ' tidallove parameters:\n' + s += '{}\n'.format(fielddisplay(self, 'h', 'tidal load Love number (deg 2)')) + s += '{}\n'.format(fielddisplay(self, 'k', 'tidal load Love number (deg 2)')) + s += '{}\n'.format(fielddisplay(self, 'k2secular', 'secular fluid Love number')) + + return s + #}}} + + def setdefaultparameters(self): # {{{ + #tidal love numbers + self.h = 0.6149 # degree 2 + self.k = 0.3055 # degree 2 + + #secular fluid love number + self.k2secular = 0.942 + + return self + #}}} + + def checkconsistency(self, md, solution, analyses): # {{{ + if not 'SealevelriseAnalysis' in analyses or solution == 'TransientSolution' and md.transient.istidallove == 0: + return md + + md = checkfield(md, 'fieldname', 'tidallove.h', 'NaN', 1, 'Inf', 1) + md = checkfield(md, 'fieldname', 'tidallove.k', 'NaN', 1, 'Inf', 1) + md = checkfield(md, 'fieldname', 'tidallove.k2secular', 'NaN', 1, 'Inf', 1) + + return md + #}}} + + def defaultoutputs(self, md): #{{{ + return [] + #}}} + + def marshall(self, prefix, md, fid): #{{{ + WriteData(fid, prefix, 'name', 'md.solidearth.love.model', 'data', 2, 'format', 'Integer') + WriteData(fid, prefix, 'object', self, 'fieldname', 'h', 'name', 'md.solidearth.love.h', 'format', 'Double', 'mattype', 1) + WriteData(fid, prefix, 'object', self, 'fieldname', 'k', 'name', 'md.solidearth.love.k', 'format', 'Double', 'mattype', 1) + WriteData(fid, prefix, 'object', self, 'fieldname', 'k2secular', 'name', 'md.solidearth.love.k2secular', 'format', 'DoubleMat', 'mattype', 1) + #}}} + + def extrude(self, md): #{{{ + return self + #}}} + Index: /issm/trunk-jpl/src/m/coordsystems/gdaltransform.py =================================================================== --- /issm/trunk-jpl/src/m/coordsystems/gdaltransform.py (revision 25124) +++ /issm/trunk-jpl/src/m/coordsystems/gdaltransform.py (revision 25125) @@ -43,5 +43,5 @@ filename_out = file_out.name - file_in.write(b'%8g %8g\n' % (x.flatten(1) y.flatten(1))) + file_in.write(b'%8g %8g\n' % (x.flatten(1), y.flatten(1))) file_in.close() Index: /issm/trunk-jpl/src/m/dev/devpath.py =================================================================== --- /issm/trunk-jpl/src/m/dev/devpath.py (revision 25124) +++ /issm/trunk-jpl/src/m/dev/devpath.py (revision 25125) @@ -11,5 +11,24 @@ raise NameError('"ISSM_DIR" environment variable is empty! You should define ISSM_DIR in your .cshrc or .bashrc!') - #Go through src/m and append any directory that contains a *.py file to PATH +#Get paths to locally installed Python libraries +tpls = [ + 'cftime', # Needed for netCDF4 + 'netCDF4' +] +for tpl in os.listdir(ISSM_DIR + '/externalpackages'): + tpl_path = ISSM_DIR + '/externalpackages/' + tpl + for dirpath, dirnames, filenames in os.walk(tpl_path): + install_lib_path = dirpath + '/install/lib' + for dirpath, dirnames, filenames in os.walk(install_lib_path): + for dirname in dirnames: + if dirname in tpls: + dir = dirpath + '/' + dirname + for file in os.listdir(dir): + if file.find(".py") != -1: + if file.find(".pyc") == -1: + if dirpath not in sys.path: + sys.path.append(dirpath) + +#Go through src/m and append any directory that contains a *.py file to PATH for root, dirs, files in os.walk(ISSM_DIR + '/src/m'): if '.svn' in dirs: @@ -21,5 +40,5 @@ sys.path.append(root) - #Also add the Nightly run directory +#Also add the Nightly run directory if ISSM_DIR + '/test/NightlyRun' not in sys.path: sys.path.append(ISSM_DIR + '/test/NightlyRun') @@ -28,5 +47,6 @@ if ISSM_DIR + '/src/wrappers/python/.libs' not in sys.path: sys.path.append(ISSM_DIR + '/src/wrappers/python/.libs') - # If using clusters, we need to have the path to the cluster settings directory + +# If using clusters, we need to have the path to the cluster settings directory if JPL_SVN is not None: jpl_path = JPL_SVN + '/usr/' + USERNAME Index: /issm/trunk-jpl/src/m/geometry/polyarea.py =================================================================== --- /issm/trunk-jpl/src/m/geometry/polyarea.py (revision 25124) +++ /issm/trunk-jpl/src/m/geometry/polyarea.py (revision 25125) @@ -4,5 +4,5 @@ -def polyarea: #{{{ +def polyarea(x, y): #{{{ ''' POLYAREA - returns the area of the 2-D polygon defined by the vertices in Index: /issm/trunk-jpl/src/m/io/loadvars.py =================================================================== --- /issm/trunk-jpl/src/m/io/loadvars.py (revision 25124) +++ /issm/trunk-jpl/src/m/io/loadvars.py (revision 25125) @@ -1,7 +1,8 @@ +from collections import OrderedDict +from netCDF4 import Dataset +import numpy as np +from re import findall import shelve -import numpy as np -from netCDF4 import Dataset -from re import findall -from collections import OrderedDict + from model import * #hack to keep python 2 compatibility @@ -15,12 +16,12 @@ def loadvars(*args, OL): - """ - LOADVARS - function to load variables to a file. - - This function loads one or more variables from a file. The names of the variables - must be supplied. If more than one variable is specified, it may be done with - a list of names or a dictionary of name as keys. The output type will correspond - to the input type. All the variables in the file may be loaded by specifying only - the file name. + ''' + LOADVARS - function to load variables from a file. + + This function loads one or more variables from a file. The names of the + variables must be supplied. If more than one variable is specified, it may + be done with a list of names or a dictionary of name as keys. The output + type will correspond to the input type. All the variables in the file may + be loaded by specifying only the file name. Usage: @@ -29,6 +30,5 @@ nvdict = loadvars('shelve.dat', {'a':None, 'b':None}) nvdict = loadvars('shelve.dat') - - """ + ''' filename = '' @@ -40,5 +40,4 @@ if not filename: filename = '/tmp/shelve.dat' - else: raise TypeError("Missing file name.") @@ -47,15 +46,11 @@ for name in args[1:]: nvdict[name] = None - elif len(args) == 2 and isinstance(args[1], list): # (filename, [names]) for name in args[1]: nvdict[name] = None - elif len(args) == 2 and isinstance(args[1], dict): # (filename, {names:values}) nvdict = args[1] - elif len(args) == 1: # (filename) pass - else: raise TypeError("Unrecognized input arguments.") @@ -78,5 +73,4 @@ print(("Variable '%s' loaded." % name)) my_shelf.close() - else: #We used netcdf for the save try: Index: /issm/trunk-jpl/src/m/materials/nye.py =================================================================== --- /issm/trunk-jpl/src/m/materials/nye.py (revision 25124) +++ /issm/trunk-jpl/src/m/materials/nye.py (revision 25125) @@ -40,5 +40,5 @@ n = 8. # Glen's exponent elif ice_type == 2: # H2O ice - A_const = 9 * 1.0e4 # s^-1 MPa + A_const = 9e4 # s^-1 MPa Q = 60000. # J mol^-1 n = 3. # Glen's exponent Index: /issm/trunk-jpl/src/m/plot/applyoptions.py =================================================================== --- /issm/trunk-jpl/src/m/plot/applyoptions.py (revision 25124) +++ /issm/trunk-jpl/src/m/plot/applyoptions.py (revision 25125) @@ -5,5 +5,5 @@ from matplotlib.ticker import MaxNLocator except ImportError: - print("could not import pylab, matplotlib has not been installed, no plotting capabilities enabled") + print("could not import pyplot, matplotlib has not been installed, no plotting capabilities enabled") import numpy as np Index: /issm/trunk-jpl/src/m/plot/plot_unit.py =================================================================== --- /issm/trunk-jpl/src/m/plot/plot_unit.py (revision 25124) +++ /issm/trunk-jpl/src/m/plot/plot_unit.py (revision 25125) @@ -7,5 +7,5 @@ from mpl_toolkits.mplot3d.art3d import Poly3DCollection except ImportError: - print("could not import pylab, matplotlib has not been installed, no plotting capabilities enabled") + print("could not import pyplot, matplotlib has not been installed, no plotting capabilities enabled") import numpy as np Index: /issm/trunk-jpl/src/m/plot/plotmodel.py =================================================================== --- /issm/trunk-jpl/src/m/plot/plotmodel.py (revision 25124) +++ /issm/trunk-jpl/src/m/plot/plotmodel.py (revision 25125) @@ -5,5 +5,5 @@ from mpl_toolkits.axes_grid1 import ImageGrid except ImportError: - print("could not import pylab, matplotlib has not been installed, no plotting capabilities enabled") + print("could not import pyplot, matplotlib has not been installed, no plotting capabilities enabled") import numpy as np Index: /issm/trunk-jpl/src/m/shp/shpread.m =================================================================== --- /issm/trunk-jpl/src/m/shp/shpread.m (revision 25124) +++ /issm/trunk-jpl/src/m/shp/shpread.m (revision 25125) @@ -21,4 +21,5 @@ %recover options options=pairoptions(varargin{:}); +invert=getfieldvalue(options,'invert',0); %some checks @@ -97,5 +98,4 @@ end -invert=getfieldvalue(options,'invert',0); if invert, for i=1:length(Struct), Index: /issm/trunk-jpl/src/m/shp/shpread.py =================================================================== --- /issm/trunk-jpl/src/m/shp/shpread.py (revision 25124) +++ /issm/trunk-jpl/src/m/shp/shpread.py (revision 25125) @@ -5,4 +5,5 @@ print("could not import shapefile, PyShp has not been installed, no shapefile reading capabilities enabled") +from helpers import OrderedStruct from pairoptions import pairoptions @@ -47,4 +48,8 @@ Sources: - https://github.com/GeospatialPython/pyshp + + TODO: + - Create class that can be used to store and pretty print shape structs + (ala OrderedStruct from src/m/qmu/helpers.py). ''' @@ -59,22 +64,14 @@ sf = shapefile.Reader(filename) - Dicts = [] - - #retrieve ID (if it exists) - fields = sf.fields - for i in range(1, len(fields)): # skip over first field, which is "DeletionFlag" - if fields[i][0] == 'id': - name = str(sf.record(i - 1)[0]) # record index is offset by one, again, because of "DeletionFlag" - break - + Structs = [] shapes = sf.shapes() for i in range(len(shapes)): - Dict = {} + Struct = OrderedStruct() shape = shapes[i] if shape.shapeType == shapefile.POINT: - Dict['x'] = shape.points[0][0] - Dict['y'] = shape.points[0][1] - Dict['density'] = 1 - Dict['Geometry'] = 'Point' + Struct.x = shape.points[0][0] + Struct.y = shape.points[0][1] + Struct.density = 1 + Struct.Geometry = 'Point' elif shape.shapeType == shapefile.POLYLINE: num_points = len(shape.points) @@ -85,11 +82,11 @@ x.append(point[0]) y.append(point[1]) - Dict['x'] = x - Dict['y'] = y - Dict['nods'] = num_points - Dict['density'] = 1 - Dict['closed'] = 1 - Dict['BoundingBox'] = shape.bbox - Dict['Geometry'] = 'Line' + Struct.x = x + Struct.y = y + Struct.nods = num_points + Struct.density = 1 + Struct.closed = 1 + Struct.BoundingBox = shape.bbox + Struct.Geometry = 'Line' elif shape.shapeType == shapefile.POLYGON: num_points = len(shape.points) @@ -100,11 +97,11 @@ x.append(point[0]) y.append(point[1]) - Dict['x'] = x - Dict['y'] = y - Dict['nods'] = num_points - Dict['density'] = 1 - Dict['closed'] = 1 - Dict['BoundingBox'] = shape.bbox - Dict['Geometry'] = 'Polygon' + Struct.x = x + Struct.y = y + Struct.nods = num_points + Struct.density = 1 + Struct.closed = 1 + Struct.BoundingBox = shape.bbox + Struct.Geometry = 'Polygon' else: # NOTE: We could do this once before looping over shapes as all @@ -115,4 +112,5 @@ # raise Exception('shpread error: geometry {} is not currently supported'.format(shape.shapeTypeName)) + name = '' fields = sf.fields @@ -121,16 +119,16 @@ # 'id' field gets special treatment if fieldname == 'id': - name = str(sf.record(j - 1)[0]) # record index is offset by one, again, because of "DeletionFlag" + name = str(sf.record(i)[j - 1]) # record index is offset by one, again, because of "DeletionFlag" else: - Dict[fieldname] = sf.record(j - 1)[0] - Dict['name'] = name - Dicts.append(Dict) + setattr(Struct, str(fieldname), sf.record(i)[j - 1]) # cast to string removes "u" from "u'fieldname'" + Struct.name = name + Structs.append(Struct) invert = options.getfieldvalue('invert', 0) if invert: - for i in range(len(Dicts)): - Dicts[i].x = np.flipud(Dicts[i].x) - Dicts[i].y = np.flipud(Dicts[i].y) + for i in range(len(Structs)): + Structs[i].x = np.flipud(Structs[i].x) + Structs[i].y = np.flipud(Structs[i].y) - return Dicts + return Structs #}}} Index: /issm/trunk-jpl/test/NightlyRun/test2001.py =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test2001.py (revision 25124) +++ /issm/trunk-jpl/test/NightlyRun/test2001.py (revision 25125) @@ -1,21 +1,30 @@ #Test Name: SquareSheetConstrainedGia2d +from socket import gethostname + import numpy as np + from model import * -from socket import gethostname -from triangle import * -from setmask import * from parameterize import * from setflowequation import * +from setmask import * from solve import * +from triangle import * -#Define a model -md = model() -md = triangle(md, '../Exp/Square.exp', 100000.) +md = triangle(model(), '../Exp/Square.exp', 100000.) md = setmask(md, '', '') md = parameterize(md, '../Par/SquareSheetConstrained.py') +#GIA +md.gia = giaivins() +md.gia.lithosphere_thickness = 100. * np.ones(md.mesh.numberofvertices) # in km +md.gia.mantle_viscosity = 1.0e21 * np.ones(md.mesh.numberofvertices) # in Pa.s +md.materials.lithosphere_shear_modulus = 6.7e10 # in Pa +md.materials.lithosphere_density = 3.32 # in g/cm^3 +md.materials.mantle_shear_modulus = 1.45e11 # in Pa +md.materials.mantle_density = 3.34 # in g/cm^3 + #Indicate what you want to compute -md.gia.cross_section_shape = 1 # for square-edged x - section +md.gia.cross_section_shape = 1 # for square-edged x-section #Define loading history (see test2001.m for the description) @@ -35,4 +44,4 @@ field_names = ['GiaW', 'GiadWdt'] field_tolerances = [1e-13, 1e-13] -field_values = [md.results.GiaSolution.GiaW, - md.results.GiaSolution.GiadWdt] +field_values = [md.results.GiaSolution.UGia, + md.results.GiaSolution.UGiaRate] Index: /issm/trunk-jpl/test/NightlyRun/test2002.py =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test2002.py (revision 25124) +++ /issm/trunk-jpl/test/NightlyRun/test2002.py (revision 25125) @@ -1,3 +1,3 @@ -#Test Name: EarthSlr +#Test Name: EarthSolidearth import numpy as np @@ -15,8 +15,8 @@ md.mesh = gmshplanet('radius', 6.371012 * 10**3, 'resolution', 700.) #500 km resolution mesh -#parameterize slr solution: -#slr loading: -md.slr.deltathickness = np.zeros((md.mesh.numberofelements)) -md.slr.sealevel = np.zeros((md.mesh.numberofvertices)) +#parameterize solidearth solution: +#solidearth loading: +md.solidearth.deltathickness = np.zeros((md.mesh.numberofelements)) +md.solidearth.sealevel = np.zeros((md.mesh.numberofvertices)) md.dsl.global_average_thermosteric_sea_level_change=np.zeros((2, )) md.dsl.sea_surface_height_change_above_geoid=np.zeros((md.mesh.numberofvertices+1, )) @@ -28,14 +28,14 @@ longe = np.sum(md.mesh.long[md.mesh.elements - 1], axis=1) / 3 pos = np.where(late < -80) -md.slr.deltathickness[pos] = -100 +md.solidearth.deltathickness[pos] = -100 #greenland pos = np.where(np.logical_and.reduce((late > 70, late < 80, longe > -60, longe < -30))) -md.slr.deltathickness[pos] = -100 +md.solidearth.deltathickness[pos] = -100 #elastic loading from love numbers: nlov = 101 -md.slr.love_h = love_numbers('h')[:nlov] -md.slr.love_k = love_numbers('k')[:nlov] -md.slr.love_l = love_numbers('l')[:nlov] +md.solidearth.love_h = love_numbers('h')[:nlov] +md.solidearth.love_k = love_numbers('k')[:nlov] +md.solidearth.love_l = love_numbers('l')[:nlov] #mask: @@ -54,5 +54,5 @@ #make sure that the elements that have loads are fully grounded: -pos = np.nonzero(md.slr.deltathickness)[0] +pos = np.nonzero(md.solidearth.deltathickness)[0] md.mask.ocean_levelset[md.mesh.elements[pos, :] - 1] = 1 @@ -61,5 +61,5 @@ md.mask.ice_levelset = icemask -md.slr.ocean_area_scaling = 0 +md.solidearth.ocean_area_scaling = 0 #geometry for the bed, arbitrary @@ -73,8 +73,8 @@ #New stuff -md.slr.spcthickness = np.nan * np.ones((md.mesh.numberofvertices, )) -md.slr.Ngia = np.zeros((md.mesh.numberofvertices, )) -md.slr.Ugia = np.zeros((md.mesh.numberofvertices, )) -md.slr.hydro_rate = np.zeros((md.mesh.numberofvertices, )) +md.solidearth.spcthickness = np.nan * np.ones((md.mesh.numberofvertices, )) +md.solidearth.Ngia = np.zeros((md.mesh.numberofvertices, )) +md.solidearth.Ugia = np.zeros((md.mesh.numberofvertices, )) +md.solidearth.hydro_rate = np.zeros((md.mesh.numberofvertices, )) #Miscellaneous @@ -82,36 +82,36 @@ #Solution parameters -md.slr.reltol = np.nan -md.slr.abstol = 1e-3 -md.slr.geodetic = 1 +md.solidearth.reltol = np.nan +md.solidearth.abstol = 1e-3 +md.solidearth.geodetic = 1 #max number of iteration reverted back to 10 (i.e., the original default value) -md.slr.maxiter = 10 +md.solidearth.maxiter = 10 #eustatic run: -md.slr.rigid = 0 -md.slr.elastic = 0 -md.slr.rotation = 0 +md.solidearth.rigid = 0 +md.solidearth.elastic = 0 +md.solidearth.rotation = 0 md = solve(md, 'Sealevelrise') Seustatic = md.results.SealevelriseSolution.Sealevel #eustatic + rigid run: -md.slr.rigid = 1 -md.slr.elastic = 0 -md.slr.rotation = 0 +md.solidearth.rigid = 1 +md.solidearth.elastic = 0 +md.solidearth.rotation = 0 md = solve(md, 'Sealevelrise') Srigid = md.results.SealevelriseSolution.Sealevel #eustatic + rigid + elastic run: -md.slr.rigid = 1 -md.slr.elastic = 1 -md.slr.rotation = 0 +md.solidearth.rigid = 1 +md.solidearth.elastic = 1 +md.solidearth.rotation = 0 md = solve(md, 'Sealevelrise') Selastic = md.results.SealevelriseSolution.Sealevel #eustatic + rigid + elastic + rotation run: -md.slr.rigid = 1 -md.slr.elastic = 1 -md.slr.rotation = 1 +md.solidearth.rigid = 1 +md.solidearth.elastic = 1 +md.solidearth.rotation = 1 md = solve(md, 'Sealevelrise') Srotation = md.results.SealevelriseSolution.Sealevel Index: /issm/trunk-jpl/test/NightlyRun/test2003.py =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test2003.py (revision 25124) +++ /issm/trunk-jpl/test/NightlyRun/test2003.py (revision 25125) @@ -1,3 +1,3 @@ -#Test Name: EarthSlr_rotationalFeedback +#Test Name: EarthSolidearth_rotationalFeedback import numpy as np from model import * @@ -13,8 +13,8 @@ md.mesh = gmshplanet('radius', 6.371012 * 1e3, 'resolution', 1000.) #500 km resolution mesh -#parameterize slr solution: -#slr loading: {{{ -md.slr.deltathickness = np.zeros((md.mesh.numberofelements, )) -md.slr.sealevel = np.zeros((md.mesh.numberofvertices, )) +#parameterize solidearth solution: +#solidearth loading: {{{ +md.solidearth.deltathickness = np.zeros((md.mesh.numberofelements, )) +md.solidearth.sealevel = np.zeros((md.mesh.numberofvertices, )) md.dsl.global_average_thermosteric_sea_level_change=np.zeros((2, )) md.dsl.sea_surface_height_change_above_geoid=np.zeros((md.mesh.numberofvertices+1, )) @@ -28,14 +28,14 @@ longe = sum(md.mesh.long[md.mesh.elements - 1], 1) / 3 pos = np.intersect1d(np.array(np.where(late < -75)), np.array(np.where(longe < 0))) -md.slr.deltathickness[pos] = -1 +md.solidearth.deltathickness[pos] = -1 #elastic loading from love numbers: nlov = 1000 -md.slr.love_h = np.array(love_numbers('h')) -md.slr.love_h = np.resize(md.slr.love_h, nlov + 1) -md.slr.love_k = np.array(love_numbers('k')) -md.slr.love_k = np.resize(md.slr.love_k, nlov + 1) -md.slr.love_l = np.array(love_numbers('l')) -md.slr.love_l = np.resize(md.slr.love_l, nlov + 1) +md.solidearth.love_h = np.array(love_numbers('h')) +md.solidearth.love_h = np.resize(md.solidearth.love_h, nlov + 1) +md.solidearth.love_k = np.array(love_numbers('k')) +md.solidearth.love_k = np.resize(md.solidearth.love_k, nlov + 1) +md.solidearth.love_l = np.array(love_numbers('l')) +md.solidearth.love_l = np.resize(md.solidearth.love_l, nlov + 1) #}}} @@ -56,5 +56,5 @@ #make sure that the elements that have loads are fully grounded: -pos = np.nonzero(md.slr.deltathickness)[0] +pos = np.nonzero(md.solidearth.deltathickness)[0] md.mask.ocean_levelset[md.mesh.elements[pos, :] - 1] = 1 @@ -65,5 +65,5 @@ # use model representation of ocea area (not the ture area) -md.slr.ocean_area_scaling = 0 +md.solidearth.ocean_area_scaling = 0 #geometry @@ -83,18 +83,18 @@ #New stuff -md.slr.spcthickness = np.nan * np.ones((md.mesh.numberofvertices, )) -md.slr.Ngia = np.zeros((md.mesh.numberofvertices, )) -md.slr.Ugia = np.zeros((md.mesh.numberofvertices, )) -md.slr.hydro_rate = np.zeros((md.mesh.numberofvertices, )) +md.solidearth.spcthickness = np.nan * np.ones((md.mesh.numberofvertices, )) +md.solidearth.Ngia = np.zeros((md.mesh.numberofvertices, )) +md.solidearth.Ugia = np.zeros((md.mesh.numberofvertices, )) +md.solidearth.hydro_rate = np.zeros((md.mesh.numberofvertices, )) #Solution parameters -md.slr.reltol = float('NaN') -md.slr.abstol = 1e-3 -md.slr.geodetic = 1 +md.solidearth.reltol = float('NaN') +md.solidearth.abstol = 1e-3 +md.solidearth.geodetic = 1 #eustatic + rigid + elastic run: -md.slr.rigid = 1 -md.slr.elastic = 1 -md.slr.rotation = 0 +md.solidearth.rigid = 1 +md.solidearth.elastic = 1 +md.solidearth.rotation = 0 md.cluster = generic('name', gethostname(), 'np', 3) #md.verbose = verbose('111111111') @@ -108,7 +108,7 @@ #eustatic + rigid + elastic + rotation run: -md.slr.rigid = 1 -md.slr.elastic = 1 -md.slr.rotation = 1 +md.solidearth.rigid = 1 +md.solidearth.elastic = 1 +md.solidearth.rotation = 1 md.cluster = generic('name', gethostname(), 'np', 3) #md.verbose = verbose('111111111') Index: /issm/trunk-jpl/test/NightlyRun/test2004.m =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test2004.m (revision 25124) +++ /issm/trunk-jpl/test/NightlyRun/test2004.m (revision 25125) @@ -1,19 +1,12 @@ %Test Name: Sea-Level-Partitions +testagainst2002=0; + %Data paths: {{{ -modeldatapath='/Users/larour/ModelData'; shppath='../Data/shp/'; gshhsshapefile=[shppath 'GSHHS_c_L1-NightlyRun.shp']; %}}} - -testagainst2002=0; - -%Data paths: {{{ -modeldatapath='/Users/larour/ModelData'; -shppath='../Data/shp/'; -gshhsshapefile=[shppath 'GSHHS_c_L1-NightlyRun.shp']; -%}}} %create sealevel model to hold our information: - sl=sealevelmodel(); +sl=sealevelmodel(); %Create basins using boundaries from shapefile: @@ -33,22 +26,22 @@ boundary('shppath',shppath,'shpfilename','SouthAntarctica','proj',proj3031)... })); - %}}} - %Ross: {{{ - sl.addbasin(basin('continent','antarctica','name','ross','proj',proj3031,'boundaries',{... - boundary('shppath',shppath,'shpfilename','SouthAntarctica','proj',proj3031),... - boundary('shppath',shppath,'shpfilename','RossIceShelf','proj',proj3031),... - boundary('shppath',shppath,'shpfilename','RossWestFront','proj',proj3031),... - boundary('shppath',shppath,'shpfilename','RossFront','proj',proj3031,'orientation','reverse')... - })); -%}}} - %HemisphereEast: {{{ - sl.addbasin(basin('continent','hemisphereeast','name','hemisphereeast','proj',laea(0,+90),'boundaries',{... %Australian projection lat,long - boundary('shppath',shppath,'shpfilename','HemisphereSplit','proj',proj4326),... - boundary('shppath',shppath,'shpfilename','SouthAntarctica','proj',proj3031),... - boundary('shppath',shppath,'shpfilename','RossFront','proj',proj3031),... - boundary('shppath',shppath,'shpfilename','RossWestFront','proj',proj3031),... - boundary('shppath',shppath,'shpfilename','EastAntarctica2','proj',proj3031,'orientation','reverse'),... - boundary('shppath',shppath,'shpfilename','NorthAntarctica','proj',proj3031)... - })); +%}}} +%Ross: {{{ +sl.addbasin(basin('continent','antarctica','name','ross','proj',proj3031,'boundaries',{... + boundary('shppath',shppath,'shpfilename','SouthAntarctica','proj',proj3031),... + boundary('shppath',shppath,'shpfilename','RossIceShelf','proj',proj3031),... + boundary('shppath',shppath,'shpfilename','RossWestFront','proj',proj3031),... + boundary('shppath',shppath,'shpfilename','RossFront','proj',proj3031,'orientation','reverse')... + })); +%}}} +%HemisphereEast: {{{ +sl.addbasin(basin('continent','hemisphereeast','name','hemisphereeast','proj',laea(0,+90),'boundaries',{... %Australian projection lat,long + boundary('shppath',shppath,'shpfilename','HemisphereSplit','proj',proj4326),... + boundary('shppath',shppath,'shpfilename','SouthAntarctica','proj',proj3031),... + boundary('shppath',shppath,'shpfilename','RossFront','proj',proj3031),... + boundary('shppath',shppath,'shpfilename','RossWestFront','proj',proj3031),... + boundary('shppath',shppath,'shpfilename','EastAntarctica2','proj',proj3031,'orientation','reverse'),... + boundary('shppath',shppath,'shpfilename','NorthAntarctica','proj',proj3031)... + })); %}}} %Antarctica excluding Ronne: {{{ @@ -59,18 +52,18 @@ boundary('shppath',shppath,'shpfilename','RossIceShelf','proj',proj3031,'orientation','reverse'),... boundary('shppath',shppath,'shpfilename','SouthAntarctica','proj',proj3031),... - boundary('shppath',shppath,'shpfilename','WestAntarctica2','proj',proj3031)... - boundary('shppath',shppath,'shpfilename','RonneWestSummit','proj',proj3031)... - boundary('shppath',shppath,'shpfilename','RonneIceShelf','proj',proj3031)... - boundary('shppath',shppath,'shpfilename','RonneEastSummit','proj',proj3031)... + boundary('shppath',shppath,'shpfilename','WestAntarctica2','proj',proj3031),... + boundary('shppath',shppath,'shpfilename','RonneWestSummit','proj',proj3031),... + boundary('shppath',shppath,'shpfilename','RonneIceShelf','proj',proj3031),... + boundary('shppath',shppath,'shpfilename','RonneEastSummit','proj',proj3031),... boundary('shppath',shppath,'shpfilename','RonneBrunt','proj',proj3031)... })); - %}}} - %Ronne: {{{ - sl.addbasin(basin('continent','antarctica','name','ronne','proj',proj3031,'boundaries',{... - boundary('shppath',shppath,'shpfilename','RonneWestSummit','proj',proj3031),... - boundary('shppath',shppath,'shpfilename','RonneIceShelf','proj',proj3031),... - boundary('shppath',shppath,'shpfilename','RonneEastSummit','proj',proj3031),... - boundary('shppath',shppath,'shpfilename','RonneFront','proj',proj3031,'orientation','reverse')... - })); +%}}} +%Ronne: {{{ +sl.addbasin(basin('continent','antarctica','name','ronne','proj',proj3031,'boundaries',{... + boundary('shppath',shppath,'shpfilename','RonneWestSummit','proj',proj3031),... + boundary('shppath',shppath,'shpfilename','RonneIceShelf','proj',proj3031),... + boundary('shppath',shppath,'shpfilename','RonneEastSummit','proj',proj3031),... + boundary('shppath',shppath,'shpfilename','RonneFront','proj',proj3031,'orientation','reverse')... + })); %}}} @@ -372,5 +365,5 @@ md.timestepping.time_step=1; -% max number of iteration reverted back to 10 (i.e., the original default value) +%max number of iterations reverted back to 10 (i.e. the original default value) md.solidearth.settings.maxiter=10; Index: /issm/trunk-jpl/test/NightlyRun/test2010.py =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test2010.py (revision 25124) +++ /issm/trunk-jpl/test/NightlyRun/test2010.py (revision 25125) @@ -14,29 +14,29 @@ md.mesh = gmshplanet('radius', rad_e, 'resolution', 1000.0) # km resolution -#parameterize slr solution: -#slr loading: {{{ +#parameterize solidearth solution: +#solidearth loading: {{{ late = sum(md.mesh.lat[md.mesh.elements - 1], 1) / 3 longe = sum(md.mesh.long[md.mesh.elements - 1], 1) / 3 -md.slr.deltathickness = np.zeros((md.mesh.numberofelements, )) +md.solidearth.deltathickness = np.zeros((md.mesh.numberofelements, )) pos = np.intersect1d(np.array(np.where(late < -75)), np.array(np.where(longe > 0))) #python does not include last element in array slices, (6:7) -> [5:7] -md.slr.deltathickness[pos[5:7]] = -1 +md.solidearth.deltathickness[pos[5:7]] = -1 -md.slr.sealevel = np.zeros((md.mesh.numberofvertices, )) +md.solidearth.sealevel = np.zeros((md.mesh.numberofvertices, )) md.dsl.global_average_thermosteric_sea_level_change=np.zeros((2, )) md.dsl.sea_surface_height_change_above_geoid=np.zeros((md.mesh.numberofvertices+1, )) md.dsl.sea_water_pressure_change_at_sea_floor=np.zeros((md.mesh.numberofvertices+1, )) -md.slr.ocean_area_scaling = 1 +md.solidearth.ocean_area_scaling = 1 #elastic loading from love numbers: nlov = 1000 -md.slr.love_h = np.array(love_numbers('h')) -md.slr.love_h = np.resize(md.slr.love_h, nlov + 1) -md.slr.love_k = np.array(love_numbers('k')) -md.slr.love_k = np.resize(md.slr.love_k, nlov + 1) -md.slr.love_l = np.array(love_numbers('l')) -md.slr.love_l = np.resize(md.slr.love_l, nlov + 1) +md.solidearth.love_h = np.array(love_numbers('h')) +md.solidearth.love_h = np.resize(md.solidearth.love_h, nlov + 1) +md.solidearth.love_k = np.array(love_numbers('k')) +md.solidearth.love_k = np.resize(md.solidearth.love_k, nlov + 1) +md.solidearth.love_l = np.array(love_numbers('l')) +md.solidearth.love_l = np.resize(md.solidearth.love_l, nlov + 1) #}}} @@ -56,5 +56,5 @@ #make sure that the elements that have loads are fully grounded: -pos = np.nonzero(md.slr.deltathickness)[0] +pos = np.nonzero(md.solidearth.deltathickness)[0] md.mask.ocean_levelset[md.mesh.elements[pos, :] - 1] = 1 @@ -79,19 +79,19 @@ # }}} #Solution parameters {{{ -md.slr.reltol = float('NaN') -md.slr.abstol = 1e-3 -md.slr.geodetic = 1 +md.solidearth.reltol = float('NaN') +md.solidearth.abstol = 1e-3 +md.solidearth.geodetic = 1 # }}} #New stuff -md.slr.spcthickness = np.nan * np.ones((md.mesh.numberofvertices, )) -md.slr.Ngia = np.zeros((md.mesh.numberofvertices, )) -md.slr.Ugia = np.zeros((md.mesh.numberofvertices, )) -md.slr.hydro_rate = np.zeros((md.mesh.numberofvertices, )) +md.solidearth.spcthickness = np.nan * np.ones((md.mesh.numberofvertices, )) +md.solidearth.Ngia = np.zeros((md.mesh.numberofvertices, )) +md.solidearth.Ugia = np.zeros((md.mesh.numberofvertices, )) +md.solidearth.hydro_rate = np.zeros((md.mesh.numberofvertices, )) #eustatic + rigid + elastic run: -md.slr.rigid = 1 -md.slr.elastic = 1 -md.slr.rotation = 1 +md.solidearth.rigid = 1 +md.solidearth.elastic = 1 +md.solidearth.rotation = 1 md.cluster = generic('name', gethostname(), 'np', 3) @@ -99,5 +99,5 @@ md = solve(md, 'Sealevelrise') eus = md.results.SealevelriseSolution.SealevelRSLEustatic -slr = md.results.SealevelriseSolution.Sealevel +solidearth = md.results.SealevelriseSolution.Sealevel moixz = md.results.SealevelriseSolution.SealevelInertiaTensorXZ moiyz = md.results.SealevelriseSolution.SealevelInertiaTensorYZ @@ -105,14 +105,14 @@ # analytical moi = > just checking FOR ICE only!!! {{{ -# ...have to mute**slr induced MOI in Tria.cpp**prior to the comparison +# ...have to mute**solidearth induced MOI in Tria.cpp**prior to the comparison #rad_e = rad_e * 1e3 # now in meters #areas = GetAreasSphericalTria(md.mesh.elements, md.mesh.lat, md.mesh.long, rad_e) #lat = late * pi / 180 lon = longe * pi / 180 -#moi_xz = sum(-md.materials.rho_freshwater. * md.slr.deltathickness. * areas. * rad_e^2. * sin(lat). * cos(lat). * cos(lon)) -#moi_yz = sum(-md.materials.rho_freshwater. * md.slr.deltathickness. * areas. * rad_e^2. * sin(lat). * cos(lat). * sin(lon)) +#moi_xz = sum(-md.materials.rho_freshwater. * md.solidearth.deltathickness. * areas. * rad_e^2. * sin(lat). * cos(lat). * cos(lon)) +#moi_yz = sum(-md.materials.rho_freshwater. * md.solidearth.deltathickness. * areas. * rad_e^2. * sin(lat). * cos(lat). * sin(lon)) # }}} #Fields and tolerances to track changes -field_names = ['eus', 'slr', 'moixz', 'moiyz', 'moizz'] +field_names = ['eus', 'solidearth', 'moixz', 'moiyz', 'moizz'] field_tolerances = [1e-13, 1e-13, 1e-13, 1e-13, 1e-13] -field_values = [eus, slr, moixz, moiyz, moizz] +field_values = [eus, solidearth, moixz, moiyz, moizz] Index: /issm/trunk-jpl/test/NightlyRun/test2424.py =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test2424.py (revision 25124) +++ /issm/trunk-jpl/test/NightlyRun/test2424.py (revision 25125) @@ -30,5 +30,5 @@ md.timestepping.time_step = .1 -md.slr.sealevel = newforcing(md.timestepping.start_time, md.timestepping.final_time, +md.solidearth.sealevel = newforcing(md.timestepping.start_time, md.timestepping.final_time, md.timestepping.time_step, -200., 200., md.mesh.numberofvertices) Index: /issm/trunk-jpl/test/NightlyRun/test2425.py =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test2425.py (revision 25124) +++ /issm/trunk-jpl/test/NightlyRun/test2425.py (revision 25125) @@ -39,5 +39,5 @@ md.geometry.bed = md.geometry.bed + 1000 md.geometry.surface = md.geometry.surface + 1000 -md.slr.sealevel = 1000 * np.ones((md.mesh.numberofvertices, )) +md.solidearth.sealevel = 1000 * np.ones((md.mesh.numberofvertices, )) md = solve(md, 'Transient', 'checkconsistency', 'no')