Index: /issm/trunk-jpl/src/m/classes/SMBgemb.py =================================================================== --- /issm/trunk-jpl/src/m/classes/SMBgemb.py (revision 25306) +++ /issm/trunk-jpl/src/m/classes/SMBgemb.py (revision 25307) @@ -1,7 +1,8 @@ import numpy as np + +from checkfield import checkfield from fielddisplay import fielddisplay -from checkfield import checkfield +from project3d import project3d from WriteData import WriteData -from project3d import project3d @@ -14,5 +15,5 @@ """ - def __init__(self): # {{{ + def __init__(self, *args): # {{{ #each one of these properties is a transient forcing to the GEMB model, loaded from meteorological data derived #from an automatic weather stations (AWS). Each property is therefore a matrix, of size (numberofvertices x number @@ -32,69 +33,69 @@ #inputs: - self.Ta = float('NaN') #2 m air temperature, in Kelvin - self.V = float('NaN') #wind speed (m/s-1) - self.dswrf = float('NaN') #downward shortwave radiation flux [W/m^2] - self.dlwrf = float('NaN') #downward longwave radiation flux [W/m^2] - self.P = float('NaN') #precipitation [mm w.e. / m^2] - self.eAir = float('NaN') #screen level vapor pressure [Pa] - self.pAir = float('NaN') #surface pressure [Pa] - self.Tmean = float('NaN') #mean annual temperature [K] - self.Vmean = float('NaN') #mean annual wind velocity [m s-1] - self.C = float('NaN') #mean annual snow accumulation [kg m-2 yr-1] - self.Tz = float('NaN') #height above ground at which temperature (T) was sampled [m] - self.Vz = float('NaN') #height above ground at which wind (V) eas sampled [m] + self.Ta = np.nan # 2 m air temperature, in Kelvin + self.V = np.nan # wind speed (m/s-1) + self.dswrf = np.nan # downward shortwave radiation flux [W/m^2] + self.dlwrf = np.nan # downward longwave radiation flux [W/m^2] + self.P = np.nan # precipitation [mm w.e. / m^2] + self.eAir = np.nan # screen level vapor pressure [Pa] + self.pAir = np.nan # surface pressure [Pa] + self.Tmean = np.nan # mean annual temperature [K] + self.Vmean = np.nan # mean annual wind velocity [m s-1] + self.C = np.nan # mean annual snow accumulation [kg m-2 yr-1] + self.Tz = np.nan # height above ground at which temperature (T) was sampled [m] + self.Vz = np.nan # height above ground at which wind (V) was sampled [m] #optional inputs: - self.aValue = float('NaN') #Albedo forcing at every element. Used only if aIdx == 0. - self.teValue = float('NaN') #Outward longwave radiation thermal emissivity forcing at every element (default in code is 1) + self.aValue = np.nan # Albedo forcing at every element. Used only if aIdx == 0. + self.teValue = np.nan # Outward longwave radiation thermal emissivity forcing at every element (default in code is 1) # Initialization of snow properties - self.Dzini = float('NaN') #cell depth (m) - self.Dini = float('NaN') #snow density (kg m-3) - self.Reini = float('NaN') #effective grain size (mm) - self.Gdnini = float('NaN') #grain dricity (0-1) - self.Gspini = float('NaN') #grain sphericity (0-1) - self.ECini = float('NaN') #evaporation/condensation (kg m-2) - self.Wini = float('NaN') #Water content (kg m-2) - self.Aini = float('NaN') #albedo (0-1) - self.Tini = float('NaN') #snow temperature (K) - self.Sizeini = float('NaN') #Number of layers + self.Dzini = np.nan # cell depth (m) + self.Dini = np.nan # snow density (kg m-3) + self.Reini = np.nan # effective grain size (mm) + self.Gdnini = np.nan # grain dricity (0-1) + self.Gspini = np.nan # grain sphericity (0-1) + self.ECini = np.nan # evaporation/condensation (kg m-2) + self.Wini = np.nan # Water content (kg m-2) + self.Aini = np.nan # albedo (0-1) + self.Tini = np.nan # snow temperature (K) + self.Sizeini = np.nan # Number of layers #settings: - self.aIdx = float('NaN') #method for calculating albedo and subsurface absorption (default is 1) - self.swIdx = float('NaN') # apply all SW to top grid cell (0) or allow SW to penetrate surface (1) (default 1) - self.denIdx = float('NaN') #densification model to use (default is 2): - self.dsnowIdx = float('NaN') #model for fresh snow accumulation density (default is 1): - self.zTop = float('NaN') # depth over which grid length is constant at the top of the snopack (default 10) [m] - self.dzTop = float('NaN') # initial top vertical grid spacing (default .05) [m] - self.dzMin = float('NaN') # initial min vertical allowable grid spacing (default dzMin/2) [m] - self.zY = float('NaN') # strech grid cells bellow top_z by a [top_dz * y ^ (cells bellow top_z)] - self.zMax = float('NaN') #initial max model depth (default is min(thickness, 250)) [m] - self.zMin = float('NaN') #initial min model depth (default is min(thickness, 130)) [m] - self.outputFreq = float('NaN') #output frequency in days (default is monthly, 30) + self.aIdx = np.nan # method for calculating albedo and subsurface absorption (default is 1) + self.swIdx = np.nan # apply all SW to top grid cell (0) or allow SW to penetrate surface (1) (default 1) + self.denIdx = np.nan # densification model to use (default is 2): + self.dsnowIdx = np.nan # model for fresh snow accumulation density (default is 1): + self.zTop = np.nan # depth over which grid length is constant at the top of the snopack (default 10) [m] + self.dzTop = np.nan # initial top vertical grid spacing (default .05) [m] + self.dzMin = np.nan # initial min vertical allowable grid spacing (default dzMin/2) [m] + self.zY = np.nan # strech grid cells bellow top_z by a [top_dz * y ^ (cells bellow top_z)] + self.zMax = np.nan # initial max model depth (default is min(thickness, 250)) [m] + self.zMin = np.nan # initial min model depth (default is min(thickness, 130)) [m] + self.outputFreq = np.nan # output frequency in days (default is monthly, 30) #specific albedo parameters: #Method 1 and 2: - self.aSnow = float('NaN') # new snow albedo (0.64 - 0.89) - self.aIce = float('NaN') # range 0.27-0.58 for old snow + self.aSnow = np.nan # new snow albedo (0.64 - 0.89) + self.aIce = np.nan # range 0.27-0.58 for old snow #Method 3: Radiation Correction Factors -> only used for met station data and Greuell & Konzelmann, 1994 albedo - self.cldFrac = float('NaN') # average cloud amount + self.cldFrac = np.nan # average cloud amount #Method 4: additonal tuning parameters albedo as a funtion of age and water content (Bougamont et al., 2005) - self.t0wet = float('NaN') # time scale for wet snow (15-21.9) - self.t0dry = float('NaN') # warm snow timescale (30) - self.K = float('NaN') # time scale temperature coef. (7) - self.adThresh = float('NaN') # Apply aIdx method to all areas with densities below this value, + self.t0wet = np.nan # time scale for wet snow (15-21.9) + self.t0dry = np.nan # warm snow timescale (30) + self.K = np.nan # time scale temperature coef. (7) + self.adThresh = np.nan # Apply aIdx method to all areas with densities below this value, # or else apply direct input value from aValue, allowing albedo to be altered. # Default value is rho water (1023 kg m-3). #densities: - self.InitDensityScaling = float('NaN') #initial scaling factor multiplying the density of ice, which describes the density of the snowpack. + self.InitDensityScaling = np.nan # initial scaling factor multiplying the density of ice, which describes the density of the snowpack. #thermo: - self.ThermoDeltaTScaling = float('NaN') #scaling factor to multiply the thermal diffusion timestep (delta t) - - self.steps_per_step = 1 - self.averaging = 0 - self.requested_outputs = [] + self.ThermoDeltaTScaling = np.nan # scaling factor to multiply the thermal diffusion timestep (delta t) + + self.steps_per_step = 1 + self.averaging = 0 + self.requested_outputs = [] #Several fields are missing from the standard GEMB model, which are capture intrinsically by ISSM. @@ -103,4 +104,12 @@ #elev: this is taken from the ISSM surface itself. + nargin = len(args) + if nargin: + if nargin == 2: + mesh = args[0] + geometry = args[1] + self.setdefaultparameters(mesh, geometry) + else: + raise Exception('constructor not supported: need mesh and geometry to set defaults') #}}} @@ -281,5 +290,5 @@ md = checkfield(md, 'fieldname', 'smb.isgraingrowth', 'values', [0, 1]) - #md = checkfield(md, 'fieldname', 'smb.issmbgradients', 'values', [0, 1]) + md = checkfield(md, 'fieldname', 'smb.issmbgradients', 'values', [0, 1]) md = checkfield(md, 'fieldname', 'smb.isalbedo', 'values', [0, 1]) md = checkfield(md, 'fieldname', 'smb.isshortwave', 'values', [0, 1]) @@ -352,5 +361,4 @@ def marshall(self, prefix, md, fid): # {{{ - yts = md.constants.yts Index: /issm/trunk-jpl/src/m/solve/WriteData.m =================================================================== --- /issm/trunk-jpl/src/m/solve/WriteData.m (revision 25306) +++ /issm/trunk-jpl/src/m/solve/WriteData.m (revision 25307) @@ -32,5 +32,5 @@ end -%Scale data if necesarry +%Scale data if necessary if strcmpi(format,'MatArray') for i=1:numel(data) Index: /issm/trunk-jpl/src/m/solve/WriteData.py =================================================================== --- /issm/trunk-jpl/src/m/solve/WriteData.py (revision 25306) +++ /issm/trunk-jpl/src/m/solve/WriteData.py (revision 25307) @@ -1,5 +1,7 @@ from copy import deepcopy -from struct import pack, error +from struct import error, pack + import numpy as np + import pairoptions @@ -47,21 +49,28 @@ # (see https://ross.ics.uci.edu/jenkins/view/All/job/Debian_Linux-Python/1036). # - if mattype != 0: - data = deepcopy(data) - - #Scale data if necesarry - if options.exist('scale'): - data = np.array(data) - scale = options.getfieldvalue('scale') - if np.size(data) > 1 and np.ndim(data) > 1 and np.size(data, 0) == timeserieslength: - data[0:-1, :] = scale * data[0:-1, :] - else: - data = scale * data - if np.size(data) > 1 and np.size(data, 0) == timeserieslength: - yts = options.getfieldvalue('yts') - if np.ndim(data) > 1: + # data = deepcopy(data) + + #Scale data if necessary + if datatype == 'MatArray': + for i in range(np.size(data)): + if options.exist('scale'): + scale = options.getfieldvalue('scale') + if np.ndim(data[i]) > 1 and data[i].shape[0] == timeserieslength: + data[i][:-2, :] = scale * data[i][:-2, :] + else: + data[i] = scale * data[i] + if np.ndim(data) > 1 and data[i].shape[0] == timeserieslength: + yts = options.getfieldvalue('yts') + data[i][-1, :] = yts * data[i][-1, :] + else: + if options.exist('scale'): + scale = options.getfieldvalue('scale') + if np.ndim(data) > 1 and data.shape[0] == timeserieslength: + data[:-2, :] = scale * data[:-2, :] + else: + data = scale * data + if np.ndim(data) > 1 and data.shape[0] == timeserieslength: + yts = options.getfieldvalue('yts') data[-1, :] = yts * data[-1, :] - else: - data[-1] = yts * data[-1] #Step 1: write the enum to identify this record uniquely Index: /issm/trunk-jpl/test/NightlyRun/test234.py =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test234.py (revision 25306) +++ /issm/trunk-jpl/test/NightlyRun/test234.py (revision 25307) @@ -10,5 +10,4 @@ from setflowequation import * from solve import * -from SMBgemb import * from partitioner import * from dmeth_params_set import * @@ -39,5 +38,5 @@ #partitioning npart = 20 -partition = partitioner(md, 'package', 'chaco', 'npart', npart, 'weighting', 'on') - 1 +partition = partitioner(md, 'package', 'chaco', 'npart', npart, 'weighting', 'on')[0] - 1 #variables Index: /issm/trunk-jpl/test/NightlyRun/test235.py =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test235.py (revision 25306) +++ /issm/trunk-jpl/test/NightlyRun/test235.py (revision 25307) @@ -10,5 +10,4 @@ from setflowequation import * from solve import * -from SMBgemb import * from partitioner import * from dmeth_params_set import * Index: /issm/trunk-jpl/test/NightlyRun/test243.py =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test243.py (revision 25306) +++ /issm/trunk-jpl/test/NightlyRun/test243.py (revision 25307) @@ -1,13 +1,15 @@ #Test Name: SquareShelfSMBGemb +from socket import gethostname +import sys + import numpy as np -import sys + from model import * -from socket import gethostname -from triangle import * -from setmask import * from parameterize import * from setflowequation import * +from setmask import * +from SMBgemb import * from solve import * -from SMBgemb import * +from triangle import * md = triangle(model(), '../Exp/Square.exp', 350000.) @@ -19,6 +21,5 @@ #Use of Gemb method for SMB computation -md.smb = SMBgemb() -md.smb.setdefaultparameters(md.mesh, md.geometry) +md.smb = SMBgemb(md.mesh, md.geometry) md.smb.dsnowIdx = 1 @@ -37,5 +38,4 @@ md.smb.eAir = np.append(np.tile(np.conjugate(inputs[b'eAir0']), (md.mesh.numberofelements, 1)), np.conjugate([inputs[b'dateN']]), axis=0) md.smb.pAir = np.append(np.tile(np.conjugate(inputs[b'pAir0']), (md.mesh.numberofelements, 1)), np.conjugate([inputs[b'dateN']]), axis=0) -md.smb.pAir = np.append(np.tile(np.conjugate(inputs[b'pAir0']), (md.mesh.numberofelements, 1)), np.conjugate([inputs[b'dateN']]), axis=0) md.smb.Vz = np.tile(np.conjugate(inputs[b'LP']['Vz']), (md.mesh.numberofelements, 1)).flatten() md.smb.Tz = np.tile(np.conjugate(inputs[b'LP']['Tz']), (md.mesh.numberofelements, 1)).flatten() @@ -44,5 +44,9 @@ #smb settings -md.smb.requested_outputs = ['SmbDz', 'SmbT', 'SmbD', 'SmbRe', 'SmbGdn', 'SmbGsp', 'SmbEC', 'SmbA', 'SmbMassBalance', 'SmbMAdd', 'SmbDzAdd', 'SmbFAC', 'SmbMeanSHF', 'SmbMeanLHF', 'SmbMeanULW', 'SmbNetLW', 'SmbNetSW'] +md.smb.requested_outputs = [ + 'SmbDz', 'SmbT', 'SmbD', 'SmbRe', 'SmbGdn', 'SmbGsp', 'SmbEC', + 'SmbA', 'SmbMassBalance', 'SmbMAdd', 'SmbDzAdd', 'SmbFAC', 'SmbMeanSHF', 'SmbMeanLHF', + 'SmbMeanULW', 'SmbNetLW', 'SmbNetSW' + ] #only run smb core: @@ -54,5 +58,5 @@ md.timestepping.start_time = 1965. md.timestepping.final_time = 1966. -md.timestepping.time_step = 1. / 365.0 +md.timestepping.time_step = 1.0 / 365 md.timestepping.interp_forcings = 0. @@ -68,20 +72,22 @@ field_tolerances = [1e-12, 2e-12, 1e-12, 1e-11, 1e-11, 2e-11, 1e-11, 1e-12, 1e-11, 1e-12, 1e-12, 1e-12, 1e-11, 2e-11, 2e-11, 1e-11, 9e-10, 2e-11] #shape is different in python solution (fixed using reshape) which can cause test failure: -field_values = [nlayers, - md.results.TransientSolution[-1].SmbDz[0, 0:nlayers].reshape(1, -1), - md.results.TransientSolution[-1].SmbT[0, 0:nlayers].reshape(1, -1), - md.results.TransientSolution[-1].SmbD[0, 0:nlayers].reshape(1, -1), - md.results.TransientSolution[-1].SmbRe[0, 0:nlayers].reshape(1, -1), - md.results.TransientSolution[-1].SmbGdn[0, 0:nlayers].reshape(1, -1), - md.results.TransientSolution[-1].SmbGsp[0, 0:nlayers].reshape(1, -1), - md.results.TransientSolution[-1].SmbA[0, 0:nlayers].reshape(1, -1), - md.results.TransientSolution[-1].SmbEC[0], - md.results.TransientSolution[-1].SmbMassBalance[0], - md.results.TransientSolution[-1].SmbMAdd[0], - md.results.TransientSolution[-1].SmbDzAdd[0], - md.results.TransientSolution[-1].SmbFAC[0], - md.results.TransientSolution[-1].SmbMeanSHF[0], - md.results.TransientSolution[-1].SmbMeanLHF[0], - md.results.TransientSolution[-1].SmbMeanULW[0], - md.results.TransientSolution[-1].SmbNetLW[0], - md.results.TransientSolution[-1].SmbNetSW[0]] +field_values = [ + nlayers, + md.results.TransientSolution[-1].SmbDz[0, 0:nlayers].reshape(1, -1), + md.results.TransientSolution[-1].SmbT[0, 0:nlayers].reshape(1, -1), + md.results.TransientSolution[-1].SmbD[0, 0:nlayers].reshape(1, -1), + md.results.TransientSolution[-1].SmbRe[0, 0:nlayers].reshape(1, -1), + md.results.TransientSolution[-1].SmbGdn[0, 0:nlayers].reshape(1, -1), + md.results.TransientSolution[-1].SmbGsp[0, 0:nlayers].reshape(1, -1), + md.results.TransientSolution[-1].SmbA[0, 0:nlayers].reshape(1, -1), + md.results.TransientSolution[-1].SmbEC[0], + md.results.TransientSolution[-1].SmbMassBalance[0], + md.results.TransientSolution[-1].SmbMAdd[0], + md.results.TransientSolution[-1].SmbDzAdd[0], + md.results.TransientSolution[-1].SmbFAC[0], + md.results.TransientSolution[-1].SmbMeanSHF[0], + md.results.TransientSolution[-1].SmbMeanLHF[0], + md.results.TransientSolution[-1].SmbMeanULW[0], + md.results.TransientSolution[-1].SmbNetLW[0], + md.results.TransientSolution[-1].SmbNetSW[0] + ] Index: /issm/trunk-jpl/test/NightlyRun/test244.py =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test244.py (revision 25306) +++ /issm/trunk-jpl/test/NightlyRun/test244.py (revision 25307) @@ -28,6 +28,5 @@ # Use of Gemb method for SMB computation -md.smb = SMBgemb() -md.smb.setdefaultparameters(md.mesh, md.geometry) +md.smb = SMBgemb(md.mesh, md.geometry) md.smb.dsnowIdx = 0 Index: /issm/trunk-jpl/test/NightlyRun/test250.py =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test250.py (revision 25306) +++ /issm/trunk-jpl/test/NightlyRun/test250.py (revision 25307) @@ -9,5 +9,4 @@ from setflowequation import * from solve import * -from SMBgemb import * from partitioner import * from dmeth_params_set import * Index: /issm/trunk-jpl/test/NightlyRun/test251.py =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test251.py (revision 25306) +++ /issm/trunk-jpl/test/NightlyRun/test251.py (revision 25307) @@ -9,5 +9,4 @@ from setflowequation import * from solve import * -from SMBgemb import * from partitioner import * from dmeth_params_set import * Index: /issm/trunk-jpl/test/NightlyRun/test252.py =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test252.py (revision 25306) +++ /issm/trunk-jpl/test/NightlyRun/test252.py (revision 25307) @@ -1,13 +1,15 @@ #Test Name: SquareShelfSMBGembClim +from socket import gethostname +import sys + import numpy as np -import sys + from model import * -from socket import gethostname -from triangle import * -from setmask import * from parameterize import * from setflowequation import * +from setmask import * +from SMBgemb import * from solve import * -from SMBgemb import * +from triangle import * md = triangle(model(), '../Exp/Square.exp', 350000.) @@ -19,6 +21,5 @@ #Use of Gemb method for SMB computation -md.smb = SMBgemb() -md.smb.setdefaultparameters(md.mesh, md.geometry) +md.smb = SMBgemb(md.mesh, md.geometry) md.smb.dsnowIdx = 4 Index: /issm/trunk-jpl/test/NightlyRun/test253.py =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test253.py (revision 25306) +++ /issm/trunk-jpl/test/NightlyRun/test253.py (revision 25307) @@ -19,6 +19,5 @@ #Use of Gemb method for SMB computation -md.smb = SMBgemb() -md.smb.setdefaultparameters(md.mesh, md.geometry) +md.smb = SMBgemb(md.mesh, md.geometry) md.smb.dsnowIdx = 3 md.smb.aIdx = 2 Index: /issm/trunk-jpl/test/NightlyRun/test418.py =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test418.py (revision 25306) +++ /issm/trunk-jpl/test/NightlyRun/test418.py (revision 25307) @@ -37,6 +37,5 @@ # - Run valgrind and fix the above # -partition, md = partitioner(md, 'package', 'chaco', 'npart', npart) -partition -= 1 +partition = partitioner(md, 'package', 'chaco', 'npart', npart)[0] - 1 vector = np.arange(1, 1 + md.mesh.numberofvertices, 1).reshape(-1, 1) Index: /issm/trunk-jpl/test/NightlyRun/test420.py =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test420.py (revision 25306) +++ /issm/trunk-jpl/test/NightlyRun/test420.py (revision 25307) @@ -18,5 +18,5 @@ #partitioning npart = 10 -partition = partitioner(md, 'package', 'chaco', 'npart', npart) - 1 +partition = partitioner(md, 'package', 'chaco', 'npart', npart)[0] - 1 md.qmu.isdakota = 1 Index: /issm/trunk-jpl/test/NightlyRun/test444.py =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test444.py (revision 25306) +++ /issm/trunk-jpl/test/NightlyRun/test444.py (revision 25307) @@ -13,5 +13,4 @@ from setflowequation import * from solve import * -from SMBgemb import * from partitioner import * from dmeth_params_set import * Index: /issm/trunk-jpl/test/Par/SquareShelf.py =================================================================== --- /issm/trunk-jpl/test/Par/SquareShelf.py (revision 25306) +++ /issm/trunk-jpl/test/Par/SquareShelf.py (revision 25307) @@ -1,10 +1,12 @@ +from arch import * +import inspect import os.path -import inspect -from arch import * + import numpy as np -from verbose import verbose + from InterpFromMeshToMesh2d import InterpFromMeshToMesh2d from paterson import paterson from SetIceShelfBC import SetIceShelfBC +from verbose import verbose #Start defining model parameters here @@ -34,8 +36,15 @@ #dbg - end +# x = x[0][:] +# y = y[0][:] +# vx = vx[0][:] +# vy = vy[0][:] +# index = index[0][:] + [md.initialization.vx] = InterpFromMeshToMesh2d(index, x, y, vx, md.mesh.x, md.mesh.y) [md.initialization.vy] = InterpFromMeshToMesh2d(index, x, y, vy, md.mesh.x, md.mesh.y) -md.initialization.vz = np.zeros((md.mesh.numberofvertices)) -md.initialization.pressure = np.zeros((md.mesh.numberofvertices)) +x = y = vx = vy = index = None +md.initialization.vz = np.zeros((md.mesh.numberofvertices, 1)) +md.initialization.pressure = np.zeros((md.mesh.numberofvertices, 1)) #dbg - begin @@ -50,19 +59,18 @@ #dbg - end +# Materials +md.initialization.temperature = (273 - 20) * np.ones((md.mesh.numberofvertices, 1)) +md.materials.rheology_B = paterson(md.initialization.temperature) +md.materials.rheology_n = 3 * np.ones((md.mesh.numberofelements, 1)) -#Materials -md.initialization.temperature = (273. - 20.) * np.ones((md.mesh.numberofvertices)) -md.materials.rheology_B = paterson(md.initialization.temperature) -md.materials.rheology_n = 3. * np.ones((md.mesh.numberofelements)) +# Friction +md.friction.coefficient = 20 * np.ones((md.mesh.numberofvertices, 1)) +md.friction.coefficient[np.where(md.mask.ocean_levelset < 0)[0]] = 0 +md.friction.p = np.ones((md.mesh.numberofelements, 1)) +md.friction.q = np.ones((md.mesh.numberofelements, 1)) -#Friction -md.friction.coefficient = 20. * np.ones((md.mesh.numberofvertices)) -md.friction.coefficient[np.nonzero(md.mask.ocean_levelset < 0.)[0]] = 0. -md.friction.p = np.ones((md.mesh.numberofelements)) -md.friction.q = np.ones((md.mesh.numberofelements)) - -#Numerical parameters -md.masstransport.stabilization = 1. -md.thermal.stabilization = 1. +# Numerical parameters +md.masstransport.stabilization = 1 +md.thermal.stabilization = 1 md.settings.waitonlock = 30 md.verbose = verbose() @@ -70,11 +78,10 @@ md.steadystate.reltol = 0.02 md.stressbalance.reltol = 0.02 -md.stressbalance.abstol = float('nan') -md.timestepping.time_step = 1. -md.timestepping.final_time = 3. +md.stressbalance.abstol = np.nan +md.timestepping.time_step = 1 +md.timestepping.final_time = 3 md.groundingline.migration = 'None' -#Boundary conditions: -# #md = SetIceShelfBC(md) +# Boundary conditions: md = SetIceShelfBC(md, '../Exp/SquareFront.exp')