Repository: simpeg/simpeg
Pull request #1299 (2023-11-08 17:17:46)

Fix issue with lengthscales in coterminal angle calculations
#### Summary
The initial implementation of coterminal angles used an average of cell sizes to determine the cell centered length. Those distances are slightly different (by 88%) than the cell gradient distances used by the `regularization_mesh.cell_gradient_` operators on octree changes. 

#### PR Checklist
* [X] If this is a work in progress PR, set as a Draft PR
* [X] Linted my code according to the [style guides](https://docs.simpeg.xyz/content/getting_started/contributing/code-style.html).
* [X] Added [tests](https://docs.simpeg.xyz/content/getting_started/practices.html#testing) to verify changes to the code.
* [N/A] Added necessary documentation to any new functions/classes following the
      expect [style](https://docs.simpeg.xyz/content/getting_started/practices.html#documentation).
* [X] Marked as ready for review (if this is was a draft PR), and converted 
      to a Pull Request
* [X] Tagged ``@simpeg/simpeg-developers`` when ready for review.

#### What does this implement/fix?
Fixes large gradients computed on the -pi/pi transition, resulting in artifacts in the model.

Current results from `examples\03-magnetics\plot_inv_mag_MVI_Sparse_TreeMesh.py.`:

Title should be L2 Spherical model...

![Figure_3](https://github.com/simpeg/simpeg/assets/55204635/74f7740b-6056-4a61-a509-cacdce037131)

After fix:

![Figure_3_new](https://github.com/simpeg/simpeg/assets/55204635/56f619cb-8242-4622-bcdf-88f432ff57c6)



Patch:
diff --git a/SimPEG/regularization/base.py b/SimPEG/regularization/base.py
index 5d092c121..656c1a257 100644
--- a/SimPEG/regularization/base.py
+++ b/SimPEG/regularization/base.py
@@ -1007,14 +1007,14 @@ class SmoothnessFirstOrder(BaseRegularization):
         .. math::
             \phi_m (\mathbf{m}) = \frac{1}{2} \Big \| \mathbf{W \, f_m} \Big \|^2
         """
-        dfm_dl = self.cell_gradient @ (self.mapping * self._delta_m(m))
+        dfm_dl = self.mapping * self._delta_m(m)
 
         if self.units is not None and self.units.lower() == "radian":
             return (
-                utils.mat_utils.coterminal(dfm_dl * self._cell_distances)
+                utils.mat_utils.coterminal(self.cell_gradient.sign() @ dfm_dl)
                 / self._cell_distances
             )
-        return dfm_dl
+        return self.cell_gradient @ dfm_dl
 
     def f_m_deriv(self, m) -> csr_matrix:
         r"""Derivative of the regularization kernel function.
@@ -1274,17 +1274,17 @@ class SmoothnessSecondOrder(SmoothnessFirstOrder):
         .. math::
             \phi_m (\mathbf{m}) = \frac{1}{2} \Big \| \mathbf{W \, f_m} \Big \|^2
         """
-        dfm_dl = self.cell_gradient @ (self.mapping * self._delta_m(m))
+        dfm_dl = self.mapping * self._delta_m(m)
 
         if self.units is not None and self.units.lower() == "radian":
-            dfm_dl = (
-                utils.mat_utils.coterminal(dfm_dl * self.length_scales)
+            return self.cell_gradient.T @ (
+                utils.mat_utils.coterminal(self.cell_gradient.sign() @ dfm_dl)
                 / self.length_scales
             )
 
-        dfm_dl2 = self.cell_gradient.T @ dfm_dl
+        dfm_dl2 = self.cell_gradient @ dfm_dl
 
-        return dfm_dl2
+        return self.cell_gradient.T @ dfm_dl2
 
     def f_m_deriv(self, m) -> csr_matrix:
         r"""Derivative of the regularization kernel function.
diff --git a/SimPEG/regularization/regularization_mesh.py b/SimPEG/regularization/regularization_mesh.py
index ae02a73ec..c07f92504 100755
--- a/SimPEG/regularization/regularization_mesh.py
+++ b/SimPEG/regularization/regularization_mesh.py
@@ -551,8 +551,10 @@ class RegularizationMesh(props.BaseSimPEG):
             Cell center distance array along the x-direction.
         """
         if getattr(self, "_cell_distances_x", None) is None:
-            Ave = self.aveCC2Fx
-            self._cell_distances_x = Ave * (self.Pac.T * self.mesh.h_gridded[:, 0])
+            self._cell_distances_x = self.cell_gradient_x.max(
+                axis=1
+            ).toarray().flatten() ** (-1.0)
+
         return self._cell_distances_x
 
     @property
@@ -565,8 +567,10 @@ class RegularizationMesh(props.BaseSimPEG):
             Cell center distance array along the y-direction.
         """
         if getattr(self, "_cell_distances_y", None) is None:
-            Ave = self.aveCC2Fy
-            self._cell_distances_y = Ave * (self.Pac.T * self.mesh.h_gridded[:, 1])
+            self._cell_distances_y = self.cell_gradient_y.max(
+                axis=1
+            ).toarray().flatten() ** (-1.0)
+
         return self._cell_distances_y
 
     @property
@@ -579,8 +583,10 @@ class RegularizationMesh(props.BaseSimPEG):
             Cell center distance array along the z-direction.
         """
         if getattr(self, "_cell_distances_z", None) is None:
-            Ave = self.aveCC2Fz
-            self._cell_distances_z = Ave * (self.Pac.T * self.mesh.h_gridded[:, 2])
+            self._cell_distances_z = self.cell_gradient_z.max(
+                axis=1
+            ).toarray().flatten() ** (-1.0)
+
         return self._cell_distances_z
 
 
diff --git a/SimPEG/utils/mat_utils.py b/SimPEG/utils/mat_utils.py
index cbf9221d3..5d128d997 100644
--- a/SimPEG/utils/mat_utils.py
+++ b/SimPEG/utils/mat_utils.py
@@ -409,7 +409,6 @@ def coterminal(theta):
         Coterminal angles
 
     """
-
     sub = theta[np.abs(theta) >= np.pi]
     sub = -np.sign(sub) * (2 * np.pi - np.abs(sub))
 


Test patch:
diff --git a/tests/base/test_regularization.py b/tests/base/test_regularization.py
index 43bcc168a..82abaca79 100644
--- a/tests/base/test_regularization.py
+++ b/tests/base/test_regularization.py
@@ -638,6 +638,22 @@ def test_cross_reg_reg_errors():
         regularization.CrossReferenceRegularization(mesh, ref_dir)
 
 
+@pytest.mark.parametrize("orientation", ("x", "y", "z"))
+def test_smoothness_first_order_coterminal_angle(orientation):
+    """
+    Test smoothness first order regularizations of angles on a treemesh
+    """
+    mesh = discretize.TreeMesh([16, 16, 16])
+    mesh.insert_cells([100, 100, 100], mesh.max_level, finalize=True)
+
+    reg = regularization.SmoothnessFirstOrder(
+        mesh, units="radian", orientation=orientation
+    )
+    angles = np.ones(mesh.n_cells) * np.pi
+    angles[5] = -np.pi
+    assert np.all(reg.f_m(angles) == 0)
+
+
 class TestParent:
     """Test parent property of regularizations."""