Changeset 25303 for issm/trunk-jpl/src/m/solve/parseresultsfromdisk.py

Timestamp:

07/27/20 12:11:11 (5 years ago)

Author:

jdquinn

Message:

BUG: Array addressing when parsing results (Python); test translation; cleanup

File:

• issm/trunk-jpl/src/m/solve/parseresultsfromdisk.py (modified) (2 diffs)

issm/trunk-jpl/src/m/solve/parseresultsfromdisk.py ¶

@@ -248 +248 @@
             rr = md.materials.radius
             rho = md.materials.density
-            rho_avg = (rho * np.diff(np.power(rr, 3), n=1, axis=0)) / np.diff(np.power(rr, 3).sum()).sum()
+            rho_avg_partial = np.diff(np.power(rr, 3), n=1, axis=0)
+            rho_avg = ((rho * rho_avg_partial) / rho_avg_partial.sum()).sum()
             temp_field = np.empty((degmax + 1, nfreq, nlayer + 1, 6))
             temp_field.fill(0.0)
@@ -254 +255 @@
                 for jj in range(nfreq):
                     for kk in range(nlayer + 1):
-                        if kk < nlayer + 1:
+                        if kk < nlayer: # NOTE: Upper bound of range is non-inclusive (compare to src/m/solve/parseresultsfromdisk.m)
                             ll = ii * (nlayer + 1) * 6 + (kk * 6 + 1) + 3
-                            temp_field[ii, jj, kk, 0] = field[ll + (1 - 1), jj] * r0        # mm = 4
-                            temp_field[ii, jj, kk, 1] = field[ll + (2 - 1), jj] * mu0       # mm = 5
-                            temp_field[ii, jj, kk, 2] = field[ll + (3 - 1), jj] * r0        # mm = 6
-                            temp_field[ii, jj, kk, 3] = field[ll + (4 - 1), jj] * mu0       # mm = 1
-                            temp_field[ii, jj, kk, 4] = field[ll + (5 - 1), jj] * r0 * g0   # mm = 2
-                            temp_field[ii, jj, kk, 5] = field[ll + (6 - 1), jj] * g0        # mm = 3
-                            print(temp_field)
+                            temp_field[ii, jj, kk, 0] = field[ll + (0 - 1), jj] * r0        # mm = 4
+                            temp_field[ii, jj, kk, 1] = field[ll + (1 - 1), jj] * mu0       # mm = 5
+                            temp_field[ii, jj, kk, 2] = field[ll + (2 - 1), jj] * r0        # mm = 6
+                            temp_field[ii, jj, kk, 3] = field[ll + (3 - 1), jj] * mu0       # mm = 1
+                            temp_field[ii, jj, kk, 4] = field[ll + (4 - 1), jj] * r0 * g0   # mm = 2
+                            temp_field[ii, jj, kk, 5] = field[ll + (5 - 1), jj] * g0        # mm = 3
                         else: # surface
-                            ll = ii * (nlayer + 1) * 6 - 2
-                            temp_field[ii, jj, kk, 0] = field[ll + (1 - 1), jj] * r0
-                            temp_field[ii, jj, kk, 2] = field[ll + (2 - 1), jj] * r0
-                            temp_field[ii, jj, kk, 4] = field[ll + (3 - 1), jj] * r0 * g0
+                            ll = (ii + 1) * (nlayer + 1) * 6 - 2
+                            temp_field[ii, jj, kk, 0] = field[ll + (0 - 1), jj] * r0
+                            temp_field[ii, jj, kk, 2] = field[ll + (1 - 1), jj] * r0
+                            temp_field[ii, jj, kk, 4] = field[ll + (2 - 1), jj] * r0 * g0
                             # surface BC
                             temp_field[ii, jj, kk, 3] = 0
                             if md.love.forcing_type == 9:
                                 temp_field[ii, jj, kk, 1] = 0
-                                temp_field[ii, jj, kk, 5] = (2 * ii - 1) / r0 - ii * field[ll + (3 - 1), jj] * g0
+                                temp_field[ii, jj, kk, 5] = (2 * (ii + 1) - 1) / r0 - (ii + 1) * field[ll + (2 - 1), jj] * g0
                             elif md.love.forcing_type == 11:
-                                temp_field[ii, jj, kk, 1] = -(2 * (ii - 1) + 1) * rho_avg / 3
-                                temp_field[ii, jj, kk, 5] = (2 * ii - 1) / r0 - ii * field[ll + (3 - 1), jj] * g0
+                                temp_field[ii, jj, kk, 1] = -(2 * ii + 1) * rho_avg / 3
+                                temp_field[ii, jj, kk, 5] = (2 * (ii + 1) - 1) / r0 - (ii + 1) * field[ll + (2 - 1), jj] * g0
             field = temp_field