Repository: simpeg/simpeg
Pull request #1038 (2021-09-02 18:54:18)

Issue-1037
Addresses issue #1037 
- Standardize sparse inversion to cell_weights normalized by volume.
- Add model_type property to base.simulation to handle sensitivity weights update.
- Minor variables cleanup


Patch:
diff --git a/SimPEG/dask/potential_fields/magnetics/simulation.py b/SimPEG/dask/potential_fields/magnetics/simulation.py
index 5137800bb..3772cb8b2 100644
--- a/SimPEG/dask/potential_fields/magnetics/simulation.py
+++ b/SimPEG/dask/potential_fields/magnetics/simulation.py
@@ -17,7 +17,7 @@ def dask_getJtJdiag(self, m, W=None):
     if getattr(self, "_gtg_diagonal", None) is None:
         if not self.is_amplitude_data:
             diag = ((W[:, None] * self.G) ** 2).sum(axis=0).compute()
-        else:  # self.modelType is amplitude
+        else:
             fieldDeriv = self.fieldDeriv
             J = (
                 fieldDeriv[0, :, None] * self.G[::3]
diff --git a/SimPEG/directives/directives.py b/SimPEG/directives/directives.py
index 6c8b94c11..a764abb50 100644
--- a/SimPEG/directives/directives.py
+++ b/SimPEG/directives/directives.py
@@ -1,5 +1,3 @@
-from __future__ import print_function
-
 import properties
 import numpy as np
 import matplotlib.pyplot as plt
@@ -8,7 +6,7 @@ import os
 import scipy.sparse as sp
 from ..data_misfit import BaseDataMisfit
 from ..objective_function import ComboObjectiveFunction
-from ..maps import SphericalSystem, ComboMap
+from ..maps import IdentityMap, Wires
 from ..regularization import (
     BaseComboRegularization,
     BaseRegularization,
@@ -1598,92 +1596,65 @@ class Update_Wj(InversionDirective):
 class UpdateSensitivityWeights(InversionDirective):
     """
     Directive to take care of re-weighting
-    the non-linear magnetic problems.
+    the non-linear problems. Assumes that the map of the regularization
+    function is either Wires or Identity.
+    Good for any problem where J is formed explicitly.
     """
-
-    mapping = None
-    JtJdiag = None
     everyIter = True
     threshold = 1e-12
-    switch = True
+    normalization: bool = True
 
     def initialize(self):
+        """
+        Calculate and update sensitivity
+        for optimization and regularization
+        """
+        for reg in self.reg.objfcts:
+            if not isinstance(getattr(reg, "mapping"), (IdentityMap, Wires)):
+                raise TypeError(
+                    f"Mapping for the regularization must be of type {IdentityMap} or {Wires}. "
+                    + f"Input mapping of type {type(reg.mapping)}."
+                )
 
-        # Calculate and update sensitivity
-        # for optimization and regularization
         self.update()
 
     def endIter(self):
-
+        """
+        Update inverse problem
+        """
         if self.everyIter:
-            # Update inverse problem
             self.update()
 
     def update(self):
-
-        # Get sum square of columns of J
-        self.getJtJdiag()
-
-        # Compute normalized weights
-        self.wr = self.getWr()
-
-        # Update the regularization
-        self.updateReg()
-
-    def getJtJdiag(self):
         """
         Compute explicitly the main diagonal of JtJ
-        Good for any problem where J is formed explicitly
+
         """
-        self.JtJdiag = []
+        jtj_diag = np.zeros_like(self.invProb.model)
         m = self.invProb.model
 
         for sim, dmisfit in zip(self.simulation, self.dmisfit.objfcts):
-
             if getattr(sim, "getJtJdiag", None) is None:
-                assert getattr(sim, "getJ", None) is not None, (
-                    "Simulation does not have a getJ attribute."
-                    + "Cannot form the sensitivity explicitly"
-                )
-                self.JtJdiag += [
-                    mkvc(np.sum((dmisfit.W * sim.getJ(m)) ** (2.0), axis=0))
-                ]
+                if getattr(sim, "getJ", None) is None:
+                    raise AttributeError(
+                        "Simulation does not have a getJ attribute."
+                        + "Cannot form the sensitivity explicitly"
+                    )
+                jtj_diag += mkvc(np.sum((dmisfit.W * sim.getJ(m)) ** 2.0, axis=0))
             else:
-                self.JtJdiag += [sim.getJtJdiag(m, W=dmisfit.W)]
-
-        return self.JtJdiag
-
-    def getWr(self):
-        """
-        Take the diagonal of JtJ and return
-        a normalized sensitivty weighting vector
-        """
+                jtj_diag += sim.getJtJdiag(m, W=dmisfit.W)
 
+        # Normalize and threshold weights
         wr = np.zeros_like(self.invProb.model)
-        if self.switch:
-            for prob_JtJ, sim, dmisfit in zip(
-                self.JtJdiag, self.simulation, self.dmisfit.objfcts
-            ):
-
-                wr += prob_JtJ + self.threshold
-
-            wr = wr ** 0.5
-            wr /= wr.max()
-        else:
-            wr += 1.0
-
-        return wr
-
-    def updateReg(self):
-        """
-        Update the cell weights with the approximated sensitivity
-        """
-
         for reg in self.reg.objfcts:
-            reg.cell_weights = reg.mapping * (self.wr)
-            if getattr(reg, "objfcts", None) is not None:
-                for obj in reg.objfcts:
-                    obj.cell_weights = obj.mapping * (self.wr)
+            wr += reg.mapping.deriv(self.invProb.model).T * (
+                    (reg.mapping * jtj_diag) / reg.objfcts[0].regmesh.vol ** 2.
+            )
+        wr /= wr.max()
+        wr += self.threshold
+        wr **= 0.5
+        for reg in self.reg.objfcts:
+            reg.cell_weights = reg.mapping * wr
 
     def validate(self, directiveList):
         # check if a beta estimator is in the list after setting the weights
diff --git a/SimPEG/potential_fields/magnetics/simulation.py b/SimPEG/potential_fields/magnetics/simulation.py
index 11b9d9f3a..3dea10b52 100644
--- a/SimPEG/potential_fields/magnetics/simulation.py
+++ b/SimPEG/potential_fields/magnetics/simulation.py
@@ -16,6 +16,7 @@ from SimPEG import Solver
 from SimPEG import props
 import properties
 from SimPEG.utils import mkvc, mat_utils, sdiag, setKwargs
+from SimPEG.utils.code_utils import deprecate_property
 
 
 class Simulation3DIntegral(BasePFSimulation):
@@ -28,16 +29,12 @@ class Simulation3DIntegral(BasePFSimulation):
         "Magnetic Susceptibility (SI)", default=1.0
     )
 
-    modelType = properties.StringChoice(
-        "Type of magnetization model",
-        choices=["susceptibility", "vector"],
-        default="susceptibility",
-    )
-
     is_amplitude_data = properties.Boolean(
         "Whether the supplied data is amplitude data", default=False
     )
 
+    _model_type: str = "scalar"
+
     def __init__(self, mesh, **kwargs):
         super().__init__(mesh, **kwargs)
         self._G = None
@@ -54,7 +51,7 @@ class Simulation3DIntegral(BasePFSimulation):
         """
         if getattr(self, "_M", None) is None:
 
-            if self.modelType == "vector":
+            if self.model_type == "vector":
                 self._M = sp.identity(self.nC) * self.survey.source_field.parameters[0]
 
             else:
@@ -80,7 +77,7 @@ class Simulation3DIntegral(BasePFSimulation):
         :parameter
         M: array (3*nC,) or (nC, 3)
         """
-        if self.modelType == "vector":
+        if self.model_type == "vector":
             self._M = sdiag(mkvc(M) * self.survey.source_field.parameters[0])
         else:
             M = M.reshape((-1, 3))
@@ -116,6 +113,31 @@ class Simulation3DIntegral(BasePFSimulation):
 
         return self._G
 
+    @property
+    def model_type(self) -> str:
+        """
+        Define the type of model. Choice of 'scalar' or 'vector' (3-components)
+        """
+        return self._model_type
+
+    @model_type.setter
+    def model_type(self, value: str):
+        if value not in ["scalar", "vector"]:
+            raise ValueError(
+                "'model_type' value should be a string: 'scalar' or 'vector'."
+                + f"Value {value} of type {type(value)} provided."
+            )
+
+        self._model_type = value
+
+    modelType = deprecate_property(
+        model_type,
+        "modelType",
+        new_name="model_type",
+        removal_version="0.16.0",
+        future_warn=True,
+    )
+
     @property
     def nD(self):
         """
diff --git a/SimPEG/regularization/sparse.py b/SimPEG/regularization/sparse.py
index 204c63b66..08d9a6b64 100755
--- a/SimPEG/regularization/sparse.py
+++ b/SimPEG/regularization/sparse.py
@@ -86,11 +86,12 @@ class SparseSmall(BaseSparse):
         if self.scale is None:
             self.scale = np.ones(self.mapping.shape[0])
 
+        weights = self.scale * self.regmesh.vol
+
         if self.cell_weights is not None:
-            return utils.sdiag((self.scale * self.cell_weights) ** 0.5) * R
+            weights *= self.cell_weights
 
-        else:
-            return utils.sdiag((self.scale * self.regmesh.vol) ** 0.5) * R
+        return utils.sdiag((weights ** 0.5)) * R
 
     def R(self, f_m):
         # if R is stashed, return that instead
@@ -180,7 +181,7 @@ class SparseDeriv(BaseSparse):
             self.scale = np.ones(self.mapping.shape[0])
 
         if self.space == "spherical":
-            Ave = getattr(self.regmesh, "aveCC2F{}".format(self.orientation))
+            ave_cc_f = getattr(self.regmesh, "aveCC2F{}".format(self.orientation))
 
             if getattr(self, "model", None) is None:
                 R = utils.speye(self.cellDiffStencil.shape[0])
@@ -189,15 +190,16 @@ class SparseDeriv(BaseSparse):
                 r = self.R(self.f_m)
                 R = utils.sdiag(r)
 
+            weights = self.scale * self.regmesh.vol
+
             if self.cell_weights is not None:
-                W = utils.sdiag((Ave * (self.scale * self.cell_weights)) ** 0.5) * R
+                weights *= self.cell_weights
 
-            else:
-                W = utils.sdiag((Ave * (self.scale * self.regmesh.vol)) ** 0.5) * R
+            W = utils.sdiag((ave_cc_f * weights ** 0.5)) * R
 
             theta = self.cellDiffStencil * (self.mapping * f_m)
-            dmdx = utils.mat_utils.coterminal(theta)
-            r = W * dmdx
+            dm_dx = utils.mat_utils.coterminal(theta)
+            r = W * dm_dx
 
         else:
             r = self.W * (self.mapping * f_m)
@@ -266,7 +268,7 @@ class SparseDeriv(BaseSparse):
             self.scale = np.ones(self.mapping.shape[0])
 
         if self.space == "spherical":
-            Ave = getattr(self.regmesh, "aveCC2F{}".format(self.orientation))
+            ave_cc_f = getattr(self.regmesh, "aveCC2F{}".format(self.orientation))
 
             if getattr(self, "model", None) is None:
                 R = utils.speye(self.cellDiffStencil.shape[0])
@@ -275,16 +277,15 @@ class SparseDeriv(BaseSparse):
                 r = self.R(self.f_m)
                 R = utils.sdiag(r)
 
-            if self.cell_weights is not None:
-                W = utils.sdiag(((Ave * (self.scale * self.cell_weights))) ** 0.5) * R
+            weights = self.scale * self.regmesh.vol
 
-            else:
-                W = utils.sdiag((Ave * (self.scale * self.regmesh.vol)) ** 0.5) * R
+            if self.cell_weights is not None:
+                weights *= self.cell_weights
 
+            W = utils.sdiag((ave_cc_f * weights ** 0.5)) * R
             theta = self.cellDiffStencil * (self.mapping * model)
-            dmdx = utils.mat_utils.coterminal(theta)
-
-            r = W * dmdx
+            dm_dx = utils.mat_utils.coterminal(theta)
+            r = W * dm_dx
 
         else:
             r = self.W * (self.mapping * model)
@@ -308,14 +309,14 @@ class SparseDeriv(BaseSparse):
 
         if self.space == "spherical":
             theta = self.cellDiffStencil * (self.mapping * f_m)
-            dmdx = utils.mat_utils.coterminal(theta)
+            dm_dx = utils.mat_utils.coterminal(theta)
 
         else:
 
             if self.gradientType == "total":
-                Ave = getattr(self.regmesh, "aveCC2F{}".format(self.orientation))
+                ave_cc_f = getattr(self.regmesh, "aveCC2F{}".format(self.orientation))
 
-                dmdx = np.abs(
+                dm_dx = np.abs(
                     self.regmesh.aveFx2CC
                     * self.regmesh.cellDiffxStencil
                     * (self.mapping * f_m)
@@ -323,7 +324,7 @@ class SparseDeriv(BaseSparse):
 
                 if self.regmesh.dim > 1:
 
-                    dmdx += np.abs(
+                    dm_dx += np.abs(
                         self.regmesh.aveFy2CC
                         * self.regmesh.cellDiffyStencil
                         * (self.mapping * f_m)
@@ -331,18 +332,18 @@ class SparseDeriv(BaseSparse):
 
                 if self.regmesh.dim > 2:
 
-                    dmdx += np.abs(
+                    dm_dx += np.abs(
                         self.regmesh.aveFz2CC
                         * self.regmesh.cellDiffzStencil
                         * (self.mapping * f_m)
                     )
 
-                dmdx = Ave * dmdx
+                dm_dx = ave_cc_f * dm_dx
 
             else:
-                dmdx = self.cellDiffStencil * (self.mapping * f_m)
+                dm_dx = self.cellDiffStencil * (self.mapping * f_m)
 
-        return dmdx
+        return dm_dx
 
     @property
     def cellDiffStencil(self):
@@ -353,7 +354,7 @@ class SparseDeriv(BaseSparse):
     @property
     def W(self):
 
-        Ave = getattr(self.regmesh, "aveCC2F{}".format(self.orientation))
+        ave_cc_f = getattr(self.regmesh, "aveCC2F{}".format(self.orientation))
 
         if getattr(self, "model", None) is None:
             R = utils.speye(self.cellDiffStencil.shape[0])
@@ -363,18 +364,13 @@ class SparseDeriv(BaseSparse):
             R = utils.sdiag(r)
         if self.scale is None:
             self.scale = np.ones(self.mapping.shape[0])
+
+        weights = self.scale * self.regmesh.vol
+
         if self.cell_weights is not None:
-            return (
-                utils.sdiag((Ave * (self.scale * self.cell_weights)) ** 0.5)
-                * R
-                * self.cellDiffStencil
-            )
-        else:
-            return (
-                utils.sdiag((Ave * (self.scale * self.regmesh.vol)) ** 0.5)
-                * R
-                * self.cellDiffStencil
-            )
+            weights *= self.cell_weights
+
+        return utils.sdiag((ave_cc_f * weights ** 0.5)) * R * self.cellDiffStencil
 
     @property
     def length_scales(self):
@@ -382,12 +378,12 @@ class SparseDeriv(BaseSparse):
         Normalized cell based weighting
 
         """
-        Ave = getattr(self.regmesh, "aveCC2F{}".format(self.orientation))
+        ave_cc_f = getattr(self.regmesh, "aveCC2F{}".format(self.orientation))
 
         if getattr(self, "_length_scales", None) is None:
             index = "xyz".index(self.orientation)
 
-            length_scales = Ave * (
+            length_scales = ave_cc_f * (
                 self.regmesh.Pac.T * self.regmesh.mesh.h_gridded[:, index]
             )
 
@@ -500,8 +496,8 @@ class Sparse(BaseComboRegularization):
 
         self.objfcts[0].norm = change["value"][:, 0]
         for i, objfct in enumerate(self.objfcts[1:]):
-            Ave = getattr(objfct.regmesh, "aveCC2F{}".format(objfct.orientation))
-            objfct.norm = Ave * change["value"][:, i + 1]
+            ave_cc_f = getattr(objfct.regmesh, "aveCC2F{}".format(objfct.orientation))
+            objfct.norm = ave_cc_f * change["value"][:, i + 1]
 
     @properties.observer("model")
     def _mirror_model_to_objfcts(self, change):
diff --git a/SimPEG/regularization/tikhonov.py b/SimPEG/regularization/tikhonov.py
index 2371f5520..d6b6ba182 100755
--- a/SimPEG/regularization/tikhonov.py
+++ b/SimPEG/regularization/tikhonov.py
@@ -54,7 +54,7 @@ class SimpleSmall(BaseRegularization):
         Weighting matrix
         """
         if self.cell_weights is not None:
-            return utils.sdiag(np.sqrt(self.cell_weights))
+            return utils.sdiag(np.sqrt(self.cell_weights * self.regmesh.vol))
         elif self._nC_residual != "*":
             return sp.eye(self._nC_residual)
         else:
@@ -122,7 +122,7 @@ class SimpleSmoothDeriv(BaseRegularization):
         )
         if self.cell_weights is not None:
 
-            W = utils.sdiag((Ave * (self.cell_weights)) ** 0.5) * W
+            W = utils.sdiag((Ave * (self.cell_weights * self.regmesh.vol)) ** 0.5) * W
         else:
             W = utils.sdiag((Ave * self.regmesh.vol) ** 0.5) * W
 
diff --git a/SimPEG/simulation.py b/SimPEG/simulation.py
index 50bc5ca59..5facecc85 100644
--- a/SimPEG/simulation.py
+++ b/SimPEG/simulation.py
@@ -166,7 +166,6 @@ class BaseSimulation(props.HasModel):
 
     ###########################################################################
     # Properties and observers
-
     @properties.observer("mesh")
     def _update_registry(self, change):
         self._REGISTRY.update(change["value"]._REGISTRY)
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 4ac36b363..28aea36aa 100644
--- a/examples/03-magnetics/plot_inv_mag_MVI_Sparse_TreeMesh.py
+++ b/examples/03-magnetics/plot_inv_mag_MVI_Sparse_TreeMesh.py
@@ -245,7 +245,7 @@ idenMap = maps.IdentityMap(nP=nC * 3)
 
 # Create the simulation
 simulation = magnetics.simulation.Simulation3DIntegral(
-    survey=survey, mesh=mesh, chiMap=idenMap, actInd=actv, modelType="vector"
+    survey=survey, mesh=mesh, chiMap=idenMap, actInd=actv, model_type="vector"
 )
 
 # Compute some data and add some random noise
diff --git a/tutorials/04-magnetics/plot_2a_magnetics_induced.py b/tutorials/04-magnetics/plot_2a_magnetics_induced.py
index 9a80c4524..df3ef2ea4 100644
--- a/tutorials/04-magnetics/plot_2a_magnetics_induced.py
+++ b/tutorials/04-magnetics/plot_2a_magnetics_induced.py
@@ -172,7 +172,7 @@ plt.show()
 simulation = magnetics.simulation.Simulation3DIntegral(
     survey=survey,
     mesh=mesh,
-    modelType="susceptibility",
+    model_type="scalar",
     chiMap=model_map,
     actInd=ind_active,
     store_sensitivities="forward_only",
diff --git a/tutorials/04-magnetics/plot_2b_magnetics_mvi.py b/tutorials/04-magnetics/plot_2b_magnetics_mvi.py
index fd4740c23..239a60a90 100644
--- a/tutorials/04-magnetics/plot_2b_magnetics_mvi.py
+++ b/tutorials/04-magnetics/plot_2b_magnetics_mvi.py
@@ -230,7 +230,7 @@ simulation = magnetics.simulation.Simulation3DIntegral(
     mesh=mesh,
     chiMap=model_map,
     actInd=ind_active,
-    modelType="vector",
+    model_type="vector",
     store_sensitivities="forward_only",
 )
 
diff --git a/tutorials/04-magnetics/plot_inv_2a_magnetics_induced.py b/tutorials/04-magnetics/plot_inv_2a_magnetics_induced.py
index dfe1dd5fe..14fbffd47 100644
--- a/tutorials/04-magnetics/plot_inv_2a_magnetics_induced.py
+++ b/tutorials/04-magnetics/plot_inv_2a_magnetics_induced.py
@@ -231,7 +231,7 @@ starting_model = background_susceptibility * np.ones(nC)
 simulation = magnetics.simulation.Simulation3DIntegral(
     survey=survey,
     mesh=mesh,
-    modelType="susceptibility",
+    model_type="scalar",
     chiMap=model_map,
     actInd=ind_active,
 )
diff --git a/tutorials/05-dcr/plot_inv_2_dcr2d.py b/tutorials/05-dcr/plot_inv_2_dcr2d.py
index 2d1901ed5..f1e5514c4 100644
--- a/tutorials/05-dcr/plot_inv_2_dcr2d.py
+++ b/tutorials/05-dcr/plot_inv_2_dcr2d.py
@@ -270,7 +270,7 @@ starting_conductivity_model = background_conductivity * np.ones(nC)
 
 # Define the problem. Define the cells below topography and the mapping
 simulation = dc.simulation_2d.Simulation2DNodal(
-    mesh, survey=survey, sigmaMap=conductivity_map, Solver=Solver
+    mesh, survey=survey, sigmaMap=conductivity_map, solver=Solver, storeJ=True
 )
 
 #######################################################################
diff --git a/tutorials/05-dcr/plot_inv_2_dcr2d_irls.py b/tutorials/05-dcr/plot_inv_2_dcr2d_irls.py
index b2890d9aa..f57697335 100644
--- a/tutorials/05-dcr/plot_inv_2_dcr2d_irls.py
+++ b/tutorials/05-dcr/plot_inv_2_dcr2d_irls.py
@@ -275,7 +275,7 @@ starting_conductivity_model = background_conductivity * np.ones(nC)
 
 # Define the problem. Define the cells below topography and the mapping
 simulation = dc.simulation_2d.Simulation2DNodal(
-    mesh, survey=survey, sigmaMap=conductivity_map, Solver=Solver
+    mesh, survey=survey, sigmaMap=conductivity_map, solver=Solver, storeJ=True
 )
 
 #######################################################################
diff --git a/tutorials/05-dcr/plot_inv_3_dcr3d.py b/tutorials/05-dcr/plot_inv_3_dcr3d.py
index 6a402030d..f5ec005de 100644
--- a/tutorials/05-dcr/plot_inv_3_dcr3d.py
+++ b/tutorials/05-dcr/plot_inv_3_dcr3d.py
@@ -276,7 +276,7 @@ starting_conductivity_model = background_conductivity * np.ones(nC)
 #
 
 dc_simulation = dc.simulation.Simulation3DNodal(
-    mesh, survey=dc_survey, sigmaMap=conductivity_map, Solver=Solver
+    mesh, survey=dc_survey, sigmaMap=conductivity_map, solver=Solver, storeJ=True
 )
 
 #################################################################
@@ -300,13 +300,7 @@ dc_data_misfit = data_misfit.L2DataMisfit(data=dc_data, simulation=dc_simulation
 
 # Define the regularization (model objective function)
 dc_regularization = regularization.Simple(
-    mesh,
-    indActive=ind_active,
-    mref=starting_conductivity_model,
-    alpha_s=1e-2,
-    alpha_x=1,
-    alpha_y=1,
-    alpha_z=1,
+    mesh, indActive=ind_active, mref=starting_conductivity_model,
 )
 
 dc_regularization.mrefInSmooth = True  # Include reference model in smoothness
@@ -468,6 +462,7 @@ cbar = mpl.colorbar.ColorbarBase(
     ax2, cmap=mpl.cm.viridis, norm=norm, orientation="vertical", format="$10^{%.1f}$"
 )
 cbar.set_label("Conductivity [S/m]", rotation=270, labelpad=15, size=12)
+plt.show()
 
 #######################################################################
 # Plotting Normalized Data Misfit or Predicted DC Data
diff --git a/tutorials/06-ip/plot_inv_2_dcip2d.py b/tutorials/06-ip/plot_inv_2_dcip2d.py
index 46314faf6..0eebb7391 100644
--- a/tutorials/06-ip/plot_inv_2_dcip2d.py
+++ b/tutorials/06-ip/plot_inv_2_dcip2d.py
@@ -287,7 +287,7 @@ starting_conductivity_model = background_conductivity * np.ones(nC)
 
 # Define the problem. Define the cells below topography and the mapping
 dc_simulation = dc.simulation_2d.Simulation2DNodal(
-    mesh, survey=dc_survey, sigmaMap=conductivity_map, Solver=Solver
+    mesh, survey=dc_survey, sigmaMap=conductivity_map, solver=Solver, storeJ=True
 )
 
 #######################################################################
@@ -525,7 +525,8 @@ ip_simulation = ip.simulation_2d.Simulation2DNodal(
     survey=ip_survey,
     etaMap=chargeability_map,
     sigma=conductivity_map * recovered_conductivity_model,
-    Solver=Solver,
+    solver=Solver,
+    storeJ=True,
 )
 
 #####################################################
diff --git a/tutorials/06-ip/plot_inv_3_dcip3d.py b/tutorials/06-ip/plot_inv_3_dcip3d.py
index d9aca8175..4e48fdc1d 100644
--- a/tutorials/06-ip/plot_inv_3_dcip3d.py
+++ b/tutorials/06-ip/plot_inv_3_dcip3d.py
@@ -325,7 +325,7 @@ starting_conductivity_model = background_conductivity * np.ones(nC)
 #
 
 dc_simulation = dc.simulation.Simulation3DNodal(
-    mesh, survey=dc_survey, sigmaMap=conductivity_map, Solver=Solver
+    mesh, survey=dc_survey, sigmaMap=conductivity_map, solver=Solver, storeJ=True
 )
 
 #################################################################
@@ -349,13 +349,7 @@ dc_data_misfit = data_misfit.L2DataMisfit(data=dc_data, simulation=dc_simulation
 
 # Define the regularization (model objective function)
 dc_regularization = regularization.Simple(
-    mesh,
-    indActive=ind_active,
-    mref=starting_conductivity_model,
-    alpha_s=0.01,
-    alpha_x=1,
-    alpha_y=1,
-    alpha_z=1,
+    mesh, indActive=ind_active, mref=starting_conductivity_model,
 )
 
 dc_regularization.mrefInSmooth = True  # Include reference model in smoothness
@@ -509,6 +503,7 @@ cbar = mpl.colorbar.ColorbarBase(
     ax2, cmap=mpl.cm.viridis, norm=norm, orientation="vertical", format="$10^{%.1f}$"
 )
 cbar.set_label("Conductivity [S/m]", rotation=270, labelpad=15, size=12)
+plt.show()
 
 #######################################################################
 # Plotting Normalized Data Misfit or Predicted DC Data
@@ -595,7 +590,8 @@ ip_simulation = ip.simulation.Simulation3DNodal(
     survey=ip_survey,
     etaMap=chargeability_map,
     sigma=conductivity_map * recovered_conductivity_model,
-    Solver=Solver,
+    solver=Solver,
+    storeJ=True,
 )
 
 #################################################
@@ -744,6 +740,7 @@ cbar = mpl.colorbar.ColorbarBase(
 )
 cbar.set_label("Intrinsic Chargeability [V/V]", rotation=270, labelpad=15, size=12)
 
+plt.show()
 ##########################################################
 # Plotting Normalized Data Misfit or Predicted IP Data
 # ----------------------------------------------------


Test patch:
diff --git a/tests/base/test_regularization.py b/tests/base/test_regularization.py
index eef0757e8..3b746b520 100644
--- a/tests/base/test_regularization.py
+++ b/tests/base/test_regularization.py
@@ -286,7 +286,7 @@ class RegularizationTests(unittest.TestCase):
         )
 
         objfct = objective_function.L2ObjectiveFunction(
-            W=utils.sdiag(np.sqrt(cell_weights)), mapping=wires.sigma
+            W=utils.sdiag(np.sqrt(cell_weights * mesh.cell_volumes)), mapping=wires.sigma
         )
 
         self.assertTrue(reg(m) == objfct(m))
diff --git a/tests/dask/test_mag_MVI_Octree.py b/tests/dask/test_mag_MVI_Octree.py
index cfc22f655..2282aca20 100644
--- a/tests/dask/test_mag_MVI_Octree.py
+++ b/tests/dask/test_mag_MVI_Octree.py
@@ -92,7 +92,7 @@ class MVIProblemTest(unittest.TestCase):
         sim = mag.Simulation3DIntegral(
             self.mesh,
             survey=survey,
-            modelType="vector",
+            model_type="vector",
             chiMap=idenMap,
             actInd=actv,
             store_sensitivities="disk",
diff --git a/tests/pf/test_mag_MVI_Octree.py b/tests/pf/test_mag_MVI_Octree.py
index 0b65b2102..d029f6574 100644
--- a/tests/pf/test_mag_MVI_Octree.py
+++ b/tests/pf/test_mag_MVI_Octree.py
@@ -91,7 +91,7 @@ class MVIProblemTest(unittest.TestCase):
         sim = mag.Simulation3DIntegral(
             self.mesh,
             survey=survey,
-            modelType="vector",
+            model_type="vector",
             chiMap=idenMap,
             actInd=actv,
             store_sensitivities="disk",