Bugfix for TDEM magnetic dipole sources

[lheagy]

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!

After
Zero -- as it should be

PR Checklist

• If this is a work in progress PR, set as a Draft PR
• Linted my code according to the style guides.
• Added tests to verify changes to the code.
• Added necessary documentation to any new functions/classes following the
  expect style.
• 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.

[codecov]

Codecov Report

All modified and coverable lines are covered by tests ✅

Project coverage is 86.44%. Comparing base (8c04621) to head (930c350).
Report is 1 commits behind head on main.

Additional details and impacted files

@@            Coverage Diff             @@
##             main    #1572      +/-   ##
==========================================
- Coverage   86.44%   86.44%   -0.01%     
==========================================
  Files         407      407              
  Lines       52475    52469       -6     
  Branches     4997     4993       -4     
==========================================
- Hits        45363    45357       -6     
- Misses       5686     5687       +1     
+ Partials     1426     1425       -1     

☔ View full report in Codecov by Sentry.
📢 Have feedback on the report? Share it here.

[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!)

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

simpeg/tests/em/tdem/test_TDEM_forward_Analytic_RawWaveform.py

Line 14 in 8c04621

def halfSpaceProblemAnaDiff(

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

simpeg/simpeg/electromagnetics/time_domain/fields.py

Lines 246 to 251 in 8c04621

	def _e(self, bSolution, source_list, tInd):
	e = self._MeSigmaI * (self._edgeCurl.T * (self._MfMui * bSolution))
	for i, src in enumerate(source_list):
	s_e = src.s_e(self.simulation, self._times[tInd])
	e[:, i] = e[:, i] - self._MeSigmaI * s_e
	return e

simpeg/simpeg/electromagnetics/time_domain/fields.py

Lines 226 to 232 in 8c04621

	def _dbdt(self, bSolution, source_list, tInd):
	# self._timeMesh.face_divergence
	dbdt = -self._edgeCurl * self._e(bSolution, source_list, tInd)
	for i, src in enumerate(source_list):
	s_m = src.s_m(self.simulation, self._times[tInd])
	dbdt[:, i] = dbdt[:, i] + s_m
	return dbdt

[lheagy]

this if-else statement wasn't necessary

[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.

[jcapriot]

Ok, thanks for sharing the notebook, was helpful! looks good to me now then.

[lheagy]

Thanks for taking the time to go through this one Joe! It was a bit of a subtle bug. I appreciate your help 😄