Index: /issm/trunk-jpl/src/m/boundaryconditions/getlovenumbers.py =================================================================== --- /issm/trunk-jpl/src/m/boundaryconditions/getlovenumbers.py (revision 25160) +++ /issm/trunk-jpl/src/m/boundaryconditions/getlovenumbers.py (revision 25161) @@ -10047,7 +10047,4 @@ #cut off series at maxdeg love_numbers = np.delete(love_numbers, range(maxdeg + 1, len(love_numbers)), axis=0) - # print(len(love_numbers)) - # print(love_numbers[-1,:]) - # exit() #retrieve right type Index: /issm/trunk-jpl/src/m/classes/friction.m =================================================================== --- /issm/trunk-jpl/src/m/classes/friction.m (revision 25160) +++ /issm/trunk-jpl/src/m/classes/friction.m (revision 25161) @@ -1,6 +1,6 @@ %FRICTION class definition % -% Usage: -% friction=friction(); +% Usage: +% friction=friction(); classdef friction Index: /issm/trunk-jpl/src/m/classes/friction.py =================================================================== --- /issm/trunk-jpl/src/m/classes/friction.py (revision 25160) +++ /issm/trunk-jpl/src/m/classes/friction.py (revision 25161) @@ -1,15 +1,15 @@ +from checkfield import checkfield from fielddisplay import fielddisplay from project3d import project3d -from checkfield import checkfield from WriteData import WriteData class friction(object): - """ + ''' FRICTION class definition - Usage: - friction = friction() - """ + Usage: + friction = friction() + ''' def __init__(self): # {{{ @@ -22,5 +22,5 @@ #set defaults self.setdefaultparameters() - self.requested_outputs = [] + #self.requested_outputs = [] #}}} @@ -34,5 +34,5 @@ string = "%s\n%s" % (string, fielddisplay(self, 'effective_pressure', 'Effective Pressure for the forcing if not coupled [Pa]')) string = "%s\n%s" % (string, fielddisplay(self, 'effective_pressure_limit', 'Neff do not allow to fall below a certain limit: effective_pressure_limit * rho_ice * g * thickness (default 0)')) - string = "%s\n%s" % (string, fielddisplay(self, 'requested_outputs', 'additional outputs requested')) + #string = "%s\n%s" % (string, fielddisplay(self, 'requested_outputs', 'additional outputs requested')) return string #}}} @@ -51,5 +51,5 @@ def setdefaultparameters(self): # {{{ - self.requested_outputs = ['default'] + #self.requested_outputs = ['default'] self.effective_pressure_limit = 0 return self @@ -62,5 +62,4 @@ def checkconsistency(self, md, solution, analyses): # {{{ - #Early return if 'StressbalanceAnalysis' not in analyses and 'ThermalAnalysis' not in analyses: @@ -76,5 +75,5 @@ elif self.coupling > 4: raise ValueError('md.friction.coupling larger than 4, not supported yet') - md = checkfield(md, 'fieldname', 'friction.requested_outputs', 'stringrow', 1) + #md = checkfield(md, 'fieldname', 'friction.requested_outputs', 'stringrow', 1) return md # }}} @@ -94,10 +93,10 @@ raise ValueError('md.friction.coupling larger than 4, not supported yet') - #process requested outputs - outputs = self.requested_outputs - indices = [i for i, x in enumerate(outputs) if x == 'default'] - if len(indices) > 0: - outputscopy = outputs[0:max(0, indices[0] - 1)] + self.defaultoutputs(md) + outputs[indices[0] + 1:] - outputs = outputscopy - WriteData(fid, prefix, 'data', outputs, 'name', 'md.friction.requested_outputs', 'format', 'StringArray') + # #process requested outputs + # outputs = self.requested_outputs + # indices = [i for i, x in enumerate(outputs) if x == 'default'] + # if len(indices) > 0: + # outputscopy = outputs[0:max(0, indices[0] - 1)] + self.defaultoutputs(md) + outputs[indices[0] + 1:] + # outputs = outputscopy + # WriteData(fid, prefix, 'data', outputs, 'name', 'md.friction.requested_outputs', 'format', 'StringArray') # }}} Index: /issm/trunk-jpl/src/m/classes/model.py =================================================================== --- /issm/trunk-jpl/src/m/classes/model.py (revision 25160) +++ /issm/trunk-jpl/src/m/classes/model.py (revision 25161) @@ -148,6 +148,6 @@ 'initialization', 'rifts', + 'dsl', 'solidearth', - 'dsl', 'debug', 'verbose', Index: /issm/trunk-jpl/src/m/classes/solidearth.py =================================================================== --- /issm/trunk-jpl/src/m/classes/solidearth.py (revision 25160) +++ /issm/trunk-jpl/src/m/classes/solidearth.py (revision 25161) @@ -88,8 +88,7 @@ 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 - + if len(pos): + outputs = np.delete(outputs, pos) #remove 'default' from outputs + outputs = np.append(outputs, self.defaultoutputs(md)) #add defaults WriteData(fid, prefix, 'data', outputs, 'name', 'md.solidearth.requested_outputs', 'format', 'StringArray') #}}} Index: /issm/trunk-jpl/src/m/coordsystems/gmtmask.py =================================================================== --- /issm/trunk-jpl/src/m/coordsystems/gmtmask.py (revision 25160) +++ /issm/trunk-jpl/src/m/coordsystems/gmtmask.py (revision 25161) @@ -1,14 +1,17 @@ -from MatlabFuncs import * -from model import * -import numpy as np from os import getenv, putenv import subprocess +import numpy as np -def gmtmask(lat, long, *varargin): - '''GMTMASK - figure out which lat, long points are on the ocean +from MatlabFuncs import * +from model import * - Usage: - mask.ocean = gmtmask(md.mesh.lat, md.mesh.long) + +def gmtmask(lat, long, *args): + ''' + GMTMASK - figure out which lat, long points are on the ocean + + Usage: + mask.ocean = gmtmask(md.mesh.lat, md.mesh.long) ''' lenlat = len(lat) @@ -16,5 +19,5 @@ #are we doing a recursive call? - if len(varargin) == 3: + if len(args) == 3: recursive = 1 else: @@ -24,8 +27,7 @@ print((' recursing: num vertices #' + str(lenlat))) else: - print(('gmtmask: num vertices #' + str(lenlat))) + print(('gmtmask: num vertices ' + str(lenlat))) - #Check lat and long size is not more than 50, 000 If so, recursively call gmtmask: - + #Check lat and long size is not more than 50,000. If so, recursively call gmtmask: if lenlat > 50000: for i in range(int(ceil(lenlat / 50000))): Index: /issm/trunk-jpl/src/m/solve/marshall.m =================================================================== --- /issm/trunk-jpl/src/m/solve/marshall.m (revision 25160) +++ /issm/trunk-jpl/src/m/solve/marshall.m (revision 25161) @@ -34,5 +34,5 @@ %Marshall current object - %disp(['marshalling ' field '...']); + %disp(['marshalling ' field '...']); %Uncomment for debugging marshall(md.(field),['md.' field],md,fid); end Index: /issm/trunk-jpl/src/m/solve/marshall.py =================================================================== --- /issm/trunk-jpl/src/m/solve/marshall.py (revision 25160) +++ /issm/trunk-jpl/src/m/solve/marshall.py (revision 25161) @@ -13,5 +13,5 @@ ''' if md.verbose.solution: - print(("marshalling file '%s.bin'." % md.miscellaneous.name)) + print("marshalling file {}.bin".format(md.miscellaneous.name)) #open file for binary writing Index: /issm/trunk-jpl/src/m/solve/solve.m =================================================================== --- /issm/trunk-jpl/src/m/solve/solve.m (revision 25160) +++ /issm/trunk-jpl/src/m/solve/solve.m (revision 25161) @@ -151,4 +151,5 @@ return; end + %wait on lock if isnan(md.settings.waitonlock), Index: /issm/trunk-jpl/test/NightlyRun/test2002.m =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test2002.m (revision 25160) +++ /issm/trunk-jpl/test/NightlyRun/test2002.m (revision 25161) @@ -15,8 +15,8 @@ late=sum(md.mesh.lat(md.mesh.elements),2)/3; longe=sum(md.mesh.long(md.mesh.elements),2)/3; -pos=find(late <-80); +pos=find(late < -80); md.solidearth.surfaceload.icethicknesschange(pos)=-100; %greenland -pos=find(late > 70 & late < 80 & longe>-60 & longe<-30); +pos=find(late>70 & late<80 & longe>-60 & longe<-30); md.solidearth.surfaceload.icethicknesschange(pos)=-100; @@ -28,6 +28,8 @@ mask=gmtmask(md.mesh.lat,md.mesh.long); icemask=ones(md.mesh.numberofvertices,1); -pos=find(mask==0); icemask(pos)=-1; -pos=find(sum(mask(md.mesh.elements),2)<3); icemask(md.mesh.elements(pos,:))=-1; +pos=find(mask==0); +icemask(pos)=-1; +pos=find(sum(mask(md.mesh.elements),2)<3); +icemask(md.mesh.elements(pos,:))=-1; md.mask.ice_levelset=icemask; md.mask.ocean_levelset=-icemask; Index: /issm/trunk-jpl/test/NightlyRun/test2003.m =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test2003.m (revision 25160) +++ /issm/trunk-jpl/test/NightlyRun/test2003.m (revision 25161) @@ -3,9 +3,8 @@ %mesh earth: md=model; - -md.mesh=gmshplanet('radius',6.371012*10^3,'resolution',1000.); %500 km resolution mesh +md.mesh=gmshplanet('radius',6.371012*10^3,'resolution',1000.); %1000 km resolution mesh %parameterize slr solution: -%slr loading: {{{ +%solidearth loading: {{{ md.solidearth.surfaceload.icethicknesschange=zeros(md.mesh.numberofelements,1); md.solidearth.sealevel=zeros(md.mesh.numberofvertices,1); @@ -27,6 +26,8 @@ mask=gmtmask(md.mesh.lat,md.mesh.long); icemask=ones(md.mesh.numberofvertices,1); -pos=find(mask==0); icemask(pos)=-1; -pos=find(sum(mask(md.mesh.elements),2)<3); icemask(md.mesh.elements(pos,:))=-1; +pos=find(mask==0); +icemask(pos)=-1; +pos=find(sum(mask(md.mesh.elements),2)<3); +icemask(md.mesh.elements(pos,:))=-1; md.mask.ice_levelset=icemask; md.mask.ocean_levelset=-icemask; @@ -41,5 +42,5 @@ % }}} -% use model representation of ocea area (not the ture area) +% use model representation of ocen area (not the true area) md.solidearth.settings.ocean_area_scaling = 0; @@ -65,5 +66,7 @@ %eustatic + rigid + elastic run: -md.solidearth.settings.rigid=1; md.solidearth.settings.elastic=1; md.solidearth.settings.rotation=0; +md.solidearth.settings.rigid=1; +md.solidearth.settings.elastic=1; +md.solidearth.settings.rotation=0; md.cluster=generic('name',oshostname(),'np',3); %md.verbose=verbose('111111111'); @@ -72,5 +75,6 @@ %eustatic + rigid + elastic + rotation run: -md.solidearth.settings.rigid=1; md.solidearth.settings.elastic=1; md.solidearth.settings.rotation=1; +md.solidearth.settings.rigid=1; +md.solidearth.settings.elastic=1; md.solidearth.settings.rotation=1; md.cluster=generic('name',oshostname(),'np',3); %md.verbose=verbose('111111111'); Index: /issm/trunk-jpl/test/NightlyRun/test2003.py =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test2003.py (revision 25160) +++ /issm/trunk-jpl/test/NightlyRun/test2003.py (revision 25161) @@ -14,9 +14,9 @@ #mesh earth: md = model() -md.mesh = gmshplanet('radius', 6.371012 * 1e3, 'resolution', 1000.) #500 km resolution mesh +md.mesh = gmshplanet('radius', 6.371012 * 1e3, 'resolution', 1000.) #1000 km resolution mesh #parameterize solidearth solution: #solidearth loading: {{{ -md.solidearth.deltathickness = np.zeros((md.mesh.numberofelements, )) +md.solidearth.surfaceload.icethicknesschange = np.zeros((md.mesh.numberofelements, )) md.solidearth.sealevel = np.zeros((md.mesh.numberofvertices, )) md.dsl.global_average_thermosteric_sea_level_change=np.zeros((2, )) @@ -30,5 +30,5 @@ 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.solidearth.deltathickness[pos] = -1 +md.solidearth.surfaceload.icethicknesschange[pos] = -1 #elastic loading from love numbers: @@ -38,33 +38,28 @@ #mask: {{{ mask = gmtmask(md.mesh.lat, md.mesh.long) -icemask = np.ones((md.mesh.numberofvertices, )) -pos = np.where(mask == 0) -#pos[0] because np.where(mask = 0) returns a 2d array, the latter parts of which are all array / s of 0s -icemask[pos[0]] = -1 -pos = np.where(sum(mask[md.mesh.elements - 1], 1) < 3) +icemask = np.ones(md.mesh.numberofvertices) +pos = np.where(mask == 0)[0] +icemask[pos] = -1 +pos = np.where(np.sum(mask[md.mesh.elements - 1], axis=1) < 3) icemask[md.mesh.elements[pos, :] - 1] = -1 md.mask.ice_levelset = icemask md.mask.ocean_levelset = -icemask -#make sure that the ice level set is all inclusive: -md.mask.land_levelset = np.zeros((md.mesh.numberofvertices, )) -md.mask.ocean_levelset = -np.ones((md.mesh.numberofvertices, )) - -#make sure that the elements that have loads are fully grounded: -pos = np.nonzero(md.solidearth.deltathickness)[0] +#make sure that the elements that have loads are fully grounded +pos = np.nonzero(md.solidearth.surfaceload.icethicknesschange)[0] md.mask.ocean_levelset[md.mesh.elements[pos, :] - 1] = 1 #make sure wherever there is an ice load, that the mask is set to ice: -icemask[md.mesh.elements[pos, :] - 1] = -1 -md.mask.ice_levelset = icemask +#pos = np.nonzero(md.solidearth.surfaceload.icethicknesschange)[0] # Do we need to do this twice? +md.mask.ice_levelset[md.mesh.elements[pos, :] - 1] = 1 # }}} -# use model representation of ocea area (not the ture area) -md.solidearth.ocean_area_scaling = 0 +# use model representation of ocen area (not the true area) +md.solidearth.settings.ocean_area_scaling = 0 #geometry di = md.materials.rho_ice / md.materials.rho_water -md.geometry.thickness = np.ones((md.mesh.numberofvertices, )) -md.geometry.surface = (1 - di) * np.zeros((md.mesh.numberofvertices, )) +md.geometry.thickness = np.ones(md.mesh.numberofvertices) +md.geometry.surface = (1 - di) * np.zeros(md.mesh.numberofvertices) md.geometry.base = md.geometry.surface - md.geometry.thickness md.geometry.bed = md.geometry.base @@ -78,33 +73,22 @@ md.miscellaneous.name = 'test2003' -#New stuff -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.solidearth.reltol = float('NaN') -md.solidearth.abstol = 1e-3 -md.solidearth.geodetic = 1 +md.solidearth.settings.reltol = np.nan +md.solidearth.settings.abstol = 1e-3 +md.solidearth.settings.computesealevelchange = 1 #eustatic + rigid + elastic run: -md.solidearth.rigid = 1 -md.solidearth.elastic = 1 -md.solidearth.rotation = 0 +md.solidearth.settings.rigid = 1 +md.solidearth.settings.elastic = 1 +md.solidearth.settings.rotation = 0 md.cluster = generic('name', gethostname(), 'np', 3) #md.verbose = verbose('111111111') -#print md.calving -#print md.gia -#print md.love -#print md.esa -#print md.autodiff md = solve(md, 'Sealevelrise') SnoRotation = md.results.SealevelriseSolution.Sealevel #eustatic + rigid + elastic + rotation run: -md.solidearth.rigid = 1 -md.solidearth.elastic = 1 -md.solidearth.rotation = 1 +md.solidearth.settings.rigid = 1 +md.solidearth.settings.elastic = 1 +md.solidearth.settings.rotation = 1 md.cluster = generic('name', gethostname(), 'np', 3) #md.verbose = verbose('111111111')