Repository: simpeg/simpeg
Pull request #1572 (2024-11-02 22:59:16)

Bugfix for TDEM magnetic dipole sources 
<!-- 
Thanks for contributing a pull request to SimPEG!
Remember to use a personal fork of SimPEG to propose changes.

Check out the stages of a pull request at
https://docs.simpeg.xyz/content/getting_started/contributing/pull-requests.html

Note that we are a team of volunteers and we appreciate your
patience during the review process.

Again, thanks for contributing!

Feel free to remove lines from this template that do not apply to you pull request.
-->

#### Summary
Bugfix so that the initial fields produced for `db/dt`, `dh/dt`, from the B field or H field simulations are zero. There was previously a bug where we set the electric source term to `Zero()` when infact it should have been evaluated. This would only effect `db/dt` (or `dh/dt`) that were measured before the first time-step, but it also impacted visualizations of the fields. 

To catch this issue, I have extended our cross-check tests for tdem simulation to also include an evaluation of the data at t=0. Running this test on the current main version of SimPEG will fail. 

As an example, using the `tdem.simulation.Simulation3DMagneticFluxDensity` for a loop source and plotting db/dt gives: 

**Before the bugfix**
Crazy db/dt! 
![image](https://github.com/user-attachments/assets/e948331a-9b60-4930-8b96-9f8dfe9ae7de)


**After**
Zero -- as it should be
![image](https://github.com/user-attachments/assets/a85f24cc-a8bd-4f10-8423-e43d49dd4321)


#### PR Checklist
* [ ] If this is a work in progress PR, set as a Draft PR
* [ ] Linted my code according to the [style guides](https://docs.simpeg.xyz/latest/content/getting_started/contributing/code-style.html).
* [ ] Added [tests](https://docs.simpeg.xyz/latest/content/getting_started/contributing/testing.html) to verify changes to the code.
* [ ] Added necessary documentation to any new functions/classes following the
      expect [style](https://docs.simpeg.xyz/latest/content/getting_started/contributing/documentation.html).
* [ ] Marked as ready for review (if this is was a draft PR), and converted 
      to a Pull Request
* [ ] Tagged ``@simpeg/simpeg-developers`` when ready for review.

<!--
Once all tests pass and the code has been reviewed and approved, it will be merged into main
-->


codecov[bot]: ## [Codecov](https://app.codecov.io/gh/simpeg/simpeg/pull/1572?dropdown=coverage&src=pr&el=h1&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=simpeg) Report
All modified and coverable lines are covered by tests :white_check_mark:
> Project coverage is 86.42%. Comparing base [(`6c5f7d5`)](https://app.codecov.io/gh/simpeg/simpeg/commit/6c5f7d5d6725579bc6fceb40e0920e528e8c299b?dropdown=coverage&el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=simpeg) to head [(`59b8545`)](https://app.codecov.io/gh/simpeg/simpeg/commit/59b8545df96e38bcde8e37d17546c92545ea71e6?dropdown=coverage&el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=simpeg).
> Report is 3 commits behind head on main.

<details><summary>Additional details and impacted files</summary>


```diff
@@            Coverage Diff             @@
##             main    #1572      +/-   ##
==========================================
+ Coverage   86.13%   86.42%   +0.28%     
==========================================
  Files         400      406       +6     
  Lines       51492    52407     +915     
  Branches     4903     4987      +84     
==========================================
+ Hits        44354    45293     +939     
+ Misses       5744     5689      -55     
- Partials     1394     1425      +31     
```

</details>

[:umbrella: View full report in Codecov by Sentry](https://app.codecov.io/gh/simpeg/simpeg/pull/1572?dropdown=coverage&src=pr&el=continue&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=simpeg).   
:loudspeaker: Have feedback on the report? [Share it here](https://about.codecov.io/codecov-pr-comment-feedback/?utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=simpeg).

jcapriot: I mean, what exactly is dB/dT at the moment the step-off happens? It depends on which direction you approach that time! So I don't know which one is "more" correct.
lheagy: Right at that moment when the transmitter shuts off, you are right that db/dt, isn't defined, but a couple of things... for the definition of initial fields, the limit is always from the left (while the transmitter is still on), so zero is a valid limit. Second, the current solution isn't physical (so isn't the limit from either side!) 

![image](https://github.com/user-attachments/assets/4f526c96-3d8f-4471-ab3f-1a4812349512)

So having the solution be zero is both physical and consistent with our definition of initial fields 
jcapriot: Thanks for the explanation there. So then is your change assuming that if the waveform has initial fields, the amplitude is 1?
lheagy: If the waveform has initial fields, those will be either `e` or `b` or some combination of both, but will be static (whatever amplitude they are). Because they are static, `dbdt` should be `Zero` initially, as there is no time variation. 
lheagy: is there any chance of getting this into a patch release relatively soon? I will be using this code in a demo this week. If we don't get it in before then, that is fine, but it would be good to have it out shortly after to avoid potential confusion
jcapriot: Ah, I see what was happening then! The source was zero before because in the formulation it makes use of the time derivative of `s_e`, so before the source was ever on it wouldn’t matter as long as it was constant. But then at the moment the source turns on, it had a big impulse on, causing those spurious dbdt’s.

Your fix at least fixes that portion, but also would have issues again if the source waveform is not 1 at the start.
For example, the piecewise linear waveform could easily have a non-unitary initial point.

From what I could tell, all of the waveforms that have initial fields support evaluating themselves with a time before its start time.
lheagy: The source was zero before, but it wasn't used in the b, h simulations (exactly as you say -- because we were using the derivative of the source term for computing in the b, h formulations). So we never caught the issue. 

However, when computing the fields it matters. This should be fine for any waveform (we test the e formulation and  it makes use of the `s_e` source term) for computing the response due to an arbitrary waveform

https://github.com/simpeg/simpeg/blob/8c046219caa74efcbfd14177d91da3f2c99ddb86/tests/em/tdem/test_TDEM_forward_Analytic_RawWaveform.py#L14

Previously, we had an if-else statement that treated b,h formulations differently for this source than for e,j, but there is no reason for that. So this pr simply gets rid of the if-else statements that were treating them differently 
lheagy: One other related note, for the figure that I showed, that is bd/dt computed at t=0 (so the moment before the transmitter shuts off) -- it wasn't an effect due to transmitter turn off, it is because we needed to have the source term included when we solve for the initial `e` fields that are used to compute `dbdt`

https://github.com/simpeg/simpeg/blob/8c046219caa74efcbfd14177d91da3f2c99ddb86/simpeg/electromagnetics/time_domain/fields.py#L246-L251 

https://github.com/simpeg/simpeg/blob/8c046219caa74efcbfd14177d91da3f2c99ddb86/simpeg/electromagnetics/time_domain/fields.py#L226-L232
lheagy: This notebook is one you can test out before / after the pr if you like: https://github.com/ubcgif/2024-aqgeo-meeting-em-tutorial/blob/notebooks/1-forward-simulation.ipynb

The widget at the bottom is the one I grabbed the screenshots from. 

Patch:
diff --git a/simpeg/electromagnetics/time_domain/sources.py b/simpeg/electromagnetics/time_domain/sources.py
index 6d18fb42d..afef51aa4 100644
--- a/simpeg/electromagnetics/time_domain/sources.py
+++ b/simpeg/electromagnetics/time_domain/sources.py
@@ -1428,11 +1428,7 @@ class MagDipole(BaseTDEMSrc):
                 self.waveform.has_initial_fields is True
                 and time < simulation.time_steps[1]
             ):
-                if simulation._fieldType == "b":
-                    return Zero()
-                elif simulation._fieldType == "e":
-                    # Compute s_e from vector potential
-                    return C.T * (MfMui * b)
+                return C.T * (MfMui * b)
             else:
                 return C.T * (MfMui * b) * self.waveform.eval(time)
 
@@ -1443,11 +1439,7 @@ class MagDipole(BaseTDEMSrc):
                 self.waveform.has_initial_fields is True
                 and time < simulation.time_steps[1]
             ):
-                if simulation._fieldType == "h":
-                    return Zero()
-                elif simulation._fieldType == "j":
-                    # Compute s_e from vector potential
-                    return C * h
+                return C * h
             else:
                 return C * h * self.waveform.eval(time)
 


Test patch:
diff --git a/tests/em/tdem/test_TDEM_crosscheck.py b/tests/em/tdem/test_TDEM_crosscheck.py
index 4c1de54eb..905b1115d 100644
--- a/tests/em/tdem/test_TDEM_crosscheck.py
+++ b/tests/em/tdem/test_TDEM_crosscheck.py
@@ -35,7 +35,7 @@ def setUp_TDEM(
     )
     mapping = maps.ExpMap(mesh) * maps.SurjectVertical1D(mesh) * activeMap
 
-    rxtimes = np.logspace(-4, -3, 20)
+    rxtimes = np.hstack([np.r_[0], np.logspace(-4, -3, 20)])
 
     if waveform.upper() == "RAW":
         t0 = 0.006