Changeset 23088

Timestamp:

08/11/18 13:06:28 (7 years ago)

Author:

kruegern

Message:

BUG: fixed various differences in the python Solid Earth (slr) tests and resulting bins

Location:

issm/trunk-jpl

Files:

• src/m/classes/slr.py (modified) (15 diffs)
• src/m/classes/stressbalance.py (modified) (1 diff)
• test/NightlyRun/test2002.py (modified) (2 diffs)
• test/NightlyRun/test2003.py (modified) (2 diffs)
• test/NightlyRun/test2010.py (modified) (1 diff)

issm/trunk-jpl/src/m/classes/slr.py ¶

@@ -11 +11 @@
         
                 Usage:
-                  slr=slr();
+                  slr=slr()
         """
         
@@ -17 +17 @@
                 self.deltathickness         = float('NaN')
                 self.sealevel               = float('NaN')
+                self.spcthickness           = float('NaN')
                 self.maxiter                = 0
                 self.reltol                 = 0
@@ -23 +24 @@
                 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 +52 @@
                         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 +64 @@
                         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 +84 @@
                 
                 #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 +95 @@
                 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 +116 @@
 
                 #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 +127 @@
                 self.transitions=[]
 
-                #default output
-                self.requested_outputs=['default']
+                #horizontal displacement?  (not by default)
+                self.horiz=0
+
                 return self
                 #}}}
@@ -137 +140 @@
                 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 +156 @@
                 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 +176 @@
                 #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 +185 @@
                 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 +192 @@
                 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 +204 @@
                 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')
         

issm/trunk-jpl/src/m/classes/stressbalance.py ¶

@@ -185 +185 @@
                 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

issm/trunk-jpl/test/NightlyRun/test2002.py ¶

@@ -46 +46 @@
 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 +109 @@
 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

issm/trunk-jpl/test/NightlyRun/test2003.py ¶

@@ -12 +12 @@
 #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 +81 @@
 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

issm/trunk-jpl/test/NightlyRun/test2010.py ¶

@@ -99 +99 @@
 # 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