Skip to content

Commit

Permalink
Waves: fix Hermitian symmetry for derivatives (#90)
Browse files Browse the repository at this point in the history 
- Correct Nyquist term in derivative amplitudes (must be zero).
- Enable tests to check derivative and displacement amplitudes are Hermitian.

Signed-off-by: Rhys Mainwaring <[email protected]>
  • Loading branch information
@srmainwaring
srmainwaring authored Nov 24, 2022
1 parent 214fd87 commit 9ee03da
Show file tree
Hide file tree
Showing 3 changed files with 210 additions and 11 deletions.
19 changes: 18 additions & 1 deletion gz-waves/src/WaveSimulationFFT.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -319,7 +319,7 @@ namespace waves
{
double kx = kx_fft_[ikx];
double kx2 = kx*kx;
for (int iky = 0; iky < ny_/2+1; ++iky)
for (int iky = 0; iky < ny_/2 + 1; ++iky)
{
double ky = ky_fft_[iky];
double ky2 = ky*ky;
Expand All @@ -333,6 +333,14 @@ namespace waves
complex hikx = hi * kx;
complex hiky = hi * ky;

// Nyquist terms for derivatives must be zero.
// For an explanation see:
// https://math.mit.edu/~stevenj/fft-deriv.pdf
if (ikx == nx_ / 2)
hikx = czero;
if (iky == ny_ / 2)
hiky = czero;

// elevation
fft_h_(ikx, iky) = h;
fft_h_ikx_(ikx, iky) = hikx;
Expand All @@ -358,6 +366,15 @@ namespace waves
complex hkyky = hok * ky2;
complex hkxky = hok * kx * ky;

if (ikx == nx_ / 2)
{
dx = czero;
}
if (iky == ny_ / 2)
{
dy = czero;
}

fft_sx_(ikx, iky) = dx;
fft_sy_(ikx, iky) = dy;
fft_h_kxkx_(ikx, iky) = hkxkx;
Expand Down
25 changes: 23 additions & 2 deletions gz-waves/src/WaveSimulationFFTRef.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -432,13 +432,21 @@ namespace waves
// height derivatives
complex hikx = hi * kx;
complex hiky = hi * ky;

// Nyquist terms for derivatives must be zero.
// For an explanation see:
// https://math.mit.edu/~stevenj/fft-deriv.pdf
if (ikx == nx_ / 2)
hikx = czero;
if (iky == ny_ / 2)
hiky = czero;

fft_h_ikx_(ikx, iky) = hikx;
fft_h_iky_(ikx, iky) = hiky;

// displacement and derivatives
if (std::abs(k) < 1.0E-8)
{
{
fft_sx_(ikx, iky) = czero;
fft_sy_(ikx, iky) = czero;
fft_h_kxkx_(ikx, iky) = czero;
Expand All @@ -452,7 +460,20 @@ namespace waves
complex hkxkx = hok * kx2;
complex hkyky = hok * ky2;
complex hkxky = hok * kx * ky;


if (ikx == nx_ / 2)
{
dx = czero;
hkxkx = complex(hkxkx.real(), 0.0);
hkxky = czero;
}
if (iky == ny_ / 2)
{
dy = czero;
hkyky = complex(hkyky.real(), 0.0);
hkxky = czero;
}

fft_sx_(ikx, iky) = dx;
fft_sy_(ikx, iky) = dy;
fft_h_kxkx_(ikx, iky) = hkxkx;
Expand Down
177 changes: 169 additions & 8 deletions gz-waves/src/WaveSimulationFFT_TEST.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -131,7 +131,6 @@ TEST_F(WaveSimulationFFTFixture, HermitianHTimeZeroReference)
}

//////////////////////////////////////////////////
#if 0
TEST_F(WaveSimulationFFTFixture, HermitianDhDxTimeZeroReference)
{
WaveSimulationFFTRefImpl model(lx_, ly_, nx_, ny_);
Expand Down Expand Up @@ -163,9 +162,8 @@ TEST_F(WaveSimulationFFTFixture, HermitianDhDxTimeZeroReference)
}
}
}
#endif

//////////////////////////////////////////////////
#if 0
TEST_F(WaveSimulationFFTFixture, HermitianDhDyTimeZeroReference)
{
WaveSimulationFFTRefImpl model(lx_, ly_, nx_, ny_);
Expand Down Expand Up @@ -197,7 +195,172 @@ TEST_F(WaveSimulationFFTFixture, HermitianDhDyTimeZeroReference)
}
}
}
#endif

//////////////////////////////////////////////////
TEST_F(WaveSimulationFFTFixture, HermitianSxTimeZeroReference)
{
WaveSimulationFFTRefImpl model(lx_, ly_, nx_, ny_);
model.ComputeBaseAmplitudes();
model.ComputeCurrentAmplitudes(0.0);

for (int ikx=0; ikx<nx_; ++ikx)
{
for (int iky=0; iky<ny_; ++iky)
{
// index for conjugate
int ckx = 0;
if (ikx != 0)
ckx = nx_ - ikx;

int cky = 0;
if (iky != 0)
cky = ny_ - iky;

// look up amplitude and conjugate
complex h = model.fft_sx_(ikx, iky);
complex hc = model.fft_sx_(ckx, cky);

// real part symmetric
EXPECT_DOUBLE_EQ(h.real(), hc.real());

// imaginary part anti-symmetric
EXPECT_DOUBLE_EQ(h.imag(), -1.0 * hc.imag());
}
}
}

//////////////////////////////////////////////////
TEST_F(WaveSimulationFFTFixture, HermitianSyTimeZeroReference)
{
WaveSimulationFFTRefImpl model(lx_, ly_, nx_, ny_);
model.ComputeBaseAmplitudes();
model.ComputeCurrentAmplitudes(0.0);

for (int ikx=0; ikx<nx_; ++ikx)
{
for (int iky=0; iky<ny_; ++iky)
{
// index for conjugate
int ckx = 0;
if (ikx != 0)
ckx = nx_ - ikx;

int cky = 0;
if (iky != 0)
cky = ny_ - iky;

// look up amplitude and conjugate
complex h = model.fft_sy_(ikx, iky);
complex hc = model.fft_sy_(ckx, cky);

// real part symmetric
EXPECT_DOUBLE_EQ(h.real(), hc.real());

// imaginary part anti-symmetric
EXPECT_DOUBLE_EQ(h.imag(), -1.0 * hc.imag());
}
}
}

//////////////////////////////////////////////////
TEST_F(WaveSimulationFFTFixture, HermitianDsxDxTimeZeroReference)
{
WaveSimulationFFTRefImpl model(lx_, ly_, nx_, ny_);
model.ComputeBaseAmplitudes();
model.ComputeCurrentAmplitudes(0.0);

for (int ikx=0; ikx<nx_; ++ikx)
{
for (int iky=0; iky<ny_; ++iky)
{
// index for conjugate
int ckx = 0;
if (ikx != 0)
ckx = nx_ - ikx;

int cky = 0;
if (iky != 0)
cky = ny_ - iky;

// look up amplitude and conjugate
complex h = model.fft_h_kxkx_(ikx, iky);
complex hc = model.fft_h_kxkx_(ckx, cky);

// real part symmetric
EXPECT_DOUBLE_EQ(h.real(), hc.real());

// imaginary part anti-symmetric
EXPECT_DOUBLE_EQ(h.imag(), -1.0 * hc.imag());
}
}
}

//////////////////////////////////////////////////
TEST_F(WaveSimulationFFTFixture, HermitianDsyDyTimeZeroReference)
{
WaveSimulationFFTRefImpl model(lx_, ly_, nx_, ny_);
model.ComputeBaseAmplitudes();
model.ComputeCurrentAmplitudes(0.0);

for (int ikx=0; ikx<nx_; ++ikx)
{
for (int iky=0; iky<ny_; ++iky)
{
// index for conjugate
int ckx = 0;
if (ikx != 0)
ckx = nx_ - ikx;

int cky = 0;
if (iky != 0)
cky = ny_ - iky;

// look up amplitude and conjugate
complex h = model.fft_h_kyky_(ikx, iky);
complex hc = model.fft_h_kyky_(ckx, cky);

// real part symmetric
EXPECT_DOUBLE_EQ(h.real(), hc.real());

// imaginary part anti-symmetric
EXPECT_DOUBLE_EQ(h.imag(), -1.0 * hc.imag());
}
}
}

//////////////////////////////////////////////////
TEST_F(WaveSimulationFFTFixture, HermitianDsxDyTimeZeroReference)
{
WaveSimulationFFTRefImpl model(lx_, ly_, nx_, ny_);
model.ComputeBaseAmplitudes();
model.ComputeCurrentAmplitudes(0.0);

for (int ikx=0; ikx<nx_; ++ikx)
{
for (int iky=0; iky<ny_; ++iky)
{
// index for conjugate
int ckx = 0;
if (ikx != 0)
ckx = nx_ - ikx;

int cky = 0;
if (iky != 0)
cky = ny_ - iky;

// look up amplitude and conjugate
complex h = model.fft_h_kxky_(ikx, iky);
complex hc = model.fft_h_kxky_(ckx, cky);

// real part symmetric
EXPECT_DOUBLE_EQ(h.real(), hc.real());

// imaginary part anti-symmetric
EXPECT_DOUBLE_EQ(h.imag(), -1.0 * hc.imag());
}
}
}

//////////////////////////////////////////////////
TEST_F(WaveSimulationFFTFixture, HermitianTimeNonZeroReference)
{
Expand Down Expand Up @@ -544,7 +707,6 @@ TEST_F(WaveSimulationFFTFixture, ElevationTimeNonZero)
}

//////////////////////////////////////////////////
#if 0
TEST_F(WaveSimulationFFTFixture, Displacement)
{
int n2 = nx_ * ny_;
Expand Down Expand Up @@ -572,11 +734,10 @@ TEST_F(WaveSimulationFFTFixture, Displacement)

for (int i=0; i<n2; ++i)
{
EXPECT_DOUBLE_EQ(sx(i, 0), ref_sx(i, 0));
EXPECT_DOUBLE_EQ(sy(i, 0), ref_sy(i, 0));
EXPECT_NEAR(sx(i, 0), ref_sx(i, 0), 1.0E-15);
EXPECT_NEAR(sy(i, 0), ref_sy(i, 0), 1.0E-15);
}
}
#endif

//////////////////////////////////////////////////
TEST_F(WaveSimulationFFTFixture, ElevationDerivatives)
Expand Down

0 comments on commit 9ee03da

Please sign in to comment.