diff --git a/examples/bubble/periodic_jet_levelset.py b/examples/bubble/periodic_jet_levelset.py
index bbb1cfe1dba86b8c9d188148cd03045fc2d5061f..1f9cdbf530a17c49cdc94dbd419d81db0e82e8e0 100644
--- a/examples/bubble/periodic_jet_levelset.py
+++ b/examples/bubble/periodic_jet_levelset.py
@@ -204,7 +204,7 @@ def compute(args):
implementation=impl, **extra_op_kwds)
### Adaptive timestep operator
- adapt_dt = AdaptiveTimeStep(dt, equivalent_CFL=True, max_dt=1e-2,
+ adapt_dt = AdaptiveTimeStep(dt, equivalent_CFL=True, max_dt=1e-4,
name='merge_dt', pretty_name='dt', )
dt_cfl = adapt_dt.push_cfl_criteria(cfl=args.cfl, Finf=min_max_U.Finf,
equivalent_CFL=True,
diff --git a/hysop/numerics/fft/fft.py b/hysop/numerics/fft/fft.py
index 798d37af6fe671be6c7e27195668756bb069cb8a..8e0fa1762a1dbf63bbde00db8b306e8a9ab83dc2 100644
--- a/hysop/numerics/fft/fft.py
+++ b/hysop/numerics/fft/fft.py
@@ -181,7 +181,7 @@ class FFTI(object):
@abstractmethod
- def irfft(self, a, out, axis=-1, build_plan=False, **kwds):
+ def irfft(self, a, out, n=None, axis=-1, build_plan=False, **kwds):
"""
Compute the one-dimensional hermitian complex to real discrete Fourier Transform scaled by 1/N.
@@ -195,24 +195,50 @@ class FFTI(object):
Last axis should match forward transform size used:
1) out.shape[axis] = 2*(a.shape[axis]-1)
2) out.shape[axis] = 2*(a.shape[axis]-1) + 1
+ n: int, optional
+ Length of the transformed axis of the output.
+ ie: n should be in [2*(a.shape[axis]-1), 2*(a.shape[axis]-1)+1]
axis: int, optional
Axis over witch to compute the transform.
Defaults to last axis.
build_plan: bool, optional
Build and return a IRFFT plan instead of applying the transform.
+ Notes
+ -----
+ To get an odd number of output points, n or out must be specified.
+
Returns
-------
- (shape0, shape1, dtype) of the output array determined from the input array.
+ (shape, dtype) of the output array determined from the input array, out and n.
"""
assert a.dtype in self.supported_ctypes
cshape = a.shape
rtype = complex_to_float_dtype(a.dtype)
- rshape = list(a.shape)
- rshape[axis] = 2*(cshape[axis]-1)
+
+ rshape_even, rshape_odd = list(a.shape), list(a.shape)
+ rshape_even[axis] = 2*(cshape[axis]-1)
+ rshape_odd[axis] = 2*(cshape[axis]-1) + 1
+
if (out is not None):
assert out.dtype in self.supported_ftypes
assert rtype == out.dtype
- assert np.arrayequal(out.shape, rshape)
+
+ ns = out.shape[axis]
+ if (n is not None):
+ assert ns==n, 'output shape mismatch with specified transformed output size.'
+ else:
+ n = ns
+
+ if (n%2==0):
+ assert np.arrayequal(out.shape, rshape_even)
+ else:
+ assert np.arrayequal(out.shape, rshape_odd)
+
+ if (n is None) or (n%2==0):
+ rshape = rshape_even
+ else:
+ rshape = rshape_odd
+
return (rshape, rtype)
diff --git a/hysop/numerics/fft/fftw_fft.py b/hysop/numerics/fft/fftw_fft.py
index 6fab55ceba4c0fa972fddb4992ce101dbe07d8bc..19eb3639162ae71240498314fdd1f931449e51d0 100644
--- a/hysop/numerics/fft/fftw_fft.py
+++ b/hysop/numerics/fft/fftw_fft.py
@@ -177,9 +177,10 @@ class FftwFFT(FFTI):
plan = FftwFFTPlan(**kwds)
return plan
- def irfft(self, a, out=None, axis=-1, **kwds):
+ def irfft(self, a, out=None, n=None, axis=-1, **kwds):
"""Planning destroys initial arrays content."""
- (shape, dtype) = super(FftwFFT, self).irfft(a=a, out=out, axis=axis, **kwds)
+ (shape, dtype) = super(FftwFFT, self).irfft(a=a, out=out, axis=axis,
+ n=n, **kwds)
out = self.allocate_output(out, shape, dtype)
if self.warn_on_misalignment:
self.check_alignment(a, out)
diff --git a/hysop/numerics/fft/numpy_fft.py b/hysop/numerics/fft/numpy_fft.py
index 6119e11a21b254e4e0e5eb95e34c1ddbd552f599..01918220a2213ef795817516fcc15405bd0c20c9 100644
--- a/hysop/numerics/fft/numpy_fft.py
+++ b/hysop/numerics/fft/numpy_fft.py
@@ -41,8 +41,8 @@ class NumpyFFTPlan(FFTPlanI):
def __call__(self, out=None, **kwds):
fn_kwds = self.kwds.copy()
fn_kwds.update(kwds)
- if (self.fn is _FFT.irfft) and (out is not None):
- fn_kwds.setdefault('n', out.shape[fn_kwds['axis']])
+ if (self.fn is _FFT.irfft):
+ assert 'n' in fn_kwds
out = first_not_None(out, self.out)
res = self.fn(**fn_kwds)
@@ -92,9 +92,9 @@ class NumpyFFT(FFTI):
**kwds)
return plan
- def irfft(self, a, out=None, axis=-1, **kwds):
- super(NumpyFFT, self).irfft(a=a, out=out, axis=axis, **kwds)
+ def irfft(self, a, out=None, n=None, axis=-1, **kwds):
+ (rshape, _) = super(NumpyFFT, self).irfft(a=a, out=out, n=n, axis=axis, **kwds)
plan = NumpyFFTPlan(fn=_FFT.irfft, a=a, out=out, axis=axis,
output_dtype=lambda x: complex_to_float_dtype(x),
- **kwds)
+ n=rshape[axis], **kwds)
return plan
diff --git a/hysop/numerics/fft/scipy_fft.py b/hysop/numerics/fft/scipy_fft.py
new file mode 100644
index 0000000000000000000000000000000000000000..2c1e977a2b49408aac04eb5bc404bc82451dae94
--- /dev/null
+++ b/hysop/numerics/fft/scipy_fft.py
@@ -0,0 +1,121 @@
+"""
+FFT iterface for fast Fourier Transforms using scipy fftpack backend.
+:class:`~hysop.numerics.ScipyFFT`
+:class:`~hysop.numerics.ScipyFFTPlan`
+"""
+
+import numpy as np
+import scipy as sp
+from scipy import fftpack as _FFT
+
+from hysop.tools.types import first_not_None
+from hysop.numerics.fft.fft import FFTPlanI, FFTI, complex_to_float_dtype, float_to_complex_dtype
+
+
+class ScipyFFTPlan(FFTPlanI):
+ """
+ Wrap a scipy fftpack call (scipy.fftpack does not offer real planning capabilities).
+ """
+
+ def __init__(self, fn, out,
+ output_dtype=lambda dtype: dtype,
+ **kwds):
+ super(ScipyFFTPlan, self).__init__()
+ assert callable(output_dtype)
+ self.fn = fn
+ self.out = out
+ self.kwds = kwds.copy()
+
+ def execute(self):
+ return self.__call__()
+
+ def __call__(self, out=None, **kwds):
+ fn = self.fn
+ fn_kwds = self.kwds.copy()
+ fn_kwds.update(kwds)
+ if (fn is _FFT.irfft):
+ assert 'n' in fn_kwds
+
+ in_ = fn_kwds['x']
+
+ if fn is _FFT.irfft:
+ axis = fn_kwds['axis']
+ assert axis in (in_.ndim-1, -1)
+ in_ = in_.view(dtype=complex_to_float_dtype(in_.dtype))
+ is_even = (fn_kwds['n']%2==0)
+ rshape = list(in_.shape)
+ rshape[-1] -= 1 + is_even
+ _in = np.empty(shape=rshape, dtype=in_.dtype)
+ _in[...,1:] = in_[...,2:in_.shape[-1]-is_even]
+ _in[...,0] = in_[...,0]
+ fn_kwds['x'] = _in
+
+ out = first_not_None(out, self.out)
+ res = fn(**fn_kwds)
+
+ if fn is _FFT.rfft:
+ axis = fn_kwds['axis']
+ assert axis in (in_.ndim-1, -1)
+ is_even = (in_.shape[axis] % 2 == 0)
+ if (out is None):
+ rshape = list(res.shape)
+ rshape[axis] += 1 + is_even
+ out = np.empty(dtype=in_.dtype, shape=rshape)
+ else:
+ rtype = complex_to_float_dtype(out.dtype)
+ out = out.view(dtype=rtype)
+ ctype = float_to_complex_dtype(out.dtype)
+ out[..., 1:out.shape[-1]-is_even] = res
+ out[..., 0] = out[..., 1]
+ out[..., 1] = 0.0
+ if is_even:
+ out[..., -1] = 0.0
+ out = out.view(dtype=ctype)
+ elif (out is not None):
+ out[...] = res
+ else:
+ out = res
+ return out
+
+
+class ScipyFFT(FFTI):
+ """
+ Interface to compute local to process FFT-like transforms using the scipy fftpack backend.
+
+ Scipy fftpack backend has some disadvantages:
+ - creation of one or two intermediate temporary buffers at each call
+ - hermitian-complex data layout is different than all other fft implementations
+ - no planning capabilities (scipy.fftpack methods are just wrapped into fake plans)
+
+ It has native float and double precision support unlike the numpy fft backend.
+ Planning won't destroy original inputs.
+ """
+
+ def __init__(self):
+ super(ScipyFFT, self).__init__()
+ self.supported_ftypes = (np.float32, np.float64,)
+ self.supported_ctypes = (np.complex64, np.complex128,)
+
+ def fft(self, a, out=None, axis=-1, **kwds):
+ super(ScipyFFT, self).fft(a=a, out=out, axis=axis, **kwds)
+ plan = ScipyFFTPlan(fn=_FFT.fft, x=a, out=out, axis=axis, **kwds)
+ return plan
+
+ def ifft(self, a, out=None, axis=-1, **kwds):
+ super(ScipyFFT, self).ifft(a=a, out=out, axis=axis, **kwds)
+ plan = ScipyFFTPlan(fn=_FFT.ifft, x=a, out=out, axis=axis, **kwds)
+ return plan
+
+ def rfft(self, a, out=None, axis=-1, **kwds):
+ super(ScipyFFT, self).rfft(a=a, out=out, axis=axis, **kwds)
+ plan = ScipyFFTPlan(fn=_FFT.rfft, x=a, out=out, axis=axis,
+ output_dtype=lambda x: float_to_complex_dtype(x),
+ **kwds)
+ return plan
+
+ def irfft(self, a, out=None, n=None, axis=-1, **kwds):
+ (rshape, _) = super(ScipyFFT, self).irfft(a=a, out=out, n=n, axis=axis, **kwds)
+ plan = ScipyFFTPlan(fn=_FFT.irfft, x=a, out=out, axis=axis,
+ output_dtype=lambda x: complex_to_float_dtype(x),
+ n=rshape[axis], **kwds)
+ return plan
diff --git a/hysop/numerics/tests/test_fft.py b/hysop/numerics/tests/test_fft.py
index 34088a72ed24b1bbdb744950c5a14d8446e06422..5f14b009a49d0178444e07e590bdbc2b181ef45b 100644
--- a/hysop/numerics/tests/test_fft.py
+++ b/hysop/numerics/tests/test_fft.py
@@ -16,6 +16,7 @@ from hysop.tools.numerics import float_to_complex_dtype
from hysop.tools.types import check_instance, first_not_None
from hysop.numerics.fft.numpy_fft import NumpyFFT
+from hysop.numerics.fft.scipy_fft import ScipyFFT
from hysop.numerics.fft.fftw_fft import FftwFFT
@@ -23,6 +24,7 @@ class TestFFT(object):
implementations = {
'numpy': NumpyFFT(),
+ 'scipy': ScipyFFT(),
'fftw': FftwFFT(warn_on_allocation=False)
}
@@ -31,9 +33,10 @@ class TestFFT(object):
msg_input_modified = 'Input array has been modified for implementation {}.'
msg_output_modified = 'Output array results are not consistent for implementation {}.'
- max_eps = 2
+ report_eps = 10
+ fail_eps = 1000
- def _test_1d(self, dtype):
+ def _test_1d(self, dtype, failures):
print
print '::Testing 1D transform, precision {}::'.format(dtype.__name__)
eps = np.finfo(dtype).eps
@@ -41,8 +44,10 @@ class TestFFT(object):
def iter_shapes():
msg = ('EVEN', 'ODD')
+ minj=(5,1)
maxj=(15,10)
- minj=(10,5)
+ #minj=(2,2)
+ #maxj=(5,5)
for i in xrange(2):
base = 2+i
print ' '+msg[i]
@@ -51,12 +56,16 @@ class TestFFT(object):
cshape = list(shape)
cshape[-1] = cshape[-1]//2 + 1
cshape = tuple(cshape)
+ rshape = list(shape)
+ rshape[-1] = (rshape[-1]//2) * 2
+ rshape = tuple(rshape)
N = np.prod(shape, dtype=np.int64)
Nc = np.prod(cshape, dtype=np.int64)
- yield (shape, cshape, N, Nc)
+ Nr = np.prod(rshape, dtype=np.int64)
+ yield (shape, cshape, rshape, N, Nc, Nr)
print ' FORWARD C2C: complex to complex forward transform'
- for (shape, cshape, N, Nc) in iter_shapes():
+ for (shape, cshape, rshape, N, Nc, Nr) in iter_shapes():
print ' FFT shape={:9s} '.format(str(shape)+':'),
Href = np.random.rand(2*N).astype(dtype).view(dtype=ctype).reshape(shape)
H = np.empty_like(Href)
@@ -74,11 +83,11 @@ class TestFFT(object):
out[...] = 0
out = plan.execute()
assert np.array_equal(out, refout), self.msg_output_modified.format(name)
- results[name] = refout
- self.check_distances(results, eps, N*self.max_eps)
+ results[name] = refout / N # forward normalization
+ self.check_distances(results, eps, self.report_eps, 'forward C2C', failures)
print ' BACKWARD C2C: complex to complex backward transform'
- for (shape, cshape, N, Nc) in iter_shapes():
+ for (shape, cshape, rshape, N, Nc, Nr) in iter_shapes():
print ' IFFT shape={:9s} '.format(str(shape)+':'),
Href = np.random.rand(2*N).astype(dtype).view(dtype=ctype).reshape(shape)
H = np.empty_like(Href)
@@ -97,10 +106,10 @@ class TestFFT(object):
out = plan.execute()
assert np.array_equal(out, refout), self.msg_output_modified.format(name)
results[name] = refout
- self.check_distances(results, eps, self.max_eps)
+ self.check_distances(results, eps, self.report_eps, 'backward C2C', failures)
print ' FORWARD R2C: real to hermitian complex transform'
- for (shape, cshape, N, Nc) in iter_shapes():
+ for (shape, cshape, rshape, N, Nc, Nr) in iter_shapes():
print ' RFFT shape={:9s} '.format(str(shape)+':'),
Href = np.random.rand(N).astype(dtype).reshape(shape)
H = np.empty_like(Href)
@@ -118,12 +127,12 @@ class TestFFT(object):
out[...] = 0
out = plan.execute()
assert np.array_equal(out, refout), self.msg_output_modified.format(name)
- results[name] = refout
- self.check_distances(results, eps, N*self.max_eps)
+ results[name] = refout / N #forward normalization
+ self.check_distances(results, eps, self.report_eps, 'R2C', failures)
print ' BACKWARD C2R: real to hermitian complex transform'
- for (shape, cshape, N, Nc) in iter_shapes():
- print ' RFFT shape={:9s} '.format(str(shape)+':'),
+ for (shape, cshape, rshape, N, Nc, Nr) in iter_shapes():
+ print ' IRFFT shape={:9s} '.format(str(rshape)+':'),
Href = np.random.rand(2*Nc).astype(dtype).view(dtype=ctype).reshape(cshape)
H = np.empty_like(Href)
H, Href = self.simd_align(H, Href)
@@ -133,7 +142,29 @@ class TestFFT(object):
plan = impl.irfft(a=H)
H[...] = Href
out = plan.execute()
- #assert out.shape == shape, self.msg_shape.format(shape, out.shape, name)
+ assert out.shape == rshape, self.msg_shape.format(rshape, out.shape, name)
+ assert out.dtype == dtype, self.msg_dtype.format(dtype, out.dtype, name)
+ assert np.array_equal(Href, H), self.msg_input_modified.format(name)
+ refout = out.copy()
+ out[...] = 0
+ out = plan.execute()
+ assert np.array_equal(out, refout), self.msg_output_modified.format(name)
+ results[name] = refout
+ self.check_distances(results, eps, self.report_eps, 'normal C2R', failures)
+
+ print ' BACKWARD FORCED C2R: real to hermitian complex transform with specified shape'
+ for (shape, cshape, rshape, N, Nc, Nr) in iter_shapes():
+ print ' IRFFT shape={:9s} '.format(str(shape)+':'),
+ Href = np.random.rand(2*Nc).astype(dtype).view(dtype=ctype).reshape(cshape)
+ H = np.empty_like(Href)
+ H, Href = self.simd_align(H, Href)
+ results = {}
+ for (name, impl) in self.implementations.iteritems():
+ if dtype in impl.supported_ftypes:
+ plan = impl.irfft(a=H, n=shape[-1])
+ H[...] = Href
+ out = plan.execute()
+ assert out.shape == shape, self.msg_shape.format(shape, out.shape, name)
assert out.dtype == dtype, self.msg_dtype.format(dtype, out.dtype, name)
assert np.array_equal(Href, H), self.msg_input_modified.format(name)
refout = out.copy()
@@ -141,24 +172,31 @@ class TestFFT(object):
out = plan.execute()
assert np.array_equal(out, refout), self.msg_output_modified.format(name)
results[name] = refout
- self.check_distances(results, eps, self.max_eps)
+ self.check_distances(results, eps, self.report_eps, 'forced C2R', failures)
- def check_distances(self, results, eps, max_eps):
+ def check_distances(self, results, eps, report_eps, tag, failures):
if len(results.keys())<2:
impl = results.keys()[0]
print 'cannot compare'
return
ss=()
- failed = False
for (r0,r1) in it.combinations(results.keys(), 2):
+ if not (results[r0].shape == results[r1].shape):
+ print
+ msg='Output shapes do not match.'
+ raise RuntimeError(msg)
E = results[r1] - results[r0]
Einf = np.max(np.abs(E))
Eeps = int(np.round(Einf/eps))
- failed |= (Eeps>max_eps)
s='|{}-{}|={}eps'.format(r0,r1,Eeps)
ss += (s,)
+ failed = (Eeps>report_eps)
+ if failed:
+ shape=results[r0].shape
+ failures.setdefault(tag, []).append((r0, r1, shape, Einf, Eeps))
+
print ', '.join(ss)
- if failed:
+ if failed and False:
print
msg='Some implementations did not agree on the result !'
raise RuntimeError(msg)
@@ -168,14 +206,54 @@ class TestFFT(object):
for array in arrays:
aligned_arrays += (pyfftw.byte_align(array),)
return aligned_arrays
+
+ def report_failures(self, failures):
+ print
+ print '== TEST FAILURES REPORT =='
+ print ' Report error has been set to {} epsilons.'.format(self.report_eps)
+ print ' Test failure has been set to {} epsilons.'.format(self.fail_eps)
+ assert self.report_eps <= self.fail_eps
+ failed = False
+ cnt = 0
+ for (dtype, fails) in failures.iteritems():
+ if not fails:
+ continue
+ print '::{} precision tests'.format(dtype)
+ for (tag, tfails) in fails.iteritems():
+ print ' |{} transform errors:'.format(tag)
+ for (r0,r1, shape, Einf, Eeps) in sorted(tfails,
+ key=lambda x: -x[-2]):
+ failed |= (Eeps >= self.fail_eps)
+ cnt+=1
+ msg=' {} vs {}:\tshape {}\t->\tE={}\t({} eps)'
+ msg=msg.format(r0, r1, shape, Einf, Eeps)
+ print msg
+ print '==========================='
+
+ if failed:
+ msg =''
+ msg+='\n************* FFT TESTS FAILED ***************'
+ msg+='\n** One or more test exceeded failure error. **'
+ msg+='\n**********************************************'
+ msg+='\n'
+ print msg
+ raise RuntimeError
+ else:
+ msg =''
+ msg+='\n*************** FFT TESTS PASSED ******************'
+ msg+='\n** Some tests may have exceeded reporting error. **'
+ msg+='\n***************************************************'
+ msg+='\n'
+ print msg
- def test_1d_transforms(self):
- self._test_1d(dtype=np.float32)
- self._test_1d(dtype=np.float64)
- #self._test_1d(dtype=np.longdouble)
def perform_tests(self):
- self.test_1d_transforms()
+ # not testing np.longdouble because only one implementation support it
+ # ie. we cannot compare results
+ failures = {}
+ for dtype in (np.float32, np.float64,):
+ self._test_1d(dtype=dtype, failures=failures.setdefault(dtype.__name__, {}))
+ self.report_failures(failures)
if __name__ == '__main__':
test = TestFFT()