diff --git a/hysop/numerics/tests/test_odesolvers.py b/hysop/numerics/tests/test_odesolvers.py
index b75b58418116b2d85df91718fbd0260984050f48..48efed05c4cc525b9c63f60b8fbc03271adcc834 100755
--- a/hysop/numerics/tests/test_odesolvers.py
+++ b/hysop/numerics/tests/test_odesolvers.py
@@ -87,8 +87,8 @@ def integrate_3d(integ, nb_steps):
npw.zeros(local_res_3d)]
# work = None
i = 1
- xref = local_res_3d[0] // 3.
- yref = local_res_3d[1] // 2.4
+ xref = local_res_3d[0] // 3
+ yref = local_res_3d[1] // 2
zref = local_res_3d[2] // 5
ref = npw.zeros((nb_steps, 4))
ref[0, 1] = y[0][xref, yref, zref]
@@ -112,30 +112,43 @@ def integrate_3d(integ, nb_steps):
def run_integ(integrator, order):
+ """apply 'integrator' (== Euler, RK3 ...)
+ and check error according to input order. 1D case
+ """
nb_steps = 100
wk_prop = integrator.get_work_properties(1, topo)['rwork']
work = []
for shape in wk_prop:
work.append(npw.zeros(shape))
- dt, err = integrate(integrator(1, topo, f=rhs, rwork=work), nb_steps)
+ dt, err = integrate(integrator(nb_components=1, topo=topo,
+ f=rhs, rwork=work), nb_steps)
assert err < dt ** order
def run_integ_3d(integrator, order):
+ """apply 'integrator' (== Euler, RK3 ...)
+ and check error according to input order. 3D case.
+
+ """
nb_steps = 50
wk_prop = integrator.get_work_properties(3, topo3d)['rwork']
work = []
for shape in wk_prop:
work.append(npw.zeros(shape))
- dt, err = integrate_3d(integrator(3, topo3d, f=rhs3d, rwork=work),
+ dt, err = integrate_3d(integrator(nb_components=3, topo=topo3d,
+ f=rhs3d, rwork=work),
nb_steps)
for e in err:
assert e < dt ** order
def run_integ_3d_no_work(integrator, order):
+ """apply 'integrator' (== Euler, RK3 ...)
+ and check error according to input order. 3D case.
+ """
nb_steps = 50
- dt, err = integrate_3d(integrator(3, topo3d, f=rhs3d), nb_steps)
+ dt, err = integrate_3d(integrator(nb_components=3, topo=topo3d,
+ f=rhs3d), nb_steps)
for e in err:
assert e < dt ** order
@@ -177,6 +190,7 @@ def test_rk4_3d():
def test_user_defined_common_workspace():
+ """One work space for all integrators"""
integrators = {Euler: 1, RK2: 2, RK3: 3, RK4: 4}
nb_steps = 50
wk_prop = {}
@@ -184,7 +198,7 @@ def test_user_defined_common_workspace():
wk_prop[integ] = integ.get_work_properties(3, topo3d)['rwork']
work = WorkSpaceTools.allocate_common_workspaces(wk_prop.values())
for integ in integrators:
- dt, err = integrate_3d(integ(3, topo3d, rhs3d, work), nb_steps)
+ dt, err = integrate_3d(integ(nb_components=3, topo=topo3d,
+ f=rhs3d, rwork=work), nb_steps)
for e in err:
assert e < dt ** integrators[integ]
-
diff --git a/hysop/operator/discrete/poisson_fft.py b/hysop/operator/discrete/poisson_fft.py
index 3027058fd38e52f41f047fc9b46a61bdfe3ac1d5..c26dfb106163d04e8fba4ea1770f0f5f48abaac7 100644
--- a/hysop/operator/discrete/poisson_fft.py
+++ b/hysop/operator/discrete/poisson_fft.py
@@ -13,7 +13,6 @@ if __FFTW_ENABLED__:
from hysop.f2hysop import fftw2py
-
class PoissonFFT(DiscreteOperator):
"""
Discretized Poisson operator based on FFTW.
diff --git a/hysop/operator/discrete/spectrum.py b/hysop/operator/discrete/spectrum.py
index defda637093bb1f59baf41de513731f898a10fb9..2bc03cbd5743b17753b18d8c1d3717523564baf7 100755
--- a/hysop/operator/discrete/spectrum.py
+++ b/hysop/operator/discrete/spectrum.py
@@ -1,12 +1,6 @@
# -*- coding: utf-8 -*-
+"""Discrete Spectrum operator using FFTW
"""
-@file spectrum.py
-Discrete Spectrum operator using FFTW (fortran)
-"""
-try:
- from hysop.f2hysop import fftw2py
-except ImportError:
- from hysop.fakef2py import fftw2py
from hysop.operator.discrete.discrete import DiscreteOperator
from hysop.constants import debug, np, HYSOP_MPI_REAL
from hysop.tools.profiler import profile
@@ -14,22 +8,28 @@ from hysop.tools.io_utils import IO
import hysop.tools.numpywrappers as npw
from hysop.mpi import MPI
import os
+from hysop import __FFTW_ENABLED__
+if __FFTW_ENABLED__:
+ from hysop.f2hysop import fftw2py
-class FFTSpectrum(DiscreteOperator):
- """
- Discretized Spectrum operator based on FFTW.
+class FFTSpectrum(DiscreteOperator):
+ """Discretized Spectrum operator based on FFTW.
"""
@debug
- def __init__(self, field, prefix=None, **kwds):
+ def __init__(self, field, **kwds):
"""
- Constructor.
- @param[in] vorticity : field to compute.
+
+ Parameters
+ ----------
+ field: :class:`~hysop.fields.discrete.DiscreteField`
+ the input field for which spectrum will be computed
"""
# Discretisation of the input field
self.field = field
-
+ msg = 'Spectrum error: implemented only for 3D problems.'
+ assert self.field.dimension == 3, msg
if self.field.nb_components > 1:
raise AttributeError("Vector case not yet implemented.")
# Base class initialisation
@@ -41,44 +41,32 @@ class FFTSpectrum(DiscreteOperator):
self._tmp = npw.zeros(((l - 1) / 2, ))
self._kx = npw.zeros(((l - 1) / 2, ))
self.res = npw.zeros(((l - 1) / 2, ))
- self._prefix = "spectrum" if prefix is None else prefix
- IO.check_dir(os.path.join(IO.default_path(),
- self._prefix + "_00000.dat"))
@debug
@profile
def apply(self, simulation=None):
assert simulation is not None, \
- "Missing dt value for diffusion computation."
+ "Missing simulation parameter."
ite = simulation.current_iteration
ghosts = self.field.topology.ghosts()
- if self.field.dimension == 3:
- fftw2py.spectrum_3d(self.field.data[0],
- self._tmp, self._kx,
- ghosts, np.min(self.domain.length))
- if self.field.topology.size == 1:
- self.res[...] = self._tmp
- else:
- self.field.topology.comm.Reduce(
- [self._tmp, self.res.shape[0], HYSOP_MPI_REAL],
- [self.res, self.res.shape[0], HYSOP_MPI_REAL],
- op=MPI.SUM, root=0)
- if self.field.topology.rank == 0:
- _file = open(os.path.join(
- IO.default_path(),
- self._prefix + "_{0:05d}.dat".format(ite)), 'w')
- for i in xrange(self.res.shape[0]):
- _file.write("{0} {1}\n".format(
- self._kx[i], self.res[i]))
- _file.close()
+ fftw2py.spectrum_3d(self.field.data[0],
+ self._tmp, self._kx,
+ ghosts, np.min(self.domain.length))
+ if self.field.topology.size == 1:
+ self.res[...] = self._tmp
else:
- raise ValueError("invalid problem dimension")
+ self.field.topology.comm.Reduce(
+ [self._tmp, self.res.shape[0], HYSOP_MPI_REAL],
+ [self.res, self.res.shape[0], HYSOP_MPI_REAL],
+ op=MPI.SUM, root=0)
+
+ if self._writer is not None and self._writer.do_write(ite):
+ nbc = self.res.shape[0]
+ self._writer.buffer[0, 0:nbc] = self._kx
+ self._writer.buffer[0, nbc:] = self.res
+ self._writer.write()
def finalize(self):
- """
- Clean memory (fftw plans and so on)
+ """Clean memory (fftw plans and so on)
"""
pass
- # TODO : fix bug that occurs when several finalize
- # of fft operators are called.
- # fftw2py.clean_fftw_solver(self.vorticity.dimension)
diff --git a/hysop/operator/spectrum.py b/hysop/operator/spectrum.py
index 1b4f5c1b05df6221985040e76850cb880dc845a6..9d244a5590ec7312b914eb1d9d707bf94797a66e 100644
--- a/hysop/operator/spectrum.py
+++ b/hysop/operator/spectrum.py
@@ -1,5 +1,4 @@
-"""
-@file spectrum.py
+"""Compute the spectrum of a field (fftw based)
"""
from hysop.operator.computational import Computational
from hysop.operator.discrete.spectrum import FFTSpectrum
@@ -12,15 +11,17 @@ class Spectrum(Computational):
Fourier spectrum computation of a scalar field.
"""
- def __init__(self, field, prefix=None, **kwds):
+ def __init__(self, field, **kwds):
"""
- Constructor for Spectrum operator
- @param[in] field : field to compute
+
+ Parameters
+ ----------
+ field: :class:`~hysop.fields.continuous.Field`
+ the input field for which spectrum will be computed
"""
super(Spectrum, self).__init__(variables=[field], **kwds)
self.field = field
self.input = [field]
- self._prefix = prefix
def discretize(self):
super(Spectrum, self)._fftw_discretize()
@@ -29,6 +30,11 @@ class Spectrum(Computational):
@opsetup
def setup(self, rwork=None, iwork=None):
self.discrete_op = FFTSpectrum(self.discreteFields[self.field],
- method=self.method,
- prefix=self._prefix)
+ method=self.method)
+ nbc = self.discrete_op.res.size
+ self._set_io('spectrum', (1, 2 * nbc))
+ self.discrete_op.set_writer(self._writer)
self._is_uptodate = True
+
+ def get_work_properties(self):
+ return {'rwork': None, 'iwork': None}
diff --git a/hysop/operator/tests/test_differential.py b/hysop/operator/tests/test_differential.py
index 2a9c983dcea199c5355c746ba0cd6ace9d16d279..f5520a4dd46cd405d9719a902f8936d22d40d702 100755
--- a/hysop/operator/tests/test_differential.py
+++ b/hysop/operator/tests/test_differential.py
@@ -134,9 +134,10 @@ def call_op(class_name, ref_formula, topo, use_work=False,
if use_work:
work_prop = op.get_work_properties()['rwork']
work = []
- for l in xrange(len(work_prop)):
- shape = work_prop[l]
- work.append(npw.zeros(shape))
+ if work_prop is not None:
+ for l in xrange(len(work_prop)):
+ shape = work_prop[l]
+ work.append(npw.zeros(shape))
op.setup(rwork=work)
check(op, ref_formula, topo, op_dim, order)
diff --git a/hysop/operator/tests/test_spectrum.py b/hysop/operator/tests/test_spectrum.py
index 7152e3f65633617b9937967645b80468e6ea26ed..e2404022e8f9f0f7117e1ef3c0cbdfa4492a63d4 100755
--- a/hysop/operator/tests/test_spectrum.py
+++ b/hysop/operator/tests/test_spectrum.py
@@ -1,23 +1,28 @@
+"""Test fftw spectrum computation
+"""
from hysop.domain.box import Box
from hysop.fields.continuous import Field
from hysop.operator.spectrum import Spectrum
from hysop.tools.parameters import Discretization
from hysop.problem.simulation import Simulation
-from hysop.operator.hdf_io import HDF_Writer
+from hysop import testsenv
import numpy as np
-from hysop.tools.io_utils import IO
pi = np.pi
sin = np.sin
cos = np.cos
def computeScal(res, x, y, z, t):
- res[0][...] = z * sin((3*pi*x)*(2*pi*y))
+ """test function, scalar input, 3D"""
+ res[0][...] = z * sin((3 * pi * x) * (2 * pi * y))
return res
+@testsenv.fftw_failed
def test_spectrum():
- dom = Box(origin=[0., 0., 0.], length=[1., 1., 1.])
+ """build/apply spectrum op
+ """
+ dom = Box()
field = Field(domain=dom, name='Field',
is_vector=False, formula=computeScal)
d3D = Discretization([257, 257, 257])