Changeset 27035 for issm/trunk/src/m/classes/fourierlove.py

Timestamp:

06/01/22 05:01:48 (3 years ago)

Author:

Mathieu Morlighem

Message:

merged trunk-jpl and trunk for revision 27033

Location:

issm/trunk

Files:

• . (modified) (1 prop)
• src (modified) (1 prop)
• src/m/classes/fourierlove.py (modified) (7 diffs)

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

@@ -22 +22 @@
         self.mu0 = 0
         self.Gravitational_Constant = 0
+        self.chandler_wobble = 0
         self.allow_layer_deletion = 0
         self.underflow_tol = 0
+        self.pw_threshold = 0
         self.integration_steps_per_layer = 0
         self.istemporal = 0
@@ -51 +53 @@
         s += '{}\n'.format(fielddisplay(self, 'allow_layer_deletion', 'allow for migration of the integration boundary with increasing spherical harmonics degree (default: 1)'))
         s += '{}\n'.format(fielddisplay(self, 'Gravitational_Constant', 'Newtonian constant of gravitation (default: 6.67259e-11 [m^3 kg^-1 s^-2])'))
+        s += '{}\n'.format(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)'))
         s += '{}\n'.format(fielddisplay(self, 'allow_layer_deletion', 'allow for migration of the integration boundary with increasing spherical harmonics degree (default: 1)'))
-        s += '{}\n'.format(fielddisplay(self, 'underflow_tol', 'threshold of deep to surface love number ratio to trigger the deletion of layers (default: 1e-16)'))
+        s += '{}\n'.format(fielddisplay(self, 'pw_threshold', 'if relative variation across frequencies is smaller than this ratio, the post-widder transform for time-dependent love numbers is bypassed (default (1e-3)'))
+        s += '{}\n'.format(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)'))
         s += '{}\n'.format(fielddisplay(self, 'integration_steps_per_layer', 'number of radial steps to propagate the yi system from the bottom to the top of each layer (default: 100)'))
         s += '{}\n'.format(fielddisplay(self, 'istemporal', ['1 for time-dependent love numbers, 0 for frequency-dependent or elastic love numbers (default: 0)', 'If 1: use fourierlove function build_frequencies_from_time to meet consistency']))
@@ -88 +92 @@
         self.mu0 = 1e11 # Pa
         self.Gravitational_Constant = 6.67259e-11 # m^3 kg^-1 s^-2
+        self.chandler_wobble = 0
         self.allow_layer_deletion = 1
         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.integration_steps_per_layer = 100
         self.istemporal = 0
@@ -113 +119 @@
         md = checkfield(md, 'fieldname', 'love.mu0', 'NaN', 1, 'Inf', 1, 'numel', 1, '>', 0)
         md = checkfield(md, 'fieldname', 'love.Gravitational_Constant', 'NaN', 1, 'Inf', 1, 'numel', 1, '>', 0)
+        md = checkfield(md, 'fieldname', 'love.chandler_wobble', 'values', [0, 1])
         md = checkfield(md, 'fieldname', 'love.allow_layer_deletion', 'values', [0, 1])
         md = checkfield(md, 'fieldname', 'love.underflow_tol', 'NaN', 1, 'Inf', 1, 'numel', 1, '>', 0)
+        md = checkfield(md, 'fieldname', 'love.pw_threshold', 'NaN', 1, 'Inf', 1, 'numel', 1, '>', 0)
         md = checkfield(md, 'fieldname', 'love.integration_steps_per_layer', 'NaN', 1, 'Inf', 1, 'numel', 1, '>', 0)
         md = checkfield(md, 'fieldname', 'love.love_kernels', 'values', [0, 1])
@@ -127 +135 @@
             raise RuntimeError('Degree 1 not supported for forcing type {}. Use sh_min >= 2 for this kind of calculation.'.format(md.love.forcing_type))
 
+        if md.love.chandler_wobble  == 1:
+            print('Warning: Chandler wobble in Love number calculator has not been validated yet')
+
         # Need 'litho' material
-        if not hasattr(md.materials, 'materials') or 'litho' not in md.materials.nature:
+        if md.materials.__class__.__name__ != 'materials' or 'litho' not in md.materials.nature:
             raise RuntimeError('Need a \'litho\' material to run a Fourier Love number analysis')
 
-        mat = np.where(np.array(nature) == 'litho')[0]
+        mat = np.where(np.array(md.materials.nature) == 'litho')[0]
         if md.love.forcing_type <= 4:
             md = checkfield(md, 'fieldname', 'love.inner_core_boundary', 'NaN', 1, 'Inf', 1, 'numel', 1, '>', 0, '<=', md.materials[mat].numlayers)
@@ -149 +160 @@
         WriteData(fid, prefix, 'object', self, 'fieldname', 'mu0', 'format', 'Double')
         WriteData(fid, prefix, 'object', self, 'fieldname', 'Gravitational_Constant', 'format', 'Double')
+        WriteData(fid, prefix, 'object', self, 'fieldname', 'chandler_wobble', 'format', 'Boolean')
         WriteData(fid, prefix, 'object', self, 'fieldname', 'allow_layer_deletion', 'format', 'Boolean')
         WriteData(fid, prefix, 'object', self, 'fieldname', 'underflow_tol', 'format', 'Double')
+        WriteData(fid, prefix, 'object', self, 'fieldname', 'pw_threshold', 'format', 'Double')
         WriteData(fid, prefix, 'object', self, 'fieldname', 'integration_steps_per_layer', 'format', 'Integer')
         WriteData(fid, prefix, 'object', self, 'fieldname', 'istemporal', 'format', 'Boolean')
@@ -171 +184 @@
         print('Temporal love numbers: Overriding md.love.nfreq and md.love.frequencies')
         self.nfreq = len(self.time) * 2 * self.n_temporal_iterations
-        self.frequencies = np.zeros((self.nfreq, 1))
+        self.frequencies = np.zeros((self.nfreq,))
         for i in range(len(self.time)):
             for j in range(2 * self.n_temporal_iterations):
-                self.frequencies[(i - 1) * 2 * self.n_temporal_iterations + j] = j * np.log(2) / self.time[i] / 2 / np.pi
+                if self.time[i] == 0:
+                    self.frequencies[(i - 1) * 2 * self.n_temporal_iterations + j] = 0 # Convention to avoid marshalling infinite numbers
+                else:
+                    self.frequencies[(i - 1) * 2 * self.n_temporal_iterations + j] = j * np.log(2) / self.time[i] / 2 / np.pi
+        return self
     #}}}