Add parameter tol_gradient to Simulation (#331)

CHANGELOG.rst

@@ -9,7 +9,22 @@ Changelog
 latest
 ------
 
-TODO: some text about empymod, scooby.
+- New hard dependency: ``empymod`` changed from being a soft dependency to
+  being a hard dependency.
+
+- Simulation
+
+  - A new parameter ``tol_gradient`` can be provided in the dict
+    ``solver_opts``; by default it is set to the value of ``tol``. ``tol`` is
+    the value used for ``compute`` (the forward), ``tol_gradient`` is used for
+    ``gradient``/``jtvec`` and ``jvec`` (the gradient). In inversions, one can
+    set the tolerance for the gradient often to a lower value, which saves
+    computation time.
+
+- Electrodes
+
+  - Fixed ``TxMagneticDipole``-representation and improved documentation of the
+    magnetic sources.
 
 - Maintenance:
 

docs/manual/cli.rst

@@ -78,6 +78,7 @@ remove the comment signs to use them.
   # linerelaxation =       # bool
   # cycle =                # string
   # tol =                  # float
+  # tol_gradient =         # float
   # verb =                 # int
   # maxit =                # int
   # nu_init =              # int

emg3d/cli/parser.py

@@ -329,7 +329,7 @@ def parse_config_file(args_dict):
                 solver[key] = cfg.get('solver_opts', key)
 
         # Check for floats.
-        for key in ['tol', ]:
+        for key in ['tol', 'tol_gradient']:
             if cfg.has_option('solver_opts', key):
                 _ = all_solver.pop(key)
                 solver[key] = float(cfg.get('solver_opts', key))

emg3d/simulations.py

@@ -149,6 +149,13 @@ class Simulation:
         parameter that is accepted by the :func:`emg3d.solver.solve` except for
         ``model``, ``sfield``, ``efield``, ``return_info``, and ``log``.
 
+        In addition to the regular parameter ``tol``, one can provide a
+        parameter ``tol_gradient``. This tolerance will be used when calling
+        ``gradient``/``jtvec`` and ``jvec``. By default, it is set to the same
+        value as ``tol``, which is used for ``compute``. However, for
+        inversions it is usually possible to relax this tolerance (e.g.,
+        ``tol=1e-6``, ``tol_gradient=1e-3``).
+
     verb : int, default: 0
         Level of verbosity. Possible options:
 
@@ -260,13 +267,16 @@ def __init__(self, survey, model, max_workers=4, gridding='single',
         self.receiver_interpolation = kwargs.pop(
                 'receiver_interpolation', 'cubic')
 
-        # Assemble solver_opts.
+        # Assemble solver_opts and store the tolerances separately.
         self.solver_opts = {
                 'verb': 1,  # Default verbosity, can be overwritten.
                 'log': -1,  # Default only log, can be overwritten.
                 **kwargs.pop('solver_opts', {}),  # User setting.
                 'return_info': True,  # return_info=True is forced.
         }
+        self.tol_forward = self.solver_opts.get('tol', 1e-6)
+        self.tol_gradient = self.solver_opts.pop(
+                'tol_gradient', self.tol_forward)
 
         # Initiate dictionaries and other values with None's.
         self._dict_grid = self._dict_initiate
@@ -425,6 +435,7 @@ def to_dict(self, what='computed', copy=False):
             raise TypeError(f"Unrecognized `what`: {what}.")
 
         # Initiate dict with input parameters.
+        self.solver_opts['tol'] = self.tol_forward
         out = {
             '__class__': self.__class__.__name__,
             'survey': self.survey.to_dict(),
@@ -440,6 +451,7 @@ def to_dict(self, what='computed', copy=False):
             'layered': self.layered,
             'layered_opts': self.layered_opts,
             'receiver_interpolation': self.receiver_interpolation,
+            'tol_gradient': self.tol_gradient,
             'file_dir': self.file_dir,
             '_input_sc2': self._input_sc2,
         }
@@ -504,6 +516,9 @@ def from_dict(cls, inp):
         cls_inp['tqdm_opts'] = inp.pop('tqdm_opts', {})
         cls_inp['layered'] = inp.pop('layered', False)
         cls_inp['layered_opts'] = inp.pop('layered_opts', {})
+        cls_inp['solver_opts'] = cls_inp['solver_opts'].copy()
+        cls_inp['solver_opts']['tol_gradient'] = inp.pop(
+                'tol_gradient', cls_inp['solver_opts'].get('tol', 1e-6))
 
         # Instantiate the class.
         out = cls(**cls_inp)
@@ -839,6 +854,7 @@ def collect_efield_inputs(inp):
                 'efield': self._dict_get('efield', source, freq),
                 'solver_opts': self.solver_opts,
             }
+            data['solver_opts']['tol'] = self.tol_forward
             return self._data_or_file('efield', source, freq, data)
 
         # Compute fields in parallel.
@@ -1195,6 +1211,7 @@ def collect_bfield_inputs(inp):
                 'efield': self._dict_get('bfield', source, freq),
                 'solver_opts': self.solver_opts
             }
+            data['solver_opts']['tol'] = self.tol_gradient
             return self._data_or_file('bfield', source, freq, data)
 
         # Compute fields in parallel.
@@ -1355,6 +1372,7 @@ def collect_gfield_inputs(inp, vector=vector):
                 'efield': None,
                 'solver_opts': self.solver_opts,
             }
+            data['solver_opts']['tol'] = self.tol_gradient
             return self._data_or_file('gfield', source, freq, data)
 
         # Compute fields (A^-1 * G * vector) in parallel.

tests/test_simulations.py

@@ -996,6 +996,35 @@ def test_jtvec_gradient(self):
         j = sim.jtvec(vector=sim.survey.data.residual*sim.survey.data.weights)
         assert_allclose(g, j)
 
+    @pytest.mark.skipif(discretize is None, reason="discretize not installed.")
+    @pytest.mark.skipif(h5py is None, reason="h5py not installed.")
+    def test_tol_gradient(self, capsys):
+        sim = simulations.Simulation(model=self.model_init, **self.sim_inp)
+        assert sim.tol_forward == 1e-6
+        assert sim.tol_gradient == 1e-6
+        sim_inp = self.sim_inp.copy()
+        sim_inp['solver_opts']['tol'] = 1e-2
+        sim_inp['solver_opts']['tol_gradient'] = 1e-1
+        sim = simulations.Simulation(model=self.model_init, **sim_inp)
+        assert sim.tol_forward == 1e-2
+        assert sim.tol_gradient == 1e-1
+        simdict = sim.to_dict()
+        sim2 = simulations.Simulation.from_dict(simdict)
+        assert sim.tol_forward == 1e-2
+        assert sim.tol_gradient == 1e-1
+        sim2.solver_opts['verb'] = 3
+        sim2.solver_opts['log'] = 1
+        _, _ = capsys.readouterr()  # empty
+        sim2.compute()
+        out, _ = capsys.readouterr()
+        assert "tol       : 0.01" in out
+        sim2.gradient
+        out, _ = capsys.readouterr()
+        assert "tol       : 0.1" in out
+        sim2.jvec(sim.model.property_x)
+        out, _ = capsys.readouterr()
+        assert "tol       : 0.1" in out
+
 
 def test_all_dir():
     assert set(simulations.__all__) == set(dir(simulations))