Index: /issm/trunk-jpl/src/m/classes/slr.py =================================================================== --- /issm/trunk-jpl/src/m/classes/slr.py (revision 23087) +++ /issm/trunk-jpl/src/m/classes/slr.py (revision 23088) @@ -11,5 +11,5 @@ Usage: - slr=slr(); + slr=slr() """ @@ -17,4 +17,5 @@ self.deltathickness = float('NaN') self.sealevel = float('NaN') + self.spcthickness = float('NaN') self.maxiter = 0 self.reltol = 0 @@ -23,6 +24,6 @@ self.love_k = 0 #ideam self.love_l = 0 #ideam - self.tide_love_h = 0 - self.tide_love_k = 0 + self.tide_love_k = 0 #ideam + self.tide_love_h = 0 #ideam self.fluid_love = 0 self.equatorial_moi = 0 @@ -51,4 +52,5 @@ string="%s\n%s"%(string,fielddisplay(self,'deltathickness','thickness change: ice height equivalent [m]')) string="%s\n%s"%(string,fielddisplay(self,'sealevel','current sea level (prior to computation) [m]')) + string="%s\n%s"%(string,fielddisplay(self,'spcthickness','thickness constraints (NaN means no constraint) [m]')) string="%s\n%s"%(string,fielddisplay(self,'reltol','sea level rise relative convergence criterion, (NaN: not applied)')) string="%s\n%s"%(string,fielddisplay(self,'abstol','sea level rise absolute convergence criterion, (default, NaN: not applied)')) @@ -62,5 +64,5 @@ string="%s\n%s"%(string,fielddisplay(self,'equatorial_moi','mean equatorial moment of inertia [kg m^2]')) string="%s\n%s"%(string,fielddisplay(self,'polar_moi','polar moment of inertia [kg m^2]')) - string="%s\n%s"%(string,fielddisplay(self,'angular_velocity','mean rotational velocity of earth [per second]')); + string="%s\n%s"%(string,fielddisplay(self,'angular_velocity','mean rotational velocity of earth [per second]')) string="%s\n%s"%(string,fielddisplay(self,'ocean_area_scaling','correction for model representation of ocean area [default: No correction]')) string="%s\n%s"%(string,fielddisplay(self,'steric_rate','rate of steric ocean expansion [mm/yr]')) @@ -82,10 +84,10 @@ #Convergence criterion: absolute, relative and residual - self.reltol=float('NaN') #default - self.abstol=0.001 #1 mm of sea level rise + self.reltol=0.01 #default + self.abstol=float('NaN') #1 mm of sea level rise #maximum of non-linear iterations. self.maxiter=5 - self.loop_increment=200; + self.loop_increment=200 #computational flags: @@ -93,20 +95,20 @@ self.rigid=1 self.elastic=1 - self.rotation=0 self.ocean_area_scaling=0 + self.rotation=1 #tidal love numbers: - self.tide_love_h=0.6149; #degree 2 - self.tide_love_k=0.3055; #degree 2 + self.tide_love_h=0.6149 #degree 2 + self.tide_love_k=0.3055 #degree 2 #secular fluid love number: - self.fluid_love=0.942; + self.fluid_love=0.942 #moment of inertia: - self.equatorial_moi=8.0077*10**37; # [kg m^2] - self.polar_moi =8.0345*10**37; # [kg m^2] + self.equatorial_moi=8.0077*10**37 # [kg m^2] + self.polar_moi =8.0345*10**37 # [kg m^2] #mean rotational velocity of earth - self.angular_velocity=7.2921*10**-5; # [s^-1] + self.angular_velocity=7.2921*10**-5 # [s^-1] #numerical discretization accuracy @@ -114,8 +116,8 @@ #steric: - self.steric_rate=0; + self.steric_rate=0 #how many time steps we skip before we run SLR solver during transient - self.geodetic_run_frequency=1; + self.geodetic_run_frequency=1 #output default: @@ -125,6 +127,7 @@ self.transitions=[] - #default output - self.requested_outputs=['default'] + #horizontal displacement? (not by default) + self.horiz=0 + return self #}}} @@ -137,4 +140,5 @@ md = checkfield(md,'fieldname','slr.deltathickness','NaN',1,'Inf',1,'size',[md.mesh.numberofelements]) md = checkfield(md,'fieldname','slr.sealevel','NaN',1,'Inf',1,'size',[md.mesh.numberofvertices]) + md = checkfield(md,'fieldname','slr.spcthickness','Inf',1,'timeseries',1) md = checkfield(md,'fieldname','slr.love_h','NaN',1,'Inf',1) md = checkfield(md,'fieldname','slr.love_k','NaN',1,'Inf',1) @@ -152,9 +156,9 @@ md = checkfield(md,'fieldname','slr.steric_rate','NaN',1,'Inf',1,'size',[md.mesh.numberofvertices]) md = checkfield(md,'fieldname','slr.degacc','size',[1,1],'>=',1e-10) + md = checkfield(md,'fieldname','slr.requested_outputs','stringrow',1) md = checkfield(md,'fieldname','slr.loop_increment','NaN',1,'Inf',1,'>=',1) md = checkfield(md,'fieldname','slr.horiz','NaN',1,'Inf',1,'values',[0,1]) md = checkfield(md,'fieldname','slr.Ngia','NaN',1,'Inf',1,'size',[md.mesh.numberofvertices]) md = checkfield(md,'fieldname','slr.Ugia','NaN',1,'Inf',1,'size',[md.mesh.numberofvertices]) - md = checkfield(md,'fieldname','slr.requested_outputs','stringrow',1) #check that love numbers are provided at the same level of accuracy: @@ -172,5 +176,5 @@ #a coupler to a planet model is provided. if self.geodetic and not md.transient.iscoupler and domaintype(md.mesh)!='mesh3dsurface': - error('model is requesting geodetic computations without being a mesh3dsurface, or being coupled to one!') + error('model is requesting geodetic computations without being a mesh3dsurface, or being coupled to one!') return md # }}} @@ -181,4 +185,5 @@ WriteData(fid,prefix,'object',self,'fieldname','deltathickness','format','DoubleMat','mattype',2) 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','spcthickness','format','DoubleMat','mattype',1,'timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts) WriteData(fid,prefix,'object',self,'fieldname','reltol','format','Double') WriteData(fid,prefix,'object',self,'fieldname','abstol','format','Double') @@ -187,10 +192,10 @@ WriteData(fid,prefix,'object',self,'fieldname','love_k','format','DoubleMat','mattype',1) WriteData(fid,prefix,'object',self,'fieldname','love_l','format','DoubleMat','mattype',1) - WriteData(fid,prefix,'object',self,'fieldname','tide_love_h','format','Double'); - WriteData(fid,prefix,'object',self,'fieldname','tide_love_k','format','Double'); - WriteData(fid,prefix,'object',self,'fieldname','fluid_love','format','Double'); - WriteData(fid,prefix,'object',self,'fieldname','equatorial_moi','format','Double'); - WriteData(fid,prefix,'object',self,'fieldname','polar_moi','format','Double'); - WriteData(fid,prefix,'object',self,'fieldname','angular_velocity','format','Double'); + WriteData(fid,prefix,'object',self,'fieldname','tide_love_h','format','Double') + WriteData(fid,prefix,'object',self,'fieldname','tide_love_k','format','Double') + WriteData(fid,prefix,'object',self,'fieldname','fluid_love','format','Double') + WriteData(fid,prefix,'object',self,'fieldname','equatorial_moi','format','Double') + WriteData(fid,prefix,'object',self,'fieldname','polar_moi','format','Double') + WriteData(fid,prefix,'object',self,'fieldname','angular_velocity','format','Double') WriteData(fid,prefix,'object',self,'fieldname','rigid','format','Boolean') WriteData(fid,prefix,'object',self,'fieldname','elastic','format','Boolean') @@ -199,10 +204,10 @@ WriteData(fid,prefix,'object',self,'fieldname','geodetic_run_frequency','format','Integer') WriteData(fid,prefix,'object',self,'fieldname','steric_rate','format','DoubleMat','mattype',1,'scale',1e-3/md.constants.yts) - WriteData(fid,prefix,'object',self,'fieldname','degacc','format','Double') - WriteData(fid,prefix,'object',self,'fieldname','loop_increment','format','Integer'); - WriteData(fid,prefix,'object',self,'fieldname','horiz','format','Integer'); WriteData(fid,prefix,'object',self,'fieldname','Ngia','format','DoubleMat','mattype',1,'scale',1e-3/md.constants.yts) WriteData(fid,prefix,'object',self,'fieldname','Ugia','format','DoubleMat','mattype',1,'scale',1e-3/md.constants.yts) + WriteData(fid,prefix,'object',self,'fieldname','degacc','format','Double') WriteData(fid,prefix,'object',self,'fieldname','transitions','format','MatArray') + WriteData(fid,prefix,'object',self,'fieldname','loop_increment','format','Integer') + WriteData(fid,prefix,'object',self,'fieldname','horiz','format','Integer') WriteData(fid,prefix,'object',self,'fieldname','geodetic','format','Integer') Index: /issm/trunk-jpl/src/m/classes/stressbalance.py =================================================================== --- /issm/trunk-jpl/src/m/classes/stressbalance.py (revision 23087) +++ /issm/trunk-jpl/src/m/classes/stressbalance.py (revision 23088) @@ -185,15 +185,9 @@ WriteData(fid,prefix,'object',self,'class','stressbalance','fieldname','rift_penalty_threshold','format','Integer') WriteData(fid,prefix,'object',self,'class','stressbalance','fieldname','referential','format','DoubleMat','mattype',1) - - if isinstance(self.loadingforce, (list, tuple, np.ndarray)): - lx=self.loadingforce[:,0]; - ly=self.loadingforce[:,1]; - lz=self.loadingforce[:,2]; - else: - lx=float('NaN'); ly=float('NaN'); lz=float('NaN'); - - WriteData(fid,prefix,'data',lx,'format','DoubleMat','mattype',1,'name','md.stressbalance.loadingforcex') - WriteData(fid,prefix,'data',ly,'format','DoubleMat','mattype',1,'name','md.stressbalance.loadingforcey') - WriteData(fid,prefix,'data',lz,'format','DoubleMat','mattype',1,'name','md.stressbalance.loadingforcez') + + if isinstance(self.loadingforce, (list, tuple, np.ndarray)) and np.size(self.loadingforce,1) == 3: + WriteData(fid,prefix,'data',self.loadingforce[:,0],'format','DoubleMat','mattype',1,'name','md.stressbalance.loadingforcex') + WriteData(fid,prefix,'data',self.loadingforce[:,1],'format','DoubleMat','mattype',1,'name','md.stressbalance.loadingforcey') + WriteData(fid,prefix,'data',self.loadingforce[:,2],'format','DoubleMat','mattype',1,'name','md.stressbalance.loadingforcez') #process requested outputs Index: /issm/trunk-jpl/test/NightlyRun/test2002.py =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test2002.py (revision 23087) +++ /issm/trunk-jpl/test/NightlyRun/test2002.py (revision 23088) @@ -46,5 +46,4 @@ icemask[md.mesh.elements[pos,:].astype(int)-1]=-1 md.mask.ice_levelset=icemask - md.mask.ocean_levelset=np.zeros((md.mesh.numberofvertices)) pos=np.where(md.mask.ice_levelset==1) @@ -110,5 +109,5 @@ md.slr.rigid=1 md.slr.elastic=1 -md.slr.rotation=1 +md.slr.rotation=0 md=solve(md,'Sealevelrise') Selastic=md.results.SealevelriseSolution.Sealevel Index: /issm/trunk-jpl/test/NightlyRun/test2003.py =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test2003.py (revision 23087) +++ /issm/trunk-jpl/test/NightlyRun/test2003.py (revision 23088) @@ -12,5 +12,5 @@ #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.) #500 km resolution mesh #parameterize slr solution: @@ -81,11 +81,11 @@ md.materials.rheology_n = 3 * np.ones((md.mesh.numberofelements,)) +#Miscellaneous +md.miscellaneous.name = 'test2003' + #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,)) - -#Miscellaneous -md.miscellaneous.name = 'test2003' #Solution parameters Index: /issm/trunk-jpl/test/NightlyRun/test2010.py =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test2010.py (revision 23087) +++ /issm/trunk-jpl/test/NightlyRun/test2010.py (revision 23088) @@ -99,5 +99,5 @@ # uncomment following 2 lines for md = solve(md,'Sealevelrise') -eus = md.results.SealevelriseSolution.SealevelEustatic +eus = md.results.SealevelriseSolution.SealevelRSLEustatic slr = md.results.SealevelriseSolution.Sealevel moixz = md.results.SealevelriseSolution.SealevelInertiaTensorXZ