diff --git a/SimPEG/directives/directives.py b/SimPEG/directives/directives.py
index 59e0fe0929..3c25cfb7df 100644
--- a/SimPEG/directives/directives.py
+++ b/SimPEG/directives/directives.py
@@ -1401,9 +1401,8 @@ def angleScale(self):
 
         for reg, var in zip(self.reg.objfcts[1:], max_s):
             for obj in reg.objfcts:
-                obj.add_set_weights(
-                    {"angle_scale": np.ones(obj.shape[0]) * max_p / var}
-                )
+                # TODO Need to make weights_shapes a public method
+                obj.set_weights(angle_scale=np.ones(obj._weights_shapes[0]) * max_p / var)
 
     def validate(self, directiveList):
         # check if a linear preconditioner is in the list, if not warn else
@@ -1587,7 +1586,7 @@ def update(self):
         for reg in self.reg.objfcts:
             if not isinstance(reg, BaseSimilarityMeasure):
                 for sub_reg in reg.objfcts:
-                    sub_reg.add_set_weights({"sensitivity": sub_reg.mapping * wr})
+                    sub_reg.set_weights(sensitivity=sub_reg.mapping * wr)
 
     def validate(self, directiveList):
         # check if a beta estimator is in the list after setting the weights
diff --git a/SimPEG/objective_function.py b/SimPEG/objective_function.py
index 09795e5352..19727db563 100644
--- a/SimPEG/objective_function.py
+++ b/SimPEG/objective_function.py
@@ -293,10 +293,6 @@ def __len__(self):
     def __getitem__(self, key):
         return self.multipliers[key], self.objfcts[key]
 
-    @property
-    def __len__(self):
-        return self.objfcts.__len__
-
     @property
     def multipliers(self):
         return self._multipliers
@@ -312,7 +308,7 @@ def multipliers(self, value):
 
         assert len(value) == len(self.objfcts), (
             "the length of multipliers should be the same as the number of"
-            " objective functions ({}), not {}".format(len(self.objfcts, len(value)))
+            " objective functions ({}), not {}".format(len(self.objfcts), len(value))
         )
 
         self._multipliers = value
diff --git a/SimPEG/regularization/__init__.py b/SimPEG/regularization/__init__.py
index 73b3f80256..a9a5a6c3b5 100644
--- a/SimPEG/regularization/__init__.py
+++ b/SimPEG/regularization/__init__.py
@@ -29,18 +29,36 @@ def __init__(self, mesh=None, **kwargs):
 
 @deprecate_class(removal_version="0.x.0", future_warn=True)
 class Simple(LeastSquaresRegularization):
-    def __init__(self, mesh=None, **kwargs):
-        super().__init__(mesh=mesh, **kwargs)
+    def __init__(self, mesh=None, alpha_x=1.0, alpha_y=1.0, alpha_z=1.0, **kwargs):
+        # These alphas are now refered to as length_scalse in the
+        # new LeastSquaresRegularization
+        super().__init__(
+            mesh=mesh,
+            length_scale_x=alpha_x,
+            length_scale_y=alpha_y,
+            length_scale_z=alpha_z,
+            **kwargs
+        )
 
 
 @deprecate_class(removal_version="0.x.0", future_warn=True)
 class Tikhonov(LeastSquaresRegularization):
-    pass
+    def __init__(
+        self, mesh=None, alpha_s=1e-6, alpha_x=1.0, alpha_y=1.0, alpha_z=1.0, **kwargs
+    ):
+        super().__init__(
+            mesh=mesh,
+            alpha_s=alpha_s,
+            alpha_x=alpha_x,
+            alpha_y=alpha_y,
+            alpha_z=alpha_z,
+            **kwargs
+        )
 
 
 @deprecate_class(removal_version="0.x.0", future_warn=True)
 class PGIwithNonlinearRelationshipsSmallness(PGIsmallness):
-    def __init__(self, gmm):
+    def __init__(self, gmm, **kwargs):
         super().__init__(gmm, non_linear_relationships=True, **kwargs)
 
 
diff --git a/SimPEG/regularization/base.py b/SimPEG/regularization/base.py
index 42ca985db6..c9a93e4ed6 100644
--- a/SimPEG/regularization/base.py
+++ b/SimPEG/regularization/base.py
@@ -2,6 +2,7 @@
 
 import numpy as np
 from discretize.base import BaseMesh
+import warnings
 
 from .. import maps
 from ..objective_function import BaseObjectiveFunction, ComboObjectiveFunction
@@ -22,16 +23,7 @@ class BaseRegularization(BaseObjectiveFunction):
 
     """
 
-    _active_cells = None
-    _mapping = None
     _model = None
-    _mesh: RegularizationMesh | None = None
-    _reference_model = None
-    _regularization_mesh = None
-    _shape = None
-    _units = None
-    _weights: {} = None
-    _W = None
 
     def __init__(
         self,
@@ -43,33 +35,55 @@ def __init__(
         weights=None,
         **kwargs,
     ):
-        self.regularization_mesh = mesh
-        self.active_cells = active_cells
+        if isinstance(mesh, BaseMesh):
+            mesh = RegularizationMesh(mesh)
+
+        if not isinstance(mesh, RegularizationMesh):
+            TypeError(
+                f"'regularization_mesh' must be of type {RegularizationMesh} or {BaseMesh}. "
+                f"Value of type {type(mesh)} provided."
+            )
+        self._regularization_mesh = mesh
+        self._weights = {}
+        if active_cells is not None:
+            self.active_cells = active_cells
         self.mapping = mapping
 
         super().__init__(**kwargs)
 
         self.reference_model = reference_model
         self.units = units
-        self.weights = weights
+        if weights is not None:
+            if not isinstance(weights, dict):
+                weights = {"user_weights": weights}
+            self.set_weights(**weights)
 
     # Properties
     @property
     def active_cells(self) -> np.ndarray:
-        """Indices of active cells in the mesh"""
-        return self._active_cells
+        """A boolean array of active cells on the regularization
+
+        Returns
+        -------
+        (n_cells, ) numpy.ndarray of bool
+
+        Notes
+        -----
+        If this is set with an array of integers, it interprets it as an array
+        listing the active cell indices.
+        """
+        return self.regularization_mesh.active_cells
 
     @active_cells.setter
     def active_cells(self, values: np.ndarray | None):
-        if self._active_cells is not None:
-            raise AttributeError("Cannot reset the 'active_cells' property.")
-
-        if (
-            isinstance(self.regularization_mesh, RegularizationMesh)
-            and getattr(self.regularization_mesh, "active_cells", None) is None
-        ):
-            self.regularization_mesh.active_cells = values
-        self._active_cells = values
+        self.regularization_mesh.active_cells = values
+        # remove any weights, and reset the volume weight if present
+        if values is not None:
+            volume_term = "volume" in self._weights
+            self._weights = {}
+            self._W = None
+            if volume_term:
+                self.set_weights(volume=self.regularization_mesh.vol)
 
     indActive = deprecate_property(
         active_cells,
@@ -96,18 +110,17 @@ def model(self, values: np.ndarray | float):
     @property
     def mapping(self) -> maps.IdentityMap:
         """Mapping applied to the model values"""
-        if getattr(self, "_mapping", None) is None:
-            self.mapping = maps.IdentityMap(nP=self._nC_residual)
         return self._mapping
 
     @mapping.setter
     def mapping(self, mapping: maps.IdentityMap):
-        if not isinstance(mapping, (maps.IdentityMap, type(None))):
+        if mapping is None:
+            mapping = maps.IdentityMap()
+        if not isinstance(mapping, maps.IdentityMap):
             raise TypeError(
                 f"'mapping' must be of type {maps.IdentityMap}. "
                 f"Value of type {type(mapping)} provided."
             )
-
         self._mapping = mapping
 
     @property
@@ -125,9 +138,13 @@ def units(self, units: str | None):
         self._units = units
 
     @property
-    def shape(self):
-        """
-        number of model parameters
+    def _weights_shapes(self):
+        """Acceptable lengths for the weights
+
+        Returns
+        -------
+        list of tuple
+            Each tuple represents accetable shapes for the weights
         """
         if (
             getattr(self, "_regularization_mesh", None) is not None
@@ -135,7 +152,7 @@ def shape(self):
         ):
             return (self.regularization_mesh.nC,)
         elif getattr(self, "_mapping", None) is not None and self.mapping.shape != "*":
-            return self.mapping.shape
+            return (self.mapping.shape[0],)
         else:
             return "*"
 
@@ -146,12 +163,12 @@ def reference_model(self) -> np.ndarray:
 
     @reference_model.setter
     def reference_model(self, values: np.ndarray | float):
+        if values is not None:
+            if isinstance(values, float):
+                values = np.ones(self._nC_residual) * values
 
-        if isinstance(values, float):
-            values = np.ones(self._nC_residual) * values
-
-        self.validate_array_type("reference_model", values, float)
-        self.validate_shape("reference_model", values, (self._nC_residual,))
+            self.validate_array_type("reference_model", values, float)
+            self.validate_shape("reference_model", values, (self._nC_residual,))
         self._reference_model = values
 
     mref = deprecate_property(
@@ -166,13 +183,6 @@ def regularization_mesh(self) -> RegularizationMesh:
         """Regularization mesh"""
         return self._regularization_mesh
 
-    @regularization_mesh.setter
-    def regularization_mesh(self, mesh: RegularizationMesh):
-        if not isinstance(mesh, RegularizationMesh):
-            mesh = RegularizationMesh(mesh)
-
-        self._regularization_mesh = mesh
-
     regmesh = deprecate_property(
         regularization_mesh,
         "regmesh",
@@ -181,41 +191,63 @@ def regularization_mesh(self, mesh: RegularizationMesh):
     )
 
     @property
-    def weights(self):
-        """Regularization weights applied to the target elements"""
-        if getattr(self, "_weights", None) is None:
-            self._weights = {}
-            self.add_set_weights({"volume": self.regularization_mesh.vol})
-
-        return self._weights
-
-    @weights.setter
-    def weights(self, weights: dict[str, np.ndarray] | np.ndarray | None):
-        self._weights = None
-        if weights is not None:
-            self.add_set_weights(weights)
-
-    cell_weights = deprecate_property(
-        weights,
-        "cell_weights",
-        "0.x.0",
-        error=False,
-    )
-
-    def add_set_weights(self, weights: dict | np.ndarray):
-        if isinstance(weights, np.ndarray):
-            weights = {"user_weights": weights}
+    def cell_weights(self):
+        warnings.warn(
+            "cell_weights are deprecated please access weights using the `set_weights`,"
+            " `get_weights`, and `remove_weights` functionality. This will be removed in 0.x.0",
+            DeprecationWarning,
+        )
+        return np.prod(list(self._weights.values()), axis=0)
+
+    @cell_weights.setter
+    def cell_weights(self, value):
+        warnings.warn(
+            "cell_weights are deprecated please access weights using the `set_weights`,"
+            " `get_weights`, and `remove_weights` functionality. This will be removed in 0.x.0",
+            DeprecationWarning,
+        )
+        self.set_weights(cell_weights=value)
 
-        if not isinstance(weights, dict):
-            raise TypeError(
-                "Weights must be provided as a dictionary or numpy.ndarray."
-            )
+    def get_weights(self, key):
+        """Returns the weights with a given key
 
+        Returns
+        -------
+        numpy.ndarray
+        """
+        return self._weights[key]
+
+    def set_weights(self, **weights):
+        """Adds (or updates) the specified weights to the regularization
+
+        Parameters:
+        -----------
+        **kwargs : key, numpy.ndarray
+            Each keyword argument is added to the weights used by the regularization.
+            They can be accessed with their keyword argument.
+
+        Examples
+        --------
+        >>> import discretize
+        >>> from SimPEG.regularization import Small
+        >>> mesh = discretize.TensorMesh([2, 3, 2])
+        >>> reg = Small(mesh)
+        >>> reg.set_weights(my_weight=np.ones(mesh.n_cells))
+        >>> reg.get_weights('my_weight')
+        array([1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.])
+        """
         for key, values in weights.items():
             self.validate_array_type("weights", values, float)
-            self.validate_shape("weights", values, (self.shape[0],))
-            self.weights[key] = values
+            self.validate_shape("weights", values, self._weights_shapes)
+            self._weights[key] = values
+        self._W = None
 
+    def remote_weights(self, key):
+        """Removes the weights with a given key"""
+        try:
+            self._weights.pop(key)
+        except KeyError:
+            raise KeyError(f"{key} is not in the weights dictionary")
         self._W = None
 
     @property
@@ -224,9 +256,8 @@ def W(self):
         Weighting matrix
         """
         if getattr(self, "_W", None) is None:
-            weights = np.prod(list(self.weights.values()), axis=0)
+            weights = np.prod(list(self._weights.values()), axis=0)
             self._W = utils.sdiag(weights ** 0.5)
-
         return self._W
 
     @property
@@ -243,7 +274,7 @@ def _nC_residual(self):
         elif nC != "*" and nC is not None:
             return self.regularization_mesh.nC
         else:
-            return self.shape[0]
+            return self._weights_shapes[0]
 
     def _delta_m(self, m):
         if self.reference_model is None:
@@ -307,7 +338,7 @@ def deriv2(self, m, v=None):
 
         .. math::
 
-            R(m) = \\frac{1}{2}\mathbf{(m-m_\\text{ref})^\\top W^\\top
+            R(m) = \\frac{1}{2}\\mathbf{(m-m_\\text{ref})^\\top W^\\top
             W(m-m_\\text{ref})}
 
         So the second derivative is straight forward:
@@ -323,6 +354,7 @@ def deriv2(self, m, v=None):
 
         return f_m_deriv.T * (self.W.T * (self.W * (f_m_deriv * v)))
 
+    # TODO move these to general validation functions in code_utils
     def validate_array_type(self, attribute, array, dtype):
         """Generic array and type validator"""
         if (
@@ -335,6 +367,7 @@ def validate_array_type(self, attribute, array, dtype):
                 f"Values of type {type(array)} provided."
             )
 
+    # TODO move these to general validation functions in code_utils
     def validate_shape(self, attribute, values, shape: tuple | tuple[tuple]):
         """Generic array shape validator"""
         if (
@@ -379,8 +412,8 @@ class Small(BaseRegularization):
     _multiplier_pair = "alpha_s"
 
     def __init__(self, mesh, **kwargs):
-
         super().__init__(mesh, **kwargs)
+        self.set_weights(volume=self.regularization_mesh.vol)
 
     def f_m(self, m):
         """
@@ -403,7 +436,6 @@ class SmoothDeriv(BaseRegularization):
     **Optional Inputs**
 
     :param discretize.base.BaseMesh mesh: SimPEG mesh
-    :param int shape: number of parameters
     :param IdentityMap mapping: regularization mapping, takes the model from model space to the space you want to regularize in
     :param numpy.ndarray reference_model: reference model
     :param numpy.ndarray active_cells: active cell indices for reducing the size of differential operators in the definition of a regularization mesh
@@ -412,35 +444,51 @@ class SmoothDeriv(BaseRegularization):
     :param numpy.ndarray weights: vector of cell weights (applied in all terms)
     """
 
-    _cell_difference = None
-    _length_scales = None
-    _orientation = None
-    _shape = None
-    _reference_model_in_smooth: bool = False
-
     def __init__(
         self, mesh, orientation="x", reference_model_in_smooth=False, **kwargs
     ):
         self.reference_model_in_smooth = reference_model_in_smooth
 
-        super().__init__(mesh=mesh, **kwargs)
+        if orientation not in ["x", "y", "z"]:
+            raise KeyError("Orientation must be 'x', 'y' or 'z'")
 
-        self.orientation = orientation
+        if orientation == "y" and mesh.dim < 2:
+            raise ValueError(
+                "Mesh must have at least 2 dimensions to regularize along the "
+                "y-direction."
+            )
+        elif orientation == "z" and mesh.dim < 3:
+            raise ValueError(
+                "Mesh must have at least 3 dimensions to regularize along the "
+                "z-direction"
+            )
+        self._orientation = orientation
+
+        super().__init__(mesh=mesh, **kwargs)
+        self.set_weights(volume=self.regularization_mesh.vol)
 
     @property
-    def shape(self):
-        return getattr(
+    def _weights_shapes(self):
+        """Acceptable lengths for the weights
+
+        Returns
+        -------
+        tuple
+            A tuple of each acceptable lengths for the weights
+        """
+        n_active_f, n_active_c = getattr(
             self.regularization_mesh, "aveCC2F{}".format(self.orientation)
         ).shape
+        return [(n_active_f,), (n_active_c,)]
 
     @property
-    def cell_difference(self):
-        """Cell difference operator"""
-        if getattr(self, "_cell_difference", None) is None:
-            self._cell_difference = getattr(
-                self.regularization_mesh, "cellDiff{}".format(self.orientation)
+    def cell_gradient(self):
+        """Cell gradient operator"""
+        if getattr(self, "_cell_gradient", None) is None:
+            self._cell_gradient = getattr(
+                self.regularization_mesh, "cell_gradient_{}".format(self.orientation)
             )
-        return self._cell_difference
+        return self._cell_gradient
 
     @property
     def reference_model_in_smooth(self) -> bool:
@@ -461,9 +509,7 @@ def reference_model_in_smooth(self, value: bool):
     def _delta_m(self, m):
         if self.reference_model is None or not self.reference_model_in_smooth:
             return m
-        return (
-            m - self.reference_model
-        )  # in case self.reference_model is Zero, returns type m
+        return m - self.reference_model
 
     @property
     def _multiplier_pair(self):
@@ -473,12 +519,12 @@ def f_m(self, m):
         """
         Model gradient
         """
-        dfm_dl = self.cell_difference @ (self.mapping * self._delta_m(m))
+        dfm_dl = self.cell_gradient @ (self.mapping * self._delta_m(m))
 
         if self.units == "radian":
             return (
-                utils.mat_utils.coterminal(dfm_dl * self.length_scales)
-                / self.length_scales
+                utils.mat_utils.coterminal(dfm_dl * self._cell_distances)
+                / self._cell_distances
             )
         return dfm_dl
 
@@ -486,49 +532,7 @@ def f_m_deriv(self, m):
         """
         Derivative of the model gradient
         """
-        return self.cell_difference @ self.mapping.deriv(self._delta_m(m))
-
-    @property
-    def weights(self):
-        """Regularization weights applied to the target elements"""
-        if getattr(self, "_weights", None) is None:
-            self._weights = {}
-            self.add_set_weights(
-                {
-                    "volume": self.regularization_mesh.vol,
-                }
-            )
-
-        return self._weights
-
-    @weights.setter
-    def weights(self, weights: dict[str, np.ndarray] | np.ndarray | None):
-        self._weights = None
-        if weights is not None:
-            self.add_set_weights(weights)
-
-    cell_weights = deprecate_property(
-        weights,
-        "cell_weights",
-        "0.x.0",
-        error=False,
-    )
-
-    def add_set_weights(self, weights: dict):
-        if isinstance(weights, np.ndarray):
-            weights = {"user_weights": weights}
-
-        if not isinstance(weights, dict):
-            raise TypeError(
-                "Weights must be provided as a dictionary or numpy.ndarray."
-            )
-
-        for key, values in weights.items():
-            self.validate_array_type("weights", values, float)
-            self.validate_shape("weights", values, ((self.shape[0],), (self.shape[1],)))
-            self.weights[key] = values
-
-        self._W = None
+        return self.cell_gradient @ self.mapping.deriv(self._delta_m(m))
 
     @property
     def W(self):
@@ -541,63 +545,24 @@ def W(self):
                 self.regularization_mesh, "aveCC2F{}".format(self.orientation)
             )
             weights = 1.0
-            for values in self.weights.values():
+            for values in self._weights.values():
                 if values.shape[0] == self.regularization_mesh.nC:
                     values = average_cell_2_face * values
                 weights *= values
             self._W = utils.sdiag(weights ** 0.5)
-
         return self._W
 
     @property
-    def length_scales(self):
+    def _cell_distances(self):
         """
-        Length scales for the cell center difference.
-        Normalized by the smallest cell dimension if 'normalized_gradient' is True.
+        Distances between cell centers for the cell center difference.
         """
-        if getattr(self, "_length_scales", None) is None:
-            Ave = getattr(
-                self.regularization_mesh, "aveCC2F{}".format(self.orientation)
-            )
-            index = "xyz".index(self.orientation)
-            self._length_scales = Ave * (
-                self.regularization_mesh.Pac.T
-                * self.regularization_mesh.mesh.h_gridded[:, index]
-            )
-
-        return self._length_scales
-
-    @length_scales.setter
-    def length_scales(self, values):
-        self.validate_array_type("length_scales", values, float)
-        self.validate_shape("length_scales", values, (self.nP,))
-        self._length_scales = values
+        return getattr(self.regularization_mesh, f"cell_distances_{self.orientation}")
 
     @property
     def orientation(self):
         return self._orientation
 
-    @orientation.setter
-    def orientation(self, value):
-        assert value in [
-            "x",
-            "y",
-            "z",
-        ], "Orientation must be 'x', 'y' or 'z'"
-
-        if value == "y":
-            assert self.regularization_mesh.dim > 1, (
-                "Mesh must have at least 2 dimensions to regularize along the "
-                "y-direction"
-            )
-
-        elif value == "z":
-            assert self.regularization_mesh.dim > 2, (
-                "Mesh must have at least 3 dimensions to regularize along the "
-                "z-direction"
-            )
-        self._orientation = value
-
 
 class SmoothDeriv2(SmoothDeriv):
     """
@@ -616,14 +581,11 @@ class SmoothDeriv2(SmoothDeriv):
     :param numpy.ndarray weights: vector of cell weights (applied in all terms)
     """
 
-    def __init__(self, mesh, orientation="x", **kwargs):
-        super().__init__(mesh=mesh, orientation=orientation, **kwargs)
-
     def f_m(self, m):
         """
         Second model derivative
         """
-        dfm_dl = self.cell_difference @ (self.mapping * self._delta_m(m))
+        dfm_dl = self.cell_gradient @ (self.mapping * self._delta_m(m))
 
         if self.units == "radian":
             dfm_dl = (
@@ -631,7 +593,7 @@ def f_m(self, m):
                 / self.length_scales
             )
 
-        dfm_dl2 = self.cell_difference.T @ dfm_dl
+        dfm_dl2 = self.cell_gradient.T @ dfm_dl
 
         return dfm_dl2
 
@@ -640,8 +602,8 @@ def f_m_deriv(self, m):
         Derivative of the second model residual
         """
         return (
-            self.cell_difference.T
-            * self.cell_difference
+            self.cell_gradient.T
+            @ self.cell_gradient
             @ self.mapping.deriv(self._delta_m(m))
         )
 
@@ -651,7 +613,7 @@ def W(self):
         Weighting matrix
         """
         if getattr(self, "_W", None) is None:
-            weights = np.prod(list(self.weights.values()), axis=0)
+            weights = np.prod(list(self._weights.values()), axis=0)
             self._W = utils.sdiag(weights ** 0.5)
 
         return self._W
@@ -669,260 +631,298 @@ def _multiplier_pair(self):
 
 
 class LeastSquaresRegularization(ComboObjectiveFunction):
-    _active_cells = None
-    _alpha_s = None
-    _alpha_x = None
-    _alpha_y = None
-    _alpha_z = None
-    _alpha_xx = None
-    _alpha_yy = None
-    _alpha_zz = None
-    _length_scale_x = None
-    _length_scale_y = None
-    _length_scale_z = None
     _model = None
-    _mapping = None
-    _normalized_gradients = True
-    _reference_model_in_smooth = False
-    _reference_model = None
-    _regularization_mesh = None
-    _units = None
-    _weights = None
 
     def __init__(
         self,
         mesh,
         active_cells=None,
-        alpha_s=None,
+        alpha_s=1.0,
         alpha_x=None,
         alpha_y=None,
         alpha_z=None,
-        alpha_xx=None,
-        alpha_yy=None,
-        alpha_zz=None,
+        alpha_xx=0.0,
+        alpha_yy=0.0,
+        alpha_zz=0.0,
         length_scale_x=None,
         length_scale_y=None,
         length_scale_z=None,
         mapping=None,
-        objfcts=None,
         reference_model=None,
         reference_model_in_smooth=False,
+        weights=None,
         **kwargs,
     ):
 
-        self.regularization_mesh = mesh
+        if isinstance(mesh, BaseMesh):
+            mesh = RegularizationMesh(mesh)
 
-        if objfcts is None:
+        if not isinstance(mesh, RegularizationMesh):
+            TypeError(
+                f"'regularization_mesh' must be of type {RegularizationMesh} or {BaseMesh}. "
+                f"Value of type {type(mesh)} provided."
+            )
+        self._regularization_mesh = mesh
+        if active_cells is not None:
+            self._regularization_mesh.active_cells = active_cells
+
+        self.alpha_s = alpha_s
+        if alpha_x is not None:
+            if length_scale_x is not None:
+                raise ValueError(
+                    "Attempted to set both alpha_x and length_scale_x at the same time. Please "
+                    "use only one of them"
+                )
+            self.alpha_x = alpha_x
+        else:
+            self.length_scale_x = length_scale_x
+
+        if alpha_y is not None:
+            if length_scale_y is not None:
+                raise ValueError(
+                    "Attempted to set both alpha_y and length_scale_y at the same time. Please "
+                    "use only one of them"
+                )
+            self.alpha_y = alpha_y
+        else:
+            self.length_scale_y = length_scale_y
+
+        if alpha_z is not None:
+            if length_scale_z is not None:
+                raise ValueError(
+                    "Attempted to set both alpha_y and length_scale_y at the same time. Please "
+                    "use only one of them"
+                )
+            self.alpha_z = alpha_z
+        else:
+            self.length_scale_z = length_scale_z
+
+        # do this to allow child classes to also pass a list of objfcts to this constructor
+        if "objfcts" not in kwargs:
             objfcts = [
                 Small(mesh=self.regularization_mesh),
                 SmoothDeriv(mesh=self.regularization_mesh, orientation="x"),
+                SmoothDeriv2(mesh=self.regularization_mesh, orientation="x"),
             ]
 
             if mesh.dim > 1:
-                objfcts.append(
-                    SmoothDeriv(mesh=self.regularization_mesh, orientation="y")
+                objfcts.extend(
+                    [
+                        SmoothDeriv(mesh=self.regularization_mesh, orientation="y"),
+                        SmoothDeriv2(mesh=self.regularization_mesh, orientation="y"),
+                    ]
                 )
 
             if mesh.dim > 2:
-                objfcts.append(
-                    SmoothDeriv(mesh=self.regularization_mesh, orientation="z")
+                objfcts.extend(
+                    [
+                        SmoothDeriv(mesh=self.regularization_mesh, orientation="z"),
+                        SmoothDeriv2(mesh=self.regularization_mesh, orientation="z"),
+                    ]
                 )
+        else:
+            objfcts = kwargs.pop("objfcts")
+        super().__init__(objfcts=objfcts, **kwargs)
+        self.mapping = mapping
+        self.reference_model = reference_model
+        self.reference_model_in_smooth = reference_model_in_smooth
+        self.alpha_xx = alpha_xx
+        self.alpha_yy = alpha_yy
+        self.alpha_zz = alpha_zz
+        if weights is not None:
+            if not isinstance(weights, dict):
+                weights = {"user_weights": weights}
+            self.set_weights(**weights)
 
-        super(LeastSquaresRegularization, self).__init__(
-            objfcts=objfcts,
-            active_cells=active_cells,
-            mapping=mapping,
-            reference_model=reference_model,
-            reference_model_in_smooth=reference_model_in_smooth,
-            alpha_s=alpha_s,
-            alpha_x=alpha_x,
-            alpha_y=alpha_y,
-            alpha_z=alpha_z,
-            alpha_xx=alpha_xx,
-            alpha_yy=alpha_yy,
-            alpha_zz=alpha_zz,
-            length_scale_x=length_scale_x,
-            length_scale_y=length_scale_y,
-            length_scale_z=length_scale_z,
-            **kwargs,
-        )
+    def set_weights(self, **weights):
+        """Update weights in children objective functions"""
+        for fct in self.objfcts:
+            fct.set_weights(**weights)
 
-    def add_set_weights(self, weights: dict):
-        """Update weights to objective functions"""
+    def remove_weights(self, key):
+        """removes weights in children objective functions"""
         for fct in self.objfcts:
-            fct.add_set_weights(weights)
+            fct.remove_weights(key)
+
+    @property
+    def cell_weights(self):
+        # All of the objective functions should have the same weights,
+        # so just grab the one from smallness here, which should also
+        # trigger the deprecation warning
+        return self.objfcts[0].cell_weights
+
+    @cell_weights.setter
+    def cell_weights(self, value):
+        warnings.warn(
+            "cell_weights are deprecated please access weights using the `set_weights`,"
+            " `get_weights`, and `remove_weights` functionality. This will be removed in 0.x.0",
+            DeprecationWarning,
+        )
+        self.set_weights(cell_weights=value)
 
     @property
     def alpha_s(self):
         """smallness weight"""
-        if getattr(self, "_alpha_s", None) is None:
-            self._alpha_s = 1.0
         return self._alpha_s
 
     @alpha_s.setter
     def alpha_s(self, value):
-        if isinstance(value, (float, int)) and value < 0:
-            raise ValueError(
-                "Input 'alpha_s' value must me of type float > 0"
-                f"Value {value} of type {type(value)} provided"
-            )
+        if value is None:
+            value = 1.0
+        try:
+            value = float(value)
+        except (ValueError, TypeError):
+            raise TypeError(f"alpha_s must be a real number, saw type{type(value)}")
+        if value < 0:
+            raise ValueError(f"alpha_s must be non-negative, not {value}")
         self._alpha_s = value
 
     @property
     def alpha_x(self):
         """weight for the first x-derivative"""
-        if getattr(self, "_alpha_x", None) is None:
-            self._alpha_x = (
-                self.length_scale_x * self.regularization_mesh.base_length
-            ) ** 2.0
         return self._alpha_x
 
     @alpha_x.setter
     def alpha_x(self, value):
-        if isinstance(value, (float, int)) and value < 0:
-            raise ValueError(
-                "Input 'alpha_x' value must me of type float > 0"
-                f"Value {value} of type {type(value)} provided"
-            )
+        try:
+            value = float(value)
+        except (ValueError, TypeError):
+            raise TypeError(f"alpha_x must be a real number, saw type{type(value)}")
+        if value < 0:
+            raise ValueError(f"alpha_x must be non-negative, not {value}")
         self._alpha_x = value
 
     @property
     def alpha_y(self):
         """weight for the first y-derivative"""
-        if getattr(self, "_alpha_y", None) is None:
-            self._alpha_y = (
-                self.length_scale_y * self.regularization_mesh.base_length
-            ) ** 2.0
         return self._alpha_y
 
     @alpha_y.setter
     def alpha_y(self, value):
-        if isinstance(value, (float, int)) and value < 0:
-            raise ValueError(
-                "Input 'alpha_y' value must me of type float > 0"
-                f"Value {value} of type {type(value)} provided"
-            )
+        try:
+            value = float(value)
+        except (ValueError, TypeError):
+            raise TypeError(f"alpha_y must be a real number, saw type{type(value)}")
+        if value < 0:
+            raise ValueError(f"alpha_y must be non-negative, not {value}")
         self._alpha_y = value
 
     @property
     def alpha_z(self):
         """weight for the first z-derivative"""
-        if getattr(self, "_alpha_z", None) is None:
-            self._alpha_z = (
-                self.length_scale_z * self.regularization_mesh.base_length
-            ) ** 2.0
         return self._alpha_z
 
     @alpha_z.setter
     def alpha_z(self, value):
-        if isinstance(value, (float, int)) and value < 0:
-            raise ValueError(
-                "Input 'alpha_z' value must me of type float > 0"
-                f"Value {value} of type {type(value)} provided"
-            )
+        try:
+            value = float(value)
+        except (ValueError, TypeError):
+            raise TypeError(f"alpha_z must be a real number, saw type{type(value)}")
+        if value < 0:
+            raise ValueError(f"alpha_z must be non-negative, not {value}")
         self._alpha_z = value
 
     @property
     def alpha_xx(self):
         """weight for the second x-derivative"""
-        if getattr(self, "_alpha_xx", None) is None:
-            self._alpha_xx = (
-                self.length_scale_x * self.regularization_mesh.base_length
-            ) ** 4.0
         return self._alpha_xx
 
     @alpha_xx.setter
     def alpha_xx(self, value):
-        if isinstance(value, (float, int)) and value < 0:
-            raise ValueError(
-                "Input 'alpha_xx' value must me of type float > 0"
-                f"Value {value} of type {type(value)} provided"
-            )
+        if value is None:
+            value = (self.length_scale_x * self.regularization_mesh.base_length) ** 4.0
+        try:
+            value = float(value)
+        except (ValueError, TypeError):
+            raise TypeError(f"alpha_xx must be a real number, saw type{type(value)}")
+        if value < 0:
+            raise ValueError(f"alpha_xx must be non-negative, not {value}")
         self._alpha_xx = value
 
     @property
     def alpha_yy(self):
         """weight for the second y-derivative"""
-        if getattr(self, "_alpha_yy", None) is None:
-            self._alpha_yy = (
-                self.length_scale_y * self.regularization_mesh.base_length
-            ) ** 4.0
         return self._alpha_yy
 
     @alpha_yy.setter
     def alpha_yy(self, value):
-        if isinstance(value, (float, int)) and value < 0:
-            raise ValueError(
-                "Input 'alpha_yy' value must me of type float > 0"
-                f"Value {value} of type {type(value)} provided"
-            )
+        if value is None:
+            value = (self.length_scale_y * self.regularization_mesh.base_length) ** 4.0
+        try:
+            value = float(value)
+        except (ValueError, TypeError):
+            raise TypeError(f"alpha_yy must be a real number, saw type{type(value)}")
+        if value < 0:
+            raise ValueError(f"alpha_yy must be non-negative, not {value}")
         self._alpha_yy = value
 
     @property
     def alpha_zz(self):
         """weight for the second z-derivative"""
-        if getattr(self, "_alpha_zz", None) is None:
-            self._alpha_zz = (
-                self.length_scale_z * self.regularization_mesh.base_length
-            ) ** 4.0
         return self._alpha_zz
 
     @alpha_zz.setter
     def alpha_zz(self, value):
-        if isinstance(value, (float, int)) and value < 0:
-            raise ValueError(
-                "Input 'alpha_zz' value must me of type float > 0"
-                f"Value {value} of type {type(value)} provided"
-            )
+        if value is None:
+            value = (self.length_scale_z * self.regularization_mesh.base_length) ** 4.0
+        try:
+            value = float(value)
+        except (ValueError, TypeError):
+            raise TypeError(f"alpha_zz must be a real number, saw type{type(value)}")
+        if value < 0:
+            raise ValueError(f"alpha_zz must be non-negative, not {value}")
         self._alpha_zz = value
 
     @property
     def length_scale_x(self):
         """Constant multiplier of the base length scale on model gradients along x."""
-        if getattr(self, "_length_scale_x", None) is None:
-            self._length_scale_x = 1.0
-        return self._length_scale_x
+        return np.sqrt(self.alpha_x / self.regularization_mesh.base_length)
 
     @length_scale_x.setter
     def length_scale_x(self, value: float):
-        if not isinstance(value, (float, int, type(None))):
-            raise ValueError(
-                "Input length scale should be of type 'float' or 'int'. "
-                f"Provided {value} of type {type(value)}."
+        if value is None:
+            value = 1.0
+        try:
+            value = float(value)
+        except (TypeError, ValueError):
+            raise TypeError(
+                f"length_scale_x must be a real number, saw type{type(value)}"
             )
-        self._length_scale_x = value
+        self.alpha_x = (value * self.regularization_mesh.base_length) ** 2
 
     @property
     def length_scale_y(self):
         """Constant multiplier of the base length scale on model gradients along y."""
-        if getattr(self, "_length_scale_y", None) is None:
-            self._length_scale_y = 1.0
-        return self._length_scale_y
+        return np.sqrt(self.alpha_y / self.regularization_mesh.base_length)
 
     @length_scale_y.setter
     def length_scale_y(self, value: float):
-        if not isinstance(value, (float, int, type(None))):
-            raise ValueError(
-                "Input length scale should be of type 'float' or 'int'. "
-                f"Provided {value} of type {type(value)}."
+        if value is None:
+            value = 1.0
+        try:
+            value = float(value)
+        except (TypeError, ValueError):
+            raise TypeError(
+                f"length_scale_y must be a real number, saw type{type(value)}"
             )
-        self._length_scale_y = value
+        self.alpha_y = (value * self.regularization_mesh.base_length) ** 2
 
     @property
     def length_scale_z(self):
         """Constant multiplier of the base length scale on model gradients along z."""
-        if getattr(self, "_length_scale_z", None) is None:
-            self._length_scale_z = 1.0
-        return self._length_scale_z
+        return np.sqrt(self.alpha_z / self.regularization_mesh.base_length)
 
     @length_scale_z.setter
     def length_scale_z(self, value: float):
-        if not isinstance(value, (float, int, type(None))):
-            raise ValueError(
-                "Input length scale should be of type 'float' or 'int'. "
-                f"Provided {value} of type {type(value)}."
+        if value is None:
+            value = 1.0
+        try:
+            value = float(value)
+        except (TypeError, ValueError):
+            raise TypeError(
+                f"length_scale_z must be a real number, saw type{type(value)}"
             )
-        self._length_scale_z = value
+        self.alpha_z = (value * self.regularization_mesh.base_length) ** 2
 
     @property
     def reference_model_in_smooth(self) -> bool:
@@ -943,23 +943,6 @@ def reference_model_in_smooth(self, value: bool):
             if getattr(fct, "reference_model_in_smooth", None) is not None:
                 fct.reference_model_in_smooth = value
 
-    @property
-    def weights(self):
-        """Fixed regularization weights applied at cell centers"""
-        return self._weights
-
-    @weights.setter
-    def weights(self, value: dict[str, np.ndarray] | np.ndarray | None):
-        for fct in self.objfcts:
-            fct.weights = value
-
-    cell_weights = deprecate_property(
-        weights,
-        "cell_weights",
-        "0.x.0",
-        error=False,
-    )
-
     # Other properties and methods
     @property
     def nP(self):
@@ -1005,38 +988,18 @@ def multipliers(self):
         """
         return [getattr(self, objfct._multiplier_pair) for objfct in self.objfcts]
 
-    @property
-    def normalized_gradients(self):
-        """
-        Pre-normalize the model gradients by the base cell size
-        """
-        return self._normalized_gradients
-
-    @normalized_gradients.setter
-    def normalized_gradients(self, value):
-        if not isinstance(value, bool):
-            raise TypeError(
-                "'normalized_gradients must be of type 'bool'. "
-                f"Value of type {type(value)} provided."
-            )
-        self._normalized_gradients = value
-        for fct in self.objfcts:
-            if hasattr(fct, "_normalized_gradients"):
-                fct.normalized_gradients = value
-
     @property
     def active_cells(self) -> np.ndarray:
         """Indices of active cells in the mesh"""
-        return self._active_cells
+        return self.regularization_mesh.active_cells
 
     @active_cells.setter
     def active_cells(self, values: np.ndarray):
         self.regularization_mesh.active_cells = values
-
+        active_cells = self.regularization_mesh.active_cells
+        # notify the objtecive functions that the active_cells changed
         for objfct in self.objfcts:
-            objfct._active_cells = values
-
-        self._active_cells = values
+            objfct.active_cells = active_cells
 
     indActive = deprecate_property(
         active_cells,
@@ -1098,39 +1061,28 @@ def units(self, units: str | None):
             )
         for fct in self.objfcts:
             fct.model = units
+        self._units = units
 
     @property
     def regularization_mesh(self) -> RegularizationMesh:
         """Regularization mesh"""
         return self._regularization_mesh
 
-    @regularization_mesh.setter
-    def regularization_mesh(self, mesh: RegularizationMesh | BaseMesh):
-        if isinstance(mesh, BaseMesh):
-            mesh = RegularizationMesh(mesh)
-
-        if not isinstance(mesh, RegularizationMesh):
-            TypeError(
-                f"'regularization_mesh' must be of type {RegularizationMesh} or {BaseMesh}. "
-                f"Value of type {type(mesh)} provided."
-            )
-        self._regularization_mesh = mesh
-
     @property
     def mapping(self) -> maps.IdentityMap:
         """Mapping applied to the model values"""
-        if getattr(self, "_mapping", None) is None:
-            self.mapping = maps.IdentityMap(nP=self._nC_residual)
         return self._mapping
 
     @mapping.setter
     def mapping(self, mapping: maps.IdentityMap):
-        if not isinstance(mapping, (maps.IdentityMap, type(None))):
+        if mapping is None:
+            mapping = maps.IdentityMap(nP=self._nC_residual)
+
+        if not isinstance(mapping, maps.IdentityMap):
             raise TypeError(
                 f"'mapping' must be of type {maps.IdentityMap}. "
                 f"Value of type {type(mapping)} provided."
             )
-
         self._mapping = mapping
 
         for fct in self.objfcts:
@@ -1157,10 +1109,9 @@ class BaseSimilarityMeasure(BaseRegularization):
     three different models.
     """
 
-    _wire_map: maps.Wires = None
-
     def __init__(self, mesh, wire_map, **kwargs):
-        super().__init__(mesh, wire_map=wire_map, **kwargs)
+        super().__init__(mesh, **kwargs)
+        self.wire_map = wire_map
 
     @property
     def wire_map(self):
diff --git a/SimPEG/regularization/cross_gradient.py b/SimPEG/regularization/cross_gradient.py
index 63e9496bd5..4900f8dc48 100644
--- a/SimPEG/regularization/cross_gradient.py
+++ b/SimPEG/regularization/cross_gradient.py
@@ -24,9 +24,6 @@ class CrossGradient(BaseSimilarityMeasure):
 
     """
 
-    # reset this here to clear out the properties attribute
-    weights = None
-
     # These are not fully implemented yet
     # grad_tol = properties.Float(
     #     "tolerance for avoiding the exteremly small gradient amplitude", default=1e-10
diff --git a/SimPEG/regularization/pgi.py b/SimPEG/regularization/pgi.py
index 46c66c6581..c4f9cea8f3 100644
--- a/SimPEG/regularization/pgi.py
+++ b/SimPEG/regularization/pgi.py
@@ -49,7 +49,6 @@ class PGIsmallness(Small):
     """
 
     _multiplier_pair = "alpha_pgi"
-    _non_linear_relationships = False
     _maplist = None
     _wiresmap = None
 
@@ -82,20 +81,20 @@ def __init__(
             )
             kwargs.pop("mapping")
 
+        weights = kwargs.pop("weights", None)
+
         super().__init__(mesh=mesh, mapping=IdentityMap(nP=self.shape[0]), **kwargs)
 
         # Save repetitive computations (see withmapping implementation)
         self._r_first_deriv = None
         self._r_second_deriv = None
 
-    def add_set_weights(self, weights: dict | np.ndarray):
-        if isinstance(weights, (np.ndarray, list)):
-            weights = {"user_weights": np.r_[weights].flatten()}
+        if weights is not None:
+            if isinstance(weights, (np.ndarray, list)):
+                weights = {"user_weights": np.r_[weights].flatten()}
+            self.set_weights(**weights)
 
-        if not isinstance(weights, dict):
-            raise TypeError(
-                "Weights must be provided as a dictionary or numpy.ndarray."
-            )
+    def set_weights(self, **weights):
 
         for key, values in weights.items():
             self.validate_array_type("weights", values, float)
@@ -105,7 +104,7 @@ def add_set_weights(self, weights: dict | np.ndarray):
 
             self.validate_shape("weights", values, (self._nC_residual,))
 
-            self.weights[key] = values
+            self._weights[key] = values
 
         self._W = None
 
@@ -571,7 +570,10 @@ def deriv2(self, m, v=None):
             else:
                 # Forming the Hessian by diagonal blocks
                 hlist = [
-                    [self._r_second_deriv[:, i, j] for i in range(len(self.wiresmap.maps))]
+                    [
+                        self._r_second_deriv[:, i, j]
+                        for i in range(len(self.wiresmap.maps))
+                    ]
                     for j in range(len(self.wiresmap.maps))
                 ]
                 Hr = sp.csc_matrix((0, 0), dtype=np.float64)
@@ -692,9 +694,9 @@ def __init__(
         alpha_x=None,
         alpha_y=None,
         alpha_z=None,
-        alpha_xx=None,
-        alpha_yy=None,
-        alpha_zz=None,
+        alpha_xx=0.0,
+        alpha_yy=0.0,
+        alpha_zz=0.0,
         gmm=None,
         wiresmap=None,
         maplist=None,
@@ -733,7 +735,7 @@ def __init__(
         for map, wire, weights in zip(self.maplist, self.wiresmap.maps, weights_list):
             objfcts += [
                 LeastSquaresRegularization(
-                    alpha_s=alpha_s,
+                    alpha_s=0.0,
                     alpha_x=alpha_x,
                     alpha_y=alpha_y,
                     alpha_z=alpha_z,
@@ -747,7 +749,7 @@ def __init__(
                 )
             ]
 
-        super(PGI, self).__init__(objfcts=objfcts)
+        super().__init__(objfcts=objfcts)
         self.reference_model_in_smooth = reference_model_in_smooth
 
     @property
diff --git a/SimPEG/regularization/regularization_mesh.py b/SimPEG/regularization/regularization_mesh.py
index 8fd191ad5a..44af364e6d 100755
--- a/SimPEG/regularization/regularization_mesh.py
+++ b/SimPEG/regularization/regularization_mesh.py
@@ -1,7 +1,6 @@
 import numpy as np
 import scipy.sparse as sp
-import warnings
-import properties
+from SimPEG.utils.code_utils import deprecate_property
 
 from .. import props
 from .. import utils
@@ -27,37 +26,88 @@ class RegularizationMesh(props.BaseSimPEG):
     """
 
     regularization_type = None  # or 'Base'
-    _active_cells = None
 
-    def __init__(self, mesh, **kwargs):
+    def __init__(self, mesh, active_cells=None, **kwargs):
         self.mesh = mesh
+        self.active_cells = active_cells
         utils.setKwargs(self, **kwargs)
 
     # active_cells = properties.Array("active indices in mesh", dtype=[bool, int])
 
     @property
     def active_cells(self):
+        """A boolean array indicating whether a cell is active
+
+        Returns
+        -------
+        (n_cells,) numpy.ndarray of bool
+
+        Notes
+        -----
+        If this is set with an array of integers, it interprets it as an array
+        listing the active cell indices.
+        """
         return self._active_cells
 
     @active_cells.setter
     def active_cells(self, values: np.ndarray):
+        if getattr(self, "_active_cells", None) is not None and not all(
+            self._active_cells == values
+        ):
+            raise AttributeError(
+                "The RegulatizationMesh already has an 'active_cells' property set."
+            )
         if values is not None:
-            if self._active_cells is not None and not all(self._active_cells == values):
-                raise AttributeError(
-                    "The RegulatizationMesh already has an 'active_cells' property set."
-                )
-
-            if (
-                not isinstance(values, np.ndarray) or values.dtype != "bool"
-            ):  # cast it to a bool otherwise
-                raise ValueError(
-                    "Input 'active_cells' must be an numpy.ndarray of type 'bool'."
-                )
+            try:
+                values = np.asarray(values)
+            except:
+                raise ValueError("Input 'active_cells' must be array_like.")
+
+            if values.dtype != bool:
+                try:
+                    tmp = np.zeros(self.mesh.nC, dtype=bool)
+                    tmp[values] = True
+                except:
+                    raise ValueError(
+                        "Values must be an array of integers or an array of bools "
+                        "indicating the active cells"
+                    )
+                if np.sum(tmp) != len(values):
+                    # This line should cause an error to be thrown if someone
+                    # accidentally passes a list of 0 & 1 integers instead of passing
+                    # it a list of booleans.
+                    raise ValueError(
+                        "Array was interpretted as a list of active indices and you "
+                        "attempted to set the same cell as active multiple times."
+                    )
+                values = tmp
 
             if values.shape != (self.mesh.nC,):
                 raise ValueError(
                     f"Input 'active_cells' must have shape {(self.mesh.nC,)}"
                 )
+            # Ensure any cached operators created when
+            # active_cells was None are deleted
+            self._vol = None
+            self._Pac = None
+            self._Pafx = None
+            self._Pafy = None
+            self._Pafz = None
+            self._aveFx2CC = None
+            self._aveFy2CC = None
+            self._aveFz2CC = None
+            self._aveCC2Fx = None
+            self._aveCC2Fy = None
+            self._aveCC2Fz = None
+            self._cell_gradient_x = None
+            self._cell_gradient_y = None
+            self._cell_gradient_z = None
+            self._faceDiffx = None
+            self._faceDiffy = None
+            self._faceDiffz = None
+            self._cell_distances_x = None
+            self._cell_distances_y = None
+            self._cell_distances_z = None
         self._active_cells = values
 
     @property
@@ -68,8 +118,10 @@ def vol(self):
         :rtype: numpy.ndarray
         :return: reduced cell volume
         """
+        if self.active_cells is None:
+            return self.mesh.cell_volumes
         if getattr(self, "_vol", None) is None:
-            self._vol = self.Pac.T * self.mesh.cell_volumes
+            self._vol = self.mesh.cell_volumes[self.active_cells]
         return self._vol
 
     @property
@@ -92,9 +144,7 @@ def dim(self):
         :rtype: int
         :return: dimension
         """
-        if getattr(self, "_dim", None) is None:
-            self._dim = self.mesh.dim
-        return self._dim
+        return self.mesh.dim
 
     @property
     def Pac(self):
@@ -156,6 +206,8 @@ def Pafy(self):
                     self.mesh.average_cell_to_total_face_y() * ind_active
                 ) >= 1
                 self._Pafy = utils.speye(self.mesh.ntFy)[:, active_cells_Fy]
+            elif self.mesh._meshType == "CYL" and self.mesh.is_symmetric:
+                return sp.csr_matrix((0, 0))
             else:
                 if self.active_cells is None:
                     self._Pafy = utils.speye(self.mesh.nFy)
@@ -212,7 +264,6 @@ def aveFx2CC(self):
                 nCinRow = utils.mkvc((self.aveCC2Fx.T).sum(1))
                 nCinRow[nCinRow > 0] = 1.0 / nCinRow[nCinRow > 0]
                 self._aveFx2CC = utils.sdiag(nCinRow) * self.aveCC2Fx.T
-
             else:
                 self._aveFx2CC = self.Pac.T * self.mesh.aveFx2CC * self.Pafx
 
@@ -250,6 +301,8 @@ def aveFy2CC(self):
                 nCinRow = utils.mkvc((self.aveCC2Fy.T).sum(1))
                 nCinRow[nCinRow > 0] = 1.0 / nCinRow[nCinRow > 0]
                 self._aveFy2CC = utils.sdiag(nCinRow) * self.aveCC2Fy.T
+            elif self.mesh._meshType == "CYL" and self.mesh.is_symmetric:
+                return sp.csr_matrix((self.nC, 0))
             else:
                 self._aveFy2CC = self.Pac.T * self.mesh.aveFy2CC * self.Pafy
 
@@ -268,6 +321,8 @@ def aveCC2Fy(self):
                 self._aveCC2Fy = (
                     self.Pafy.T * self.mesh.average_cell_to_total_face_y() * self.Pac
                 )
+            elif self.mesh._meshType == "CYL" and self.mesh.is_symmetric:
+                return sp.csr_matrix((0, self.nC))
             else:
                 self._aveCC2Fy = (
                     utils.sdiag(1.0 / (self.aveFy2CC.T).sum(1)) * self.aveFy2CC.T
@@ -315,29 +370,31 @@ def aveCC2Fz(self):
     def base_length(self):
         """The smallest core cell size"""
         if getattr(self, "_base_length", None) is None:
-            self._base_length = self.mesh.h_gridded.min()
+            self._base_length = self.mesh.edge_lengths.min()
         return self._base_length
 
     @property
     def cell_gradient(self):
         if self.dim == 1:
-            return self.cellDiffx
+            return self.cell_gradient_x
         elif self.dim == 2:
-            return sp.vstack([self.cellDiffx, self.cellDiffy])
+            return sp.vstack([self.cell_gradient_x, self.cell_gradient_y])
         else:
-            return sp.vstack([self.cellDiffx, self.cellDiffy, self.cellDiffz])
+            return sp.vstack(
+                [self.cell_gradient_x, self.cell_gradient_y, self.cell_gradient_z]
+            )
 
     @property
-    def cellDiffx(self):
+    def cell_gradient_x(self):
         """
-        cell centered difference in the x-direction
+        cell centered gradient in the x-direction
 
         :rtype: scipy.sparse.csr_matrix
         :return: differencing matrix for active cells in the x-direction
         """
-        if getattr(self, "_cellDiffx", None) is None:
+        if getattr(self, "_cell_gradient_x", None) is None:
             if self.mesh._meshType == "TREE":
-                self._cellDiffx = (
+                self._cell_gradient_x = (
                     self.Pafx.T
                     * utils.sdiag(
                         self.mesh.average_cell_to_total_face_x()
@@ -347,20 +404,22 @@ def cellDiffx(self):
                     * self.Pac
                 )
             else:
-                self._cellDiffx = self.Pafx.T * self.mesh.cellGradx * self.Pac
-        return self._cellDiffx
+                self._cell_gradient_x = (
+                    self.Pafx.T * self.mesh.cell_gradient_x * self.Pac
+                )
+        return self._cell_gradient_x
 
     @property
-    def cellDiffy(self):
+    def cell_gradient_y(self):
         """
-        cell centered difference in the y-direction
+        cell centered gradient in the y-direction
 
         :rtype: scipy.sparse.csr_matrix
         :return: differencing matrix for active cells in the y-direction
         """
-        if getattr(self, "_cellDiffy", None) is None:
+        if getattr(self, "_cell_gradient_y", None) is None:
             if self.mesh._meshType == "TREE":
-                self._cellDiffy = (
+                self._cell_gradient_y = (
                     self.Pafy.T
                     * utils.sdiag(
                         self.mesh.average_cell_to_total_face_y()
@@ -370,20 +429,22 @@ def cellDiffy(self):
                     * self.Pac
                 )
             else:
-                self._cellDiffy = self.Pafy.T * self.mesh.cellGrady * self.Pac
-        return self._cellDiffy
+                self._cell_gradient_y = (
+                    self.Pafy.T * self.mesh.cell_gradient_y * self.Pac
+                )
+        return self._cell_gradient_y
 
     @property
-    def cellDiffz(self):
+    def cell_gradient_z(self):
         """
         cell centered difference in the z-direction
 
         :rtype: scipy.sparse.csr_matrix
         :return: differencing matrix for active cells in the z-direction
         """
-        if getattr(self, "_cellDiffz", None) is None:
+        if getattr(self, "_cell_gradient_z", None) is None:
             if self.mesh._meshType == "TREE":
-                self._cellDiffz = (
+                self._cell_gradient_z = (
                     self.Pafz.T
                     * utils.sdiag(
                         self.mesh.average_cell_to_total_face_z()
@@ -393,8 +454,41 @@ def cellDiffz(self):
                     * self.Pac
                 )
             else:
-                self._cellDiffz = self.Pafz.T * self.mesh.cellGradz * self.Pac
-        return self._cellDiffz
+                self._cell_gradient_z = (
+                    self.Pafz.T * self.mesh.cell_gradient_z * self.Pac
+                )
+        return self._cell_gradient_z
+
+    cellDiffx = deprecate_property(
+        cell_gradient_x, "cellDiffx", "0.x.0", error=False, future_warn=False
+    )
+    cellDiffy = deprecate_property(
+        cell_gradient_y, "cellDiffy", "0.x.0", error=False, future_warn=False
+    )
+    cellDiffz = deprecate_property(
+        cell_gradient_z, "cellDiffz", "0.x.0", error=False, future_warn=False
+    )
+
+    @property
+    def cell_distances_x(self):
+        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])
+        return self._cell_distances_x
+
+    @property
+    def cell_distances_y(self):
+        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])
+        return self._cell_distances_x
+
+    @property
+    def cell_distances_z(self):
+        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])
+        return self._cell_distances_x
 
     @property
     def faceDiffx(self):
diff --git a/SimPEG/regularization/sparse.py b/SimPEG/regularization/sparse.py
index 70d2acbd64..7e8cff1a1a 100644
--- a/SimPEG/regularization/sparse.py
+++ b/SimPEG/regularization/sparse.py
@@ -2,6 +2,8 @@
 
 import numpy as np
 
+from discretize.base import BaseMesh
+
 from .base import (
     BaseRegularization,
     LeastSquaresRegularization,
@@ -17,12 +19,11 @@ class BaseSparse(BaseRegularization):
     Base class for building up the components of the Sparse Regularization
     """
 
-    _irls_scaled = True
-    _irls_threshold = 1e-8
-    _norm = 2.0
-
-    def __init__(self, mesh, **kwargs):
+    def __init__(self, mesh, norm=2.0, irls_scaled=True, irls_threshold=1e-8, **kwargs):
         super().__init__(mesh=mesh, **kwargs)
+        self.norm = norm
+        self.irls_scaled = irls_scaled
+        self.irls_threshold = irls_threshold
 
     @property
     def irls_scaled(self) -> bool:
@@ -45,15 +46,13 @@ def irls_threshold(self):
         """
         Constant added to the denominator of the IRLS weights for stability.
         """
-        if self._irls_threshold is None:
-            raise AttributeError("'irls_threhsold' must be set on the regularization.")
         return self._irls_threshold
 
     @irls_threshold.setter
     def irls_threshold(self, value):
+        value = float(value)
         if value <= 0:
             raise ValueError("Value of 'irls_threshold' should be larger than 0.")
-
         self._irls_threshold = value
 
     @property
@@ -61,15 +60,15 @@ def norm(self):
         """
         Value of the norm
         """
-        if getattr(self, "_norm", None) is None:
-            self.norm = 2.0
         return self._norm
 
     @norm.setter
     def norm(self, value: float | np.ndarray | None):
-        if value is not None:
+        if value is None:
+            value = 2.0
+        else:
             if isinstance(value, (float, int)):
-                value = np.ones(self.shape[0]) * value
+                value = np.ones(self._weights_shapes[0]) * value
 
             if np.any(value < 0) or np.any(value > 2):
                 raise ValueError(
@@ -112,15 +111,12 @@ class SparseSmall(BaseSparse, Small):
 
     _multiplier_pair = "alpha_s"
 
-    def __init__(self, mesh, **kwargs):
-        super().__init__(mesh=mesh, **kwargs)
-
     def update_weights(self, m):
         """
         Compute and store the irls weights.
         """
         f_m = self.f_m(m)
-        self.add_set_weights({"irls": self.get_lp_weights(f_m)})
+        self.set_weights(irls=self.get_lp_weights(f_m))
 
 
 class SparseDeriv(BaseSparse, SmoothDeriv):
@@ -128,9 +124,11 @@ class SparseDeriv(BaseSparse, SmoothDeriv):
     Base Class for sparse regularization on first spatial derivatives
     """
 
-    _gradient_type = "total"
-
-    def __init__(self, mesh, orientation="x", **kwargs):
+    def __init__(self, mesh, orientation="x", gradient_type="total", **kwargs):
+        if "gradientType" in kwargs:
+            self.gradientType = kwargs.pop("gradientType")
+        else:
+            self.gradient_type = gradient_type
         super().__init__(mesh=mesh, orientation=orientation, **kwargs)
 
     def update_weights(self, m):
@@ -151,7 +149,10 @@ def update_weights(self, m):
                         self.regularization_mesh.Pac.T
                         * self.regularization_mesh.mesh.h_gridded[:, ii]
                     )
-                    dm = getattr(self.regularization_mesh, f"cellDiff{comp}") * delta_m
+                    dm = (
+                        getattr(self.regularization_mesh, f"cell_gradient_{comp}")
+                        * delta_m
+                    )
 
                     if self.units == "radian":
                         dm = utils.mat_utils.coterminal(dm)
@@ -167,7 +168,7 @@ def update_weights(self, m):
         else:
             f_m = self.f_m(m)
 
-        self.add_set_weights({"irls": self.get_lp_weights(f_m)})
+        self.set_weights(irls=self.get_lp_weights(f_m))
 
     @property
     def gradient_type(self) -> str:
@@ -196,35 +197,47 @@ class Sparse(LeastSquaresRegularization):
 
     .. math::
 
-        R(m) = \\frac{1}{2}\mathbf{(m-m_\\text{ref})^\\top W^\\top R^\\top R
+        R(m) = \\frac{1}{2}\\mathbf{(m-m_\\text{ref})^\\top W^\\top R^\\top R
         W(m-m_\\text{ref})}
 
     where the IRLS weight
 
     .. math::
 
-        R = \eta TO FINISH LATER!!!
+        R = \\eta TO FINISH LATER!!!
 
     So the derivative is straight forward:
 
     .. math::
 
-        R(m) = \mathbf{W^\\top R^\\top R W (m-m_\\text{ref})}
+        R(m) = \\mathbf{W^\\top R^\\top R W (m-m_\\text{ref})}
 
     The IRLS weights are recomputed after each beta solves.
     It is strongly recommended to do a few Gauss-Newton iterations
     before updating.
     """
 
-    _irls_scaled = True
-    _irls_threshold = 1e-8
-    _gradient_type = "total"
-    _norms = None
-
     def __init__(
-        self, mesh, active_cells=None, norms=None, gradient_type="total", **kwargs
+        self,
+        mesh,
+        active_cells=None,
+        norms=None,
+        gradient_type="total",
+        irls_scaled=True,
+        irls_threshold=1e-8,
+        **kwargs,
     ):
-        self.regularization_mesh = RegularizationMesh(mesh)
+        if not isinstance(mesh, RegularizationMesh):
+            mesh = RegularizationMesh(mesh)
+
+        if not isinstance(mesh, RegularizationMesh):
+            TypeError(
+                f"'regularization_mesh' must be of type {RegularizationMesh} or {BaseMesh}. "
+                f"Value of type {type(mesh)} provided."
+            )
+        self._regularization_mesh = mesh
+        if active_cells is not None:
+            self._regularization_mesh.active_cells = active_cells
 
         objfcts = [
             SparseSmall(mesh=self.regularization_mesh),
@@ -237,14 +250,24 @@ def __init__(
         if mesh.dim > 2:
             objfcts.append(SparseDeriv(mesh=self.regularization_mesh, orientation="z"))
 
+        gradientType = kwargs.pop("gradientType", None)
         super().__init__(
             self.regularization_mesh,
             objfcts=objfcts,
-            active_cells=active_cells,
-            norms=norms,
-            gradient_type=gradient_type,
             **kwargs,
         )
+        if norms is None:
+            norms = [1] * (mesh.dim + 1)
+        self.norms = norms
+
+        if gradientType is not None:
+            # Trigger deprecation warning
+            self.gradientType = gradientType
+        else:
+            self.gradient_type = gradient_type
+
+        self.irls_scaled = irls_scaled
+        self.irls_threshold = irls_threshold
 
     @property
     def gradient_type(self) -> str:
@@ -256,8 +279,8 @@ def gradient_type(self) -> str:
     @gradient_type.setter
     def gradient_type(self, value: str):
         for fct in self.objfcts:
-            if hasattr(fct, "_gradient_type"):
-                fct._gradient_type = value
+            if hasattr(fct, "gradient_type"):
+                fct.gradient_type = value
 
         self._gradient_type = value
 
@@ -300,10 +323,9 @@ def irls_scaled(self, value: bool):
                 "'irls_scaled must be of type 'bool'. "
                 f"Value of type {type(value)} provided."
             )
-        self._irls_scaled = value
-
         for fct in self.objfcts:
             fct.irls_scaled = value
+        self._irls_scaled = value
 
     @property
     def irls_threshold(self):
@@ -314,6 +336,7 @@ def irls_threshold(self):
 
     @irls_threshold.setter
     def irls_threshold(self, value):
+        value = float(value)
         if value <= 0:
             raise ValueError("Value of 'irls_threshold' should be greater than 0.")
 
diff --git a/examples/01-maps/plot_sumMap.py b/examples/01-maps/plot_sumMap.py
index 5341a5daf1..e60a257525 100644
--- a/examples/01-maps/plot_sumMap.py
+++ b/examples/01-maps/plot_sumMap.py
@@ -122,9 +122,7 @@ def run(plotIt=True):
 
     # Take the cell number out of the scaling.
     # Want to keep high sens for large volumes
-    scale = utils.sdiag(
-        np.r_[utils.mkvc(1.0 / homogMap.P.sum(axis=0)), np.ones(nC)]
-    )
+    scale = utils.sdiag(np.r_[utils.mkvc(1.0 / homogMap.P.sum(axis=0)), np.ones(nC)])
 
     # for ii in range(survey.nD):
     #     wr += (
@@ -132,9 +130,14 @@ def run(plotIt=True):
     #         / data.standard_deviation[ii]
     #     ) ** 2.0 / np.r_[homogMap.P.T * mesh.cell_volumes[actv], mesh.cell_volumes[actv]] **2.
 
-    wr = prob.getJtJdiag(np.ones(sumMap.shape[1])) / np.r_[homogMap.P.T * mesh.cell_volumes[actv], mesh.cell_volumes[actv]] **2.
+    wr = (
+        prob.getJtJdiag(np.ones(sumMap.shape[1]))
+        / np.r_[homogMap.P.T * mesh.cell_volumes[actv], mesh.cell_volumes[actv]] ** 2.0
+    )
     # Scale the model spaces independently
-    wr[wires.homo.index] /= np.max((wires.homo * wr)) * utils.mkvc(homogMap.P.sum(axis=0).flatten())
+    wr[wires.homo.index] /= np.max((wires.homo * wr)) * utils.mkvc(
+        homogMap.P.sum(axis=0).flatten()
+    )
     wr[wires.hetero.index] /= np.max(wires.hetero * wr)
     wr = wr ** 0.5
 
@@ -148,7 +151,9 @@ def run(plotIt=True):
     reg_m1.mref = np.zeros(sumMap.shape[1])
 
     # Regularization for the voxel model
-    reg_m2 = regularization.Sparse(mesh, active_cells=actv, mapping=wires.hetero, gradient_type="components")
+    reg_m2 = regularization.Sparse(
+        mesh, active_cells=actv, mapping=wires.hetero, gradient_type="components"
+    )
     reg_m2.cell_weights = wires.hetero * wr
     reg_m2.norms = [0, 0, 0, 0]
     reg_m2.mref = np.zeros(sumMap.shape[1])
diff --git a/examples/03-magnetics/plot_inv_mag_MVI_Sparse_TreeMesh.py b/examples/03-magnetics/plot_inv_mag_MVI_Sparse_TreeMesh.py
index 0299a3f3a0..d2d2038b10 100644
--- a/examples/03-magnetics/plot_inv_mag_MVI_Sparse_TreeMesh.py
+++ b/examples/03-magnetics/plot_inv_mag_MVI_Sparse_TreeMesh.py
@@ -380,18 +380,24 @@ def plotVectorSectionsOctree(
 
 # Create a Combo Regularization
 # Regularize the amplitude of the vectors
-reg_a = regularization.Sparse(mesh, gradient_type="components", active_cells=actv, mapping=wires.amp)
+reg_a = regularization.Sparse(
+    mesh, gradient_type="components", active_cells=actv, mapping=wires.amp
+)
 reg_a.norms = [0.0, 0.0, 0.0, 0.0]  # Sparse on the model and its gradients
 reg_a.reference_model = np.zeros(3 * nC)
 
 # Regularize the vertical angle of the vectors
-reg_t = regularization.Sparse(mesh, gradient_type="components", active_cells=actv, mapping=wires.theta)
+reg_t = regularization.Sparse(
+    mesh, gradient_type="components", active_cells=actv, mapping=wires.theta
+)
 reg_t.alpha_s = 0.0  # No reference angle
 reg_t.space = "spherical"
 reg_t.norms = [0.0, 0.0, 0.0, 0.0]  # Only norm on gradients used
 
 # Regularize the horizontal angle of the vectors
-reg_p = regularization.Sparse(mesh, gradient_type="components", active_cells=actv, mapping=wires.phi)
+reg_p = regularization.Sparse(
+    mesh, gradient_type="components", active_cells=actv, mapping=wires.phi
+)
 reg_p.alpha_s = 0.0  # No reference angle
 reg_p.space = "spherical"
 reg_p.norms = [0.0, 0.0, 0.0, 0.0]  # Only norm on gradients used
diff --git a/examples/10-pgi/plot_inv_1_PGI_Linear_1D_joint_WithRelationships.py b/examples/10-pgi/plot_inv_1_PGI_Linear_1D_joint_WithRelationships.py
index 4dede30e74..3a57cefe40 100644
--- a/examples/10-pgi/plot_inv_1_PGI_Linear_1D_joint_WithRelationships.py
+++ b/examples/10-pgi/plot_inv_1_PGI_Linear_1D_joint_WithRelationships.py
@@ -159,9 +159,7 @@ def g(k):
     chi0_ratio=np.r_[1.0, 1.0], verbose=True, n_pw_iter=10
 )
 scaling_schedule = directives.JointScalingSchedule(verbose=True)
-alpha0_ratio = np.r_[
-    1e+6, 1e+4
-]
+alpha0_ratio = np.r_[1e6, 1e4]
 alphas = directives.AlphasSmoothEstimate_ByEig(
     alpha0_ratio=alpha0_ratio, n_pw_iter=10, verbose=True
 )
@@ -236,9 +234,13 @@ def g(k):
 
 # LeastSquaresRegularization Inversion
 
-reg1 = regularization.LeastSquaresRegularization(mesh, alpha_s=1.0, alpha_x=1.0, mapping=wires.m1)
+reg1 = regularization.LeastSquaresRegularization(
+    mesh, alpha_s=1.0, alpha_x=1.0, mapping=wires.m1
+)
 reg1.cell_weights = wr1
-reg2 = regularization.LeastSquaresRegularization(mesh, alpha_s=1.0, alpha_x=1.0, mapping=wires.m2)
+reg2 = regularization.LeastSquaresRegularization(
+    mesh, alpha_s=1.0, alpha_x=1.0, mapping=wires.m2
+)
 reg2.cell_weights = wr2
 reg = reg1 + reg2
 
diff --git a/examples/20-published/plot_heagyetal2017_casing.py b/examples/20-published/plot_heagyetal2017_casing.py
index cc543593a5..5563e26a64 100644
--- a/examples/20-published/plot_heagyetal2017_casing.py
+++ b/examples/20-published/plot_heagyetal2017_casing.py
@@ -1479,5 +1479,6 @@ def run(plotIt=True, runTests=False, reRun=False, saveFig=False):
         except PermissionError:
             pass
 
+
 if __name__ == "__main__":
     run(plotIt=True, runTests=False, reRun=False, saveFig=False)
diff --git a/examples/20-published/plot_laguna_del_maule_inversion.py b/examples/20-published/plot_laguna_del_maule_inversion.py
index 8c2c307347..503314238a 100644
--- a/examples/20-published/plot_laguna_del_maule_inversion.py
+++ b/examples/20-published/plot_laguna_del_maule_inversion.py
@@ -241,7 +241,8 @@ def run(plotIt=True, cleanAfterRun=True):
         plt.ylabel("Elevation")
         plt.gca().set_aspect("equal", adjustable="box")
         cb = plt.colorbar(
-            dat1[0], ax=ax,
+            dat1[0],
+            ax=ax,
             orientation="vertical",
             ticks=np.linspace(vmin, vmax, 4),
             cmap="bwr",
diff --git a/tests/base/test_maps.py b/tests/base/test_maps.py
index cebc0e5ca2..32a5251ba2 100644
--- a/tests/base/test_maps.py
+++ b/tests/base/test_maps.py
@@ -572,7 +572,7 @@ def test_sphericalInclusions(self):
         mesh = discretize.TensorMesh([4, 5, 3])
         mapping = maps.SelfConsistentEffectiveMedium(mesh, sigma0=1e-1, sigma1=1.0)
         m = np.abs(np.random.rand(mesh.nC))
-        mapping.test(m=m, dx=0.05, num=3)
+        mapping.test(m=m, dx=0.05 * np.ones(mesh.n_cells), num=3)
 
     def test_spheroidalInclusions(self):
         mesh = discretize.TensorMesh([4, 3, 2])
@@ -580,7 +580,7 @@ def test_spheroidalInclusions(self):
             mesh, sigma0=1e-1, sigma1=1.0, alpha0=0.8, alpha1=0.9, rel_tol=1e-8
         )
         m = np.abs(np.random.rand(mesh.nC))
-        mapping.test(m=m, dx=0.05, num=3)
+        mapping.test(m=m, dx=0.05 * np.ones(mesh.n_cells), num=3)
 
 
 if __name__ == "__main__":
diff --git a/tests/base/test_regularization.py b/tests/base/test_regularization.py
index b7c4446a0f..2e67f2fbcb 100644
--- a/tests/base/test_regularization.py
+++ b/tests/base/test_regularization.py
@@ -118,7 +118,9 @@ def test_regularization_ActiveCells(self):
                         continue
 
                     nP = int(active_cells.sum())
-                    reg = r(mesh, active_cells=active_cells, mapping=maps.IdentityMap(nP=nP))
+                    reg = r(
+                        mesh, active_cells=active_cells, mapping=maps.IdentityMap(nP=nP)
+                    )
                     m = np.random.rand(mesh.nC)[active_cells]
                     mref = np.ones_like(m) * np.mean(m)
                     reg.reference_model = mref
@@ -155,7 +157,9 @@ def test_regularizationMesh(self):
                         + 0.5
                     )
 
-                regularization_mesh = regularization.RegularizationMesh(mesh, active_cells=indAct)
+                regularization_mesh = regularization.RegularizationMesh(
+                    mesh, active_cells=indAct
+                )
 
                 assert (regularization_mesh.vol == mesh.vol[indAct]).all()
 
@@ -168,13 +172,16 @@ def test_property_mirroring(self):
 
             self.assertTrue(reg.nP == reg.regularization_mesh.nC)
 
-            [self.assertTrue(np.all(fct.active_cells == active_cells)) for fct in reg.objfcts]
+            [
+                self.assertTrue(np.all(fct.active_cells == active_cells))
+                for fct in reg.objfcts
+            ]
 
             # test assignment of cell weights
             cell_weights = np.random.rand(active_cells.sum())
-            reg.add_set_weights(cell_weights)
+            reg.set_weights(user_weights=cell_weights)
             [
-                self.assertTrue(np.all(fct.weights["user_weights"] == cell_weights))
+                self.assertTrue(np.all(fct.get_weights("user_weights") == cell_weights))
                 for fct in reg.objfcts
             ]
 
@@ -267,9 +274,7 @@ def test_mappings_and_cell_weights(self):
 
         wires = maps.Wires(("sigma", mesh.nC), ("mu", mesh.nC))
 
-        reg = regularization.Small(
-            mesh, mapping=wires.sigma, weights=cell_weights
-        )
+        reg = regularization.Small(mesh, mapping=wires.sigma, weights=cell_weights)
 
         objfct = objective_function.L2ObjectiveFunction(
             W=utils.sdiag(np.sqrt(cell_weights * mesh.cell_volumes)),
@@ -293,7 +298,8 @@ def test_update_of_sparse_norms(self):
             np.all(
                 reg.norms
                 == np.kron(
-                    np.ones((reg.regularization_mesh.Pac.shape[1], 1)), np.c_[2.0, 2.0, 2.0, 2.0]
+                    np.ones((reg.regularization_mesh.Pac.shape[1], 1)),
+                    np.c_[2.0, 2.0, 2.0, 2.0],
                 )
             )
         )
@@ -308,7 +314,8 @@ def test_update_of_sparse_norms(self):
             np.all(
                 reg.norms
                 == np.kron(
-                    np.ones((reg.regularization_mesh.Pac.shape[1], 1)), np.c_[0.0, 1.0, 1.0, 1.0]
+                    np.ones((reg.regularization_mesh.Pac.shape[1], 1)),
+                    np.c_[0.0, 1.0, 1.0, 1.0],
                 )
             )
         )
@@ -321,15 +328,21 @@ def test_linked_properties(self):
         mesh = discretize.TensorMesh([8, 7, 6])
         reg = regularization.LeastSquaresRegularization(mesh)
 
-        [self.assertTrue(reg.regularization_mesh is fct.regularization_mesh) for fct in reg.objfcts]
+        [
+            self.assertTrue(reg.regularization_mesh is fct.regularization_mesh)
+            for fct in reg.objfcts
+        ]
         [self.assertTrue(reg.mapping is fct.mapping) for fct in reg.objfcts]
 
         D = reg.regularization_mesh.cellDiffx
-        reg.regularization_mesh._cellDiffx = 4 * D
+        reg.regularization_mesh._cell_gradient_x = 4 * D
         v = np.random.rand(D.shape[1])
         [
             self.assertTrue(
-                np.all(reg.regularization_mesh._cellDiffx * v == fct.regularization_mesh.cellDiffx * v)
+                np.all(
+                    reg.regularization_mesh._cell_gradient_x * v
+                    == fct.regularization_mesh.cellDiffx * v
+                )
             )
             for fct in reg.objfcts
         ]
@@ -337,10 +350,15 @@ def test_linked_properties(self):
         active_cells = mesh.gridCC[:, 2] < 0.4
         reg.active_cells = active_cells
         self.assertTrue(np.all(reg.regularization_mesh.active_cells == active_cells))
-        [self.assertTrue(np.all(reg.active_cells == fct.active_cells)) for fct in reg.objfcts]
+        [
+            self.assertTrue(np.all(reg.active_cells == fct.active_cells))
+            for fct in reg.objfcts
+        ]
 
         [
-            self.assertTrue(np.all(reg.active_cells == fct.regularization_mesh.active_cells))
+            self.assertTrue(
+                np.all(reg.active_cells == fct.regularization_mesh.active_cells)
+            )
             for fct in reg.objfcts
         ]
 
@@ -378,7 +396,7 @@ def test_active_cells_nc_residual(self):
         temp = np.logspace(np.log10(1.0), np.log10(12.0), 19)
         temp_pad = temp[-1] * 1.3 ** np.arange(npad)
         hz = np.r_[temp_pad[::-1], temp[::-1], temp, temp_pad]
-        mesh = discretize.CylMesh([hx, 1, hz], "00C")
+        mesh = discretize.CylMesh([hx, 3, hz], "00C")
         active = mesh.cell_centers[:, 2] < 0.0
 
         reg = regularization.LeastSquaresRegularization(mesh, active_cells=active)
diff --git a/tests/dask/test_grav_inversion_linear.py b/tests/dask/test_grav_inversion_linear.py
index 855324b50e..28c23fed8b 100644
--- a/tests/dask/test_grav_inversion_linear.py
+++ b/tests/dask/test_grav_inversion_linear.py
@@ -67,7 +67,7 @@ def setUp(self):
         # We can now create a density model and generate data
         # Here a simple block in half-space
         model = np.zeros(self.mesh.shape_cells)
-        model[(midx - 2): (midx + 2), (midy - 2): (midy + 2), -6:-2] = 0.5
+        model[(midx - 2) : (midx + 2), (midy - 2) : (midy + 2), -6:-2] = 0.5
         model = utils.mkvc(model)
         self.model = model[actv]
 
diff --git a/tests/dask/test_mag_inversion_linear_Octree.py b/tests/dask/test_mag_inversion_linear_Octree.py
index 64e21c1d82..1ba7236a8d 100644
--- a/tests/dask/test_mag_inversion_linear_Octree.py
+++ b/tests/dask/test_mag_inversion_linear_Octree.py
@@ -115,7 +115,9 @@ def setUp(self):
         )
 
         # Create a regularization
-        reg = regularization.Sparse(self.mesh, indActive=actv, mapping=idenMap, gradient_type="components")
+        reg = regularization.Sparse(
+            self.mesh, indActive=actv, mapping=idenMap, gradient_type="components"
+        )
         reg.norms = [0, 0, 0, 0]
         reg.reference_model = np.zeros(nC)
 
diff --git a/tests/em/vrm/test_vrminv.py b/tests/em/vrm/test_vrminv.py
index 27e48bda53..d2f906b4d6 100644
--- a/tests/em/vrm/test_vrminv.py
+++ b/tests/em/vrm/test_vrminv.py
@@ -61,7 +61,12 @@ def test_basic_inversion(self):
         Survey.noise_floor = 1e-11
 
         dmis = data_misfit.L2DataMisfit(data=dobs, simulation=Problem)
-        W = mkvc((np.sum(np.array(Problem.A) ** 2, axis=0))/meshObj.cell_volumes**2.) ** 0.25
+        W = (
+            mkvc(
+                (np.sum(np.array(Problem.A) ** 2, axis=0)) / meshObj.cell_volumes ** 2.0
+            )
+            ** 0.25
+        )
         reg = regularization.LeastSquaresRegularization(
             meshObj, alpha_s=0.01, alpha_x=1.0, alpha_y=1.0, alpha_z=1.0, weights=W
         )
diff --git a/tests/pf/test_grav_inversion_linear.py b/tests/pf/test_grav_inversion_linear.py
index 0f60530697..861f4b0d39 100644
--- a/tests/pf/test_grav_inversion_linear.py
+++ b/tests/pf/test_grav_inversion_linear.py
@@ -20,14 +20,11 @@
 
 class GravInvLinProblemTest(unittest.TestCase):
     def setUp(self):
-
         # Create a self.mesh
         dx = 5.0
-
         hxind = [(dx, 5, -1.3), (dx, 5), (dx, 5, 1.3)]
         hyind = [(dx, 5, -1.3), (dx, 5), (dx, 5, 1.3)]
         hzind = [(dx, 5, -1.3), (dx, 6)]
-
         self.mesh = discretize.TensorMesh([hxind, hyind, hzind], "CCC")
 
         # Get index of the center
@@ -87,8 +84,7 @@ def setUp(self):
         # Create a regularization
         reg = regularization.Sparse(self.mesh, active_cells=actv, mapping=idenMap)
         reg.norms = [0, 0, 0, 0]
-        reg.gradientType = "component"
-        # reg.eps_p, reg.eps_q = 5e-2, 1e-2
+        reg.gradientType = "components"
 
         # Data misfit function
         dmis = data_misfit.L2DataMisfit(simulation=sim, data=data)
@@ -113,7 +109,6 @@ def test_grav_inverse(self):
         # Run the inversion
         mrec = self.inv.run(self.model)
         residual = np.linalg.norm(mrec - self.model) / np.linalg.norm(self.model)
-        print(residual)
 
         # import matplotlib.pyplot as plt
         # plt.figure()
diff --git a/tests/pf/test_mag_inversion_linear.py b/tests/pf/test_mag_inversion_linear.py
index 85ef268c3e..2cedef6998 100644
--- a/tests/pf/test_mag_inversion_linear.py
+++ b/tests/pf/test_mag_inversion_linear.py
@@ -29,11 +29,9 @@ def setUp(self):
 
         # Create a mesh
         dx = 5.0
-
         hxind = [(dx, 5, -1.3), (dx, 5), (dx, 5, 1.3)]
         hyind = [(dx, 5, -1.3), (dx, 5), (dx, 5, 1.3)]
         hzind = [(dx, 5, -1.3), (dx, 6)]
-
         self.mesh = discretize.TensorMesh([hxind, hyind, hzind], "CCC")
 
         # Get index of the center
@@ -97,12 +95,10 @@ def setUp(self):
         # Create a regularization
         reg = regularization.Sparse(self.mesh, indActive=actv, mapping=idenMap)
         reg.norms = [0, 0, 0, 0]
-        reg.gradientType = "component"
-        # reg.eps_p, reg.eps_q = 1e-3, 1e-3
+        reg.gradientType = "components"
 
         # Data misfit function
         dmis = data_misfit.L2DataMisfit(simulation=sim, data=data)
-        # dmis.W = 1/wd
 
         # Add directives to the inversion
         opt = optimization.ProjectedGNCG(
@@ -124,9 +120,7 @@ def test_mag_inverse(self):
 
         # Run the inversion
         mrec = self.inv.run(self.model)
-
         residual = np.linalg.norm(mrec - self.model) / np.linalg.norm(self.model)
-        print(residual)
 
         # plt.figure()
         # ax = plt.subplot(1, 2, 1)
@@ -141,7 +135,6 @@ def test_mag_inverse(self):
         # plt.show()
 
         self.assertTrue(residual < 0.05)
-        # self.assertTrue(residual < 0.05)
 
     def tearDown(self):
         # Clean up the working directory
diff --git a/tutorials/13-joint_inversion/plot_inv_3_cross_gradient_pf.py b/tutorials/13-joint_inversion/plot_inv_3_cross_gradient_pf.py
index cc6693720e..a9ec0ef017 100755
--- a/tutorials/13-joint_inversion/plot_inv_3_cross_gradient_pf.py
+++ b/tutorials/13-joint_inversion/plot_inv_3_cross_gradient_pf.py
@@ -307,7 +307,9 @@
 dmis_mag = data_misfit.L2DataMisfit(data=data_object_mag, simulation=simulation_mag)
 
 # Define the regularization (model objective function).
-reg_grav = regularization.LeastSquaresRegularization(mesh, indActive=ind_active, mapping=wires.density)
+reg_grav = regularization.LeastSquaresRegularization(
+    mesh, indActive=ind_active, mapping=wires.density
+)
 reg_mag = regularization.LeastSquaresRegularization(
     mesh, indActive=ind_active, mapping=wires.susceptibility
 )