diff --git a/CMakeLists.txt b/CMakeLists.txt
index ccbee338cd6143184fad5cb0d3fb1cf21310c92b..77936e8f3bc4c8c7beaaab9567e8ab1daae9ba78 100755
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -546,7 +546,7 @@ if(WITH_TESTS)
if(NOT USE_MPI)
set(WITH_MPI_TESTS "OFF")
endif()
- include(HySoPTests)
+ include(hysop_tests)
endif(WITH_TESTS)
# ============= Documentation setup =============
diff --git a/cmake/HySoPTests.cmake b/cmake/hysop_tests.cmake
similarity index 97%
rename from cmake/HySoPTests.cmake
rename to cmake/hysop_tests.cmake
index bf9948b98835380afc4f2e231d62c8e5fe77aef1..09f38645bebed72e8ae2487eb7417d2c40a4e983 100755
--- a/cmake/HySoPTests.cmake
+++ b/cmake/hysop_tests.cmake
@@ -101,7 +101,7 @@ endforeach()
# Doctest are run for every line which contains '>>>'
set(py_doctest_files)
foreach(testdir ${py_src_dirs})
- file(GLOB testfiles hysop/${testdir}/[a-zA-Z]*.py)
+ file(GLOB testfiles RELATIVE ${CMAKE_SOURCE_DIR} hysop/${testdir}/[a-zA-Z]*.py)
foreach(testfile ${testfiles})
file(STRINGS ${testfile} test_doctest REGEX ">>>")
if(NOT "${test_doctest}" STREQUAL "")
@@ -110,6 +110,7 @@ foreach(testdir ${py_src_dirs})
endforeach()
endforeach()
+
# === Create tests from py_test_files ===
foreach(testfile ${py_test_files})
get_filename_component(testName ${testfile} NAME_WE)
@@ -121,6 +122,7 @@ endforeach()
# Add files containing doctests
foreach(testfile ${py_doctest_files})
get_filename_component(testName ${testfile} NAME_WE)
- add_test("doctest_${testName}" py.test -v --doctest-modules ${testfile})
+ set(exename ${testDir}/${testfile})
+ add_test("doctest_${testName}" py.test -v --doctest-modules ${exename})
endforeach()
diff --git a/docs/sphinx/users_guide/poisson.rst b/docs/sphinx/users_guide/poisson.rst
index 853c660647b86011f81ce9e79748ed15ae079dfd..edba19f18f365ddb1892755faeab3e3e1e12a7f9 100644
--- a/docs/sphinx/users_guide/poisson.rst
+++ b/docs/sphinx/users_guide/poisson.rst
@@ -10,7 +10,7 @@ Poisson Operator
The class :class:`Poisson` is used to model and solve the following problem :
.. math::
- \nabla \psi = -\omega \\
+ \Delta \psi = -\omega \\
v = \nabla \times \psi
with
diff --git a/hysop/domain/tests/test_box.py b/hysop/domain/tests/test_box.py
index 24a12f79a52c1a6222ea2fd973dce23fb8b16332..9665ae861dc081b3c4f9f3291118129bd56a801e 100755
--- a/hysop/domain/tests/test_box.py
+++ b/hysop/domain/tests/test_box.py
@@ -122,6 +122,8 @@ def test_topo_plane():
def test_topo_from_mesh():
+ """Build a topology from a previoulsy defined mesh.
+ """
# e.g. for fftw
dom = Box(proc_tasks=proc_tasks)
from hysop import __FFTW_ENABLED__
diff --git a/hysop/fields/tests/func_for_tests.py b/hysop/fields/tests/func_for_tests.py
index 3e2eaf562a1fdc3560b039510e0a6274851d4a2d..5dc679504773b72ebb1588e2ababe349358acd45 100755
--- a/hysop/fields/tests/func_for_tests.py
+++ b/hysop/fields/tests/func_for_tests.py
@@ -121,3 +121,17 @@ def w_TG(res, x, y, z, t):
res[1][...] = - sin(x) * cos(y) * sin(z)
res[2][...] = 2. * sin(x) * sin(y) * cos(z)
return res
+
+
+def v_TG_2d(res, x, y, t):
+ """Taylor-Green-like 2d velocity"""
+ res[0][...] = sin(x) * cos(y)
+ res[1][...] = - cos(x) * sin(y)
+ return res
+
+
+# Function to compute reference vorticity
+def w_TG_2d(res, x, y, t):
+ """Taylor-Green-like vorticity, 2d domain"""
+ res[0][...] = - cos(x) * sin(y)
+ return res
diff --git a/hysop/operator/computational.py b/hysop/operator/computational.py
index d02936b5a6c03fbcdf1bd21a212209a256e39a05..b2bb201b9f6400059f48f5bb0ed598054d74ab78 100755
--- a/hysop/operator/computational.py
+++ b/hysop/operator/computational.py
@@ -6,7 +6,6 @@ from hysop.constants import debug
from hysop.operator.continuous import Operator, opapply
from hysop.mpi.topology import Cartesian
from hysop.tools.parameters import Discretization
-from hysop.tools.profiler import profile
import hysop.tools.numpywrappers as npw
from hysop import __FFTW_ENABLED__
if __FFTW_ENABLED__:
diff --git a/hysop/operator/curlAndDiffusion.py b/hysop/operator/curlAndDiffusion.py
deleted file mode 100644
index 8b86d2bff462b80e59862a6aa7746417db7d271b..0000000000000000000000000000000000000000
--- a/hysop/operator/curlAndDiffusion.py
+++ /dev/null
@@ -1,79 +0,0 @@
-# -*- coding: utf-8 -*-
-"""
-@file diffusion.py
-
-Operator for diffusion problem.
-
-"""
-from hysop.operator.continuous import Operator
-try:
- from hysop.f2hysop import fftw2py
-except ImportError:
- from hysop.fakef2py import fftw2py
-from hysop.operator.discrete.diffusion_fft import DiffusionFFT
-from hysop.constants import debug
-from hysop.operator.continuous import opsetup
-
-
-class CurlDiffusion(Operator):
- """
- Diffusion operator
- \f{eqnarray*}
- \omega = Op(\omega)
- \f} with :
- \f{eqnarray*}
- \frac{\partial \omega}{\partial t} &=& \nu\Delta\omega
- \f}
- """
-
- @debug
- def __init__(self, velocity, vorticity, **kwds):
- """
- Constructor.
- Create a Diffusion operator using FFT.
-
- @param velocity ContinuousVectorField : velocity variable.
- @param vorticity ContinuousVectorField : vorticity variable.
- @param viscosity : viscosity of the considered medium.
- """
- super(CurlDiffusion, self).__init__(variables=[velocity, vorticity], **kwds)
- self.velocity = velocity
- self.vorticity = vorticity
- raise ValueError("This operator is obsolete and must be reviewed.\
- Do not use it.")
-
- @debug
- @opsetup
- def setup(self):
- """
- Diffusion operator discretization method.
- Create a discrete Diffusion operator from given specifications.
- """
- if self._comm is None:
- from hysop.mpi import main_comm as comm
- else:
- comm = self._comm
-
- localres, localoffset = fftw2py.init_fftw_solver(
- self.resolutions[self.vorticity],
- self.domain.length, comm=comm.py2f())
-
- topodims = self.resolutions[self.vorticity] / localres
- print ('topodims DIFFUSION', topodims)
- # variables discretization
-
- for v in self.variables:
- topo = self.domain.getOrCreateTopology(self.domain.dimension,
- self.resolutions[v],
- topodims,
- comm=comm)
- vd = v.discretize(topo)
- self.discreteFields[v] = vd
-
- self.discrete_op =\
- DiffusionFFT(self.discreteFields[self.velocity],
- self.discreteFields[self.vorticity],
- self.method, **self.config)
-
- self.discrete_op.setup()
- self._is_uptodate = True
diff --git a/hysop/operator/diffusion.py b/hysop/operator/diffusion.py
index bbdf1239fe7947767626e6c5e6291028d29eeba6..a7d81be9b21302658ce74375f30820690265871a 100644
--- a/hysop/operator/diffusion.py
+++ b/hysop/operator/diffusion.py
@@ -6,59 +6,61 @@ See :ref:`diffusion` in HySoP user guide.
"""
from hysop.operator.computational import Computational
-from hysop.operator.discrete.diffusion_fft import DiffusionFFT
+from hysop.operator.discrete.diffusion_fft import DiffusionFFT,\
+ CurlAndDiffusionFFT
from hysop.constants import debug
from hysop.operator.continuous import opsetup
-from hysop.methods_keys import SpaceDiscretisation
-import hysop.default_methods as default
-from hysop import __FFTW_ENABLED__
+from hysop.methods_keys import SpaceDiscretisation, GhostUpdate
+from hysop import __FFTW_ENABLED__, __GPU_ENABLED__
+if __GPU_ENABLED__:
+ from hysop.gpu.gpu_diffusion import GPUDiffusion
class Diffusion(Computational):
- """ Diffusion of a field.
-
-
- \f{eqnarray*}
- \omega = Op(\omega)
- \f} with :
- \f{eqnarray*}
- \frac{\partial \omega}{\partial t} &=& \nu\Delta\omega
- \f}
+ """Diffusion of a field.
"""
+ _authorized_methods = ['fftw', 'on_gpu']
+
@debug
- def __init__(self, viscosity, vorticity=None, **kwds):
- """
- Constructor for the diffusion operator.
- @param[in,out] vorticity : vorticity field. If None, it must be passed
- through variables argument
- @param[in] viscosity : viscosity of the considered medium.
- """
- if vorticity is not None:
- super(Diffusion, self).__init__(variables=[vorticity], **kwds)
- else:
- super(Diffusion, self).__init__(**kwds)
+ def __init__(self, viscosity, vorticity, **kwds):
+ """Diffusion operator.
+ See :ref:`diffusion` in HySoP user guide.
- # The only available method at the time is fftw
+ Parameters
+ ----------
+ viscosity : double
+ constant viscosity value
+ vorticity : :class:`~hysop.fields.continuous.Field`
+ vorticity field, in/out parameter.
+ kwds : base class parameters.
+ """
+ assert 'variables' not in kwds, 'variables parameter is useless.'
+ super(Diffusion, self).__init__(variables=[vorticity], **kwds)
+ msg = 'Diffusion : unknown method for space discretization'
if self.method is None:
+ import hysop.default_methods as default
self.method = default.DIFFUSION
- ## input/output field, solution of the problem
- if vorticity is not None:
- self.vorticity = vorticity
- else:
- self.vorticity = self.variables.keys()[0]
- ## viscosity
+ assert self.method[SpaceDiscretisation] in self._authorized_methods,\
+ msg
+ msg = 'Diffusion : on_gpu resolution is not available on your system.'
+ if self.method[SpaceDiscretisation] is 'on_gpu':
+ assert __GPU_ENABLED__, msg
+ # input/output field, solution of the problem
+ self.vorticity = vorticity
+ # viscosity
self.viscosity = viscosity
-
- self.kwds = kwds
-
+ # kwds required for gpu discrete operator
+ self.kwds = kwds.copy()
+ if 'discretization' in self.kwds:
+ self.kwds.pop('discretization')
self.input = [self.vorticity]
self.output = [self.vorticity]
def discretize(self):
if self.method[SpaceDiscretisation] is 'fftw' and __FFTW_ENABLED__:
super(Diffusion, self)._fftw_discretize()
- elif self.method[SpaceDiscretisation] is 'fd':
+ elif self.method[SpaceDiscretisation] is 'on_gpu':
super(Diffusion, self)._standard_discretize()
else:
raise AttributeError("Method not yet implemented.")
@@ -67,18 +69,75 @@ class Diffusion(Computational):
@opsetup
def setup(self, rwork=None, iwork=None):
if self.method[SpaceDiscretisation] is 'fftw':
- self.discrete_op = DiffusionFFT(
- self.discreteFields[self.vorticity], self.viscosity,
- method=self.method)
- elif self.method[SpaceDiscretisation] is 'fd':
- from hysop.gpu.gpu_diffusion import GPUDiffusion
- kw = self.kwds.copy()
- if 'discretization' in kw.keys():
- kw.pop('discretization')
- self.discrete_op = GPUDiffusion(
- self.discreteFields[self.vorticity],
- viscosity=self.viscosity,
- **kw)
+ dop = DiffusionFFT
+ elif self.method[SpaceDiscretisation] is 'on_gpu':
+ dop = GPUDiffusion
+
+ self.discrete_op = dop(
+ viscosity=self.viscosity,
+ vorticity=self.discreteFields[self.vorticity],
+ method=self.method, **self.kwds)
+ self._is_uptodate = True
+
+ def get_work_properties(self):
+ return {'rwork': None, 'iwork': None}
+
+
+class CurlAndDiffusion(Computational):
+ """Solve curl of velocity and its diffusion in one shot
+ (in Fourier domain).
+ """
+
+ _authorized_methods = ['fftw']
+
+ @debug
+ def __init__(self, viscosity, velocity, vorticity, **kwds):
+ """Diffusion operator.
+ See :ref:`diffusion` in HySoP user guide.
+
+ Parameters
+ ----------
+ viscosity : double
+ constant viscosity value
+ velocity : :class:`~hysop.fields.continuous.Field
+ vorticity field, in/out parameter.
+ vorticity : :class:`~hysop.fields.continuous.Field
+ vorticity field, in/out parameter.
+ kwds : base class parameters.
+
+ Notes:
+ * vorticity parameter is optional since the field
+ can be provided in the usual way for operators using
+ variables parameter::
+
+ op = Diffusion(variables={vorticity: topo}, ...)
+
+ """
+ self.viscosity = viscosity
+ self.velocity = velocity
+ self.vorticity = vorticity
+ msg = 'fftw required for CurlAndDiffusion. '
+ msg += 'Try to recompile with WITH_FFTW=ON'
+ assert __FFTW_ENABLED__, msg
+ msg = 'CurlAndDiffusion : unknown method for space discretization'
+ self.method = {SpaceDiscretisation: 'fftw', GhostUpdate: True}
+ self.input = [self.velocity]
+ self.output = [self.vorticity]
+ assert 'variables' not in kwds, 'variables parameter is useless.'
+ super(CurlAndDiffusion, self).__init__(
+ variables=[velocity, vorticity], **kwds)
+
+ def discretize(self):
+ super(CurlAndDiffusion, self)._fftw_discretize()
+
+ @debug
+ @opsetup
+ def setup(self, rwork=None, iwork=None):
+ self.discrete_op = CurlAndDiffusionFFT(
+ viscosity=self.viscosity,
+ velocity=self.discreteFields[self.velocity],
+ vorticity=self.discreteFields[self.vorticity],
+ method=self.method)
self._is_uptodate = True
def get_work_properties(self):
diff --git a/hysop/operator/discrete/curlAndDiffusion_fft.py b/hysop/operator/discrete/curlAndDiffusion_fft.py
deleted file mode 100644
index d8828c7fd1e5c6769737938358f7c09c5a3352f4..0000000000000000000000000000000000000000
--- a/hysop/operator/discrete/curlAndDiffusion_fft.py
+++ /dev/null
@@ -1,117 +0,0 @@
-# -*- coding: utf-8 -*-
-"""
-@file diffusion_fft.py
-Discrete Diffusion 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
-from hysop.tools.profiler import profile
-
-
-class DiffusionFFT(DiscreteOperator):
- """
- Discretized Poisson operator based on FFTW.
- See details in hysop.operator.diffusion.
-
- """
- @debug
- def __init__(self, vorticity, viscosity, method=None):
- """
- Constructor.
- @param[in,out] vorticity : discretisation of the field \f$ \omega \f$.
- @param[in] viscosity : \f$\nu\f$, viscosity of the considered medium.
- """
-
- DiscreteOperator.__init__(self, [vorticity], method,
- name="Diffusion FFT")
- ## Velocity.
- self.velocity = velocity
- ## Vorticity.
- self.vorticity = vorticity
- ## Viscosity.
- self.viscosity = viscosity
- ## Boolean : pure vort diffusion or curl(u) + vort diffusion.
- self.with_curl = with_curl
-
- @debug
- @profile
- def apply(self, simulation):
-
- if (self.vorticity.dimension == 2):
-
- if self.with_curl:
- raise ValueError("Not yet implemented")
-
- else:
- self.vorticity.data = \
- fftw2py.solve_diffusion_2d(self.viscosity * dt,
- self.vorticity.data)
-
- elif (self.vorticity.dimension == 3):
-
- if self.with_curl:
- self.vorticity.data[0], self.vorticity.data[1],\
- self.vorticity.data[2] = \
- fftw2py.solve_curl_diffusion_3d(self.viscosity * dt,
- self.velocity.data[0],
- self.velocity.data[1],
- self.velocity.data[2],
- self.vorticity.data[0],
- self.vorticity.data[1],
- self.vorticity.data[2])
-
- else:
- self.vorticity.data[0], self.vorticity.data[1],\
- self.vorticity.data[2] = \
- fftw2py.solve_diffusion_3d(self.viscosity * dt,
- self.vorticity.data[0],
- self.vorticity.data[1],
- self.vorticity.data[2])
-
- else:
- raise ValueError("invalid problem dimension")
-# ind0a = self.topology.mesh.local_start[0]
-# ind0b = self.topology.mesh.local_end[0] + 1
-# ind1a = self.topology.mesh.local_start[1]
-# ind1b = self.topology.mesh.local_end[1] + 1
-# ind2a = self.topology.mesh.local_start[2]
-# ind2b = self.topology.mesh.local_end[2] + 1
-
-# vorticityNoG = [npw.zeros((self.resolution - 2 * self.ghosts))
-# for d in xrange(self.dim)]
-# velocityNoG = [nwp.zeros((self.resolution - 2 * self.ghosts))
-# for d in xrange(self.dim)]
-# for i in xrange(self.dim):
-# vorticityNoG[i][...] = self.vorticity[i][ind0a:ind0b,
-# ind1a:ind1b, ind2a:ind2b]
-# velocityNoG[i][...] = self.velocity[i][ind0a:ind0b,
-# ind1a:ind1b, ind2a:ind2b]
-
-# # Curl + Vorticity diffusion
-## vorticityNoG[0][...], vorticityNoG[1][...], vorticityNoG[2][...] = \
-## fftw2py.solve_curl_diffusion_3d(self.viscosity * dt,
-## velocityNoG[0][...],
-## velocityNoG[1][...],
-## velocityNoG[2][...],
-## vorticityNoG[0][...],
-## vorticityNoG[1][...],
-## vorticityNoG[2][...])
-
-# # Pure Vorticity diffusion
-# vorticityNoG[0][...], vorticityNoG[1][...], vorticityNoG[2][...] = \
-# fftw2py.solve_diffusion_3d(self.viscosity * dt,
-# vorticityNoG[0][...],
-# vorticityNoG[1][...],
-# vorticityNoG[2][...])
-
-# for i in xrange(self.dim):
-# self.vorticity[i][ind0a:ind0b, ind1a:ind1b, ind2a:ind2b] = \
-# vorticityNoG[i][...]
-
- def __str__(self):
- s = "Diffusion_d (DiscreteOperator). " + DiscreteOperator.__str__(self)
- return s
diff --git a/hysop/operator/discrete/diffusion_fft.py b/hysop/operator/discrete/diffusion_fft.py
index f93ce20d36bd9e26af9ebc1ea8e82e14d32229cd..461a74c78e9c2f48924f01f7d6bb202df531b873 100644
--- a/hysop/operator/discrete/diffusion_fft.py
+++ b/hysop/operator/discrete/diffusion_fft.py
@@ -1,5 +1,8 @@
# -*- coding: utf-8 -*-
"""Discrete Diffusion operator using FFTW (fortran)
+
+See :ref:`diffusion` in HySoP user guide.
+
"""
try:
from hysop.f2hysop import fftw2py
@@ -19,14 +22,15 @@ class DiffusionFFT(DiscreteOperator):
"""
@debug
- def __init__(self, vorticity, viscosity, **kwds):
+ def __init__(self, viscosity, vorticity, **kwds):
"""Discrete diffusion operator, based on fftw solver.
Parameters
----------
+ viscosity : double
+ constant viscosity value
vorticity : :class:`~hysop.fields.discrete.DiscreteField`
vorticity field, in/out parameter
- viscosity : double
kwds : base class arguments
"""
# Discretisation of the solution field
@@ -79,3 +83,52 @@ class DiffusionFFT(DiscreteOperator):
# TODO : fix bug that occurs when several finalize
# of fft operators are called.
# fftw2py.clean_fftw_solver(self.vorticity.dimension)
+
+
+class CurlAndDiffusionFFT(DiscreteOperator):
+ """Discretized Curl/Diffusion operator based on FFTW.
+ See details in hysop.operator.diffusion.
+
+ """
+ @debug
+ def __init__(self, viscosity, velocity, vorticity, **kwds):
+ """Solve diffusion problem in Fourier
+ domain (velocity curl + diffusion in one shot)
+
+ Parameters
+ ----------
+ viscosity : double
+ constant viscosity value
+ velocity : :class:`~hysop.fields.discrete.DiscreteField`
+ velocity field, input parameter
+ vorticity : :class:`~hysop.fields.discrete.DiscreteField`
+ velocity field, output parameter
+ kwds : base class arguments (vorticity, viscosity ...)
+ """
+ self.velocity = velocity
+ self.vorticity = vorticity
+ self.viscosity = viscosity
+ assert 'variables' not in kwds, 'variables parameter is useless.'
+ super(CurlAndDiffusionFFT, self).__init__(
+ variables=[velocity, vorticity], **kwds)
+ self.input = [self.velocity]
+ self.output = [self.vorticity]
+ msge = 'CurlAndDiffusion: only implemented for 3d domain.'
+ assert self.velocity.dimension == 3, msge
+
+ @debug
+ @profile
+ def apply(self, simulation=None):
+ dt = simulation.time_step
+ gh_velo = self.velocity.topology.ghosts()
+ gh_vorti = self.vorticity.topology.ghosts()
+ self.vorticity.data[0], self.vorticity.data[1],\
+ self.vorticity.data[2] = \
+ fftw2py.solve_curl_diffusion_3d(self.viscosity * dt,
+ self.velocity.data[0],
+ self.velocity.data[1],
+ self.velocity.data[2],
+ self.vorticity.data[0],
+ self.vorticity.data[1],
+ self.vorticity.data[2],
+ gh_velo, gh_vorti)
diff --git a/hysop/operator/discrete/poisson_fft.py b/hysop/operator/discrete/poisson_fft.py
index c26dfb106163d04e8fba4ea1770f0f5f48abaac7..5d3a0fb8f011aae5249ba01b73192cd8dccbe0d3 100644
--- a/hysop/operator/discrete/poisson_fft.py
+++ b/hysop/operator/discrete/poisson_fft.py
@@ -1,5 +1,8 @@
# -*- coding: utf-8 -*-
"""Discrete operator for Poisson problem (fftw based)
+
+.. currentmodule hysop.operator.discrete
+
"""
import hysop.tools.numpywrappers as npw
@@ -7,35 +10,39 @@ from hysop.operator.discrete.discrete import DiscreteOperator
from hysop.operator.discrete.reprojection import Reprojection
from hysop.constants import debug
from hysop.tools.profiler import profile
-from hysop import __FFTW_ENABLED__
-
-if __FFTW_ENABLED__:
+try:
from hysop.f2hysop import fftw2py
+except ImportError:
+ msg = 'fftw package not available for your hysop install.'
+ msg += 'Try to recompile with WITH_FFTW=ON'
+ raise ImportError(msg)
class PoissonFFT(DiscreteOperator):
- """
- Discretized Poisson operator based on FFTW.
+ """Discretized Poisson operator based on FFTW.
See details in hysop.operator.poisson
"""
@debug
def __init__(self, output_field, input_field, projection=None,
- filterSize=None, correction=None, formulation=None, **kwds):
- """
- Constructor.
- @param[out] output_field : discretization of the solution field
- @param[in] input_field : discretization of the RHS (mind the minus rhs!)
- @param projection : if None, no projection. Else:
- - either the value of the frequency of reprojection, never updated.
- - or Reprojection discrete operator. In that case, a criterion
- depending on the input_field will be computed at each time step, if
- criterion > threshold, then frequency projection is active.
- @param filterSize :
- @param correction : operator used to shift output_field according
- to a given input (fixed) flowrate.
- See hysop.operator.velocity_correction.
- Default = None.
+ filter_size=None, correction=None, formulation=None, **kwds):
+ """Poisson operator, incompressible flow.
+
+ Parameters
+ ----------
+ output_field : :class:`~hysop.fields.discrete.DiscreteField
+ solution field
+ input_field : :class:`~hysop.fields.discrete.DiscreteField
+ right-hand side
+ projection : double or tuple, optional
+ projection criterion, see notes below.
+ filter_size :
+ correction : :class:`~velocity_correction.VelocityCorrection_D
+ operator used to shift output_field according
+ to a given input (fixed) flowrate.
+ See hysop.operator.velocity_correction.
+ formulation :
+ kwds : base class parameters.
"""
# Base class initialisation
assert 'variables' not in kwds, 'variables parameter is useless.'
@@ -48,11 +55,13 @@ class PoissonFFT(DiscreteOperator):
# Solenoidal projection of input_field ?
self.projection = projection
# Filter size array = domainLength/(CoarseRes-1)
- self.filterSize = filterSize
+ self.filter_size = filter_size
# If 2D problem, input_field must be a scalar
- self.dim = self.output_field.domain.dimension
- if self.dim == 2:
+ self._dim = self.output_field.domain.dimension
+ if self._dim == 2:
assert self.input_field.nb_components == 1
+ assert (self._dim >= 2),\
+ "Wrong problem dimension: only 2D and 3D cases are implemented."
self.correction = correction
self.formulation = formulation
self.input = [self.input_field]
@@ -66,26 +75,24 @@ class PoissonFFT(DiscreteOperator):
self._select_solve()
def _select_solve(self):
-
"""
TODO : add pressure solver selection
f(output.nbComponents) + pb type 'pressure-poisson'
"""
-
- # Multiresolution ?
- multires = self.output_field.topology.mesh != self.input_field.topology.mesh
-
+ # Multiresolution?
+ multires = \
+ self.output_field.topology.mesh != self.input_field.topology.mesh
# connexion to the required apply function
- if self.dim == 2:
- self._solve = self._solve2D
- elif self.dim == 3:
+ if self._dim == 2:
+ self._solve = self._solve_2d
+ elif self._dim == 3:
# If there is a projection, input_field is also an output
if self.projection is not None:
self.output.append(self.input_field)
if multires:
- self._solve = self._solve3D_proj_multires
+ self._solve = self._solve_3d_proj_multires
else:
- self._solve = self._solve3D_proj
+ self._solve = self._solve_3d_proj
if isinstance(self.projection, Reprojection):
self.do_projection = self.do_projection_with_op
@@ -94,42 +101,34 @@ class PoissonFFT(DiscreteOperator):
else:
if multires:
- self._solve = self._solve3D_multires
+ self._solve = self._solve_3d_multires
elif self.formulation is not None:
self._solve = self._solve_3d_scalar_fd
else:
- self._solve = self._solve3D
- else:
- raise AttributeError('Not implemented for 1D problems.')
+ self._solve = self._solve_3d
- # Operator to shift output_field according to an input required flowrate
+ # swtich to the proper solving method (with or without correction)
if self.correction is not None:
self.solve = self._solve_and_correct
else:
self.solve = self._solve
-
- if not __FFTW_ENABLED__:
- self.solve = self._solve_error
-
- def _solve_error():
- """Error message when fftw fortran interface is not available.
- """
- msg = "Warning: fftw fortran interface is not available,"
- msg += "the fft solver does nothing."
- raise NotImplementedError('msg')
def do_projection_with_op(self, simu):
+ """Compute projection criterion and return
+ true if projection is needed.
+ """
self.projection.apply(simu)
ite = simu.current_iteration
return self.projection.do_projection(ite)
def do_projection_no_op(self, simu):
+ """return true if projection is needed.
+ """
ite = simu.current_iteration
return ite % self.projection == 0
- def _solve2D(self, simu=None):
- """
- Solve 2D poisson problem
+ def _solve_2d(self, simu=None):
+ """ Solve 2D poisson problem, no projection, no correction.
"""
ghosts_v = self.output_field.topology.ghosts()
ghosts_w = self.input_field.topology.ghosts()
@@ -140,75 +139,75 @@ class PoissonFFT(DiscreteOperator):
ghosts_w, ghosts_v)
def _project(self):
- """
- apply projection onto input_field
+ """apply projection onto input_field
"""
ghosts_w = self.input_field.topology.ghosts()
self.input_field.data[0], self.input_field.data[1], \
self.input_field.data[2] = \
- fftw2py.projection_om_3d(self.input_field.data[0],
- self.input_field.data[1],
- self.input_field.data[2], ghosts_w)
+ fftw2py.projection_om_3d(self.input_field.data[0],
+ self.input_field.data[1],
+ self.input_field.data[2], ghosts_w)
- def _solve3D_multires(self, simu=None):
- """
- 3D, multiresolution
+ def _solve_3d_multires(self, simu=None):
+ """3D, multiresolution
"""
# Projects input_field values from fine to coarse grid
# in frequencies space by nullifying the smallest modes
vortFilter = npw.copy(self.input_field.data)
vortFilter[0], vortFilter[1], vortFilter[2] = \
- fftw2py.multires_om_3d(self.filterSize[0], self.filterSize[1],
- self.filterSize[2], self.input_field.data[0],
- self.input_field.data[1],
- self.input_field.data[2])
+ fftw2py.multires_om_3d(self.filter_size[0], self.filter_size[1],
+ self.filter_size[2],
+ self.input_field.data[0],
+ self.input_field.data[1],
+ self.input_field.data[2])
# Solves Poisson equation using filter input_field
ghosts_v = self.output_field.topology.ghosts()
ghosts_w = self.input_field.topology.ghosts()
- self.output_field.data[0], self.output_field.data[1], self.output_field.data[2] = \
+ self.output_field.data[0], self.output_field.data[1],\
+ self.output_field.data[2] = \
fftw2py.solve_poisson_3d(vortFilter[0], vortFilter[1],
vortFilter[2], self.output_field.data[0],
self.output_field.data[1],
self.output_field.data[2],
ghosts_w, ghosts_v)
- def _solve3D_proj_multires(self, simu):
- """
- 3D, multiresolution, with projection
+ def _solve_3d_proj_multires(self, simu):
+ """3D, multiresolution, with projection
"""
if self.do_projection(simu):
self._project()
- self._solve3D_multires()
+ self._solve_3d_multires()
- def _solve3D_proj(self, simu):
- """
- 3D, with projection
+ def _solve_3d_proj(self, simu):
+ """3D, with projection
"""
if self.do_projection(simu):
self._project()
- self._solve3D()
+ self._solve_3d()
- def _solve3D(self,simu=None):
- """
- Basic solve
+ def _solve_3d(self,simu=None):
+ """Basic solve
"""
# Solves Poisson equation using usual input_field
ghosts_v = self.output_field.topology.ghosts()
ghosts_w = self.input_field.topology.ghosts()
- self.output_field.data[0], self.output_field.data[1], self.output_field.data[2] =\
+ self.output_field.data[0], self.output_field.data[1],\
+ self.output_field.data[2] =\
fftw2py.solve_poisson_3d(self.input_field.data[0],
self.input_field.data[1],
self.input_field.data[2],
self.output_field.data[0],
self.output_field.data[1],
- self.output_field.data[2], ghosts_w, ghosts_v)
+ self.output_field.data[2],
+ ghosts_w, ghosts_v)
def _solve_and_correct(self, simu):
+ """Solve Poisson problem and apply correction on velocity.
+ """
self._solve(simu.current_iteration)
self.correction.apply(simu)
-
def _solve_3d_scalar_fd(self, simu=None):
"""solve poisson-pressure like problem
input = 3D vector field
@@ -222,7 +221,7 @@ class PoissonFFT(DiscreteOperator):
ghosts = self.output_field.topology.ghosts()
self.output_field.data[0] = fftw2py.pressure_3d(
self.input_field.data[0], ghosts)
-
+
def _solve_3d_scalar(self, simu=None):
"""solve poisson-pressure like problem
input = 3D vector field
diff --git a/hysop/operator/discrete/reprojection.py b/hysop/operator/discrete/reprojection.py
index 86506809712df2b94b312f4d3e7aa85d334a86a7..0e680bb4bff94768f7736da23157d00e52f89189 100644
--- a/hysop/operator/discrete/reprojection.py
+++ b/hysop/operator/discrete/reprojection.py
@@ -1,10 +1,8 @@
# -*- coding: utf-8 -*-
-"""
-@file reprojection.py
-Compute reprojection criterion and divergence maximum
+"""Compute reprojection criterion and divergence maximum
"""
import numpy as np
-from hysop.constants import debug, HYSOP_MPI_REAL
+from hysop.constants import debug
from hysop.methods_keys import SpaceDiscretisation
from hysop.operator.discrete.discrete import DiscreteOperator
from hysop.numerics.finite_differences import FDC4
@@ -16,37 +14,40 @@ from hysop.tools.profiler import profile
class Reprojection(DiscreteOperator):
- """
- Update the reprojection frequency, according to the current
+ """Update the reprojection frequency, according to the current
value of the vorticity field.
"""
def __init__(self, vorticity, threshold, frequency, **kwds):
"""
- Constructor.
- @param vorticity: discretization of the vorticity field
- @param threshold : update frequency when criterion is greater than
- this threshold
- @param frequency : set frequency of execution of the reprojection
+
+ Parameters
+ ----------
+ vorticity : :class:`~hysop.fields.discrete.DiscreteField
+ vorticity field
+ threshold : double
+ update frequency when criterion is greater than this threshold
+ frequency : double
+ frequency of execution of the reprojection
"""
+ # ensure that space discretisation method is set
if 'method' in kwds and kwds['method'] is None:
kwds['method'] = {SpaceDiscretisation: FDC4}
- ## vorticity field
+ # vorticity field
self.vorticity = vorticity
assert 'variables' not in kwds, 'variables parameter is useless.'
super(Reprojection, self).__init__(variables=[vorticity], **kwds)
- ## Frequency for reprojection
+ # Frequency for reprojection
self.frequency = frequency
- ## The initial value will be used as default during
+ # The initial value will be used as default during
# simulation
self._default_frequency = frequency
- # constant defining the reprojection criterion :
+ # constant threshold defining the reprojection criterion :
# if the latter is greater than this constant, then a reprojection
# is needed
self.threshold = threshold
- # local counter
self._counter = 0
- ## Numerical methods for space discretization
+ # Numerical methods for space discretization
assert SpaceDiscretisation in self.method
self.method = self.method[SpaceDiscretisation]
self.input = [vorticity]
@@ -123,7 +124,6 @@ class Reprojection(DiscreteOperator):
self._writer.write()
def do_projection(self, ite):
- """
- True if projection must be done
+ """True if projection must be done
"""
return ite % self.frequency == 0
diff --git a/hysop/operator/discrete/velocity_correction.py b/hysop/operator/discrete/velocity_correction.py
index c09104af5f9a0f59a617f738bada36ef3e5f3a1a..8a187025207e3a97bc8c4a4db83f1202d5779b49 100755
--- a/hysop/operator/discrete/velocity_correction.py
+++ b/hysop/operator/discrete/velocity_correction.py
@@ -1,10 +1,11 @@
# -*- coding: utf-8 -*-
-"""
-@file operator/discrete/velocity_correction.py
+"""Operator used to shift velocity value
+ to fit with a required input flowrate.
-Correction of the velocity field.
-"""
+Check details in :ref:`velocity_correction` in HySoP user guide.
+
+"""
from hysop.constants import debug
from hysop.operator.discrete.discrete import DiscreteOperator
from hysop.fields.variable_parameter import VariableParameter
@@ -14,44 +15,52 @@ from hysop.constants import XDIR, YDIR, ZDIR
class VelocityCorrection_D(DiscreteOperator):
- """
- The velocity field is corrected after solving the
+ """The velocity field is corrected after solving the
Poisson equation. For more details about calculations,
- see the "velocity_correction.pdf" explanations document
- in Docs.
+ see :ref:`velocity_correction` in HySoP user guide.
"""
@debug
def __init__(self, velocity, vorticity, req_flowrate, cb, **kwds):
- """
- @param[in, out] velocity field to be corrected
- @param[in] vorticity field used to compute correction
- @param[in] req_flowrate : required value for the
- flowrate (VariableParameter object)
- @param[in] surf : surface (hysop.domain.obstacle.planes.SubPlane)
- used to compute reference flow rates. Default = surface at x_origin,
- normal to x-dir.
+ """Update velocity field (solution of Poisson equation)
+ in order to prescribe proper mean flow and ensure
+ the desired inlet flowrate.
+
+ Parameters
+ ----------
+ velocity : :class:`~hysop.fields.discrete.DiscreteField`
+ in/out velocity vector field.
+ vorticity : :class:`~hysop.fields.discrete.DiscreteField`
+ input vorticity vector field.
+ req_flowrate : a
+ :class:`~hysop.fields.variable_parameter.VariableParameter`
+ the desired inlet flowrate value
+ cb : :class:`~hysop.domain.control_box.ControlBox`
+ volume of control used to compute a reference for correction.
+ kwds : base class parameters
"""
assert 'variables' not in kwds, 'variables parameter is useless.'
super(VelocityCorrection_D, self).__init__(
variables=[velocity, vorticity], **kwds)
- ## velocity discrete field
+ # velocity discrete field
self.velocity = velocity
- ## vorticity discrete field
+ # vorticity discrete field
self.vorticity = vorticity
- ## domain dimension
- self.dim = self.velocity.domain.dimension
# If 2D problem, vorticity must be a scalar
- if self.dim == 2:
+ if self._dim == 2:
assert self.vorticity.nb_components == 1
- assert (self.dim >= 2),\
+ self._update_velocity_components = self._update_velocity_2d
+ elif self._dim == 3:
+ self._update_velocity_components = self._update_velocity_3d
+
+ assert (self._dim >= 2),\
"Wrong problem dimension: only 2D and 3D cases are implemented."
self.input = self.variables
self.output = [self.velocity]
- ## A reference topology
+ # A reference topology
self.topo = self.velocity.topology
- ## Volume of control
+ # Volume of control
self.cb = cb
self.cb.discretize(self.topo)
# A reference surface, i.e. input surface for flow in x direction
@@ -62,13 +71,13 @@ class VelocityCorrection_D(DiscreteOperator):
cb_length = self.cb.real_length[self.topo]
self._inv_ds = 1. / npw.prod(cb_length[sdirs])
self._inv_dvol = 1. / npw.prod(cb_length)
- ## Expected value for the flow rate through self.surfRef
+ # Expected value for the flow rate through self.surfRef
self.req_flowrate = req_flowrate
assert isinstance(self.req_flowrate, VariableParameter),\
"the required flowrate must be a VariableParameter object."
- ## The correction that must be applied on each
- ## component of the velocity.
- self.velocity_shift = npw.zeros(self.dim)
+ # The correction that must be applied on each
+ # component of the velocity.
+ self.velocity_shift = npw.zeros(self._dim)
nbf = self.velocity.nb_components + self.vorticity.nb_components
# temp buffer, used to save flow rates and mean
# values of vorticity
@@ -83,20 +92,19 @@ class VelocityCorrection_D(DiscreteOperator):
# surface for the flow.
self.x_coord = self.topo.mesh.coords[XDIR] - x0
- def computeCorrection(self):
- """
- Compute the required correction for the current state
+ def compute_correction(self):
+ """Compute the required correction for the current state
but do not apply it onto velocity.
"""
- ## Computation of the flowrates evaluated from
- ## current (ie non corrected) velocity
+ # Computation of the flowrates evaluated from
+ # current (ie non corrected) velocity
nbf = self.velocity.nb_components + self.vorticity.nb_components
localrates = npw.zeros((nbf))
for i in xrange(self.velocity.nb_components):
localrates[i] = self._in_surf.integrate_dfield_on_proc(
self.velocity, component=i)
start = self.velocity.nb_components
- ## Integrate vorticity over the whole domain
+ # Integrate vorticity over the whole domain
for i in xrange(self.vorticity.nb_components):
localrates[start + i] = self.cb.integrate_dfield_on_proc(
self.vorticity, component=i)
@@ -110,7 +118,7 @@ class VelocityCorrection_D(DiscreteOperator):
self.rates[:start] *= self._inv_ds
self.rates[start:] *= self._inv_dvol
-
+
# Set velocity_shift == [Vx_shift, vort_mean[Y], vort_mean[Z]]
# or (in 2D) velocity_shift == [Vx_shift, vort_mean]
# Velocity shift for main dir component
@@ -119,6 +127,20 @@ class VelocityCorrection_D(DiscreteOperator):
# Shifts in other directions depend on x coord
# and will be computed during apply.
+ def _update_velocity_2d(self, vort_mean):
+ """update velocity value, 2d case
+ """
+ self.velocity[YDIR][...] += self.req_flowrate_val[YDIR] +\
+ vort_mean[XDIR] * self.x_coord - self.rates[YDIR]
+
+ def _update_velocity_3d(self, vort_mean):
+ """update velocity value, 3d case
+ """
+ self.velocity[YDIR][...] += self.req_flowrate_val[YDIR] + \
+ vort_mean[ZDIR] * self.x_coord - self.rates[YDIR]
+ self.velocity[ZDIR][...] += self.req_flowrate_val[ZDIR] - \
+ vort_mean[YDIR] * self.x_coord - self.rates[ZDIR]
+
@debug
@profile
def apply(self, simulation=None):
@@ -129,10 +151,10 @@ class VelocityCorrection_D(DiscreteOperator):
self.req_flowrate_val = self.req_flowrate.data * self._inv_ds
# Computation of the required velocity shift
# for the current state
- self.computeCorrection()
+ self.compute_correction()
compute_index = self.topo.mesh.compute_index
- # Apply shift to velocity
+ # Apply shift to velocity (x component)
self.velocity[XDIR][compute_index] += self.velocity_shift[XDIR]
start = self.velocity.nb_components
# reminder : self.rates =[vx_shift, flowrates[Y], flowrate[Z],
@@ -146,14 +168,4 @@ class VelocityCorrection_D(DiscreteOperator):
self._writer.buffer[0, 2:] = vort_mean[...]
self._writer.write()
- if self.dim == 2:
- # Correction of the Y-velocity component
- self.velocity[YDIR][...] += self.req_flowrate_val[YDIR] + \
- vort_mean[XDIR] * self.x_coord - self.rates[YDIR]
-
- elif self.dim == 3:
- # Correction of the Y and Z-velocity components
- self.velocity[YDIR][...] += self.req_flowrate_val[YDIR] + \
- vort_mean[ZDIR] * self.x_coord - self.rates[YDIR]
- self.velocity[ZDIR][...] += self.req_flowrate_val[ZDIR] - \
- vort_mean[YDIR] * self.x_coord - self.rates[ZDIR]
+ self._update_velocity_components(vort_mean)
diff --git a/hysop/operator/poisson.py b/hysop/operator/poisson.py
index 324ed1168d074be0c8756d0d46475c910af00daa..da9328ae7ae043546db0f2bd609c65fad210d999 100644
--- a/hysop/operator/poisson.py
+++ b/hysop/operator/poisson.py
@@ -1,5 +1,7 @@
# -*- coding: utf-8 -*-
-"""Poisson problem.
+"""Operator to solve Poisson problem.
+
+See :ref:`poisson` in HySoP user guide.
"""
from hysop.operator.computational import Computational
from hysop.operator.discrete.poisson_fft import PoissonFFT
@@ -8,61 +10,68 @@ from hysop.operator.velocity_correction import VelocityCorrection
from hysop.operator.reprojection import Reprojection
from hysop.methods_keys import SpaceDiscretisation, Formulation
from hysop.operator.continuous import opsetup
-import hysop.default_methods as default
from hysop import __FFTW_ENABLED__
class Poisson(Computational):
+ """Solve Poisson problem (in: vorticity, out velocity)
+ for incompressible flow.
"""
- \f{eqnarray*}
- v = Op(\omega)
- \f} with :
- \f{eqnarray*}
- \Delta \phi &=& -\omega \\
- v &=& \nabla \times \phi
- \f}
- """
+
+ _authorized_methods = ['fftw']
@debug
def __init__(self, output_field, input_field, flowrate=None,
projection=None, **kwds):
- """
- Constructor for the Poisson problem.
-
- @param[out] output_field : solution field
- @param[in] input_field : rhs field
- @param[in] flowrate : a flow rate value (through input_field surf,
- normal to xdir) used to compute a correction of the solution field.
- Default = 0 (no correction). See hysop.operator.output_field_correction.
- @param projection : if None, no projection. Else:
- - either the value of the frequency of reprojection, never update.
- - or a tuple = (frequency, threshold).
- In that case, a criterion
- depending on the input_field will be computed at each time step, if
- criterion > threshold, then frequency projection is active.
-
- Note about method:
- - SpaceDiscretisation == fftw
- - Formulation = 'velocity' or 'pressure'
- velocity : laplacian(phi) = -w and v = nabla X psi, in = vorticity, out = velo
- pressure : laplacian(p) = -nabla.(u.nabla u, in = velo, out = pressure
+ """Poisson operator, incompressible flow.
+
+ Parameters
+ ----------
+ output_field : :class:`~hysop.fields.continuous.Field
+ solution field
+ input_field : :class:`~hysop.fields.continuous.Field`
+ right-hand side
+ flowrate: :class:`~hysop.fields.variable_parameter.VariableParameter`
+ or double, optional
+ flow rate value through input surface (normal to xdir),
+ used to calibrate solution, default=0.
+ projection : double or tuple, optional
+ projection criterion, see notes below.
+ kwds : base class parameters.
+
+
+ Notes:
+ * projection might be:
+ * None: no projection
+ * the value of the frequency of reprojection (constant)
+ * a tuple (frequency, threshold). In that case, a criterion
+ depending on the input_field will be computed at each time step,
+ and if criterion > threshold., then frequency projection is active.
+ * About method parameter:
+ - SpaceDiscretisation == fftw
+ - Formulation = 'velocity' or 'pressure'
+ velocity : laplacian(phi) = -w and v = nabla X psi, in = vorticity, out = velo
+ pressure : laplacian(p) = -nabla.(u.nabla u, in = velo, out = pressure
"""
# Warning : for fftw all variables must have
# the same resolution.
assert 'variables' not in kwds, 'variables parameter is useless.'
super(Poisson, self).__init__(variables=[output_field, input_field],
**kwds)
- ## solution of the problem
+ # solution of the problem
self.output_field = output_field
- ## -(right-hand side)
+ # -(right-hand side)
self.input_field = input_field
if self.method is None:
+ import hysop.default_methods as default
self.method = default.POISSON
+ if self.method[SpaceDiscretisation] is 'fftw':
+ assert __FFTW_ENABLED__
+ msg = 'Poisson : unknown method for space discretization'
+ assert self.method[SpaceDiscretisation] in self._authorized_methods,\
+ msg
- if self.method[SpaceDiscretisation] is not 'fftw':
- raise AttributeError("Method not yet implemented.")
-
- # Deterlination of the Poisson equation formulation :
+ # Set Poisson equation formulation :
# Velo Poisson eq or Pressure Poisson eq
self.formulation = None
if self.method[Formulation] is not 'velocity':
@@ -70,6 +79,7 @@ class Poisson(Computational):
self.input = [self.input_field]
self.output = [self.output_field]
+ # Enable correction if required
if flowrate is not None:
self.withCorrection = True
self._flowrate = flowrate
@@ -84,8 +94,9 @@ class Poisson(Computational):
def discretize(self):
# Poisson solver based on fftw
- if self.method[SpaceDiscretisation] is 'fftw' and __FFTW_ENABLED__:
+ if self.method[SpaceDiscretisation] is 'fftw':
super(Poisson, self)._fftw_discretize()
+ # prepare correction and projection, if needed
if self.withCorrection:
toporef = self.discreteFields[self.output_field].topology
if 'discretization' in self._config:
@@ -102,8 +113,6 @@ class Poisson(Computational):
threshold, freq,
**self._config)
self.projection.discretize()
- else:
- raise AttributeError("Method not implemented.")
@debug
@opsetup
diff --git a/hysop/operator/tests/test_diffusion.py b/hysop/operator/tests/test_diffusion.py
index 5f87c4b335a991956ca4449fc206ac0e55f73eb1..effb9ebe3dcc4d9786560246ad2c53d122fdee46 100755
--- a/hysop/operator/tests/test_diffusion.py
+++ b/hysop/operator/tests/test_diffusion.py
@@ -1,13 +1,16 @@
+"""Tests for diffusion and curl_diffusion operators.
+"""
# -*- coding: utf-8 -*-
import hysop as pp
-from hysop.operator.diffusion import Diffusion
+from hysop.operator.diffusion import Diffusion, CurlAndDiffusion
from hysop.problem.simulation import Simulation
from hysop.tools.parameters import Discretization
import numpy as np
import hysop.tools.numpywrappers as npw
import math
from hysop import testsenv
+from hysop.fields.tests.func_for_tests import v_TG, w_TG
pi = math.pi
sin = np.sin
cos = np.cos
@@ -19,13 +22,6 @@ d3D = Discretization([33, 33, 33])
d2D = Discretization([33, 33])
-def computeVort(res, x, y, z, t):
- res[0][...] = sin(x * cc[0]) * sin(y * cc[1]) * cos(z * cc[2])
- res[1][...] = cos(x * cc[0]) * sin(y * cc[1]) * sin(z * cc[2])
- res[2][...] = cos(x * cc[0]) * cos(y * cc[1]) * sin(z * cc[2])
- return res
-
-
def computeVort2D(res, x, y, t):
# todo ...
res[0][...] = 4 * pi ** 2 * (cos(x * cc[0]) * sin(y * cc[1])) * \
@@ -37,44 +33,62 @@ def computeVort2D(res, x, y, t):
def test_diffusion_3d():
"""Vector field diffusion, 3d domain"""
dom = pp.Box(length=LL)
- vorticity = pp.Field(domain=dom, formula=computeVort,
+ vorticity = pp.Field(domain=dom, formula=w_TG,
name='Vorticity', is_vector=True)
diff = Diffusion(viscosity=0.3, vorticity=vorticity, discretization=d3D)
diff.discretize()
diff.setup()
topo = diff.discreteFields[vorticity].topology
- simu = Simulation(nb_iter=10)
+ simu = Simulation(end=0.1, time_step=0.01)
vorticity.initialize(topo=topo)
diff.apply(simu)
diff.finalize()
@testsenv.fftw_failed
-def test_diffusion_3d_2():
- """Vector field diffusion, 3d domain, """
- dom = pp.Box(length=LL)
- vorticity = pp.Field(domain=dom, formula=computeVort,
- name='Vorticity', is_vector=True)
- diff = Diffusion(viscosity=0.3, variables={vorticity: d3D})
+def test_diffusion_2d():
+ """Vector field diffusion 2d domain"""
+ dom = pp.Box(length=LL[:2])
+ vorticity = pp.Field(domain=dom, formula=computeVort2D, name='Vorticity')
+ diff = Diffusion(viscosity=0.3, vorticity=vorticity, discretization=d2D)
diff.discretize()
diff.setup()
topo = diff.discreteFields[vorticity].topology
- simu = Simulation(nb_iter=10)
+ simu = Simulation(end=0.1, time_step=0.01)
vorticity.initialize(topo=topo)
diff.apply(simu)
diff.finalize()
@testsenv.fftw_failed
-def test_diffusion_2d():
- """Vector field diffusion 2d domain"""
- dom = pp.Box(length=LL[:2])
- vorticity = pp.Field(domain=dom, formula=computeVort2D, name='Vorticity')
- diff = Diffusion(viscosity=0.3, vorticity=vorticity, discretization=d2D)
+def test_curl_and_diffusion_3d():
+ """Vector field diffusion, 3d domain"""
+ dom = pp.Box(length=LL)
+ vorticity_ref = pp.Field(domain=dom, name='Wref', formula=w_TG,
+ is_vector=True)
+ vorticity = pp.Field(domain=dom, name='Vorticity', is_vector=True)
+ velocity = pp.Field(domain=dom, name='Velocity', formula=v_TG,
+ is_vector=True)
+
+ diff_0 = Diffusion(viscosity=0.3, vorticity=vorticity_ref,
+ discretization=d3D)
+ diff_0.discretize()
+ diff_0.setup()
+ topo = diff_0.discreteFields[vorticity_ref].topology
+ diff = CurlAndDiffusion(viscosity=0.3, velocity=velocity,
+ vorticity=vorticity, discretization=d3D)
diff.discretize()
diff.setup()
- topo = diff.discreteFields[vorticity].topology
- simu = Simulation(nb_iter=10)
- vorticity.initialize(topo=topo)
+
+ simu = Simulation(end=0.1, time_step=0.01)
+ vorticity_ref.initialize(topo=topo)
+ velocity.initialize(topo=topo)
+ simu.initialize()
+ wd = vorticity.discretize(topo).data
+ wd_ref = vorticity_ref.discretize(topo).data
+ diff_0.apply(simu)
diff.apply(simu)
- diff.finalize()
+ for i in xrange(3):
+ assert np.allclose(wd[i], wd_ref[i])
+
+ diff_0.finalize()
diff --git a/hysop/operator/tests/test_velocity_correction.py b/hysop/operator/tests/test_velocity_correction.py
index 7c883152045cf49aac2e72ee5bb5563b646294a7..e5196bada8e0435e2b084d726feb659905caca79 100755
--- a/hysop/operator/tests/test_velocity_correction.py
+++ b/hysop/operator/tests/test_velocity_correction.py
@@ -1,11 +1,12 @@
-# -*- coding: utf-8 -*-
-
+"""Test correction operator
+"""
import hysop as pp
from hysop.operator.velocity_correction import VelocityCorrection
from hysop.problem.simulation import Simulation
import numpy as np
import hysop.tools.numpywrappers as npw
from hysop.tools.parameters import Discretization
+from hysop.fields.tests.func_for_tests import v_TG, v_TG_2d, w_TG, w_TG_2d
pi = np.pi
cos = np.cos
sin = np.sin
@@ -14,47 +15,21 @@ uinf = 1.0
tol = 1e-12
-## Function to compute TG velocity
-def computeVel(res, x, y, z, t):
- res[0][...] = sin(x) * cos(y) * cos(z)
- res[1][...] = - cos(x) * sin(y) * cos(z)
- res[2][...] = 0.
- return res
-
-
-## Function to compute reference vorticity
-def computeVort(res, x, y, z, t):
- res[0][...] = - cos(x) * sin(y) * sin(z)
- res[1][...] = - sin(x) * cos(y) * sin(z)
- res[2][...] = 2. * sin(x) * sin(y) * cos(z)
- return res
-
-
-## Function to compute TG velocity
-def computeVel2D(res, x, y, t):
- res[0][...] = sin(x) * cos(y)
- res[1][...] = - cos(x) * sin(y)
- return res
-
-
-## Function to compute reference vorticity
-def computeVort2D(res, x, y, t):
- res[0][...] = - cos(x) * sin(y)
- return res
-
-## Global resolution
+# Global resolution
g = 0
d2D = Discretization([33, 33], [g, g])
d3D = Discretization([33, 33, 33], [g, g, g])
-def test_velocity_correction_3D():
+def test_velocity_correction_3d():
+ """Apply velocity correction, 3d domain
+ """
# Domain
box = pp.Box(length=[2.0 * pi, pi, pi])
# Vector Fields
- velo = pp.Field(domain=box, formula=computeVel,
+ velo = pp.Field(domain=box, formula=v_TG,
name='Velocity', is_vector=True)
- vorti = pp.Field(domain=box, formula=computeVort,
+ vorti = pp.Field(domain=box, formula=w_TG,
name='Vorticity', is_vector=True)
# Usual Cartesian topology definition
@@ -80,17 +55,19 @@ def test_velocity_correction_3D():
assert (np.abs(flowrate - ref_rate) < tol).all()
-def test_velocity_correction_2D():
- ## Domain
+def test_velocity_correction_2d():
+ """Apply velocity correction, 2d domain
+ """
+ # Domain
box = pp.Box(length=[2.0 * pi, pi], origin=[0., 0.])
- ## Vector Fields
- velo = pp.Field(domain=box, formula=computeVel2D,
+ # Vector Fields
+ velo = pp.Field(domain=box, formula=v_TG_2d,
name='Velocity', is_vector=True)
- vorti = pp.Field(domain=box, formula=computeVort2D,
+ vorti = pp.Field(domain=box, formula=w_TG_2d,
name='Vorticity', is_vector=False)
- ## Usual Cartesian topology definition
+ # Usual Cartesian topology definition
topo = box.create_topology(discretization=d2D)
ref_rate = npw.zeros(2)
diff --git a/hysop/operator/velocity_correction.py b/hysop/operator/velocity_correction.py
index 9c863f19a7b0d8055930755eab411810a97f0222..7d3a5dc9196422b0a7c4e4a933f24e0212f9c176 100755
--- a/hysop/operator/velocity_correction.py
+++ b/hysop/operator/velocity_correction.py
@@ -1,8 +1,9 @@
# -*- coding: utf-8 -*-
-"""
-@file operator/velocity_correction.py
+"""Operator used to shift velocity value
+ to fit with a required input flowrate.
+
+Check details in :ref:`velocity_correction` in HySoP user guide.
-Operator to shift velocity to fit with a required input flowrate.
"""
from hysop.constants import debug
@@ -13,44 +14,43 @@ from hysop.operator.continuous import opsetup
class VelocityCorrection(Computational):
- """
- The velocity field is corrected after solving the
+ """The velocity field is corrected after solving the
Poisson equation. For more details about calculations,
- see the "velocity_correction.pdf" explanations document
- in Docs.
+ see :ref:`velocity_correction` in HySoP user guide.
"""
@debug
def __init__(self, velocity, vorticity, req_flowrate, **kwds):
- """
- Corrects the values of the velocity field after
- solving Poisson equation in order to prescribe proper
- mean flow and ensure the desired inlet flowrate.
+ """Update velocity field (solution of Poisson equation)
+ in order to prescribe proper mean flow and ensure
+ the desired inlet flowrate.
- @param[in, out] velocity field to be corrected
- @param[in] vorticity field used to compute correction
- @param resolutions : grid resolutions of velocity and vorticity
- @param[in] req_flowrate : required value for the flowrate
- @param topo : a predefined topology to discretize velocity/vorticity
+ Parameters
+ ----------
+ velocity : :class:`~hysop.fields.continuous.Field`
+ in/out velocity continuous vector field.
+ vorticity : :class:`~hysop.fields.continuous.Field`
+ input vorticity vector field.
+ req_flowrate : a
+ :class:`~hysop.fields.variable_parameter.VariableParameter`
+ the desired inlet flowrate value
+ kwds : base class parameters
"""
assert 'variables' not in kwds, 'variables parameter is useless.'
- super(VelocityCorrection, self).__init__(variables=[velocity,
- vorticity], **kwds)
- ## velocity variable (vector)
+ super(VelocityCorrection, self).__init__(
+ variables=[velocity, vorticity], **kwds)
+ # velocity field
self.velocity = velocity
- ## velocity variable
+ # vorticity field
self.vorticity = vorticity
-
self.input = [self.velocity, self.vorticity]
self.output = [self.velocity]
- ## Expected value for the flow rate through input surface
+ # Expected value for the flow rate through input surface
self.req_flowrate = req_flowrate
dom = self.velocity.domain
+ # volume of control used to compute a reference for correction.
self.cb = ControlBox(origin=dom.origin, length=dom.length,
parent=dom)
- ## Extra parameters that may be required for discrete operator
- ## (at the time, only io_params)
- self.config = kwds
def discretize(self):
super(VelocityCorrection, self)._standard_discretize()
@@ -70,12 +70,11 @@ class VelocityCorrection(Computational):
self.discrete_op.set_writer(self._writer)
self._is_uptodate = True
- def computeCorrection(self):
- """
- Compute the required correction for the current state
+ def compute_correction(self):
+ """Compute the required correction for the current state
but do not apply it onto velocity.
"""
- self.discrete_op.computeCorrection()
+ self.discrete_op.compute_correction()
def get_work_properties(self):
return {'rwork': None, 'iwork': None}