Changeset 26744 for issm/trunk/test/NightlyRun/test2004.py

Timestamp:

12/22/21 10:39:44 (3 years ago)

Author:

Mathieu Morlighem

Message:

merged trunk-jpl and trunk for revision 26742

Location:

issm/trunk

Files:

• . (modified) (1 prop)
• test (modified) (1 prop)
• test/NightlyRun/test2004.py (modified) (14 diffs)

issm/trunk/test/NightlyRun/test2004.py

@@ -114 +114 @@
 hmin = hmin * 1000
 hmax = hmax * 1000
-tolerance = 100
+tolerance = 100 # tolerance of 100m on Earth position when merging 3d meshes
 threshold = 5
 defaultoptions = [
@@ -145 +145 @@
     # Vertex connectivity
     md.mesh.vertexconnectivity = NodeConnectivity(md.mesh.elements, md.mesh.numberofvertices)
+
     # Add model to sl icecaps
     sl.addicecap(md)
@@ -178 +179 @@
         print('      reading bedrock')
         md.geometry.bed = -np.ones((md.mesh.numberofvertices, ))
-    #}}}
-
-    # SLR #{{{
+        md.geometry.base = md.geometry.bed
+        md.geometry.thickness = 1000 * np.ones((md.mesh.numberofvertices, ))
+        md.geometry.surface = md.geometry.bed + md.geometry.thickness
+    #}}}
+
+    # SLC #{{{
     if bas.iscontinentany('antarctica'):
         if testagainst2002:
             # TODO: Check if the following works as expected: 'pos' is empty, so nothing is assigned to 'md.solidearth.surfaceload.icethicknesschange[pos]'
-            md.solidearth.surfaceload.icethicknesschange = np.zeros((md.mesh.numberofelements, ))
+            md.masstransport.spcthickness = np.zeros((md.mesh.numberofvertices, ))
             # Antarctica
             late = np.sum(md.mesh.lat[md.mesh.elements - 1], axis=1) / 3
@@ -190 +194 @@
             pos = np.where(late < -85)[0]
             ratio = 0.225314032985172 / 0.193045366574523
-            md.solidearth.surfaceload.icethicknesschange[pos] = -100 * ratio
+            md.masstransport.spcthickness[md.mesh.elements[pos]] = md.masstransport.spcthickness[md.mesh.elements[pos]] - 100 * ratio
         else:
             delH = np.loadtxt('../Data/AIS_delH_trend.txt')
@@ -201 +205 @@
             pos = np.where(longe > 360)[0]
             longe[pos] = longe[pos] - 360
-            delHAIS = InterpFromMesh2d(index, longAIS, latAIS, delHAIS, longe, late)  # NOTE: Compare to corresponding output under MATLAB to understand why we offset triangle indices by 1 (only caught because Triangle.cpp was producing triangles with negative areas)
+            delHAIS = InterpFromMesh2d(index, longAIS, latAIS, delHAIS, longe, late) # NOTE: Compare to corresponding output under MATLAB to understand why we offset triangle indices by 1 (only caught because Triangle.cpp was producing triangles with negative areas)
             northpole = find_point(md.mesh.long, md.mesh.lat, 0, 90)
             delHAIS[northpole] = 0
-            md.solidearth.surfaceload.icethicknesschange = np.mean(delHAIS[md.mesh.elements - 1], axis=1) / 100
-
-        md.solidearth.initialsealevel = np.zeros((md.mesh.numberofvertices, ))
-
-        md.dsl.global_average_thermostatic_sea_level_change = np.zeros((2, 1))
-        md.dsl.sea_surface_height_change_above_geoid = np.zeros((md.mesh.numberofvertices + 1, 1))
-        md.dsl.sea_water_pressure_change_at_sea_floor = np.zeros((md.mesh.numberofvertices + 1, 1))
+            md.masstransport.spcthickness = delHAIS / 100
+
+        md.initialization.sealevel = np.zeros((md.mesh.numberofvertices, ))
+
+        md.dsl.global_average_thermosteric_sea_level = np.zeros((2, 1))
+        md.dsl.sea_surface_height_above_geoid = np.zeros((md.mesh.numberofvertices + 1, 1))
+        md.dsl.sea_water_pressure_at_sea_floor = np.zeros((md.mesh.numberofvertices + 1, 1))
     #}}}
 
@@ -226 +230 @@
 # Parameterize continents #{{{
 for ind in sl.basinindx('continent', ['hemisphereeast', 'hemispherewest']):
-    print("Masks for basin {}".format(sl.icecaps[ind].miscellaneous.name))
+    print('Masks for basin {}'.format(sl.icecaps[ind].miscellaneous.name))
     md = sl.icecaps[ind]
     bas = sl.basins[ind]
@@ -283 +287 @@
     #}}}
 
-    # SLR loading/calibration #{{{
-    md.solidearth.surfaceload.icethicknesschange = np.zeros((md.mesh.numberofelements, ))
+    # SLC loading/calibration #{{{
+    md.masstransport.spcthickness = np.zeros((md.mesh.numberofvertices, ))
 
     if testagainst2002:
@@ -293 +297 @@
         pos = np.where(np.logical_and.reduce((late > 70, late < 80, longe > -60, longe < -30)))[0]
         ratio = .3823 / .262344
+        md.masstransport.spcthickness[md.mesh.elements[pos]] = md.masstransport.spcthickness[md.mesh.elements[pos]] - 100 * ratio
 
         # Correct mask
@@ -307 +312 @@
         pos = np.where(longe > 360)[0]
         longe[pos] = longe[pos] - 360
-        delHGIS = InterpFromMeshToMesh2d(index, longGIS, latGIS, delHGIS, longe, late) # NOTE: Compare to corresponding ouptut under MATLAB to understand why we offset triangle indices by 1 (only caught because Triangle.cpp was producing triangles with negative areas)
-        delHGISe = np.mean(delHGIS[md.mesh.elements - 1], axis=1).flatten()
+        delHGIS = InterpFromMeshToMesh2d(index, longGIS, latGIS, delHGIS, longe, late)
 
         delH = np.loadtxt('../Data/GLA_delH_trend_15regions.txt')
@@ -319 +323 @@
         pos = np.where(longe > 360)[0]
         longe[pos] = longe[pos] - 360
-        delHGLA = InterpFromMeshToMesh2d(index, longGLA, latGLA, delHGLA, longe, late) # NOTE: Compare to corresponding ouptut under MATLAB to understand why we offset triangle indices by 1 (only caught because Triangle.cpp was producing triangles with negative areas)
-        delHGLAe = np.mean(delHGLA[md.mesh.elements - 1], axis=1).flatten()
-
-        pos = np.nonzero(delHGISe)[0]
-        md.solidearth.surfaceload.icethicknesschange[pos] = delHGISe[pos] / 100
-        pos = np.nonzero(delHGLAe)[0]
-        md.solidearth.surfaceload.icethicknesschange[pos] = delHGLAe[pos] / 100
+        delHGLA = InterpFromMeshToMesh2d(index, longGLA, latGLA, delHGLA, longe, late)
+
+        # NOTE: For some reason, cannot apply pos to multiple arrays in a 
+        #       singlelike we might do in MATLAB. Instead, we iterate over 
+        #       elements of pos.
+        #
+        pos = np.nonzero(delHGIS)[0]
+        for p in pos:
+            md.masstransport.spcthickness[p] = md.masstransport.spcthickness[p] - delHGIS[p] / 100
+        pos = np.nonzero(delHGLA)[0]
+        for p in pos:
+            md.masstransport.spcthickness[p] = md.masstransport.spcthickness[p] - delHGLA[p] / 100
 
         # Adjust mask accordingly
-        pos = np.nonzero(md.solidearth.surfaceload.icethicknesschange)[0]
-        flags = np.zeros((md.mesh.numberofvertices, ))
-        flags[md.mesh.elements[pos, :] - 1] = 1
-        pos = np.nonzero(flags)[0]
+        pos = np.nonzero(md.masstransport.spcthickness)[0]
         md.mask.ice_levelset[pos] = -1
         md.mask.ocean_levelset[pos] = 1
 
-    md.solidearth.initialsealevel = np.zeros((md.mesh.numberofvertices, ))
-
-    md.dsl.global_average_thermostatic_sea_level_change = np.zeros((2, 1))
-    md.dsl.sea_surface_height_change_above_geoid = np.zeros((md.mesh.numberofvertices + 1, 1))
-    md.dsl.sea_water_pressure_change_at_sea_floor = np.zeros((md.mesh.numberofvertices + 1, 1))
-    #}}}
+    md.initialization.sealevel = np.zeros((md.mesh.numberofvertices, ))
+
+    md.dsl.global_average_thermosteric_sea_level = np.zeros((2, 1))
+    md.dsl.sea_surface_height_above_geoid = np.zeros((md.mesh.numberofvertices + 1, 1))
+    md.dsl.sea_water_pressure_at_sea_floor = np.zeros((md.mesh.numberofvertices + 1, 1))
+    #}}}
+
     # Geometry #{{{
     di = md.materials.rho_ice / md.materials.rho_water
     md.geometry.bed = -np.ones((md.mesh.numberofvertices, ))
+    md.geometry.base = md.geometry.bed
+    md.geometry.thickness = 1000 * np.ones((md.mesh.numberofvertices, ))
+    md.geometry.surface = md.geometry.bed + md.geometry.thickness
     #}}}
     # Materials #{{{
@@ -375 +385 @@
 sl.transfer('mask.ocean_levelset')
 sl.transfer('geometry.bed')
+sl.transfer('geometry.surface')
+sl.transfer('geometry.thickness')
+sl.transfer('geometry.base')
 sl.transfer('mesh.lat')
 sl.transfer('mesh.long')
-sl.transfer('solidearth.surfaceload.icethicknesschange') #
-sl.transfer('solidearth.initialsealevel')
-sl.transfer('dsl.sea_surface_height_change_above_geoid')
-sl.transfer('dsl.sea_water_pressure_change_at_sea_floor')
+sl.transfer('masstransport.spcthickness') #
+sl.transfer('initialization.sealevel')
+sl.transfer('dsl.sea_surface_height_above_geoid')
+sl.transfer('dsl.sea_water_pressure_at_sea_floor')
 
 # Radius
@@ -398 +411 @@
 
 # Solve Sea-level equation on Earth only #{{{
-md = sl.earth #we don't do computations on ice sheets or land
+md = sl.earth # we don't do computations on ice sheets or land
 
 #Materials
-md.materials=materials('hydro')
+md.materials = materials('hydro')
 
 # Elastic loading from love numbers
@@ -411 +424 @@
 
 # New stuff
-md.dsl.global_average_thermosteric_sea_level_change = np.array([[(1.1 + .38)], [0]]) # steric + water storage AR5
+md.dsl.global_average_thermosteric_sea_level = np.array([[(1.1 + .38)], [0]]) # steric + water storage AR5
 
 # Solutuion parameters
 md.solidearth.settings.reltol = np.nan
 md.solidearth.settings.abstol = 1e-3
-md.solidearth.settings.computesealevelchange = 1
+md.solidearth.settings.sealevelloading = 1
+md.solidearth.settings.isgrd = 1
+md.solidearth.settings.ocean_area_scaling = 1
+md.solidearth.settings.grdmodel = 1
 md.timestepping.time_step = 1
+
+# Physics
+md.transient.issmb = 0
+md.transient.isstressbalance = 0
+md.transient.isthermal = 0
+md.transient.ismasstransport = 1
+md.transient.isslc = 1
+
+# Initializations
+md.basalforcings.groundedice_melting_rate = np.zeros((md.mesh.numberofvertices,))
+md.basalforcings.floatingice_melting_rate = np.zeros((md.mesh.numberofvertices,))
+md.initialization.vx = np.zeros((md.mesh.numberofvertices,))
+md.initialization.vy = np.zeros((md.mesh.numberofvertices,))
+md.initialization.sealevel = np.zeros((md.mesh.numberofvertices,))
+md.initialization.bottompressure = np.zeros((md.mesh.numberofvertices,))
+md.initialization.dsl = np.zeros((md.mesh.numberofvertices,))
+md.initialization.str = 0
+md.smb.mass_balance = np.zeros((md.mesh.numberofvertices,))
 
 # Max number of iterations reverted back to 10 (i.e. the original default value)
@@ -423 +457 @@
 
 # Eustatic run:
-md.solidearth.settings.rigid = 0
+md.solidearth.settings.selfattraction = 0
 md.solidearth.settings.elastic = 0
 md.solidearth.settings.rotation = 0
+md.solidearth.settings.viscous = 0
 md.solidearth.requested_outputs = [
     'default',
-    'SurfaceloadIceThicknessChange',
+    'DeltaIceThickness',
     'Sealevel',
-    'SealevelRSLRate',
-    'SealevelriseCumDeltathickness',
-    'SealevelNEsaRate',
-    'SealevelUEsaRate',
-    'NGiaRate',
-    'UGiaRate',
-    'SealevelEustaticMask',
-    'SealevelEustaticOceanMask'
+    'SealevelUGrd',
+    'SealevelchangeBarystaticMask',
+    'SealevelchangeBarystaticOceanMask',
 ]
-md = solve(md, 'Sealevelrise')
-Seustatic = md.results.SealevelriseSolution.Sealevel
-
-# Eustatic + rigid run
-md.solidearth.settings.rigid = 1
+md = solve(md, 'Transient')
+Seustatic = md.results.TransientSolution.Sealevel
+
+# Eustatic + selfattraction run
+md.solidearth.settings.selfattraction = 1
 md.solidearth.settings.elastic = 0
 md.solidearth.settings.rotation = 0
-md = solve(md, 'Sealevelrise')
-Srigid = md.results.SealevelriseSolution.Sealevel
-
-# Eustatic + rigid + elastic run
-md.solidearth.settings.rigid = 1
+md.solidearth.settings.viscous = 0
+md = solve(md, 'Transient')
+Sselfattraction = md.results.TransientSolution.Sealevel
+
+# Eustatic + selfattraction + elastic run
+md.solidearth.settings.selfattraction = 1
 md.solidearth.settings.elastic = 1
 md.solidearth.settings.rotation = 0
-md = solve(md, 'Sealevelrise')
-Selastic = md.results.SealevelriseSolution.Sealevel
-
-# Eustatic + rigid + elastic + rotation run
-md.solidearth.settings.rigid = 1
+md.solidearth.settings.viscous = 0
+md = solve(md, 'Transient')
+Selastic = md.results.TransientSolution.Sealevel
+
+# Eustatic + selfattraction + elastic + rotation run
+md.solidearth.settings.selfattraction = 1
 md.solidearth.settings.elastic = 1
 md.solidearth.settings.rotation = 1
-md = solve(md, 'Sealevelrise')
-Srotation = md.results.SealevelriseSolution.Sealevel
+md.solidearth.settings.viscous = 0
+md = solve(md, 'Transient')
+Srotation = md.results.TransientSolution.Sealevel
 
 #Fields and tolerances to track changes
 field_names = ['Eustatic', 'Rigid', 'Elastic', 'Rotation']
 field_tolerances = [1e-13, 1e-13, 1e-13, 1e-13]
-field_values = [Seustatic, Srigid, Selastic, Srotation]
+field_values = [Seustatic, Sselfattraction, Selastic, Srotation]