Changeset 27327

Timestamp:

10/24/22 20:23:59 (2 years ago)

Author:

jdquinn

Message:

CHG: Python translations; cleanup

Location:

issm/trunk-jpl

Files:

• src/m/classes/fourierlove.py (modified) (5 diffs)
• src/m/classes/love.m (modified) (4 diffs)
• src/m/classes/love.py (added)
• src/m/classes/model.py (modified) (1 diff)
• src/m/classes/sealevelmodel.py (modified) (6 diffs)
• src/m/solve/solveslm.py (modified) (1 diff)
• test/NightlyRun/test2084.py (modified) (1 diff)

issm/trunk-jpl/src/m/classes/fourierlove.py

@@ -1 +1 @@
 import numpy as np
 
+from checkfield import checkfield
 from fielddisplay import fielddisplay
-from checkfield import checkfield
 from WriteData import WriteData
 
@@ -43 +43 @@
         # - Correct display to match MATLAB
         #
-        s = '   Fourier Love class:\n'
-        s += '{}\n'.format(fielddisplay(self, 'nfreq', 'number of frequencies sampled (default: 1, elastic) [Hz]'))
+        s = '{}\n'.format(fielddisplay(self, 'nfreq', 'number of frequencies sampled (default: 1, elastic) [Hz]'))
         s += '{}\n'.format(fielddisplay(self, 'frequencies', 'frequencies sampled (convention defaults to 0 for the elastic case) [Hz]'))
         s += '{}\n'.format(fielddisplay(self, 'sh_nmax', 'maximum spherical harmonic degree (default: 256, .35 deg, or 40 km at equator)'))
@@ -64 +63 @@
         s += '{}\n'.format(fielddisplay(self, 'forcing_type', 'integer indicating the nature and depth of the forcing for the Love number calculation (default: 11):'))
         s += '{}\n'.format('                                                     1:  Inner core boundary -- Volumic Potential')
-        s += '{}\n'.format('                                                     2:  Inner core boundary --  Pressure')
-        s += '{}\n'.format('                                                     3:  Inner core boundary --  Loading')
-        s += '{}\n'.format('                                                     4:  Inner core boundary --  Tangential traction')
-        s += '{}\n'.format('                                                     5:  Core mantle boundary --  Volumic Potential')
-        s += '{}\n'.format('                                                     6:  Core mantle boundary --  Pressure')
-        s += '{}\n'.format('                                                     7:  Core mantle boundary --  Loading')
-        s += '{}\n'.format('                                                     8:  Core mantle boundary --  Tangential traction')
-        s += '{}\n'.format('                                                     9:  Surface--  Volumic Potential')
-        s += '{}\n'.format('                                                     10: Surface--  Pressure')
-        s += '{}\n'.format('                                                     11: Surface--  Loading')
-        s += '{}\n'.format('                                                     12: Surface--  Tangential traction ')
+        s += '{}\n'.format('                                                     2:  Inner core boundary -- Pressure')
+        s += '{}\n'.format('                                                     3:  Inner core boundary -- Loading')
+        s += '{}\n'.format('                                                     4:  Inner core boundary -- Tangential traction')
+        s += '{}\n'.format('                                                     5:  Core mantle boundary -- Volumic Potential')
+        s += '{}\n'.format('                                                     6:  Core mantle boundary -- Pressure')
+        s += '{}\n'.format('                                                     7:  Core mantle boundary -- Loading')
+        s += '{}\n'.format('                                                     8:  Core mantle boundary -- Tangential traction')
+        s += '{}\n'.format('                                                     9:  Surface -- Volumic Potential')
+        s += '{}\n'.format('                                                     10: Surface -- Pressure')
+        s += '{}\n'.format('                                                     11: Surface -- Loading')
+        s += '{}\n'.format('                                                     12: Surface -- Tangential traction ')
         s += '{}\n'.format(fielddisplay(self, 'inner_core_boundary', 'interface index in materials.radius locating forcing. Only used for forcing_type 1--4 (default: 1)'))
         s += '{}\n'.format(fielddisplay(self, 'core_mantle_boundary', 'interface index in materials.radius locating forcing. Only used for forcing_type 5--8 (default: 2)'))
@@ -87 +86 @@
         self.sh_nmax = 256 # .35 degree, 40 km at the equator
         self.sh_nmin = 1
-        # Work on matlab script for computing g0 for given Earth's structure
+        # Work on Python script for computing g0 for given Earth's structure
         self.g0 = 9.81 # m/s^2
         self.r0 = 6371 * 1e3 # m
@@ -129 +128 @@
 
         md = checkfield(md, 'fieldname', 'love.istemporal', 'values', [0, 1])
+
         if md.love.istemporal:
             md = checkfield(md, 'fieldname', 'love.n_temporal_iterations', 'NaN', 1, 'Inf', 1, 'numel', 1, '>', 0)
             md = checkfield(md, 'fieldname', 'love.time', 'NaN', 1, 'Inf', 1, 'numel', md.love.nfreq / 2 / md.love.n_temporal_iterations)
-        if md.love.sh_nmin <= 1 and (md.love.forcing_type == 9 or md.love.forcing_type == 5 or md.love.forcing_type == 1):
+        if md.love.sh_nmin <= 1 and (md.love.forcing_type == 1 or md.love.forcing_type == 5 or md.love.forcing_type == 9):
             raise RuntimeError('Degree 1 not supported for forcing type {}. Use sh_min >= 2 for this kind of calculation.'.format(md.love.forcing_type))
 

issm/trunk-jpl/src/m/classes/love.m

@@ -33 +33 @@
                 hypergeom_table1=0;
                 hypergeom_table2=0;
-                hypergeom_nalpha=0;
+                hypergeom_table1=0;
                 hypergeom_nz=0;
                 hypergeom_z=0;
@@ -68 +68 @@
                         self.underflow_tol=1e-16; %threshold of deep to surface love number ratio to trigger the deletion of layer 
                         self.pw_threshold=1e-3; %if relative variation across frequencies is smaller than this ratio, the post-widder transform for time-dependent love numbers is bypassed 
-                        self.min_integration_steps      = 50;
-                        self.max_integration_dr = 10e3;
+                        self.min_integration_steps=50;
+                        self.max_integration_dr=1e4;
                         self.integration_scheme=1;
                         self.istemporal=0;
@@ -96 +96 @@
                         fielddisplay(self,'mu0','adimensioning constant for stress (default: 10^11) [Pa]');
                         fielddisplay(self,'Gravitational_Constant','Newtonian constant of gravitation (default: 6.67259e-11 [m^3 kg^-1 s^-2])');
-                        fielddisplay(self,'chandler_wobble','includes the inertial terms for the chandler wobble in the rotational feedback love numbers, only for forcing_type=11 (default: 0) (/!\ 1 is untested yet)');
+                        fielddisplay(self,'chandler_wobble','includes the inertial terms for the chandler wobble in the rotational feedback love numbers, only for forcing_type=11 (default: 0) (/!\ 1 has not been validated yet)');
                         fielddisplay(self,'allow_layer_deletion','allow for migration of the integration boundary with increasing spherical harmonics degree (default: 1)');                    
                         fielddisplay(self,'underflow_tol','threshold of deep to surface love number ratio to trigger the deletion of layers (default: 1e-16)');
@@ -111 +111 @@
                         fielddisplay(self,'core_mantle_boundary','interface index in materials.radius locating forcing. Only used for forcing_type 5--8 (default: 2)'); 
                         fielddisplay(self,'complex_computation','return love numbers as 0: real (useful for elastic or temporal forms), 1: complex numbers (useful for Fourier spectral form) (default: 0)'); 
-                        fielddisplay(self,'quad_precision','toogle computation love numbers and post-widder transform with 32 digit precision, useful for temporal form (default: 1)'); 
+                        fielddisplay(self,'quad_precision','toggle computation love numbers and post-widder transform with 32 digit precision, useful for temporal form (default: 1)'); 
                         fielddisplay(self,'debug','outputs yi system matrix prior to solving (default: 0)'); 
                         fielddisplay(self,'hypergeom_table1','table 1 for hypergeometric function, only for EBM rheology (default: [1])'); 

issm/trunk-jpl/src/m/classes/model.py

@@ -271 +271 @@
         self.calving = calving()
         self.frontalforcings = frontalforcings()
-        self.love = fourierlove()
+        self.love = love()
         self.esa = esa()
         self.sampling = sampling()

issm/trunk-jpl/src/m/classes/sealevelmodel.py

@@ -88 +88 @@
     def checkconsistency(slm, solutiontype):  # {{{
         # Is the coupler turned on?
-        for i in range(len(slm.icecaps)):
-            if not slm.icecaps[i].transient.iscoupler:
-                print('Warning: sealevelmodel.py::checkconsistency: icecap model {} should have the transient coupler option turned on!'.format(slm.icecaps[i].miscellaneous.name))
-
-        if not slm.earth.transient.iscoupler:
-            print('Warning: sealevelmodel.py::checkconsistency: earth model should have the transient coupler option turned on!')
+        #for i in range(len(slm.icecaps)):
+        #    if not slm.icecaps[i].transient.iscoupler:
+        #        print('Warning: sealevelmodel.py::checkconsistency: icecap model {} should have the transient coupler option turned on!'.format(slm.icecaps[i].miscellaneous.name))
+
+        #if not slm.earth.transient.iscoupler:
+        #    print('Warning: sealevelmodel.py::checkconsistency: earth model should have the transient coupler option turned on!')
 
         # Check that the transition vectors have the right size
-        for i in range(len(slm.icecaps)):
-            if slm.icecaps[i].mesh.numberofvertices != len(slm.earth.slr.transitions[i]):
-                raise Exception('sealevelmodel.py::checkconsistency: issue with size of transition vector for ice cap: {} name: {}'.format(i, slm.icecaps[i].miscellaneous.name))
+        if slm.earth.mesh.numberofvertices != len(slm.earth.solidearth.transfercount):
+            raise Exception('sealevelmodel.py::checkconsistency: earth.solidearth.transfercount should be of size earth.mesh.numberofvertices')
 
         # Check that run frequency is the same everywhere
@@ -108 +107 @@
         for i in range(len(slm.icecaps)):
             md = slm.icecaps[i]
-            if np.nonzero(md.dsl.steric_rate - slm.earth.dsl.steric_rate[slm.earth.dsl.transitions[i]]) != []:
+            if np.nonzero(md.dsl.sea_surface_height_above_geoid - slm.earth.dsl.sea_surface_height_above_geoid[slm.transitions[i]]) != []:
                 raise Exception('sealevelmodel.py::checkconsistency: steric rate on ice cap {} is not the same as for the earth'.format(md.miscellaneous.name))
 
@@ -114 +113 @@
         for i in range(len(slm.icecaps)):
             md = slm.icecaps[i]
-            if md.solidearthsettings.isgrd != slm.earth.solidearthsettings.isgrd:
+            if md.solidearth.settings.isgrd != slm.earth.solidearth.settings.isgrd:
                 raise RuntimeError('sealevelmodel.py::checkconsistency: isgrd on ice cap {} is not the same as for the earth\n'.format(md.miscellaneous.name))
 
@@ -221 +220 @@
         self.transitions = []
         self.eltransitions = []
+        self.earth.solidearth.transfercount = np.zeros(self.earth.mesh.numberofvertices)
 
         # For elements
@@ -240 +240 @@
             self.transitions.append(meshintersect3d(self.earth.mesh.x, self.earth.mesh.y, self.earth.mesh.z, mdi.mesh.x, mdi.mesh.y, mdi.mesh.z, 'force', force))
             self.eltransitions.append(meshintersect3d(xe, ye, ze, xei, yei, zei, 'force', force))
+
+            self.earth.solidearth.transfercount[self.transitions[i]] = self.earth.solidearth/transfercount[self.transitions[i]] + 1
+
+        for i in range(len(self.icecaps)):
+            self.icecaps[i].solidearth.transfercount = self.earth.solidearth.transfercount[self.transitions[i]]
     # }}}
 
@@ -246 +251 @@
         for i in range(len(self.basins)):
             flags[self.transitions[i]] = i
-        plotmodel(self.earth, 'data', flags, 'coastline', 'on')
+        plotmodel(self.earth, 'data', flags, 'coastlines', 'on')
     # }}}
 

issm/trunk-jpl/src/m/solve/solveslm.py

@@ -77 +77 @@
         nps.append(slm.earth.cluster.np)
 
-    BuildQueueScriptMultipleModels(cluster, slm.private.runtimename, slm.miscellaneous.name, slm.private.solution, valgrind, privateruntimenames, miscellaneousnames, nps)
+    BuildQueueScriptMultipleModels(cluster, slm.private.runtimename, slm.miscellaneous.name, slm.private.solution, privateruntimenames, miscellaneousnames, nps)
 
     # Upload all required files, given that each individual solution for icecaps and earth model already did

issm/trunk-jpl/test/NightlyRun/test2084.py

@@ -56 +56 @@
 md.love.pw_threshold = 1e-3
 md.love.Gravitational_Constant = 6.6732e-11
-md.love.integration_steps_per_layer = 100
 md.love.allow_layer_deletion = 1
 md.love.forcing_type = 11