Index: ../trunk-jpl/test/Par/GiaIvinsBenchmarksAB.py =================================================================== --- ../trunk-jpl/test/Par/GiaIvinsBenchmarksAB.py (revision 25171) +++ ../trunk-jpl/test/Par/GiaIvinsBenchmarksAB.py (revision 25172) @@ -14,29 +14,29 @@ rad = 600000. nv = md.mesh.numberofvertices if (np.isnan(md.geometry.thickness)): - md.geometry.thickness = np.zeros((md.mesh.numberofvertices, )) + md.geometry.thickness = np.zeros((md.mesh.numberofvertices, 1)) for i in range(nv): dist = np.sqrt(md.mesh.x[i]**2 + md.mesh.y[i]**2) if (dist <= rad): md.geometry.thickness[i] = 2000.0 else: - md.geometry.thickness[i] = 1.0 # non - zero thickness + md.geometry.thickness[i] = 1.0 # non-zero thickness md.geometry.thickness = md.geometry.thickness.reshape(-1, 1) -md.geometry.base = np.zeros((md.mesh.numberofvertices, )) +md.geometry.base = np.zeros((md.mesh.numberofvertices, 1)) md.geometry.surface = md.geometry.thickness + md.geometry.base.reshape(-1, 1) #would otherwise create a 91x91 matrix -#Ice density used for benchmarking, not 917 kg / m^3 -md.materials.rho_ice = 1000 #kg m^3 +#Ice density used for benchmarking, not 917 kg/m^3 +md.materials.rho_ice = 1000 #kg/m^3 -#GIA parameters specific to Experiments A and B +#GIA parameters specific to Experiments A and B md.gia=giaivins(); -md.gia.mantle_viscosity = 1e21 * np.ones((md.mesh.numberofvertices, )) #in Pa.s -md.gia.lithosphere_thickness = 100 * np.ones((md.mesh.numberofvertices, )) #in km +md.gia.mantle_viscosity = 1e21 * np.ones((md.mesh.numberofvertices, 1)) #in Pa.s +md.gia.lithosphere_thickness = 100 * np.ones((md.mesh.numberofvertices, 1)) #in km md.materials.lithosphere_shear_modulus = 6.7 * 1e10 #in Pa -md.materials.lithosphere_density = 3.36 #in g / cm^3 +md.materials.lithosphere_density = 3.36 #in g/cm^3 md.materials.mantle_shear_modulus = 1.45 * 1e11 #in Pa -md.materials.mantle_density = 3.38 #in g / cm^3 +md.materials.mantle_density = 3.38 #in g/cm^3 #Initial velocity x = archread('../Data/SquareSheetConstrained.arch', 'x') @@ -52,19 +52,19 @@ x = None y = None index = None -md.initialization.vz = np.zeros((md.mesh.numberofvertices, )) -md.initialization.pressure = np.zeros((md.mesh.numberofvertices, )) +md.initialization.vz = np.zeros((md.mesh.numberofvertices, 1)) +md.initialization.pressure = np.zeros((md.mesh.numberofvertices, 1)) #Materials -md.initialization.temperature = (273. - 20.) * np.ones((md.mesh.numberofvertices, )) +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, )) +md.materials.rheology_n = 3. * np.ones((md.mesh.numberofelements, 1)) #Friction -md.friction.coefficient = 20. * np.ones((md.mesh.numberofvertices, )) +md.friction.coefficient = 20. * np.ones((md.mesh.numberofvertices, 1)) md.friction.coefficient[np.where(md.mask.ocean_levelset < 0.)] = 0. -md.friction.p = np.ones((md.mesh.numberofelements, )) -md.friction.q = np.ones((md.mesh.numberofelements, )) +md.friction.p = np.ones((md.mesh.numberofelements, 1)) +md.friction.q = np.ones((md.mesh.numberofelements, 1)) #Numerical parameters md.groundingline.migration = 'None' @@ -75,7 +75,7 @@ md.stressbalance.restol = 0.05 md.steadystate.reltol = 0.05 md.stressbalance.reltol = 0.05 -md.stressbalance.abstol = float('NaN') +md.stressbalance.abstol = np.nan md.timestepping.time_step = 1. md.timestepping.final_time = 3. Index: ../trunk-jpl/test/Par/GiaIvinsBenchmarksCD.py =================================================================== --- ../trunk-jpl/test/Par/GiaIvinsBenchmarksCD.py (revision 25171) +++ ../trunk-jpl/test/Par/GiaIvinsBenchmarksCD.py (revision 25172) @@ -1,18 +1,22 @@ #Geometry specific to Experiments C and D +import inspect import os.path -import inspect + import numpy as np + from arch import * +from giaivins import giaivins from InterpFromMeshToMesh2d import * from paterson import * from verbose import * from SetIceSheetBC import * + rad = 800000 nv = md.mesh.numberofvertices if (np.isnan(md.geometry.thickness)): - md.geometry.thickness = np.zeros((md.mesh.numberofvertices, )) + md.geometry.thickness = np.zeros((md.mesh.numberofvertices, 1)) for i in range(nv): dist = np.sqrt(md.mesh.x[i]**2 + md.mesh.y[i]**2) if (dist <= rad): @@ -21,15 +25,16 @@ md.geometry.thickness[i] = 1.0 # non - zero thickness md.geometry.thickness = md.geometry.thickness.reshape(-1, 1) -md.geometry.base = np.zeros((md.mesh.numberofvertices, )) +md.geometry.base = np.zeros((md.mesh.numberofvertices, 1)) md.geometry.surface = md.geometry.thickness + md.geometry.base.reshape(-1, 1) #would otherwise create a 91x91 matrix #Ice density used for benchmarking, not 917 kg / m^3 -md.materials.rho_ice = 1000 #kg m^3 +md.materials.rho_ice = 1000 #kg m^3 -#GIA parameters specific to Experiments A and B -md.gia.mantle_viscosity = 1e21 * np.ones((md.mesh.numberofvertices, )) #in Pa.s -md.gia.lithosphere_thickness = 100 * np.ones((md.mesh.numberofvertices, )) #in km +#GIA parameters specific to Experiments A and B +md.gia=giaivins(); +md.gia.mantle_viscosity = 1e21 * np.ones((md.mesh.numberofvertices, 1)) #in Pa.s +md.gia.lithosphere_thickness = 100 * np.ones((md.mesh.numberofvertices, 1)) #in km 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 @@ -49,19 +54,19 @@ x = None y = None index = None -md.initialization.vz = np.zeros((md.mesh.numberofvertices, )) -md.initialization.pressure = np.zeros((md.mesh.numberofvertices, )) +md.initialization.vz = np.zeros((md.mesh.numberofvertices, 1)) +md.initialization.pressure = np.zeros((md.mesh.numberofvertices, 1)) #Materials -md.initialization.temperature = (273. - 20.) * np.ones((md.mesh.numberofvertices, )) +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, )) +md.materials.rheology_n = 3. * np.ones((md.mesh.numberofelements, 1)) #Friction -md.friction.coefficient = 20. * np.ones((md.mesh.numberofvertices, )) +md.friction.coefficient = 20. * np.ones((md.mesh.numberofvertices, 1)) md.friction.coefficient[np.where(md.mask.ocean_levelset < 0.)] = 0. -md.friction.p = np.ones((md.mesh.numberofelements, )) -md.friction.q = np.ones((md.mesh.numberofelements, )) +md.friction.p = np.ones((md.mesh.numberofelements, 1)) +md.friction.q = np.ones((md.mesh.numberofelements, 1)) #Numerical parameters md.groundingline.migration = 'None' @@ -72,7 +77,7 @@ md.stressbalance.restol = 0.05 md.steadystate.reltol = 0.05 md.stressbalance.reltol = 0.05 -md.stressbalance.abstol = float('NaN') +md.stressbalance.abstol = np.nan md.timestepping.time_step = 1. md.timestepping.final_time = 3. Index: ../trunk-jpl/test/NightlyRun/test2071.py =================================================================== --- ../trunk-jpl/test/NightlyRun/test2071.py (revision 25171) +++ ../trunk-jpl/test/NightlyRun/test2071.py (revision 25172) @@ -16,11 +16,11 @@ md = parameterize(md, '../Par/GiaIvinsBenchmarksCD.py') #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 md.timestepping.start_time = 0.3 # for t \approx 0 kyr : to get eleastic response! -md.timestepping.final_time = 2500000 # 2, 500 kyr +md.timestepping.final_time = 2500000 # 2,500 kyr md.geometry.thickness = np.array([np.append(md.geometry.thickness * 0.0, 0.0), np.append(md.geometry.thickness / 2.0, 0.1), np.append(md.geometry.thickness, 0.2), Index: ../trunk-jpl/src/m/classes/model.m =================================================================== --- ../trunk-jpl/src/m/classes/model.m (revision 25171) +++ ../trunk-jpl/src/m/classes/model.m (revision 25172) @@ -1263,7 +1263,6 @@ md.calving = calving(); md.frontalforcings = frontalforcings(); md.gia = giamme(); - md.esa = esa(); md.love = fourierlove(); md.esa = esa(); md.autodiff = autodiff(); Index: ../trunk-jpl/src/m/classes/model.py =================================================================== --- ../trunk-jpl/src/m/classes/model.py (revision 25171) +++ ../trunk-jpl/src/m/classes/model.py (revision 25172) @@ -86,29 +86,27 @@ ''' self.mesh = mesh2d() self.mask = mask() + self.constants = constants() self.geometry = geometry() - self.constants = constants() + self.initialization = initialization() self.smb = SMBforcing() self.basalforcings = basalforcings() + self.friction = friction() + self.rifts = rifts() + self.solidearth = solidearth() + self.dsl = dsl() + self.timestepping = timestepping() + self.groundingline = groundingline() self.materials = matice() self.damage = damage() - self.friction = friction() self.flowequation = flowequation() - self.timestepping = timestepping() - self.initialization = initialization() - self.rifts = rifts() - self.dsl = dsl() - self.solidearth = solidearth() - self.debug = debug() self.verbose = verbose() self.settings = issmsettings() self.toolkits = toolkits() self.cluster = generic() - self.balancethickness = balancethickness() self.stressbalance = stressbalance() - self.groundingline = groundingline() self.hydrology = hydrologyshreve() self.masstransport = masstransport() self.thermal = thermal() @@ -124,10 +122,9 @@ self.inversion = inversion() self.qmu = qmu() self.amr = amr() - + self.radaroverlay = radaroverlay() self.results = results() self.outputdefinition = outputdefinition() - self.radaroverlay = radaroverlay() self.miscellaneous = miscellaneous() self.private = private() #}}} @@ -136,19 +133,20 @@ # ordered list of properties since vars(self) is random return ['mesh', 'mask', + 'constants', 'geometry', - 'constants', + 'initialization', 'smb', 'basalforcings', + 'friction', + 'rifts', + 'solidearth', + 'dsl', + 'timestepping', + 'groundingline', 'materials', 'damage', - 'friction', 'flowequation', - 'timestepping', - 'initialization', - 'rifts', - 'dsl', - 'solidearth', 'debug', 'verbose', 'settings', @@ -156,7 +154,6 @@ 'cluster', 'balancethickness', 'stressbalance', - 'groundingline', 'hydrology', 'masstransport', 'thermal', @@ -165,16 +162,16 @@ 'levelset', 'calving', 'frontalforcings', + 'gia', 'love', - 'gia', 'esa', 'autodiff', 'inversion', 'qmu', 'amr', + 'radaroverlay', + 'results', 'outputdefinition', - 'results', - 'radaroverlay', 'miscellaneous', 'private'] # }}}