diff --git a/docs/sphinx/conf.py.in b/docs/sphinx/conf.py.in
index 31fe2749ff73432ec2b93c543c466235a800adf2..8399b4a984cfc674a93c06e90e22f278927f4e77 100644
--- a/docs/sphinx/conf.py.in
+++ b/docs/sphinx/conf.py.in
@@ -185,10 +185,10 @@ html_static_path = ['_static']
# '**': ['menu.html',],# 'searchbox.html'],
# 'using/windows': ['windowssidebar.html', 'searchbox.html'],
#}
-html_sidebars = {
- '**': ['menu.html',],# 'searchbox.html'],
-# 'using/windows': ['windowssidebar.html', 'searchbox.html'],
-}
+# html_sidebars = {
+# '**': ['menu.html', ],# 'searchbox.html'],
+# 'using/windows': ['windowssidebar.html', 'searchbox.html'],
+# }
# Additional templates that should be rendered to pages, maps page names to
# template names.
diff --git a/docs/sphinx/hysop.bib b/docs/sphinx/hysop.bib
index 5ba5c1351f9b04a32d262f71e1ef471f682db7a9..63a5a0eae611353ecdb131dcf8a714f378faf3d6 100644
--- a/docs/sphinx/hysop.bib
+++ b/docs/sphinx/hysop.bib
@@ -42,3 +42,10 @@
Title = {A comparison of methods for evaluating time-dependent fluid dynamic forces on bodies, using only velocity fields and their derivatives},
Volume = {13},
Year = {1999}}
+
+@misc{cottet2000vortex,
+ title={Vortex methods: theory and practice},
+ author={Cottet, Georges-Henri and Koumoutsakos, Petros D},
+ year={2000},
+ publisher={Cambridge University Press}
+}
\ No newline at end of file
diff --git a/docs/sphinx/install_guide/index.rst b/docs/sphinx/install_guide/index.rst
new file mode 100644
index 0000000000000000000000000000000000000000..4e6ed9939f2779d0577feba7632894921c216f82
--- /dev/null
+++ b/docs/sphinx/install_guide/index.rst
@@ -0,0 +1,80 @@
+.. hysop_install_guide:
+
+HySoP installation guide
+========================
+
+
+HySoP is written in Python (main user interface and high level functionnalities) and Fortran.
+
+Quick start and install
+-----------------------
+
+::
+
+ cd path_to_build
+ cmake path_to_sources
+ make -j N
+ make install
+
+Note:
+
+* -jN for make command runs the build process on N processes. Choose a value of N that fits with your machine.
+* depending on your system you may need to set PYTHONPATH to where hysop has been installed::
+
+ export PYTHONPATH=hysop_install_path:${PYTHONPATH}
+
+ Check the end of cmake output, where the proper command will be indicated.
+
+
+
+
+Then just run python interactively and ::
+
+
+ >>> import hysop
+ >>> dir(hysop)
+
+
+or run python on an hysop example file::
+
+
+ mpirun -np 4 python your_file.py
+
+
+.. _hysop_config:
+
+Configuration options
+---------------------
+
+cmake configuration may be customized by user, using::
+
+ cmake -DOPTION_NAME=option_value path_to_source
+
+'OPTION_NAME' being one of the options described below.
+
+* FFTW_DIR : where to find fftw if it's not in a "standard" place.
+* WITH_SCALES=ON/OFF : to compile an HySoP version including scales (default = on)
+* WITH_TESTS=ON/OFF: enable testing (i.e. prepare target "make test", default = off)
+
+
+.. _hysop_dependencies:
+
+Dependencies
+------------
+
+* python > 2.7
+* a proper mpi implementation including a fortran compiler
+* fftw
+* cmake > 2.8
+* numpy, mpi4py, h5py
+
+
+.. _hysop_install:
+
+What will be installed
+----------------------
+
+* python package in PREFIX/lib/pythonX.Y/site-packages/hysop/
+* library in PREFIX/lib/pythonX.Y/site-packages/hysop/lib
+* fortran modules in PREFIX/lib/pythonX.Y/site-packages/hysop/include/Modules
+* CMake configuration files in PREFIX/lib/pythonX.Y/site-packages/hysop/share/CMake
diff --git a/docs/sphinx/users_guide/absorption_bc.rst b/docs/sphinx/users_guide/absorption_bc.rst
index c50a907f0e27d5b2e1b7ef64a3b730b945acda2f..1f16a5490d04810b789d19757ff1a67c2d593ebf 100644
--- a/docs/sphinx/users_guide/absorption_bc.rst
+++ b/docs/sphinx/users_guide/absorption_bc.rst
@@ -17,11 +17,11 @@ where f is the smoothing function defined as follows :
.. math::
- f(x) = \begin{cases}
+ f(x) = \begin{cases}
1, \text{ if } x < x_b\\
- \frac{tanh(\epsilon(x - x_c)) - tanh(\epsilon(x_e - x_c))}, \text{ if } x_b < x < x_e \\
+ \frac{tanh(\epsilon(x - x_c)) - tanh(\epsilon(x_e - x_c))}{tanh(\epsilon(x_b - x_c)) - tanh(\epsilon(x_e - x_c))}, \text{ if } x_b < x < x_e \\
0 \text{ if } x > x_e
- \end{cases}
+ \end{cases}
with :math:`x_c = \frac{x_b + x_e}{2}` (we refer to the figure below for a schematic description of the absorption process).
@@ -56,3 +56,36 @@ not the case when the vorticity is absorbed as follows :
Indeed, defining the absorbed vorticity as the curl of the absorbed velocity to project the vorticity in the space of divergence free vector fields
(see :ref:`vorticity_projection`). In practice, vorticity absorption based on equation 6 is responsible for reflexion effects and produces spurious oscillations (see drag plots for instance).
From an algorithmic point of view, the vorticity absorption is applied before solving the Poisson equation.
+
+
+The corresponding operator is :class:AbsorptionBC`, build with the standard parameters of all operators, plus:
+
+* a float, req_flow_rate, to set required value of the flow rate at the outlet.
+* a 'x_range' vector, to set the coordinates delimitating the absorption domain, x_range=[xb, xe]
+* an optional 'filter_func' parameter, list of two python functions, used to set the filter function and its derivative
+
+
+Example::
+
+ d3d = Discretization([33, 33, 33])
+ dom = Box([1., ] * dim)
+ v = Field(dom, name='velo', is_vector=True)
+ w = Field(dom, name='vorti', is_vector=True)
+ # -- Build and apply absorption operator on w --
+ # set flowrate
+ flowrate = VariableParameter(data=10.)
+ # set layer size
+ x_range = [0.9, 1.]
+ op = AbsorptionBC(v, w, flowrate, x_range,
+ discretization=d3d)
+ op.discretize()
+ op.setup()
+
+ simu = Simulation(nb_iter=100)
+ simu.initialize()
+
+ op.apply(simu)
+
+
+One may also choose user-defined functions to apply in the absorption area (in place of the tanh-like function shown above).
+See the example in test_absorption_bc.py file.
diff --git a/docs/sphinx/users_guide/custom.rst b/docs/sphinx/users_guide/custom.rst
new file mode 100644
index 0000000000000000000000000000000000000000..74499c519d2aa96d0e3463d3b909efcf0746b926
--- /dev/null
+++ b/docs/sphinx/users_guide/custom.rst
@@ -0,0 +1,44 @@
+.. _custom:
+
+User-defined operator
+=====================
+
+A specific operator class, :class:`~hysop.operator.custom.Custom` is available in hysop to allow users to
+define their own operator.apply, thanks to a python function.
+
+
+.. code::
+
+ # v1, v2, v3 some continuous fields
+ iop = IOParams('somefile.dat', fileformat=IO.ASCII)
+ my_op = Custom(function=my_function, in_fields=[v1, v2], out_fields=[v3],
+ diagnostics_shape=(1,2), io_params=iop, ...)
+ my_op.discretize()
+ my_op.setup()
+
+ simu = Simulation(nb_iter=10)
+ my_op.apply(simu)
+
+
+A call to my_op.apply(simu) will run 'my_func(simu, [v1, v2], [v3], buffer)'
+and write 'buffer' into a file (defined as usual with iop).
+
+* my_func arguments list MUST match the following::
+
+ def my_func(sim, in_fields, out_fields=None, diagnostics=None):
+ # optionally compute out_fields depending on in_fields
+ out_fields[0].data[1] = 2 * in_fields[0].data[0]
+ # optionally compute some diagnostics
+ diagnostics[0,0] = sim.time
+ diagnostics[0,1] = out_fields[0].data[0].max()
+
+* diagnostics_shape arg must corresponds to the shape of the buffer written in the function
+ and buffer must be a 2D array (i.e. diagnostics_shape a tuple like (nx, ny)).
+
+* if diagnostics_shape is None, no file output will be done and obviously in that case,
+ diagnostics must not be filled in my_func
+
+* as usual real use-case are available in operator/test/test_custom.py
+
+
+
diff --git a/docs/sphinx/users_guide/operators.rst b/docs/sphinx/users_guide/operators.rst
index a9107f4226a2cffeb9d3e54cddd356477342797d..861af66d70722a7d53c6996b431b894d10872465 100644
--- a/docs/sphinx/users_guide/operators.rst
+++ b/docs/sphinx/users_guide/operators.rst
@@ -116,3 +116,4 @@ Review
advection
stretching
energy_enstrophy
+ custom
diff --git a/docs/sphinx/users_guide/stretching.rst b/docs/sphinx/users_guide/stretching.rst
index 79d00d7f01a686a53192459cfc7bcc1b0a86e7da..c8a37b0e951aa5177e595d804734a51f00e21a25 100644
--- a/docs/sphinx/users_guide/stretching.rst
+++ b/docs/sphinx/users_guide/stretching.rst
@@ -4,6 +4,19 @@ Stretching
==========
.. currentmodule:: hysop.operator.stretching
+
+The stretching equation, :math:`\frac{\partial\omega}{\partial t} = (\omega.\nabla)u`
+which models the lengthening of vortices in a three-dimensional flow,
+may be given under several formulations :cite:`Cottet2000vortex`:
+
+
+.. math::
+
+ \frac{\partial\omega}{\partial t} &=& [\nabla u] . \omega \\
+ &=& [\nabla u]^T . \omega \\
+ &=& \nabla.(\omega:u)
+
+
Solve
@@ -20,7 +33,7 @@ Two formulations are proposed :
.. math::
- \frac{\partial \omega}{\partial t} = (\omega .\nabla)u
+ \frac{\partial \omega}{\partial t} = [\nabla u]^T . \omega
The default is 'conservative'.
@@ -46,10 +59,35 @@ TO BE REVIEWED
Linearized Stretching
---------------------
-Specific formulation to solve the linearized stretching equation,
+In the context of a global linear stability study of the Navier-Stokes equation, one may have to formulate and solve
+a "linearized stretching". See details in :ref:`global_linear_stability`.
+
+Let us decompose the vorticity and the velocity fields into the sum of the steady state solution, :math:`(.)_b`, and a perturbation,
+:math:`(.)'`:
.. math::
- \frac{\partial \omega}{\partial t} &=& (\omega \cdot \nabla)u_b + (\om_b \cdot \nabla)u
+ (\omega,u) = (\omega' + \omega_b, u' + u_b)
+
+then the stretching part of the governing equations writes:
+
+.. math::
+
+ \frac{\partial\omega'}{\partial t} = \underbrace{(\omega'\cdot\nabla)u_b}_{A} + \underbrace{(\omega_b\cdot\nabla)u'}_{B}
+
+
+The corresponding operator is :class:`StretchingLinearized`, built with the same parameters as for the standard stretching
+plus two fields corresponding to the steady state part of the velocity and the vorticity ('BF' below stands for 'Base Flow'):
+
+
+.. code::
+
+ stretch_lin = StretchingLinearized(velocity_BF=v_b, vorticity_BF=w_b,
+ velocity=v, vorticity=w, ...)
+
+
+Notice that time-integrator for part A of the equation above is user-defined but set to Euler scheme for part B.
+
+
diff --git a/hysop/constants.py.in b/hysop/constants.py.in
new file mode 100644
index 0000000000000000000000000000000000000000..2dda2f611000d144d3f26f9a12a1505e2edeb4a6
--- /dev/null
+++ b/hysop/constants.py.in
@@ -0,0 +1,135 @@
+"""Default values for constant parameters used
+in hysop package.
+
+Mostly deals with float/int numbers precision and
+some labels for numerical methods and operators.
+
+This file is generated from constants.py.in
+and it's probably a bad idea to modify it.
+
+"""
+from hysop import __DEBUG__, __PROFILE__
+import numpy as np
+import math
+__MPI_ENABLED__ = "@USE_MPI@" is "ON"
+if __MPI_ENABLED__:
+ from hysop.mpi import MPI
+
+PI = math.pi
+
+HYSOP_REAL = @PYTHON_PREC@
+"""Set default type for real and integer numbers"""
+
+SIZEOF_HYSOP_REAL = int(HYSOP_REAL(1.).nbytes)
+"""Size in memory of hysop real default type"""
+
+HYSOP_INDEX = np.uint32
+"""type for array indices"""
+
+HYSOP_INTEGER = np.int32
+"""type for integers"""
+
+HYSOP_DIM = np.int16
+"""integer used for arrays dimensions"""
+
+if __MPI_ENABLED__:
+ HYSOP_MPI_REAL = @MPI_PYTHON_PREC@
+ """float type used in MPI"""
+
+ HYSOP_MPI_INTEGER = MPI.INT
+ """integer type used in MPI"""
+
+ ORDERMPI = @MPI_DATA_LAYOUT@
+ """Default array layout for MPI"""
+
+
+ORDER = @DATA_LAYOUT@
+"""default array layout (fortran or C convention)"""
+
+EPS = np.finfo(HYSOP_REAL).eps
+"""machine epsilon for HYSOP_REAL"""
+
+# to check array ordering with :
+# assert tab.flags.f_contiguous is CHECK_F_CONT
+if ORDER is 'F':
+ CHECK_F_CONT = True
+ """True if array layout is fortran"""
+else:
+ CHECK_F_CONT = False
+
+XDIR = 0
+"""label for x direction"""
+
+YDIR = 1
+"""label for y direction"""
+
+ZDIR = 2
+"""label for z direction"""
+
+PERIODIC = 99
+"""Tag for periodic boundary conditions"""
+
+S_DIR = ["_X", "_Y", "_Z"]
+"""Directions string"""
+
+WITH_GUESS = 1
+"""
+Optimisation level for time integrators :
+no need to recompute the right-hand side, an initial guess
+must be given in input arguments.
+ """
+
+NOALIAS = 2
+"""Optimisation level for time integrators :
+no need to recompute the right-hand side, an initial guess
+must be given in input arguments and we ensure that
+y is different from result arg.
+"""
+
+DEFAULT_TASK_ID = 999
+"""Default value for task id (mpi task)"""
+
+
+def _debugdecorator(f):
+ if __DEBUG__:
+ # function f is being decorated
+ def decorator(*args, **kw):
+ # Print informations on decorated function
+ if f.__name__ is '__new__':
+ fullclassname = args[0].__mro__[0].__module__ + '.'
+ fullclassname += args[0].__mro__[0].__name__
+ print ('Instanciate :', fullclassname,)
+ print (' (Inherits from : ',)
+ print ([c.__name__ for c in args[0].__mro__[1:]], ')')
+ else:
+ print ('Call : ', f.__name__, 'in ', f.__code__.co_filename,)
+ print (f.__code__.co_firstlineno)
+ # Calling f
+ r = f(*args, **kw)
+ if f.__name__ is '__new__':
+ print (' |-> : ', repr(r))
+ return r
+ return decorator
+ else:
+ #define empty debug decorator:
+ return f
+
+debug = _debugdecorator
+"""Debug decorator
+
+Usage: @debug before function definition you want to debug
+"""
+
+# redefine profile decorator
+if __PROFILE__:
+ from memory_profiler import profile
+ prof = profile
+ """Profiling decorator
+
+ Usage: @profile before function you want to profile
+ """
+
+else:
+ def prof(f):
+ # Nothing ...
+ return f
diff --git a/hysop/default_methods.py b/hysop/default_methods.py
index df41a328e367e5f79221d04f65ae1954c4646eb2..ad08c9a7434d4ce5a960de1101afa6b403dc74a8 100644
--- a/hysop/default_methods.py
+++ b/hysop/default_methods.py
@@ -9,8 +9,7 @@ from hysop.numerics.interpolation import Linear
from hysop.numerics.remeshing import L2_1, Linear as Rmsh_Linear
from hysop.numerics.finite_differences import FDC4, FDC2
from hysop.problem.simulation import O2
-
-#from hysop.operator.discrete.stretching import Conservative
+from hysop.operator.discrete.stretching import Conservative
ADVECTION = {TimeIntegrator: RK2, Interpolation: Linear,
@@ -44,7 +43,7 @@ POISSON = {SpaceDiscretisation: 'fftw', GhostUpdate: True,
"""Poisson operators default setup
"""
-STRETCHING = {TimeIntegrator: RK3, Formulation: "Conservative",
+STRETCHING = {TimeIntegrator: RK3, Formulation: Conservative,
SpaceDiscretisation: FDC4}
"""Stretching operators default setup
"""
diff --git a/hysop/f2hysop.pyf.in b/hysop/f2hysop.pyf.in
new file mode 100644
index 0000000000000000000000000000000000000000..a11194e7df7e5f7f02bb55b202ff01bf3faa85e9
--- /dev/null
+++ b/hysop/f2hysop.pyf.in
@@ -0,0 +1,20 @@
+! -*- f90 -*-
+
+! Generated file - Do not edit.
+
+! Note: the context of this file is case sensitive.
+
+python module f2hysop ! in
+
+ interface
+ ! Example
+ include '@CMAKE_SOURCE_DIR@/hysop/fortran/template.pyf'
+ ! precision
+ include '@CMAKE_SOURCE_DIR@/hysop/f2py/parameters.pyf'
+ ! fftw
+ include '@CMAKE_SOURCE_DIR@/hysop/f2py/fftw2py.pyf'
+ ! scales
+ include '@CMAKE_SOURCE_DIR@/hysop/f2py/scales2py.pyf'
+ end interface
+
+end python module f2hysop
\ No newline at end of file
diff --git a/hysop/methods.py b/hysop/methods.py
index 65e2e9b0d27c63fb24d3b907552c59348027e5c1..3437bc3c028e446e9998f3f2f82c91a593c8145d 100644
--- a/hysop/methods.py
+++ b/hysop/methods.py
@@ -23,6 +23,11 @@ that can be imported by final user.
"""
import hysop.numerics.odesolvers as odesolvers
+import hysop.numerics.interpolation as interpolation
+import hysop.numerics.remeshing as remesh
+import hysop.numerics.finite_differences as fd
+import hysop.operator.discrete.stretching as strd
+
RK2 = odesolvers.RK2
"""Runge-Kutta 2 integration"""
@@ -36,7 +41,6 @@ Euler = odesolvers.Euler
"""Euler integration"""
# Remesh
-import hysop.numerics.remeshing as remesh
L2_1 = remesh.L2_1
L2_2 = remesh.L2_2
L2_3 = remesh.L2_3
@@ -54,16 +58,13 @@ Rmsh_Linear = remesh.Linear
# A completer ...
# Interpolation
-import hysop.numerics.interpolation as interpolation
Linear = interpolation.Linear
# Finite differences
-import hysop.numerics.finite_differences as fd
FDC4 = fd.FDC4
FDC2 = fd.FDC2
FD2C2 = fd.FD2C2
# Stretching formulations
-import hysop.operator.discrete.stretching as strd
Conservative = strd.Conservative
GradUW = strd.GradUW
diff --git a/hysop/operator/absorption_bc.py b/hysop/operator/absorption_bc.py
index 8e9c46f9f07f7953bb8bbaedf69e3ba46cc2613d..7927106886d195c24400774e63981f62bc062eab 100755
--- a/hysop/operator/absorption_bc.py
+++ b/hysop/operator/absorption_bc.py
@@ -3,6 +3,8 @@
by vorticity absorption in order to set the far field
velocity to u_inf at the inlet)
+See :ref:`absorption_bc`.
+
"""
from hysop.constants import debug
from hysop.operator.discrete.absorption_bc import AbsorptionBC as Dab
@@ -119,3 +121,6 @@ class AbsorptionBC(Computational):
self._set_io('absorption_BC', (1, 2 + self.domain.dimension))
self.discrete_op.set_writer(self._writer)
self._is_uptodate = True
+
+ def wait(self):
+ print "TEMP WAIT FOR TEST"
diff --git a/hysop/operator/analytic.py b/hysop/operator/analytic.py
index 9e3cf0c179ac01ca4c7b59090299efe19597c9ff..46cbe93fbba2e327fe1bc49490d940536cca0956 100644
--- a/hysop/operator/analytic.py
+++ b/hysop/operator/analytic.py
@@ -1,6 +1,4 @@
-"""
-@file operator/analytic.py
-Initialize fields on a grid, with a user-defined function
+"""Initialize fields on a grid, with a user-defined function
"""
from hysop.constants import debug
from hysop.operator.continuous import opsetup, opapply
@@ -15,18 +13,27 @@ class Analytic(Computational):
@debug
def __init__(self, formula=None, vectorize_formula=False, **kwds):
- """
- Operator to apply a user-defined formula onto a list of fields.
- @param formula : the formula to be applied
- @param vectorize_formula : true if formula must be vectorized (numpy),
- default = false.
+ """ Apply a user-defined formula onto a list of fields.
+
+ Parameters
+ ----------
+ formula : python function
+ the formula to be applied
+ vectorize_formula : boolean, optional
+ true if formula must be vectorized (numpy), default = false.
+
+ Notes
+ -----
+ see :ref:`analytical_operator` for details on
+ the authorized signature for input formula or check
+ test_analytic.py
"""
super(Analytic, self).__init__(**kwds)
isGPU = False
if 'method' in kwds.keys() and Support in kwds['method'].keys():
isGPU = kwds['method'][Support].find('gpu') >= 0
if formula is not None:
- ## A formula applied to all variables of this operator
+ # A formula applied to all variables of this operator
self.formula = formula
for v in self.variables:
v.set_formula(formula, vectorize_formula)
diff --git a/hysop/operator/baroclinic.py b/hysop/operator/baroclinic.py
index 7937b55a558377a6420c858b106630bfb1aed6c5..83ff508e7f5e7d960a1b2026b398698135497955 100644
--- a/hysop/operator/baroclinic.py
+++ b/hysop/operator/baroclinic.py
@@ -75,6 +75,8 @@ class Baroclinic(Computational):
self._is_uptodate = True
-
def initialize_velocity(self):
self.discrete_op.initialize_velocity()
+
+ def get_work_properties(self):
+ return {'rwork': None, 'iwork': None}
diff --git a/hysop/operator/continuous.py b/hysop/operator/continuous.py
index 65623755a9f7b67f2493164d63cc20dba0f61754..055497d796bde3676e6885ca0317f47efafdd8cb 100755
--- a/hysop/operator/continuous.py
+++ b/hysop/operator/continuous.py
@@ -266,11 +266,8 @@ class Operator(object):
# default values for iop
self.io_params = IOParams(filename, fileformat=IO.ASCII)
else:
- if self.io_params.fileformat is not IO.ASCII:
- self.io_params.fileformat = IO.ASCII
- msg = 'Warning, wrong file format for operator output.'
- msg += 'This will be reset to ASCII.'
- print msg
+ msg = 'Error, wrong file format for operator output.'
+ assert self.io_params.fileformat is IO.ASCII, msg
self._writer = io.Writer(io_params=self.io_params,
mpi_params=self._mpis,
buffshape=buffshape)
@@ -322,8 +319,8 @@ def opapply(f):
"""decorate 'apply' method"""
# get 'continuous' base class and run its apply function
# --> call wait function of ops in wait_list
- t0 = ftime()
name = inspect.getmro(args[0].apply.im_class)
+ t0 = ftime()
name[-2].apply(args[0])
#t0 = ftime()
res = f(*args, **kwargs)
diff --git a/hysop/operator/custom.py b/hysop/operator/custom.py
index 5af161e895e8daab3d2fd387bbd7a37ce7a821d2..8fdc409f89646021533725e6b5474e7a23dad8bf 100644
--- a/hysop/operator/custom.py
+++ b/hysop/operator/custom.py
@@ -1,45 +1,65 @@
-"""
-@file custom.py
+"""Interface to set a user-defined operator
"""
from hysop.operator.computational import Computational
-from hysop.operator.discrete.custom import CustomMonitor as CM
-from hysop.operator.discrete.custom import CustomOp as CO
+from hysop.operator.discrete.custom import Custom as CO
from hysop.operator.continuous import opsetup
-class CustomOp(Computational):
- def __init__(self, in_fields, out_fields, function, **kwds):
- super(CustomOp, self).__init__(**kwds)
+class Custom(Computational):
+ """User-defined generic operator
+ """
+ def __init__(self, function, in_fields, out_fields=None,
+ diagnostics_shape=None, **kwds):
+ """
+
+ Parameters
+ ----------
+ in_fields: list of :class:`~hysop.fields.discrete.DiscreteField`
+ input fields args. for function, see notes below
+ out_fields: list of :class:`~hysop.fields.discrete.DiscreteField`
+ output fields args for function, see notes below
+ function: python function
+ a user defined function, called by this op.apply method.
+ diagnostics_shape: tuple, optional
+ shape of the data expected to be written into file.
+
+ Notes
+ -----
+ A function is used to set the behavior of the current operator,
+ during apply call.
+ This function must look like::
+
+ def some_func(simulation, in_fields, out_fields, diagnostics=None):
+ # do things ...
+
+ and compute out_fields values and optionnaly some diagnostics.
+ See :ref:`custom`
+ """
+ super(Custom, self).__init__(**kwds)
self.function = function
self.input = in_fields
- self.output = out_fields
+ if out_fields is not None:
+ self.output = out_fields
+ self._diagnostics_shape = diagnostics_shape
+
+ def discretize(self):
+ super(Custom, self)._standard_discretize()
@opsetup
def setup(self, rwork=None, iwork=None):
if not self._is_uptodate:
- self.discretize()
self.discrete_op = CO(
- [self.discreteFields[f] for f in self.input],
- [self.discreteFields[f] for f in self.output],
- self.function,
+ in_fields=[self.discreteFields[f] for f in self.input],
+ out_fields=[self.discreteFields[f] for f in self.output],
+ function=self.function,
variables=self.discreteFields.values())
- self._is_uptodate = True
+ if self._diagnostics_shape is not None:
+ assert isinstance(self._diagnostics_shape, tuple)
+ assert len(self._diagnostics_shape) == 2
+ self._set_io(self.function.__name__, self._diagnostics_shape)
+ self.discrete_op.set_writer(self._writer)
-
-class CustomMonitor(Computational):
- def __init__(self, function, res_shape=1, **kwds):
- super(CustomMonitor, self).__init__(**kwds)
- self.function = function
- self.res_shape = res_shape
- self.input = self.variables
-
- @opsetup
- def setup(self, rwork=None, iwork=None):
- if not self._is_uptodate:
- self.discretize()
- self.discrete_op = CM(self.function, self.res_shape,
- variables=self.discreteFields.values())
- # Output setup
- self._set_io(self.function.__name__, (1, 1 + self.res_shape))
- 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/discrete/__init__.py b/hysop/operator/discrete/__init__.py
index 70f15434bdf7ca63593af802292c2fbb54364f1c..e69de29bb2d1d6434b8b29ae775ad8c2e48c5391 100644
--- a/hysop/operator/discrete/__init__.py
+++ b/hysop/operator/discrete/__init__.py
@@ -1,24 +0,0 @@
-## @package hysop.operator.discrete
-# Discrete operators classes.
-#
-# A DiscreteOperator is an object that represents the discretisation of a
-# continuous operator for a specific given method and of its variables for
-# some resolutions.
-#
-# Example: if we want to perform the advection of a scalar at
-# velocity v using scales with M4 remesh, we define the following continuous
-# operator :
-#\code
-# nbElem = [65, 65, 65]
-# advec = Advection(velo, scal,
-# resolutions={velo: nbElem,
-# scal: nbElem},
-# method = 'scales, p_M4',)
-# ...
-# advec.setup()
-# ...
-# advec.apply()
-# \endcode
-#
-# setup call will result in the creation of a ScalesAdvection operator
-# and apply will perform a call to scale's solver.
diff --git a/hysop/operator/discrete/custom.py b/hysop/operator/discrete/custom.py
index 20b6aa4daf32ae48f394566a99509914fca27823..85929894c8827ed5a30e5fd182085b68abe2c715 100644
--- a/hysop/operator/discrete/custom.py
+++ b/hysop/operator/discrete/custom.py
@@ -1,28 +1,55 @@
+"""User-defined discrete operator"""
+
from hysop.operator.discrete.discrete import DiscreteOperator
-class CustomOp(DiscreteOperator):
- def __init__(self, in_fields, out_fields, function, **kwds):
- super(CustomOp, self).__init__(**kwds)
- self.function = function
- self._in_fields = in_fields
- self._out_fields = out_fields
+class Custom(DiscreteOperator):
+ """User-defined operator: action defined by an external function.
+ """
+ def __init__(self, function, in_fields, out_fields=None, **kwds):
+ """
- def apply(self, simulation):
- self.function(simulation, self._in_fields, self._out_fields)
+ Parameters
+ ----------
+ in_fields: list of :class:`~hysop.fields.discrete.DiscreteField`
+ input fields args. for function, see notes below
+ out_fields: list of :class:`~hysop.fields.discrete.DiscreteField`
+ output fields args for function, see notes below
+ function: python function
+ a user defined function, called by this op.apply method.
+ Notes
+ -----
+ A function is used to set the behavior of the current operator,
+ during apply call.
+ This function must look like::
-class CustomMonitor(DiscreteOperator):
- def __init__(self, function, res_shape=1, **kwds):
- super(CustomMonitor, self).__init__(**kwds)
+ def some_func(simulation, in_fields, out_fields=None, diag=None):
+ # do things ...
+
+ and compute out_fields values.
+ """
+ # callback for apply function
self.function = function
- self.res_shape = res_shape
+ super(Custom, self).__init__(**kwds)
+ # in/out fields must obviously belong to variables
+ self.input = in_fields
+ if out_fields is not None:
+ self.input += out_fields
+ self.output = out_fields
+ msg = 'Custom: all in/out fields must belong to op variables.'
+ assert set(self.input).intersection(set(self.variables)) == \
+ set(self.input), msg
+ self._in_fields = in_fields
+ self._out_fields = out_fields
def apply(self, simulation=None):
+ if self._writer is not None:
+ diagnostics = self._writer.buffer
+ else:
+ diagnostics = None
+ self.function(simulation, self._in_fields, self._out_fields,
+ diagnostics)
ite = simulation.current_iteration
- values = self.function(simulation,
- self.variables)
if self._writer is not None and self._writer.do_write(ite):
- self._writer.buffer[0, 0] = simulation.time
- self._writer.buffer[0, 1:] = values
self._writer.write()
diff --git a/hysop/operator/discrete/profiles.py b/hysop/operator/discrete/profiles.py
index 97e1a3682ab610f0e5e47577e42649eb3067b8a2..68373ef77a6365fc46cc7bf3366247cc953dba46 100644
--- a/hysop/operator/discrete/profiles.py
+++ b/hysop/operator/discrete/profiles.py
@@ -87,17 +87,6 @@ class Profiles(DiscreteOperator):
assert self._rwork[0].shape == shape_v
assert self._rwork[1].shape == shape_w
- def get_work_properties(self):
-
- v_ind = self.velocity.topology.mesh.compute_index
- w_ind = self.vorticity.topology.mesh.compute_index
- shape_v = self.velocity.data[0][v_ind].shape
- shape_w = self.velocity.data[0][w_ind].shape
- if shape_v == shape_w:
- return {'rwork': [shape_v], 'iwork': None}
- else:
- return {'rwork': [shape_v, shape_w], 'iwork': None}
-
@debug
@profile
def apply(self, simulation=None):
diff --git a/hysop/operator/discrete/stretching.py b/hysop/operator/discrete/stretching.py
index f35e2d07ba9d8d62abc9a869f9be13a68c7be3a8..4aebd947b92b94f861d514cbcb56c98116731de8 100755
--- a/hysop/operator/discrete/stretching.py
+++ b/hysop/operator/discrete/stretching.py
@@ -11,7 +11,7 @@ Formulations:
from hysop.constants import debug, WITH_GUESS
from hysop.methods_keys import TimeIntegrator, SpaceDiscretisation
from hysop.operator.discrete.discrete import DiscreteOperator
-import hysop.numerics.differential_operations as diff_op
+from hysop.numerics.differential_operations import GradVxW, DivWV
import hysop.tools.numpywrappers as npw
from hysop.numerics.update_ghosts import UpdateGhosts
from hysop.tools.profiler import profile
@@ -30,13 +30,19 @@ class Stretching(DiscreteOperator):
__metaclass__ = ABCMeta
@debug
- def __init__(self, velocity, vorticity, formulation, rhs, **kwds):
+ def __init__(self, velocity, vorticity, formulation,
+ rhs, extra_ghosts_fields=None, **kwds):
"""Abstract interface for stretching operator
Parameters
----------
velocity, vorticity : :class:`~hysop.fields.discrete.DiscreteField`
formulation : one of the discrete stretching classes.
+ rhs : python function
+ right-hand side for the time integrator
+ extra_ghosts_fields : int, optional
+ number of additional numpy arrays to be synchronized (mpi),
+ default = None, i.e. synchro of vorticity and velocity components.
**kwds : extra parameters for base class
"""
@@ -69,9 +75,12 @@ class Stretching(DiscreteOperator):
self.nb_components = 3 # Stretching only in 3D and for vector fields.
# prepare ghost points synchro for velocity and vorticity
+ nb_ghosts_fields = self.velocity.nb_components + \
+ self.vorticity.nb_components
+ if extra_ghosts_fields is not None:
+ nb_ghosts_fields += extra_ghosts_fields
self._synchronize = UpdateGhosts(self.velocity.topology,
- self.velocity.nb_components +
- self.vorticity.nb_components)
+ nb_ghosts_fields)
# A function to compute the gradient of the velocity.
# Work vector is provided in input.
@@ -161,8 +170,8 @@ class Conservative(Stretching):
"""
return self.strFunc(y, self.velocity.data, result)
- super(Conservative, self).__init__(formulation=diff_op.DivWV,
- rhs=rhs, **kwds)
+ super(Conservative, self).__init__(
+ formulation=DivWV, rhs=rhs, **kwds)
@profile
def _compute(self, dt, t):
@@ -183,8 +192,7 @@ class GradUW(Stretching):
self.strFunc(self.velocity.data, y, result, self.diagnostics)
return result
- super(GradUW, self).__init__(formulation=diff_op.GradVxW,
- rhs=rhs, **kwds)
+ super(GradUW, self).__init__(formulation=GradVxW, rhs=rhs, **kwds)
# stability constant
# Depends on time integration method
@@ -236,52 +244,61 @@ class StretchingLinearized(Stretching):
UpdateGhosts(self.vorticity_BF.topology,
self.vorticity_BF.nb_components)
self.usual_op = usual_op
+ # boolean used to switch between two rhs forms:
+ # either rhs(t, y) = div(y:u) (if true)
+ # or rhs(t,y) = div(y:w_bf) (if false)
+ self._divwu = True
- def rhs(t, y, result, form):
+ def rhs(t, y, result):
"""rhs used in time integrator
"""
- if form == "div(w:u)":
+ if self._divwu:
result = self.strFunc(y, self.velocity.data, result)
else:
result = self.strFunc(y, self.vorticity_BF.data, result)
return result
- super(StretchingLinearized, self).__init__(formulation=diff_op.DivWV,
- rhs=rhs, **kwds)
+ super(StretchingLinearized, self).__init__(
+ formulation=DivWV, rhs=rhs,
+ extra_ghosts_fields=self.vorticity_BF.nb_components,
+ **kwds)
- def _integrate(self, dt, t):
+ def _compute(self, dt, t):
# - Call time integrator (1st term over 3) -
# Init workspace with a first evaluation of the div(wb:u') term in the
# rhs of the integrator
+ self._divwu = True
self._ti_work[:self.nb_components] = \
self.timeIntegrator.f(t, self.vorticity_BF.data,
- self._ti_work[:self.nb_components],
- "div(w:u)")
+ self._ti_work[:self.nb_components])
# perform integration and save result in-place
self.vorticity.data = self.timeIntegrator(t, self.vorticity.data, dt,
result=self.vorticity.data)
# - Call time integrator (2nd term over 3) -
# Init workspace with a first evaluation of the div(u':wb) term in the
# rhs of the integrator
+ self._divwu = False
self._ti_work[:self.nb_components] = \
self.timeIntegrator.f(t, self.velocity.data,
- self._ti_work[:self.nb_components],
- "div(u:w)")
+ self._ti_work[:self.nb_components])
# perform integration and save result in-place
self.vorticity.data = self.timeIntegrator(t, self.vorticity.data, dt,
result=self.vorticity.data)
- def _compute(self, dt, t):
- # No subcycling for this formulation
- self._integrate(dt, t)
+ def apply(self, simulation=None):
+ assert simulation is not None, \
+ "Missing simulation value for computation."
+
+ # time step
+ dt = simulation.time_step
+ # current time
+ t = simulation.time
- def _apply(self, t, dt):
# Synchronize ghost points
- self._synchronize(self.velocity.data + self.vorticity.data)
- self._synchronize_vort_BF(self.vorticity_BF.data)
+ self._synchronize(self.velocity.data + self.vorticity.data +
+ self.vorticity_BF.data)
# Compute the 2 first "stretching" terms (div(wb:u') and div(u':wb))
# and update vorticity for each of them
self._compute(dt, t)
# Compute the 3rd stretching term (div(w':ub)) and update vorticity
- self.usual_op._apply(t, dt)
-
+ self.usual_op.apply(simulation)
diff --git a/hysop/operator/profiles.py b/hysop/operator/profiles.py
index 3f74bbb2949f5415c8120e1b016876df0470421d..3ce3cd9ce21cbb480712a88d278cd54189a0bc1e 100644
--- a/hysop/operator/profiles.py
+++ b/hysop/operator/profiles.py
@@ -1,7 +1,4 @@
-# -*- coding: utf-8 -*-
-"""
-@file profiles.py
-Compute and print velo/vorti profiles
+"""Compute and print velo/vorti profiles
"""
from hysop.operator.discrete.profiles import Profiles as ProfD
from hysop.operator.computational import Computational
@@ -13,7 +10,7 @@ class Profiles(Computational):
Compute and print velo/vorti profiles
"""
- def __init__(self, velocity, vorticity, prof_coords,
+ def __init__(self, velocity, vorticity, prof_coords,
direction, beginMeanComput, **kwds):
"""
Constructor.
@@ -30,24 +27,21 @@ class Profiles(Computational):
assert 'variables' not in kwds, 'variables parameter is useless.'
super(Profiles, self).__init__(variables=[velocity, vorticity],
**kwds)
- ## velocity field
+ # velocity field
self.velocity = velocity
- ## vorticity field
+ # vorticity field
self.vorticity = vorticity
- ## X and Y coordinates of the profile
+ # X and Y coordinates of the profile
self.prof_coords = prof_coords
- ## profile direction (0, 1 or 2)
+ # profile direction (0, 1 or 2)
self.direction = direction
- ## time at which the computation of mean profile must begin
+ # time at which the computation of mean profile must begin
self.beginMeanComput = beginMeanComput
self.input = [velocity, vorticity]
self.output = []
def get_work_properties(self):
- if not self._is_discretized:
- msg = 'The operator must be discretized '
- msg += 'before any call to this function.'
- raise RuntimeError(msg)
+ super(Profiles, self).get_work_properties()
vd = self.discreteFields[self.velocity]
wd = self.discreteFields[self.vorticity]
v_ind = vd.topology.mesh.compute_index
@@ -72,4 +66,3 @@ class Profiles(Computational):
self._set_io('profile', (1, 9))
self.discrete_op.set_writer(self._writer)
self._is_uptodate = True
-
diff --git a/hysop/operator/stretching.py b/hysop/operator/stretching.py
index c0e242cbd3d7fa48b379448e8aaf12359528c65b..80095bb777ec5d782f36fa38fa702c34d2ff6ab1 100755
--- a/hysop/operator/stretching.py
+++ b/hysop/operator/stretching.py
@@ -16,7 +16,7 @@ from hysop.operator.computational import Computational
from hysop.operator.continuous import opsetup
from hysop.operator.discrete.stretching import Conservative, GradUW
from hysop.operator.discrete.stretching import StretchingLinearized as SLD
-import hysop.numerics.differential_operations as diff_op
+from hysop.numerics.differential_operations import GradVxW, DivWV
from hysop.numerics.odesolvers import Euler
@@ -25,6 +25,7 @@ class Stretching(Computational):
"""
_authorized_methods = [FDC4]
+ _authorized_formulations = [Conservative, GradUW]
@debug
def __init__(self, velocity, vorticity, **kwds):
@@ -61,11 +62,9 @@ class Stretching(Computational):
msg
# Formulation used for the stretching equation.
# Default = conservative form.
- if self.method[Formulation] == "GradUW":
- self.formulation = GradUW
- else:
- self.formulation = Conservative
-
+ self.formulation = self.method[Formulation]
+ msg = 'Stretching error : unknown formulation.'
+ assert self.formulation in self._authorized_formulations, msg
self.input = [self.velocity, self.vorticity]
self.output = [self.vorticity]
@@ -81,9 +80,10 @@ class Stretching(Computational):
res = ti.get_work_properties(nbc, topo)
# ---> differential operator work space
if self.formulation is GradUW:
- dop = diff_op.GradVxW
+ dop = GradVxW
+
elif self.formulation is Conservative:
- dop = diff_op.DivWV
+ dop = DivWV
res['rwork'] += dop.get_work_properties(topo)['rwork']
return res
@@ -103,12 +103,10 @@ class Stretching(Computational):
class StretchingLinearized(Stretching):
- """Solve the linearized stretching equation, i.e:
+ """Solve the linearized stretching equation
+
+ See details in :ref:`stretching`.
- \f{eqnarray*}
- \frac{\partial \omega}{\partial t} &=& (\om \cdot \nabla)u_b +
- (\om_b \cdot \nabla)u
- \f}
"""
@debug
def __init__(self, velocity_BF, vorticity_BF, **kwds):
diff --git a/hysop/operator/tests/test_absorption_bc.py b/hysop/operator/tests/test_absorption_bc.py
index 173e207043e400d21257c10e24b1f20f7d3434d6..8f97dd6c244ee6f982e361783404a2fa5724e12f 100644
--- a/hysop/operator/tests/test_absorption_bc.py
+++ b/hysop/operator/tests/test_absorption_bc.py
@@ -28,7 +28,7 @@ def init_test(dim, discr, ext_filter, filter_func=None):
flowrate = VariableParameter(data=uinf)
x_range = [0.9, 1.]
op = AbsorptionBC(v, w, flowrate, x_range,
- discretization=d3d, filter_func=filter_func)
+ discretization=discr, filter_func=filter_func)
op.discretize()
op.setup()
wk_prop = op.get_work_properties()['rwork']
diff --git a/hysop/operator/tests/test_custom.py b/hysop/operator/tests/test_custom.py
new file mode 100644
index 0000000000000000000000000000000000000000..6f92c54ff6500b5c80c762ad72d6ef948ea10f1b
--- /dev/null
+++ b/hysop/operator/tests/test_custom.py
@@ -0,0 +1,84 @@
+"""Test user-defined (custom) operator
+"""
+
+from hysop.operator.custom import Custom
+from hysop import Field, Discretization, Box, Simulation, IOParams, IO
+import numpy as np
+import os
+
+d2d = Discretization([33, 65], [2, 0])
+d3d = Discretization([33, 65, 33], [2, 0, 1])
+
+
+def func3d(sim, f_in, f_out, d=None):
+ """A test function, for a 3d domain,
+ no diagnostics
+ """
+ nbc = f_in[0].nb_components
+ for d in xrange(nbc):
+ f_out[0][d] = d * f_in[0][d] + + f_in[1][0] + np.cos(sim.time)
+
+
+def func3d_with_diag(sim, f_in, f_out, diagnostics):
+ """A test function, for a 3d domain,
+ no diagnostics
+ """
+ nbc = f_in[0].nb_components
+ for d in xrange(nbc):
+ f_out[0][d] = d * f_in[0][d] + f_in[1][0] + np.cos(sim.time)
+ diagnostics[0, 0] = sim.time
+ diagnostics[0, 1] = f_out[0].data[0].min()
+ diagnostics[0, 1] = f_out[0].data[0].max()
+
+
+def init_custom(dim, discr, func, do_write=False):
+ """Test build and apply, without
+ writer.
+ """
+ dom = Box(length=[1., ] * dim)
+ v1 = Field(name='v1', is_vector=True, domain=dom)
+ v2 = Field(name='v2', is_vector=False, domain=dom)
+ v3 = Field(name='v3', is_vector=True, domain=dom)
+ topo = dom.create_topology(discr)
+ v1.randomize(topo)
+ v2.randomize(topo)
+ if do_write:
+ iop = IOParams('/tmp/hysop_custom_test.dat', fileformat=IO.ASCII)
+ d_shape = (1, 3)
+ else:
+ iop = None
+ d_shape = None
+ op = Custom(in_fields=[v1, v2], out_fields=[v3],
+ variables=[v1, v2, v3], discretization=discr,
+ function=func, io_params=iop,
+ diagnostics_shape=d_shape)
+ op.discretize()
+ op.setup()
+ sim = Simulation(nb_iter=10)
+ sim.initialize()
+ nbc = v1.nb_components
+ vd1 = v1.discretize(topo)
+ vd2 = v2.discretize(topo)
+ vd3 = v3.discretize(topo)
+ while not sim.is_over:
+ op.apply(sim)
+ for d in xrange(nbc):
+ tmp1 = d * vd1[d] + + vd2[0] + np.cos(sim.time)
+ assert np.allclose(vd3[d], tmp1)
+ sim.advance()
+ if do_write:
+ assert os.path.isfile(iop.filename)
+
+
+def test_custom_op_1():
+ """Test build and apply, without
+ writer.
+ """
+ init_custom(3, d3d, func3d)
+
+
+def test_custom_op_2():
+ """Test build and apply, with
+ output in a file
+ """
+ init_custom(3, d3d, func3d_with_diag, True)
diff --git a/hysop/operator/tests/test_diff_poisson_3D.py b/hysop/operator/tests/test_diff_poisson_3D.py
index 9dbea2debb24c60df61232202dd76a34a856b649..0268a81acc2bac6be3205253dcccdc8364a51210 100755
--- a/hysop/operator/tests/test_diff_poisson_3D.py
+++ b/hysop/operator/tests/test_diff_poisson_3D.py
@@ -4,6 +4,7 @@ from hysop.operator.poisson import Poisson
from hysop.operator.diffusion import Diffusion
from math import sqrt, pi, exp
from hysop.problem.simulation import Simulation
+from hysop import testsenv
def computeVel(x, y, z):
@@ -36,6 +37,7 @@ def computeVort(x, y, z):
return wx, wy, wz
+@testsenv.fftw_failed
def test_Diff_Poisson():
# Parameters
nb = 33
@@ -44,16 +46,16 @@ def test_Diff_Poisson():
from hysop.tools.parameters import Discretization
d3D = Discretization([nb, nb, nb])
- ## Domain
+ # Domain
box = pp.Box(length=boxLength, origin=boxMin)
- ## Fields
+ # Fields
velo = pp.Field(domain=box, formula=computeVel,
name='Velocity', is_vector=True)
vorti = pp.Field(domain=box, formula=computeVort,
name='Vorticity', is_vector=True)
- ## FFT Diffusion operators and FFT Poisson solver
+ # FFT Diffusion operators and FFT Poisson solver
diffusion = Diffusion(variables={vorti: d3D}, viscosity=0.002)
poisson = Poisson(velo, vorti, discretization=d3D)
diff --git a/hysop/operator/tests/test_differential.py b/hysop/operator/tests/test_differential.py
index 9eaa3895b537fd371434e01bcad3d3fd734697d3..2a9c983dcea199c5355c746ba0cd6ace9d16d279 100755
--- a/hysop/operator/tests/test_differential.py
+++ b/hysop/operator/tests/test_differential.py
@@ -9,12 +9,11 @@ from hysop.methods_keys import SpaceDiscretisation
from hysop.methods import FDC4, FDC2
from hysop.operator.differential import Curl, Grad, DivAdvection
from hysop.tools.parameters import Discretization
-
+from hysop import testsenv
# Domain and topologies definitions
-
-nb = 65
import math
+nb = 65
Lx = Ly = Lz = 2. * math.pi
box1_3d = Box(length=[Lx, Ly, Lz], origin=[0., 0., 0.])
box1_2d = Box(length=[Lx, Ly], origin=[0., 0.])
@@ -91,6 +90,8 @@ def grad_velo_2d(res, x, y, t):
def check(op, ref_formula, topo, op_dim=3, order=4):
+ """Apply operator 'op' and compare its results with some references.
+ """
# Reference field
ref = Field(domain=topo.domain, formula=ref_formula, nb_components=op_dim,
name='reference')
@@ -116,6 +117,10 @@ def check(op, ref_formula, topo, op_dim=3, order=4):
def call_op(class_name, ref_formula, topo, use_work=False,
op_dim=3, method=None, order=4, vform=velocity_f):
+ """init and set an operator of type 'class_name'
+ and call check function on this operator.
+ """
+
# Velocity and result fields
velo = Field(domain=topo.domain, formula=vform, is_vector=True,
name='velo')
@@ -139,6 +144,10 @@ def call_op(class_name, ref_formula, topo, use_work=False,
def call_op_fft(class_name, ref_formula, dom, discr,
op_dim=3, method=None, order=4, vform=velocity_f):
+ """init and set an operator of type 'class_name'
+ and call check function on this operator. FFTW method.
+ """
+
# Velocity and result fields
velo = Field(domain=dom, formula=vform, is_vector=True,
name='velo')
@@ -180,6 +189,12 @@ def test_curl_fd_3_2d():
method=method, vform=velocity_f2d)
+def test_curl_fd_work():
+ method = {SpaceDiscretisation: FDC4}
+ call_op(Curl, vorticity_f, topo3_3d, method=method, use_work=True)
+
+
+@testsenv.fftw_failed
def test_curl_fft_1():
method = {SpaceDiscretisation: 'fftw'}
d1_3d_nog = Discretization([nb, nb, nb])
@@ -228,17 +243,17 @@ def test_grad_3():
call_op(Grad, grad_velo, topo3_3d, op_dim=9, method=method, order=2)
+def test_grad_3_work():
+ method = {SpaceDiscretisation: FDC2}
+ call_op(Grad, grad_velo, topo3_3d, op_dim=9, method=method, order=2,
+ use_work=True)
+
+
def test_grad_3_2d():
method = {SpaceDiscretisation: FDC4}
call_op(Grad, grad_velo_2d, topo3_2d, op_dim=4, method=method,
vform=velocity_f2d)
-
-def test_curl_fd_work():
- method = {SpaceDiscretisation: FDC4}
- call_op(Curl, vorticity_f, topo3_3d, use_work=True, method=method)
-
-
def divadvection_func(res, x, y, z, t):
res[0][...] = - cos(z) * cos(z) * (cos(x) * cos(x) - sin(x) * sin(x)) - \
cos(z) * cos(z) * (cos(y) * cos(y) - sin(y) * sin(y))
diff --git a/hysop/operator/tests/test_diffusion.py b/hysop/operator/tests/test_diffusion.py
index 6473a8b5d9ba226f8e7ba1717ff72fa2cf2335a6..cc3e80f86a461c55a3e2b4bf2a8f3df1681a005a 100755
--- a/hysop/operator/tests/test_diffusion.py
+++ b/hysop/operator/tests/test_diffusion.py
@@ -7,10 +7,11 @@ from hysop.tools.parameters import Discretization
import numpy as np
import hysop.tools.numpywrappers as npw
import math
+from hysop import testsenv
pi = math.pi
sin = np.sin
cos = np.cos
-## Physical Domain description
+# Physical Domain description
dim = 3
LL = 2 * pi * npw.ones((dim))
cc = 2 * pi / LL
@@ -32,6 +33,7 @@ def computeVort2D(res, x, y, t):
return res
+@testsenv.fftw_failed
def test_Diffusion3D():
dom = pp.Box(length=LL)
@@ -50,6 +52,7 @@ def test_Diffusion3D():
diff.finalize()
+@testsenv.fftw_failed
def test_Diffusion3D_2():
dom = pp.Box(length=LL)
@@ -68,6 +71,7 @@ def test_Diffusion3D_2():
diff.finalize()
+@testsenv.fftw_failed
def test_Diffusion2D():
dom = pp.Box(length=LL[:2])
diff --git a/hysop/operator/tests/test_operators.py b/hysop/operator/tests/test_operators.py
index 1b184f29bd99d3faa2e2df79e48e613ebdbdb16e..f49e10dc2c667895b52c7546330781d4cb6bda12 100644
--- a/hysop/operator/tests/test_operators.py
+++ b/hysop/operator/tests/test_operators.py
@@ -1,7 +1,5 @@
# -*- coding: utf-8 -*-
-"""
-@file test_operators.py
-tests for operators general interface
+"""tests for operators general interface
"""
from hysop.mpi.tests.utils import create_multitask_context, CPU, GPU, OTHER
diff --git a/hysop/operator/tests/test_Stretching.py b/hysop/operator/tests/test_stretching.py
similarity index 65%
rename from hysop/operator/tests/test_Stretching.py
rename to hysop/operator/tests/test_stretching.py
index 3a0f89b2fed9091b517e316a9d4d795279cbb653..97dae7246c14c4e86a8bd398d49e3707760d3fdf 100755
--- a/hysop/operator/tests/test_Stretching.py
+++ b/hysop/operator/tests/test_stretching.py
@@ -3,7 +3,7 @@
"""
from hysop import Field, Box
import numpy as np
-from hysop.operator.stretching import Stretching, StretchingLinearized
+from hysop.operator.stretching import Stretching, StretchingLinearized, GradUW
from hysop.problem.simulation import Simulation
from hysop.methods_keys import TimeIntegrator, Formulation,\
SpaceDiscretisation
@@ -15,46 +15,13 @@ cos = np.cos
sin = np.sin
-def compute_vel(res, x, y, z, t):
- amodul = cos(pi * 1. / 3.)
- pix = pi * x
- piy = pi * y
- piz = pi * z
- pi2x = 2. * pix
- pi2y = 2. * piy
- pi2z = 2. * piz
- res[0][...] = 2. * sin(pix) * sin(pix) \
- * sin(pi2y) * sin(pi2z) * amodul
- res[1][...] = - sin(pi2x) * sin(piy) \
- * sin(piy) * sin(pi2z) * amodul
- res[2][...] = - sin(pi2x) * sin(piz) \
- * sin(piz) * sin(pi2y) * amodul
- return res
-
-
-def compute_vort(res, x, y, z, t):
- amodul = cos(pi * 1. / 3.)
- pix = pi * x
- piy = pi * y
- piz = pi * z
- pi2x = 2. * pix
- pi2y = 2. * piy
- pi2z = 2. * piz
- res[0][...] = 2. * pi * sin(pi2x) * amodul *\
- (- cos(pi2y) * sin(piz) * sin(piz) +
- sin(piy) * sin(piy) * cos(pi2z))
-
- res[1][...] = 2. * pi * sin(pi2y) * amodul *\
- (2. * cos(pi2z) * sin(pix) * sin(pix) +
- sin(piz) * sin(piz) * cos(pi2x))
-
- res[2][...] = -2. * pi * sin(pi2z) * amodul *\
- (cos(pi2x) * sin(piy) * sin(piy) +
- sin(pix) * sin(pix) * cos(pi2y))
+# 3d discretisation with size 2 ghosts layer
+d3d = Discretization([65, ] * 3, [2, ] * 3)
- return res
-def computeVelBF(res, x, y, z, t):
+def compute_vel(res, x, y, z, t):
+ """3d vector field
+ """
amodul = cos(pi * 1. / 3.)
pix = pi * x
piy = pi * y
@@ -71,7 +38,9 @@ def computeVelBF(res, x, y, z, t):
return res
-def computeVortBF(res, x, y, z, t):
+def compute_vort(res, x, y, z, t):
+ """3d vector field
+ """
amodul = cos(pi * 1. / 3.)
pix = pi * x
piy = pi * y
@@ -97,13 +66,8 @@ def computeVortBF(res, x, y, z, t):
def init(method=None, work=False):
"""Build, init, setup for operator
"""
- nb = 33
- box_length = [1., 1., 1.]
- box_min = [0., 0., 0.]
- nb_elem = Discretization([nb, nb, nb], [2, 2, 2])
-
# Domain
- box = Box(length=box_length, origin=box_min)
+ box = Box()
# Fields
velo = Field(
@@ -112,17 +76,22 @@ def init(method=None, work=False):
vorti = Field(
domain=box, formula=compute_vort,
name='Vorticity', is_vector=True)
- op = Stretching(velo, vorti, discretization=nb_elem, method=method)
+ # Stretching operator
+ op = Stretching(velo, vorti, discretization=d3d, method=method)
op.discretize()
rwork = None
if work:
+ # Find required size for internal work vector, if needed
wkp = op.get_work_properties()['rwork']
+ # Allocate work space
rwork = WorkSpaceTools.check_work_array(len(wkp), wkp[0])
topo = op.discreteFields[velo].topology
- velo.initialize(topo=topo)
- vorti.initialize(topo=topo)
op.setup(rwork=rwork)
simulation = Simulation(start=0, end=1., time_step=0.05)
+ # initialize fields and simu
+ velo.initialize(topo=topo)
+ vorti.initialize(topo=topo)
+ simulation.initialize()
return op, simulation
@@ -136,7 +105,7 @@ def test_stretching():
def test_stretching_external_work():
- """User-defined work arrays
+ """Default setup but with user-defined work arrays
"""
op, simu = init(work=True)
op.apply(simu)
@@ -145,43 +114,58 @@ def test_stretching_external_work():
def test_stretching_graduw():
"""GradUW formulation
"""
- method = {TimeIntegrator: RK3, Formulation: "GradUW",
+ method = {TimeIntegrator: RK3, Formulation: GradUW,
SpaceDiscretisation: FDC4}
op, simu = init(method=method)
op.apply(simu)
+# def test_compare_stretching():
+# """Run conservative and graduv form,
+# check if results are close enough
+# """
+# op1, simu = init()
+# op1.apply(simu)
+# method = {TimeIntegrator: RK3, Formulation: GradUW,
+# SpaceDiscretisation: FDC4}
+# op2, simu = init(method=method)
+# op2.apply(simu)
+# w1 = op1.discrete_op.vorticity
+# w2 = op2.discrete_op.vorticity
+# h = w1.topology.mesh.space_step.max()
+# for d in xrange(op1.domain.dimension):
+# print np.abs((w1[d] - w2[d])).max()
+# assert np.allclose(w1[d], w2[d], atol=h ** 2)
+
+
def init_linearized(method=None, work=False):
- """Build, init, setup for operator
+ """Build, init, setup for linearized stretching
"""
- nb = 33
- box_length = [1., 1., 1.]
- box_min = [0., 0., 0.]
- nb_elem = Discretization([nb, nb, nb], [2, 2, 2])
-
# Domain
- box = Box(length=box_length, origin=box_min)
+ box = Box()
- # Fields
+ # Base fields:
+ velo_bf = Field(
+ domain=box, formula=compute_vel,
+ name='VelocityBF', is_vector=True)
+ vorti_bf = Field(
+ domain=box, formula=compute_vort,
+ name='VorticityBF', is_vector=True)
+ # Perturbations
velo = Field(
domain=box, formula=compute_vel,
name='Velocity', is_vector=True)
vorti = Field(
domain=box, formula=compute_vort,
name='Vorticity', is_vector=True)
- velo_bf = Field(
- domain=box, formula=computeVelBF,
- name='VelocityBF', is_vector=True)
- vorti_bf = Field(
- domain=box, formula=computeVortBF,
- name='VorticityBF', is_vector=True)
# Usual stretching operator
- stretch1 = Stretching(velo, vorti, discretization=nb_elem)
+ stretch1 = Stretching(velo, vorti, discretization=d3d)
+ # Linearized stretching
stretch2 = StretchingLinearized(velocity=velo, vorticity=vorti,
- velocity_BF=velo_bf,
- vorticity_BF=vorti_bf,
- discretization=nb_elem, method=method)
+ velocity_BF=velo_bf,
+ vorticity_BF=vorti_bf,
+ discretization=d3d, method=method)
stretch1.discretize()
stretch2.discretize()
rwork = None
@@ -189,6 +173,7 @@ def init_linearized(method=None, work=False):
wkp = stretch2.get_work_properties()['rwork']
rwork = WorkSpaceTools.check_work_array(len(wkp), wkp[0])
topo = stretch1.discreteFields[velo].topology
+ # initialize all fields
velo.initialize(topo=topo)
vorti.initialize(topo=topo)
velo_bf.initialize(topo=topo)
@@ -196,6 +181,7 @@ def init_linearized(method=None, work=False):
stretch1.setup()
stretch2.setup(rwork=rwork)
simulation = Simulation(start=0, end=1., time_step=0.05)
+ simulation.initialize()
return stretch1, stretch2, simulation
@@ -210,13 +196,7 @@ def test_stretching_linearized():
def test_stretching_external_work_graduv():
"""User-defined work arrays for GradUW formulation
"""
- method = {TimeIntegrator: RK3, Formulation: "GradUW",
+ method = {TimeIntegrator: RK3, Formulation: GradUW,
SpaceDiscretisation: FDC4}
op, simu = init(work=True, method=method)
op.apply(simu)
-
-if __name__ == "__main__":
- test_stretching()
- test_stretching_external_work()
- test_stretching_graduw()
- test_stretching_external_work_graduv()
diff --git a/trashed_examples/LevelSet3D/levelSet3D.py b/trashed_examples/LevelSet3D/levelSet3D.py
index 2ee37c14589160681b89cace46f1a8f2b7c050b7..ef4b982ed43346ca34f9816921c242ebeb1750cb 100644
--- a/trashed_examples/LevelSet3D/levelSet3D.py
+++ b/trashed_examples/LevelSet3D/levelSet3D.py
@@ -12,7 +12,7 @@ from hysop.numerics.remeshing import L6_4 as remesh_formula
from hysop.numerics.remeshing import Linear
from hysop.operator.advection import Advection
from hysop.operator.analytic import Analytic
-from hysop.operator.custom import CustomMonitor
+from hysop.operators import Custom
from hysop.operator.hdf_io import HDF_Writer
from hysop.gpu.gpu_transfer import DataTransfer
from hysop.problem.problem import Problem
@@ -34,7 +34,7 @@ def vitesse(res, x, y, z, t=0.):
return res
-def volume(_simu, var):
+def volume(_simu, var, **kwds):
v_loc = np.sum(var[0].data[0] > 0.5) * np.prod(
var[0].topology.mesh.space_step)
return main_comm.allreduce(sendobj=v_loc, op=MPI.SUM)
@@ -96,9 +96,9 @@ topo_s = advec.advec_dir[0].discreteFields[scal].topology
velocity = Analytic(variables={velo: topo_v},
method={Support: 'gpu'})
-volume_m = CustomMonitor(function=volume, res_shape=1,
- variables={scal: topo_s},
- io_params=IOParams(filename="volume.dat"))
+volume_m = Custom(function=volume, diagnostics_shape=(1, 2),
+ in_fields=[scal], variables={scal: topo_s},
+ io_params=IOParams(filename="volume.dat"))
p = HDF_Writer(variables={scal: topo_s},
io_params=IOParams(frequency=outputModulo,
filename="levelset",
diff --git a/trashed_examples/LevelSet3D/levelSet3D_Scales.py b/trashed_examples/LevelSet3D/levelSet3D_Scales.py
index 251b45a43f73837ee812b1ddaaffdb9b5e5f4ff1..4572efca65c05ea6a821d5ec7c7fc02ed89eef70 100644
--- a/trashed_examples/LevelSet3D/levelSet3D_Scales.py
+++ b/trashed_examples/LevelSet3D/levelSet3D_Scales.py
@@ -6,7 +6,7 @@ from hysop.tools.parameters import Discretization, IOParams
from hysop.methods_keys import Scales, MultiScale
from hysop.operator.advection import Advection
from hysop.operator.analytic import Analytic
-from hysop.operator.custom import CustomMonitor
+from hysop.operators import Custom
from hysop.operator.hdf_io import HDF_Writer
from hysop.problem.problem import Problem
from hysop.problem.simulation import Simulation
@@ -27,7 +27,7 @@ def vitesse(res, x, y, z, t=0.):
return res
-def volume(_simu, var):
+def volume(_simu, var, **kwds):
v_loc = np.sum(var[0].data[0] > 0.5) * np.prod(
var[0].topology.mesh.space_step)
return main_comm.allreduce(sendobj=v_loc, op=MPI.SUM)
@@ -68,9 +68,9 @@ topo_s = advec.discreteFields[scal].topology
velocity = Analytic(variables={velo: topo_v})
-volume_m = CustomMonitor(function=volume, res_shape=1,
- variables={scal: topo_s},
- io_params=IOParams(filename="volume.dat"))
+volume_m = Custom(function=volume, res_shape=1,
+ in_fields=[scal], discretization=topo_s,
+ io_params=IOParams(filename="volume.dat"))
p = HDF_Writer(variables={scal: topo_s},
io_params=IOParams(frequency=outputModulo,
filename="levelset",
diff --git a/trashed_examples/LevelSet3D/levelSet3D_Scales_MultiScale.py b/trashed_examples/LevelSet3D/levelSet3D_Scales_MultiScale.py
index b50c7baed82c2fc85b038266d226360af67f4ae3..8b970f960006b1dade145fbb5d8d4bb85a06f456 100644
--- a/trashed_examples/LevelSet3D/levelSet3D_Scales_MultiScale.py
+++ b/trashed_examples/LevelSet3D/levelSet3D_Scales_MultiScale.py
@@ -6,7 +6,7 @@ from hysop.tools.parameters import Discretization, IOParams
from hysop.methods_keys import Scales, MultiScale
from hysop.operator.advection import Advection
from hysop.operator.analytic import Analytic
-from hysop.operator.custom import CustomMonitor
+from hysop.operators import Custom
from hysop.operator.hdf_io import HDF_Writer
from hysop.problem.problem import Problem
from hysop.problem.simulation import Simulation
@@ -27,7 +27,7 @@ def vitesse(res, x, y, z, t=0.):
return res
-def volume(_simu, var):
+def volume(_simu, var, **kwds):
v_loc = np.sum(var[0].data[0] > 0.5) * np.prod(
var[0].topology.mesh.space_step)
return main_comm.allreduce(sendobj=v_loc, op=MPI.SUM)
@@ -69,9 +69,9 @@ advec.discretize()
velocity = Analytic(variables={velo: topo_v})
-volume_m = CustomMonitor(function=volume, res_shape=1,
- variables={scal: topo_s},
- io_params=IOParams(filename="volume.dat"))
+volume_m = Custom(function=volume, diagnostics_shape=(1, 2),
+ in_fields=[scal], variables={scal: topo_s},
+ io_params=IOParams(filename="volume.dat"))
p = HDF_Writer(variables={scal: topo_s},
io_params=IOParams(frequency=outputModulo,
filename="levelset",
diff --git a/trashed_examples/LevelSet3D/levelSet3D_gpu.py b/trashed_examples/LevelSet3D/levelSet3D_gpu.py
index b12d157daac86fcc9af693a5d49003d7552b6fb3..e51895da6366910a6085b1fd090d6a696582473c 100755
--- a/trashed_examples/LevelSet3D/levelSet3D_gpu.py
+++ b/trashed_examples/LevelSet3D/levelSet3D_gpu.py
@@ -12,7 +12,7 @@ from hysop.numerics.remeshing import L6_4 as remesh_formula
from hysop.numerics.remeshing import Linear
from hysop.operator.advection import Advection
from hysop.operator.analytic import Analytic
-from hysop.operator.custom import CustomMonitor
+from hysop.operators import Custom
from hysop.operator.hdf_io import HDF_Writer
from hysop.gpu.gpu_transfer import DataTransfer
from hysop.problem.problem import Problem
@@ -34,7 +34,7 @@ def vitesse(res, x, y, z, t=0.):
return res
-def volume(_simu, var):
+def volume(_simu, var, **kwds):
v_loc = np.sum(var[0].data[0] > 0.5) * np.prod(
var[0].topology.mesh.space_step)
return main_comm.allreduce(sendobj=v_loc, op=MPI.SUM)
@@ -88,9 +88,9 @@ topo_s = advec.advec_dir[0].discreteFields[scal].topology
velocity = Analytic(variables={velo: topo_v})
-volume_m = CustomMonitor(function=volume, res_shape=1,
- variables={scal: topo_s},
- io_params=IOParams(filename="volume.dat"))
+volume_m = Custom(function=volume, diagnostics_shape=1,
+ in_fields=[scal], variables={scal: topo_s},
+ io_params=IOParams(filename="volume.dat"))
p = HDF_Writer(variables={scal: topo_s},
io_params=IOParams(frequency=outputModulo,
filename="levelset",
diff --git a/trashed_examples/LevelSet3D/levelSet3D_gpu_MultiScale.py b/trashed_examples/LevelSet3D/levelSet3D_gpu_MultiScale.py
index 2bfef8a41982257b679f9180b8b25a79b9b833bb..3d05e5019213b42d16b5135878ac931f47d1840d 100644
--- a/trashed_examples/LevelSet3D/levelSet3D_gpu_MultiScale.py
+++ b/trashed_examples/LevelSet3D/levelSet3D_gpu_MultiScale.py
@@ -12,7 +12,7 @@ from hysop.numerics.remeshing import L6_4 as remesh_formula
from hysop.numerics.remeshing import Linear
from hysop.operator.advection import Advection
from hysop.operator.analytic import Analytic
-from hysop.operator.custom import CustomMonitor
+from hysop.operators import Custom
from hysop.operator.hdf_io import HDF_Writer
from hysop.gpu.gpu_transfer import DataTransfer
from hysop.problem.problem import Problem
@@ -34,7 +34,7 @@ def vitesse(res, x, y, z, t=0.):
return res
-def volume(_simu, var):
+def volume(_simu, var, **kwds):
v_loc = np.sum(var[0].data[0] > 0.5) * np.prod(
var[0].topology.mesh.space_step)
return main_comm.allreduce(sendobj=v_loc, op=MPI.SUM)
@@ -89,9 +89,9 @@ topo_s = advec.advec_dir[0].discreteFields[scal].topology
velocity = Analytic(variables={velo: topo_v})
-volume_m = CustomMonitor(function=volume, res_shape=1,
- variables={scal: topo_s},
- io_params=IOParams(filename="volume.dat"))
+volume_m = Custom(function=volume, diagnostics_shape=(1, 2),
+ in_fields=[scal], variables={scal: topo_s},
+ io_params=IOParams(filename="volume.dat"))
p = HDF_Writer(variables={scal: topo_s},
io_params=IOParams(frequency=outputModulo,
filename="levelset",
diff --git a/trashed_examples/LevelSet3D/levelSet3D_only_gpu.py b/trashed_examples/LevelSet3D/levelSet3D_only_gpu.py
index 5a87c8e07f120f1c90c351b78b58cc290b46d794..a163a5fd25e6216c551015af7378c953c4672cce 100644
--- a/trashed_examples/LevelSet3D/levelSet3D_only_gpu.py
+++ b/trashed_examples/LevelSet3D/levelSet3D_only_gpu.py
@@ -12,7 +12,7 @@ from hysop.numerics.remeshing import L6_4 as remesh_formula
from hysop.numerics.remeshing import Linear
from hysop.operator.advection import Advection
from hysop.operator.analytic import Analytic
-from hysop.operator.custom import CustomMonitor
+from hysop.operators import Custom
from hysop.operator.hdf_io import HDF_Writer
from hysop.gpu.gpu_transfer import DataTransfer
from hysop.problem.problem import Problem
@@ -20,7 +20,7 @@ from hysop.problem.simulation import Simulation
sin, cos, pi, sqrt = np.sin, np.cos, np.pi, np.sqrt
-def volume(_simu, var):
+def volume(_simu, var, **kwds):
v_loc = np.sum(var[0].data[0] > 0.5) * np.prod(
var[0].topology.mesh.space_step)
return main_comm.allreduce(sendobj=v_loc, op=MPI.SUM)
@@ -78,9 +78,9 @@ topo_s = advec.advec_dir[0].discreteFields[scal].topology
velocity = Analytic(variables={velo: topo_v},
method={Support: 'gpu'})
-volume_m = CustomMonitor(function=volume, res_shape=1,
- variables={scal: topo_s},
- io_params=IOParams(filename="volume.dat"))
+volume_m = Custom(function=volume, diagnostics_shape=(1, 2),
+ in_fields=[scal], variables={scal: topo_s},
+ io_params=IOParams(filename="volume.dat"))
p = HDF_Writer(variables={scal: topo_s},
io_params=IOParams(frequency=outputModulo,
filename="levelset",
diff --git a/trashed_examples/LevelSet3D/levelSet3D_only_gpu_MultiScale.py b/trashed_examples/LevelSet3D/levelSet3D_only_gpu_MultiScale.py
index 5a689da0c2cf7752957fa921bd2ddb76b9a5fc2a..0726d2328d3ba3fdbb281acce52108d4e5c5b79d 100644
--- a/trashed_examples/LevelSet3D/levelSet3D_only_gpu_MultiScale.py
+++ b/trashed_examples/LevelSet3D/levelSet3D_only_gpu_MultiScale.py
@@ -12,7 +12,7 @@ from hysop.numerics.remeshing import L6_4 as remesh_formula
from hysop.numerics.remeshing import Linear
from hysop.operator.advection import Advection
from hysop.operator.analytic import Analytic
-from hysop.operator.custom import CustomMonitor
+from hysop.operators import Custom
from hysop.operator.hdf_io import HDF_Writer
from hysop.gpu.gpu_transfer import DataTransfer
from hysop.problem.problem import Problem
@@ -20,7 +20,7 @@ from hysop.problem.simulation import Simulation
sin, cos, pi, sqrt = np.sin, np.cos, np.pi, np.sqrt
-def volume(_simu, var):
+def volume(simu, var, **kwds):
v_loc = np.sum(var[0].data[0] > 0.5) * np.prod(
var[0].topology.mesh.space_step)
return main_comm.allreduce(sendobj=v_loc, op=MPI.SUM)
@@ -80,9 +80,9 @@ topo_s = advec.advec_dir[0].discreteFields[scal].topology
velocity = Analytic(variables={velo: topo_v},
method={Support: 'gpu'})
-volume_m = CustomMonitor(function=volume, res_shape=1,
- variables={scal: topo_s},
- io_params=IOParams(filename="volume.dat"))
+volume_m = Custom(function=volume, diagnostics_shape=(1, 2),
+ in_fields=[scal], variables={scal: topo_s},
+ io_params=IOParams(filename="volume.dat"))
p = HDF_Writer(variables={scal: topo_s},
io_params=IOParams(frequency=outputModulo,
filename="levelset",
diff --git a/trashed_examples/LevelSet3D/levelSet3D_python.py b/trashed_examples/LevelSet3D/levelSet3D_python.py
index fc845ee45bc9242367b87fec1b67ce2a2322a848..3a3d071d9243ca0eb4ce620e73978329100a8b63 100644
--- a/trashed_examples/LevelSet3D/levelSet3D_python.py
+++ b/trashed_examples/LevelSet3D/levelSet3D_python.py
@@ -12,7 +12,7 @@ from hysop.numerics.remeshing import L2_1 as remesh_formula
from hysop.numerics.interpolation import Linear
from hysop.operator.advection import Advection
from hysop.operator.analytic import Analytic
-from hysop.operator.custom import CustomMonitor
+from hysop.operators import Custom
from hysop.operator.hdf_io import HDF_Writer
from hysop.gpu.gpu_transfer import DataTransfer
from hysop.problem.problem import Problem
@@ -34,7 +34,7 @@ def vitesse(res, x, y, z, t=0.):
return res
-def volume(_simu, var):
+def volume(_simu, var, **kwds):
v_loc = np.sum(var[0].data[0] > 0.5) * np.prod(
var[0].topology.mesh.space_step)
return main_comm.allreduce(sendobj=v_loc, op=MPI.SUM)
@@ -100,9 +100,9 @@ topo_s = advec.advec_dir[0].discreteFields[scal].topology
velocity = Analytic(variables={velo: topo_v})
#method={Support: 'gpu'})
-volume_m = CustomMonitor(function=volume, res_shape=1,
- variables={scal: topo_s},
- io_params=IOParams(filename="volume.dat"))
+volume_m = Custom(function=volume, diagnostics_shape=(1, 2),
+ in_fields=[scal], variables={scal: topo_s},
+ io_params=IOParams(filename="volume.dat"))
p = HDF_Writer(variables={scal: topo_s},
io_params=IOParams(frequency=outputModulo,
filename="levelset",
diff --git a/trashed_examples/Multiphase/NS_planeJet_hybrid_MS_MP.py b/trashed_examples/Multiphase/NS_planeJet_hybrid_MS_MP.py
index 0529eebf6eb0a9f50727fd60fe460482381f0a36..bcf74505fc3fb023a20fe7a8f10f2a043b913dad 100644
--- a/trashed_examples/Multiphase/NS_planeJet_hybrid_MS_MP.py
+++ b/trashed_examples/Multiphase/NS_planeJet_hybrid_MS_MP.py
@@ -39,7 +39,7 @@ from hysop.operator.stretching import Stretching
from hysop.operator.poisson import Poisson
from hysop.operator.differential import Curl
from hysop.operator.adapt_timestep import AdaptTimeStep
-from hysop.operator.custom import CustomMonitor
+from hysop.operators import Custom
from hysop.operator.redistribute_inter import RedistributeInter
from hysop.operator.redistribute_intra import RedistributeIntra
from hysop.gpu.gpu_transfer import DataTransfer
@@ -372,13 +372,14 @@ if IS_OUTPUT:
io_params=IOParams(frequency=1,
filename='energy.dat',
fileformat=ASCII))
- maxvelo = CustomMonitor(variables={velo: topo_F_0g_1d},
- function=calc_maxvelo,
- res_shape=3,
- mpi_params=mpi_params_UW,
- io_params=IOParams(frequency=1,
- filename='maxvelo.dat',
- fileformat=ASCII))
+ maxvelo = Custom(in_fields=[velo],
+ variables={velo: topo_F_0g_1d},
+ function=calc_maxvelo,
+ diagnostics_shape=(1, 4),
+ mpi_params=mpi_params_UW,
+ io_params=IOParams(frequency=1,
+ filename='maxvelo.dat',
+ fileformat=ASCII))
if box.is_on_task(TASK_UW):
for op in (p_velo, p_velo_xy, p_velo_xz, energy, maxvelo):
op.discretize()
diff --git a/trashed_examples/Multiphase/RTI.py b/trashed_examples/Multiphase/RTI.py
index 4534e86e61e2c3e6a5e131084012d1e86a212df0..2e23829a6b00e086a4fbe07a68629b54074ce483 100644
--- a/trashed_examples/Multiphase/RTI.py
+++ b/trashed_examples/Multiphase/RTI.py
@@ -42,7 +42,7 @@ from hysop.operator.multiresolution_filter import MultiresolutionFilter
from hysop.operator.multiphase_gradp import MultiphaseGradP
from hysop.operator.multiphase_baroclinic_rhs import MultiphaseBaroclinicRHS
from hysop.operator.adapt_timestep import AdaptTimeStep
-from hysop.operator.custom import CustomMonitor
+from hysop.operators import Custom
from hysop.operator.redistribute_inter import RedistributeInter
from hysop.operator.redistribute_intra import RedistributeIntra
from hysop.gpu.gpu_transfer import DataTransfer
@@ -338,13 +338,14 @@ if IS_OUTPUT:
io_params=IOParams(frequency=1,
filename='energy.dat',
fileformat=ASCII))
- maxvelo = CustomMonitor(variables={velo: topo_C_1d},
- function=calc_maxvelo,
- res_shape=3,
- mpi_params=mpi_params_UW,
- io_params=IOParams(frequency=1,
- filename='maxvelo.dat',
- fileformat=ASCII))
+ maxvelo = Custom(in_fields=[velo],
+ variables={velo: topo_C_1d},
+ function=calc_maxvelo,
+ diagnostics_shape=(1, 4),
+ mpi_params=mpi_params_UW,
+ io_params=IOParams(frequency=1,
+ filename='maxvelo.dat',
+ fileformat=ASCII))
if box.is_on_task(TASK_UW):
for op in (p_velo, energy, maxvelo):
op.discretize()
diff --git a/trashed_examples/Plane_jet/NS_planeJet_hybrid_MS.py b/trashed_examples/Plane_jet/NS_planeJet_hybrid_MS.py
index bee1729fe1b86ccafb94b021338fe88ea2b9c978..6f65a4f347897448a95269ecbcb04f6a8ae99cc9 100644
--- a/trashed_examples/Plane_jet/NS_planeJet_hybrid_MS.py
+++ b/trashed_examples/Plane_jet/NS_planeJet_hybrid_MS.py
@@ -33,7 +33,7 @@ from hysop.operator.stretching import Stretching
from hysop.operator.poisson import Poisson
from hysop.operator.differential import Curl
from hysop.operator.adapt_timestep import AdaptTimeStep
-from hysop.operator.custom import CustomMonitor
+from hysop.operators import Custom
from hysop.operator.redistribute_inter import RedistributeInter
from hysop.operator.redistribute_intra import RedistributeIntra
from hysop.gpu.gpu_transfer import DataTransfer
@@ -299,13 +299,14 @@ if IS_OUTPUT:
io_params=IOParams(frequency=1,
filename='energy.dat',
fileformat=ASCII))
- maxvelo = CustomMonitor(variables={velo: topo_CPU_velo_scales},
- function=calc_maxvelo,
- res_shape=3,
- mpi_params=mpi_params_UW,
- io_params=IOParams(frequency=1,
- filename='maxvelo.dat',
- fileformat=ASCII))
+ maxvelo = Custom(in_fields=[velo],
+ variables={velo: topo_CPU_velo_scales},
+ function=calc_maxvelo,
+ diagnostics_shape=(1, 4),
+ mpi_params=mpi_params_UW,
+ io_params=IOParams(frequency=1,
+ filename='maxvelo.dat',
+ fileformat=ASCII))
if box.isOnTask(TASK_UW):
for op in (p_velo, p_velo_xy, p_velo_xz, energy, maxvelo):
op.discretize()
diff --git a/trashed_examples/RMI/RMI_hybrid.py b/trashed_examples/RMI/RMI_hybrid.py
index 7f2076058b223e14f49e3867aed4564430df7c81..e895a25898d893256719293fff857570657810c0 100644
--- a/trashed_examples/RMI/RMI_hybrid.py
+++ b/trashed_examples/RMI/RMI_hybrid.py
@@ -40,7 +40,7 @@ from hysop.operator.baroclinic import Baroclinic
from hysop.operator.penalization import PenalizeVorticity
from hysop.operator.poisson import Poisson
from hysop.operator.adapt_timestep import AdaptTimeStep
-from hysop.operator.custom import CustomMonitor, CustomOp
+from hysop.operators import Custom
from hysop.operator.redistribute_inter import RedistributeInter
from hysop.operator.redistribute_intra import RedistributeIntra
from hysop.gpu.gpu_transfer import DataTransfer
@@ -129,17 +129,15 @@ def func_scal_to_rho(simu, f_in, f_out):
temp_maxvelo = npw.zeros((3, ))
-maxvelo_values = npw.zeros((3, ))
-def calc_maxvelo(simu, v):
+def calc_maxvelo(simu, v, w=None, maxvelo_values):
temp_maxvelo[0] = np.max(np.abs(v[0].data[0]))
temp_maxvelo[1] = np.max(np.abs(v[0].data[1]))
temp_maxvelo[2] = np.max(np.abs(v[0].data[2]))
v[0].topology.comm.Allreduce(sendbuf=[temp_maxvelo, 3, HYSOP_MPI_REAL],
recvbuf=[maxvelo_values, 3, HYSOP_MPI_REAL],
op=MPI.MAX)
- return maxvelo_values
ctime = MPI.Wtime()
# Domain
@@ -222,10 +220,11 @@ advec = Advection(velo,
method={Scales: 'p_64', MultiScale: 'L4_4'})
stretch = Stretching(velo, vorti, discretization=d_uw_ghosts,
mpi_params=mpi_params_UW)
-scal_to_rho = CustomOp([scal, ], [rho, ], func_scal_to_rho,
- variables={scal: d_uw_ghosts,
- rho: d_uw_ghosts},
- mpi_params=mpi_params_UW)
+scal_to_rho = Custom(in_fields=[scal, ], out_fields=[rho, ],
+ function=func_scal_to_rho,
+ variables={scal: d_uw_ghosts,
+ rho: d_uw_ghosts},
+ mpi_params=mpi_params_UW)
penal = PenalizeVorticity(velocity=velo, vorticity=vorti,
discretization=d_uw_ghosts,
obstacles=[bc_t, bc_b], coeff=1e8,
@@ -291,13 +290,14 @@ if IS_OUTPUT:
io_params=IOParams(frequency=1,
filename='energy.dat',
fileformat=ASCII))
- maxvelo = CustomMonitor(variables={velo: topo_CPU_velo_scales},
- function=calc_maxvelo,
- res_shape=3,
- mpi_params=mpi_params_UW,
- io_params=IOParams(frequency=1,
- filename='maxvelo.dat',
- fileformat=ASCII))
+ maxvelo = Custom(in_fields=[velo],
+ variables={velo: topo_CPU_velo_scales},
+ function=calc_maxvelo,
+ diagnostics_shape=(1, 4),
+ mpi_params=mpi_params_UW,
+ io_params=IOParams(frequency=1,
+ filename='maxvelo.dat',
+ fileformat=ASCII))
if box.is_on_task(TASK_UW):
for op in (p_velo, p_rho, energy, maxvelo):
op.discretize()