Skip to content

FDTD: a sine excitation repeated its signal one timestep early - #249

Open
MrHighVoltage wants to merge 1 commit into
thliebig:masterfrom
MrHighVoltage:excitation-period-rounding
Open

MrHighVoltage wants to merge 1 commit into
thliebig:masterfrom
MrHighVoltage:excitation-period-rounding

Conversation

@MrHighVoltage

@MrHighVoltage MrHighVoltage commented Nov 5, 2026 •

Copy link
Copy Markdown
Contributor

AI: generated with Claude Opus 5.5 (Claude Code), reviewed and tested by me.

A sine excitation repeats its signal every int(T/dT) timesteps. T/dT is a whole number by construction, but in floating point it can come out just below it, e.g. 432.99999999999994 instead of 433, and int() then cuts off a whole timestep. At 2.4 GHz this happens for 554 of the timestep counts from 2 to 10000. The sine then skips a sample every period, and the steady-state detection never ends the run.

Why the period is a whole number of timesteps: Operator::CalcECOperator() shrinks the timestep after the CFL limit, a forced timestep and the timestep factor (operator.cpp#L1125-L1131):

unsigned int TS = ceil(m_Exc->GetSignalPeriod()/dT);
double new_dT = m_Exc->GetSignalPeriod()/TS;

Only the sine has a signal period (excitation.cpp#L73). In exact arithmetic T/dT = TS. In double, dT = T/TS is rounded, and T/dT comes back as TS plus or minus about 1e-16 relative (at 2.4 GHz at most 9e-13 timesteps, for TS up to 10000). Rounding gives TS exactly; there is no real fraction of a timestep to round away.

How the period is used: the excitation plays a table of sin(2π f0 n dT), two periods long (excitation.cpp#L293), at index (numTS - delay) % p (engine_ext_excitation.cpp#L45-L55). With p = TS the wrap is seamless. With p = TS-1:

  • the last sample of each period is skipped: the source runs at TS/(TS-1) of f0 (+0.23 % at TS 433), with a phase step of 360°/TS (0.83°) at each wrap;
  • the steady-state detection, which openEMS adds to every sine run, rounds the period (operator_ext_steadystate.cpp#L50) and compares windows of TS timesteps. The signal moves by one sample against them every period, so they never match, and the run goes on to the max. number of timesteps;
  • the period printed at the start is one timestep short (openems.cpp#L1457).

The same int() is in the TF/SF extension (engine_ext_tfsf.cpp#L72) and in HIP and Metal (hip_ext_excitation.hip#L139, hip_ext_tfsf.hip#L163, metal_ext_excitation.mm#L140, metal_ext_tfsf.mm#L167). The host extensions serve the basic, multithreaded, CPU and gpu-reference engines.

Fix: Excitation::GetSignalPeriodTS() returns lround(T/dT), and all eight places use it.

Test: the new python/Tests/Sinusoidal_Period.py drives a TEM wave in a 1D channel (0.5 mm mesh, PEC and PMC walls, PML at both ends) at 2.4 GHz, 433 timesteps per period, with a 200-period cap. On master the period is printed as 432 timesteps, the probe phasor turns by 0.83° per period, and the run goes on to the cap (86.6k timesteps). Fixed, the steady-state detection ends the run after 4 periods (1.8k timesteps), and the probe phasor changes by 0.005 % between the last two periods. The same on the default (multithreaded), basic, cpu, gpu-reference and HIP (RX 6800) engines, each of which runs to the cap on master. SteadyState.py, PlaneWave_TFSF.py, Engine_Compare.py (--engine=cpu and --engine=gpu) and the C++ unit tests still pass. Metal was not compiled.

Note for the gpu-cellmodel branch: gpu_ext.cpp#L143 (excitation) and #L448 (TF/SF) copy the int() for the reference, HIP and Metal backends, and replace the HIP and Metal hunks here:

-			s.op.period = int(m_Op->m_Exc->GetSignalPeriod()/m_Op->m_Exc->GetTimestep());
+			s.op.period = m_Op->m_Exc->GetSignalPeriodTS();
...
-				P.op.period = int(op->m_Exc->GetSignalPeriod()/op->m_Exc->GetTimestep());
+				P.op.period = op->m_Exc->GetSignalPeriodTS();

With that and the helper, Sinusoidal_Period.py passes on gpu and gpu-reference there (on a local branch on fc20a33).

The operator shrinks the timestep so that a period T of a periodic
excitation is a whole number TS of timesteps, dT = T/TS. The excitation
and the TF/SF, on the host and in the HIP and Metal kernels, then
repeated their signal every int(T/dT) timesteps, but T/(T/TS) can come
out just below TS, e.g. TS 433 at 2.4 GHz on the mesh of the new test;
at 2.4 GHz for 554 of the TS from 2 to 10000. The signal then
skipped a sample every period and ran at TS/(TS-1) times its frequency,
and the steady-state detection, which rounds the period, never saw two
equal periods: the run went on to the max. number of timesteps.

Excitation::GetSignalPeriodTS() rounds, and every user of the period in
timesteps takes it from there: the host extensions (the basic,
multithreaded, CPU and gpu-reference engines), the HIP and Metal kernels
and the period openEMS prints. Nothing changes where int() was right.

python/Tests/Sinusoidal_Period.py runs a 2.4 GHz sine on such a mesh: it
must end by the steady-state detection, with a periodic probe signal.

Generated-by: Claude Opus 5.5 (Claude Code)
Signed-off-by: Georg Zachl <it@zachlge.org>
@MrHighVoltage
MrHighVoltage force-pushed the excitation-period-rounding branch from f87b02a to bceabe9 Compare November 9, 2026 22:40

@thliebig thliebig left a comment

Copy link
Copy Markdown
Owner

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The python test is very extensive again for this little rounding error? Not sure we really need it?

Please rebase/adept to current master. Otherwise looks good, thanks

Comment thread FDTD/excitation.h
Comment on lines +86 to +89
/*!
Rounded, not truncated: the operator shrinks the timestep to a whole number of timesteps
per period, but the period divided by that timestep can come out just below it.
*/

Copy link
Copy Markdown
Owner

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

We really need not to document the original issue in the API doc...

Comment thread CHANGELOG.md
Comment on lines +76 to +78
- A sinusoidal excitation could repeat its signal one timestep early, e.g. at
2.4 GHz on some meshes: it ran slightly above its frequency, and the
steady-state detection never ended the simulation.

Copy link
Copy Markdown
Owner

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This was "just" a small issue, not sure it is worth 3 lines in the change log if any?

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

2 participants