From 4f898569eb9add2d6f9439cb515360b212b0c0d9 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Franck=20P=C3=A9rignon?= <franck.perignon@imag.fr>
Date: Tue, 24 May 2016 14:26:40 +0200
Subject: [PATCH] Work arrays management in ops, fix tests
---
CMakeLists.txt | 2 +-
docs/sphinx/contacts.rst | 14 ++
hysop/default_methods.py | 27 +++-
hysop/domain/box.py | 6 +-
hysop/domain/tests/test_box.py | 10 +-
hysop/domain/tests/test_subset.py | 3 +-
hysop/numerics/odesolvers.py | 2 +-
hysop/numerics/tests/test_odesolvers.py | 18 ++-
hysop/operator/absorption_BC.py | 2 +-
hysop/operator/adapt_timestep.py | 2 +-
hysop/operator/custom.py | 2 +-
hysop/operator/discrete/discrete.py | 30 +----
hysop/operator/discrete/energy_enstrophy.py | 86 +++++--------
hysop/operator/discrete/stretching.py | 136 ++++++++------------
hysop/operator/drag_and_lift.py | 4 +-
hysop/operator/energy_enstrophy.py | 76 +++++------
hysop/operator/monitoringPoints.py | 2 +-
hysop/operator/poisson.py | 8 +-
hysop/operator/profiles.py | 2 +-
hysop/operator/reprojection.py | 2 +-
hysop/operator/residual.py | 2 +-
hysop/operator/stretching.py | 74 +++--------
hysop/operator/velocity_correction.py | 2 +-
hysop/operators.py | 29 +++++
hysop/tools/misc.py | 6 +-
hysop/tools/tests/test_work_arrays.py | 71 +++++++++-
26 files changed, 312 insertions(+), 306 deletions(-)
create mode 100644 docs/sphinx/contacts.rst
create mode 100644 hysop/operators.py
diff --git a/CMakeLists.txt b/CMakeLists.txt
index 3b5fb6a5d..e2750a842 100755
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -469,7 +469,7 @@ add_custom_target(pyclean COMMAND rm -f ${PYCFILES}
# ============= Tests =============
if(WITH_TESTS)
if(NOT ENABLE_LONG_TESTS)
- set(TEST_TIMEOUT 10 CACHE STRING "Default value for test runtime limit (in seconds)")
+ set(TEST_TIMEOUT 15 CACHE STRING "Default value for test runtime limit (in seconds)")
else()
set(TEST_TIMEOUT 500 CACHE STRING "Default value for test runtime limit (in seconds)")
endif()
diff --git a/docs/sphinx/contacts.rst b/docs/sphinx/contacts.rst
new file mode 100644
index 000000000..f8646da86
--- /dev/null
+++ b/docs/sphinx/contacts.rst
@@ -0,0 +1,14 @@
+.. _contacts:
+
+
+========
+Contacts
+========
+
+.. topic:: HySoP development team
+
+ `Laboratoire Jean Kuntzmann`_, Grenoble.
+
+ **E-Mail** : hysop-members@lists.forge.imag.fr
+
+.. _Laboratoire Jean Kuntzmann: http://ljk.imag.fr
diff --git a/hysop/default_methods.py b/hysop/default_methods.py
index 425cafd5f..d1e20b5c0 100644
--- a/hysop/default_methods.py
+++ b/hysop/default_methods.py
@@ -4,39 +4,54 @@ from hysop.methods_keys import TimeIntegrator, Interpolation, GhostUpdate,\
Remesh, Support, Splitting, MultiScale, Formulation, SpaceDiscretisation, \
dtCrit, Precision
from hysop.constants import HYSOP_REAL
-from hysop.numerics.integrators.runge_kutta2 import RK2
-from hysop.numerics.integrators.runge_kutta3 import RK3
-from hysop.numerics.finite_differences import FDC4, FDC2
+from hysop.numerics.odesolvers import RK2, RK3
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.operator.discrete.stretching import Conservative
ADVECTION = {TimeIntegrator: RK2, Interpolation: Linear,
Remesh: L2_1, Support: '', Splitting: 'o2', MultiScale: L2_1,
Precision: HYSOP_REAL}
-
-from hysop.numerics.finite_differences import FDC4, FDC2
+"""Advection operators default setup
+"""
DIFFERENTIAL = {SpaceDiscretisation: FDC4, GhostUpdate: True}
+"""Differential operators (finite differencies based) default setup
+"""
ADAPT_TIME_STEP = {TimeIntegrator: RK3, SpaceDiscretisation: FDC4,
dtCrit: 'vort'}
+"""Adaptative time step default setup
+"""
BAROCLINIC = {SpaceDiscretisation: FDC4}
+"""Baroclinic operators default setup
+"""
MULTIPHASEBAROCLINIC = {SpaceDiscretisation: FDC4}
MULTIPHASEGRADP = {SpaceDiscretisation: FDC4}
DIFFUSION = {SpaceDiscretisation: 'fftw', GhostUpdate: True}
+"""Differential operators (finite differencies based) default setup
+"""
POISSON = {SpaceDiscretisation: 'fftw', GhostUpdate: True,
Formulation: 'velocity'}
+"""Poisson operators default setup
+"""
STRETCHING = {TimeIntegrator: RK3, Formulation: "Conservative",
SpaceDiscretisation: FDC4}
+"""Stretching operators default setup
+"""
FORCES = {SpaceDiscretisation: FDC2}
+"""Drag/lift computation operators default setup
+"""
-MULTIRESOLUTION_FILER = {Remesh: Rmsh_Linear, }
+MULTIRESOLUTION_FILTER = {Remesh: Rmsh_Linear, }
+"""filters operators default setup
+"""
diff --git a/hysop/domain/box.py b/hysop/domain/box.py
index cbf6317fc..7a45d2e81 100644
--- a/hysop/domain/box.py
+++ b/hysop/domain/box.py
@@ -40,16 +40,16 @@ class Box(Domain):
super(Box, self).__init__(**kwds)
- ## Box length.
+ # Box length.
if length is None:
length = [1.0] * self.dimension
if origin is None:
origin = [0.] * self.dimension
self.length = npw.const_realarray(length)
- ## Box origin
+ # Box origin
self.origin = npw.const_realarray(origin)
- ## Maximum Box position. max = origin + length
+ # Maximum Box position. max = origin + length
self.end = self.origin + self.length
# set periodic boundary conditions
if self.boundaries is None:
diff --git a/hysop/domain/tests/test_box.py b/hysop/domain/tests/test_box.py
index d9ab03cd0..24a12f79a 100755
--- a/hysop/domain/tests/test_box.py
+++ b/hysop/domain/tests/test_box.py
@@ -124,20 +124,16 @@ def test_topo_plane():
def test_topo_from_mesh():
# e.g. for fftw
dom = Box(proc_tasks=proc_tasks)
- from hysop.f2hysop import fftw2py
- if dom.is_on_task(CPU):
+ from hysop import __FFTW_ENABLED__
+ if dom.is_on_task(CPU) and __FFTW_ENABLED__:
+ from hysop.f2hysop import fftw2py
localres, global_start = fftw2py.init_fftw_solver(
r3D.resolution, dom.length, comm=comm_s.py2f())
print localres, global_start
topo = dom.create_plane_topology_from_mesh(localres=localres,
global_start=global_start,
discretization=r3D)
- elif dom.is_on_task(GPU):
- topo = dom.create_topology(discretization=r3DGh, dim=2)
- if dom.is_on_task(CPU):
assert (topo.mesh.resolution == localres).all()
assert (topo.mesh.start() == global_start).all()
assert topo.dimension == 1
assert (topo.shape == [1, 1, comm_s.Get_size()]).all()
- elif dom.is_on_task(GPU):
- assert topo.size == 1
diff --git a/hysop/domain/tests/test_subset.py b/hysop/domain/tests/test_subset.py
index 8501fc3b4..8f3de686d 100755
--- a/hysop/domain/tests/test_subset.py
+++ b/hysop/domain/tests/test_subset.py
@@ -183,7 +183,8 @@ def init_subtract_lists(discr, fileref):
sd = Field(domain=topo.domain, is_vector=False, name='s2').discretize(topo)
indices = Subset.subtract_list_of_sets(box, obs, topo)
sd[0][indices] = 2.
- return np.allclose(sd[0][topo.mesh.compute_index], vd[0][topo.mesh.compute_index])
+ return np.allclose(sd[0][topo.mesh.compute_index],
+ vd[0][topo.mesh.compute_index])
def test_subtract_list_3d():
diff --git a/hysop/numerics/odesolvers.py b/hysop/numerics/odesolvers.py
index 00d55837f..1ef90454c 100755
--- a/hysop/numerics/odesolvers.py
+++ b/hysop/numerics/odesolvers.py
@@ -101,7 +101,7 @@ class ODESolver(object):
self.in_indices[i].start
for i in xrange(len(self.in_indices))])
else:
- subsize = wk_prop['rwork']
+ subsize = topo.mesh.resolution
return WorkSpaceTools.check_work_array(len(wk_prop), subsize, rwork)
diff --git a/hysop/numerics/tests/test_odesolvers.py b/hysop/numerics/tests/test_odesolvers.py
index f1122cda8..b75b58418 100755
--- a/hysop/numerics/tests/test_odesolvers.py
+++ b/hysop/numerics/tests/test_odesolvers.py
@@ -44,8 +44,12 @@ def rhs3d(t, y, work):
return work
-# -- 1D cases --
def integrate(integ, nb_steps):
+ """Execute odesolver 'integ' nb_steps times.
+
+ * integ must be an instance of an odesolver (i.e. Euler, RK2 ...).
+ * the system must be 1D.
+ """
t = start
time_points = np.linspace(start, end, nb_steps)
dtt = time_points[1] - time_points[0]
@@ -67,8 +71,12 @@ def integrate(integ, nb_steps):
return dtt, err
-# -- 3D cases --
def integrate_3d(integ, nb_steps):
+ """Execute odesolver 'integ' nb_steps times.
+
+ * integ must be an instance of an odesolver (i.e. Euler, RK2 ...).
+ * 3D system
+ """
t = start
time_points = np.linspace(start, end, nb_steps)
dtt = time_points[1] - time_points[0]
@@ -114,7 +122,7 @@ def run_integ(integrator, order):
def run_integ_3d(integrator, order):
- nb_steps = 100
+ nb_steps = 50
wk_prop = integrator.get_work_properties(3, topo3d)['rwork']
work = []
for shape in wk_prop:
@@ -126,7 +134,7 @@ def run_integ_3d(integrator, order):
def run_integ_3d_no_work(integrator, order):
- nb_steps = 100
+ nb_steps = 50
dt, err = integrate_3d(integrator(3, topo3d, f=rhs3d), nb_steps)
for e in err:
assert e < dt ** order
@@ -170,7 +178,7 @@ def test_rk4_3d():
def test_user_defined_common_workspace():
integrators = {Euler: 1, RK2: 2, RK3: 3, RK4: 4}
- nb_steps = 100
+ nb_steps = 50
wk_prop = {}
for integ in integrators:
wk_prop[integ] = integ.get_work_properties(3, topo3d)['rwork']
diff --git a/hysop/operator/absorption_BC.py b/hysop/operator/absorption_BC.py
index 303ac78ae..4d27323be 100755
--- a/hysop/operator/absorption_BC.py
+++ b/hysop/operator/absorption_BC.py
@@ -69,5 +69,5 @@ class AbsorptionBC(Computational):
cb=self.cb, rwork=rwork, iwork=iwork)
# Output setup
self._set_io('absorption_BC', (1, 2 + self.domain.dimension))
- self.discrete_op.setWriter(self._writer)
+ self.discrete_op.set_writer(self._writer)
self._is_uptodate = True
diff --git a/hysop/operator/adapt_timestep.py b/hysop/operator/adapt_timestep.py
index b298ea720..1e0750355 100755
--- a/hysop/operator/adapt_timestep.py
+++ b/hysop/operator/adapt_timestep.py
@@ -123,7 +123,7 @@ class AdaptTimeStep(Computational):
rwork=rwork, iwork=iwork)
# Output setup
self._set_io('dt_adapt', (1, 7))
- self.discrete_op.setWriter(self._writer)
+ self.discrete_op.set_writer(self._writer)
self._is_uptodate = True
def wait(self):
diff --git a/hysop/operator/custom.py b/hysop/operator/custom.py
index d7d8d9124..5af161e89 100644
--- a/hysop/operator/custom.py
+++ b/hysop/operator/custom.py
@@ -41,5 +41,5 @@ class CustomMonitor(Computational):
variables=self.discreteFields.values())
# Output setup
self._set_io(self.function.__name__, (1, 1 + self.res_shape))
- self.discrete_op.setWriter(self._writer)
+ self.discrete_op.set_writer(self._writer)
self._is_uptodate = True
diff --git a/hysop/operator/discrete/discrete.py b/hysop/operator/discrete/discrete.py
index 2069227a8..94e5418a3 100755
--- a/hysop/operator/discrete/discrete.py
+++ b/hysop/operator/discrete/discrete.py
@@ -99,29 +99,7 @@ class DiscreteOperator(object):
for v in self.variables:
assert v.topology.is_consistent_with(toporef)
- def get_work_properties(self):
- """
- Get properties of internal work arrays.
- Returns
- -------
- dictionnary
- keys = 'rwork' and 'iwork', values = list of shapes of internal
- arrays required by this operator (real arrays for rwork, integer
- arrays for iwork).
-
- Example
- -------
- >>> works_prop = op.get_work_properties()
- >>> print works_prop
- {'rwork': [(12, 12), (45, 12, 33)], 'iwork': None}
-
- means that the operator requires two real arrays of
- shape (12,12) and (45, 12, 33) and no integer arrays
-
- """
- return {'rwork': None, 'iwork': None}
-
- def _set_work_arrays(self, rwork, iwork):
+ def _set_work_arrays(self, rwork=None, iwork=None):
"""
To set the internal work arrays used by this operator.
Parameters
@@ -134,7 +112,7 @@ class DiscreteOperator(object):
"""
pass
- def setWriter(self, writer):
+ def set_writer(self, writer):
"""
Assign a writer to the current operator
"""
@@ -159,8 +137,8 @@ class DiscreteOperator(object):
def __str__(self):
"""Common printings for discrete operators."""
- shortName = str(self.__class__).rpartition('.')[-1][0:-2]
- s = shortName + " discrete operator. \n"
+ short_name = str(self.__class__).rpartition('.')[-1][0:-2]
+ s = short_name + " discrete operator. \n"
if self.input is not None:
s += "Input fields : \n"
for f in self.input:
diff --git a/hysop/operator/discrete/energy_enstrophy.py b/hysop/operator/discrete/energy_enstrophy.py
index 974cf1665..06c772270 100644
--- a/hysop/operator/discrete/energy_enstrophy.py
+++ b/hysop/operator/discrete/energy_enstrophy.py
@@ -1,38 +1,41 @@
# -*- coding: utf-8 -*-
-"""
-@file energy_enstrophy.py
-Compute Energy and Enstrophy
+"""Discrete operator to compute Energy and Enstrophy
"""
from hysop.constants import debug
from hysop.tools.profiler import profile
import hysop.tools.numpywrappers as npw
from hysop.operator.discrete.discrete import DiscreteOperator
+from hysop.tools.misc import WorkSpaceTools
class EnergyEnstrophy(DiscreteOperator):
- """
- Discretization of the energy/enstrophy computation process.
+ """Discretization of the energy/enstrophy computation process.
"""
def __init__(self, velocity, vorticity, is_normalized=True, **kwds):
"""
- Constructor.
- @param velocity : discretization of the velocity field
- @param vorticity : discretization of the vorticity field
- @param coeffs : dict of coefficients
+ Parameters
+ ----------
+ velocity : :class:`~hysop.operator.field.discrete.DiscreteField`
+ velocity discrete fields
+ vorticity : :class:`~hysop.operator.field.discrete.DiscreteField`
+ vorticity discrete fields
+ isNormalized : boolean
+ true if enstrophy and energy values have to be normalized
+ by the domain lengths.
"""
- ## velocity field
+ # velocity field
self.velocity = velocity
- ## vorticity field
+ # vorticity field
self.vorticity = vorticity
assert 'variables' not in kwds, 'variables parameter is useless.'
super(EnergyEnstrophy, self).__init__(variables=[velocity, vorticity],
**kwds)
- ## Coeffs for integration
+ # Coeffs for integration
self.coeff = {}
- ## Global energy
+ # Global energy
self.energy = 0.0
- ## Global enstrophy
+ # Global enstrophy
self.enstrophy = 0.0
topo_w = self.vorticity.topology
topo_v = self.velocity.topology
@@ -52,37 +55,11 @@ class EnergyEnstrophy(DiscreteOperator):
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
- # setup for rwork, iwork is useless.
- if rwork is None:
- # --- Local allocation ---
- if shape_v == shape_w:
- self._rwork = [npw.zeros(shape_v)]
- else:
- self._rwork = [npw.zeros(shape_v), npw.zeros(shape_w)]
- else:
- assert isinstance(rwork, list), 'rwork must be a list.'
- # --- External rwork ---
- self._rwork = rwork
- if shape_v == shape_w:
- assert len(self._rwork) == 1
- assert self._rwork[0].shape == shape_v
- else:
- assert len(self._rwork) == 2
- 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}
+ size_v = self.velocity.data[0][v_ind].size
+ size_w = self.vorticity.data[0][w_ind].size
+ size_work = max(size_w, size_v)
+ lwork = 1
+ self._rwork = WorkSpaceTools.check_work_array(lwork, size_work, rwork)
@debug
@profile
@@ -94,23 +71,24 @@ class EnergyEnstrophy(DiscreteOperator):
vd = self.velocity
# get the list of computation points (no ghosts)
nbc = vd.nb_components
- v_ind = self.velocity.topology.mesh.compute_index
+ ic = self.velocity.topology.mesh.compute_index
# Integrate (locally) velocity ** 2
local_energy = 0.
+ size_v = vd[0][ic].size
for i in xrange(nbc):
- self._rwork[0][...] = vd[i][v_ind] ** 2
- local_energy += npw.real_sum(self._rwork[0])
+ self._rwork[0][:size_v] = (vd[i][ic] ** 2).flat
+ local_energy += npw.real_sum(self._rwork[0][:size_v])
# --- Enstrophy computation ---
- vortd = self.vorticity
- nbc = vortd.nb_components
- w_ind = self.vorticity.topology.mesh.compute_index
+ wd = self.vorticity
+ nbc = wd.nb_components
+ ic = self.vorticity.topology.mesh.compute_index
+ size_w = wd[0][ic].size
# Integrate (locally) vorticity ** 2
- work = self._rwork[-1]
local_enstrophy = 0.
for i in xrange(nbc):
- work[...] = vortd[i][w_ind] ** 2
- local_enstrophy += npw.real_sum(work)
+ self._rwork[0][:size_w] = (wd[i][ic] ** 2).flat
+ local_enstrophy += npw.real_sum(self._rwork[0][:size_w])
# --- Reduce energy and enstrophy values overs all proc ---
# two ways : numpy or classical. Todo : check perf and comm
diff --git a/hysop/operator/discrete/stretching.py b/hysop/operator/discrete/stretching.py
index 4f1527605..8e3d2e16d 100755
--- a/hysop/operator/discrete/stretching.py
+++ b/hysop/operator/discrete/stretching.py
@@ -5,33 +5,26 @@ Formulations:
* :class:`~hysop.operator.discrete.stretching.Conservative`
* :class:`~hysop.operator.discrete.stretching.GradUW`
-* :class:`~hysop.operator.discrete.stretching.Symmetric`
"""
from hysop.constants import debug, WITH_GUESS
from hysop.methods_keys import TimeIntegrator, SpaceDiscretisation
from hysop.operator.discrete.discrete import DiscreteOperator
-from hysop.numerics.integrators.euler import Euler
-from hysop.numerics.integrators.runge_kutta2 import RK2
-from hysop.numerics.integrators.runge_kutta3 import RK3
-from hysop.numerics.integrators.runge_kutta4 import RK4
import hysop.numerics.differential_operations as diff_op
import hysop.tools.numpywrappers as npw
from hysop.numerics.update_ghosts import UpdateGhosts
from hysop.tools.profiler import profile
+from hysop.tools.misc import WorkSpaceTools
from abc import ABCMeta, abstractmethod
import math
+import numpy as np
ceil = math.ceil
class Stretching(DiscreteOperator):
- """
- Abstract interface to stretching discrete operators.
- Three formulations are available:
- - conservative, see operator.discrete.stretching.Conservative
- - 'Grad(UxW), see operator.discrete.stretching.GradUW
- - Symmetric, see operator.discrete.stretching.Symmetric
+ """Abstract interface to stretching discrete operators.
+ Check user guide to see the available formulations.
"""
__metaclass__ = ABCMeta
@@ -61,6 +54,8 @@ class Stretching(DiscreteOperator):
self._ti_work = None
# Work vector used by numerical diff operator.
self._str_work = None
+ # all works vectors are set in base class init,
+ # thanks to _set_work_arrays function.
super(Stretching, self).__init__(variables=[self.velocity,
self.vorticity], **kwds)
@@ -75,94 +70,74 @@ class Stretching(DiscreteOperator):
# prepare ghost points synchro for velocity and vorticity
self._synchronize = UpdateGhosts(self.velocity.topology,
- self.velocity.nb_components
- + self.vorticity.nb_components)
+ self.velocity.nb_components +
+ self.vorticity.nb_components)
- # A function to compute the gradient of a vector field
+ # A function to compute the gradient of the velocity.
# Work vector is provided in input.
self.strFunc = \
self.formulation(topo=self.velocity.topology,
work=self._str_work,
method=self.method[SpaceDiscretisation])
-
# Time integrator
self.timeIntegrator = \
self.method[TimeIntegrator](self.nb_components,
- self._ti_work,
self.velocity.topology,
+ rwork=self._ti_work,
f=rhs,
optim=WITH_GUESS)
+ # In self._integrate ti_work is used as in/out parameters in the rhs,
+ # of the time-integrator so it needs a reshape.
+ # Warning : this is done without any memory realloc.
+ for i in xrange(len(self._ti_work)):
+ self._ti_work[i] = self.timeIntegrator.rwork[i]
def _set_work_arrays(self, rwork=None, iwork=None):
- shape_v = self.velocity.data[0][...].shape
ti = self.method[TimeIntegrator]
- # work list length for time-integrator
- work_length_ti = ti.getWorkLengths(3)
- rwork_length = work_length_ti + self.formulation.get_work_length()
- # setup for rwork, iwork is useless.
- if rwork is None:
- # --- Local allocation ---
- self._rwork = []
- for _ in xrange(rwork_length):
- self._rwork.append(npw.zeros(shape_v))
- else:
- assert isinstance(rwork, list), 'rwork must be a list.'
- # --- External rwork ---
- self._rwork = rwork
- assert len(self._rwork) == rwork_length
- for wk in rwork:
- assert wk.shape == shape_v
- self._ti_work = self._rwork[:work_length_ti]
- self._str_work = self._rwork[work_length_ti:]
+ topo = self.velocity.topology
+ nbc = self.velocity.nb_components
+ # properties of working arrays required for time-integrator
+ wk_prop = ti.get_work_properties(nbc, topo)['rwork']
+ lenti = len(wk_prop)
+ # properties of working arrays required for differential operations,
+ wk_prop += self.formulation.get_work_properties(topo)['rwork']
+ subshape = np.prod(topo.mesh.resolution)
+ self._rwork = WorkSpaceTools.check_work_array(len(wk_prop),
+ subshape, rwork)
+ self._ti_work = self._rwork[:lenti]
+ self._str_work = self._rwork[lenti:]
@profile
def update_ghosts(self):
+ """Ghost points synchronization
"""
- Update ghost points values
- """
- self._synchronize(self.velocity.data + self.vorticity.data)
-
- def _apply(self, t, dt):
- # Synchronize ghost points
self._synchronize(self.velocity.data + self.vorticity.data)
- # Compute stretching term (rhs) and update vorticity
- self._compute(dt, t)
def apply(self, simulation=None):
+ """Solve stretching equation and update vorticity
"""
- """
assert simulation is not None, \
"Missing simulation value for computation."
- # current time
- t = simulation.time
+
# time step
dt = simulation.time_step
- self._apply(t, dt)
-
-# 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
-#
-# # Synchronize ghost points of velocity
-# self._synchronize(self.velocity.data + self.vorticity.data)
-# self._compute(dt, t)
+ # current time
+ t = simulation.time
+
+ # Synchronize ghost points of velocity
+ self._synchronize(self.velocity.data + self.vorticity.data)
+ self._compute(dt, t)
@abstractmethod
def _compute(self, dt, t):
- """
+ """Abstract interface to local stretching resolution
"""
@profile
def _integrate(self, dt, t):
- # - Call time integrator -
+ """Apply time integrator
+ """
# Init workspace with a first evaluation of the
# rhs of the integrator
self._ti_work[:self.nb_components] = \
@@ -182,8 +157,9 @@ class Conservative(Stretching):
# Right-hand side for time integration
def rhs(t, y, result):
- result = self.strFunc(y, self.velocity.data, result)
- return result
+ """rhs used in time integrator
+ """
+ return self.strFunc(y, self.velocity.data, result)
super(Conservative, self).__init__(formulation=diff_op.DivWV,
rhs=rhs, **kwds)
@@ -199,10 +175,10 @@ class GradUW(Stretching):
"""
def __init__(self, **kwds):
- # a vector to save diagnostics computed from GradVxW (max div ...)
- self.diagnostics = npw.ones(2)
def rhs(t, y, result):
+ """rhs used in time integrator
+ """
result, self.diagnostics =\
self.strFunc(self.velocity.data, y, result, self.diagnostics)
return result
@@ -211,22 +187,15 @@ class GradUW(Stretching):
rhs=rhs, **kwds)
# stability constant
- self.cststretch = 0.
# Depends on time integration method
- timeint = self.method[TimeIntegrator]
- classtype = timeint.mro()[0]
- if classtype is Euler:
- self.cststretch = 2.0
- elif classtype is RK2:
- self.cststretch = 2.0
- elif classtype is RK3:
- self.cststretch = 2.5127
- elif classtype is RK4:
- self.cststretch = 2.7853
+ self.cststretch = self.method[TimeIntegrator].stability_coeff()
+ # a vector to save diagnostics computed from GradVxW (max div ...)
+ self.diagnostics = npw.zeros(2)
+ self.diagnostics[1] = self.cststretch
@debug
@profile
- def _compute(self, t, dt):
+ def _compute(self, dt, t):
# Compute the number of required subcycles
ndt, subdt = self._check_stability(dt)
assert sum(subdt) == dt
@@ -255,12 +224,11 @@ class GradUW(Stretching):
return nb_cycles, subdt
-class StretchingLinearized(Stretching):
- """By default: Conservative formulation
+class StretchingLinearized(Conservative):
+ """Conservative formulation of the linearized stretching.
"""
-
+
def __init__(self, vorticity_BF, usual_op, **kwds):
-
# vorticity of the base flow (steady solution)
self.vorticity_BF = vorticity_BF
# prepare ghost points synchro for vorticity_BF
diff --git a/hysop/operator/drag_and_lift.py b/hysop/operator/drag_and_lift.py
index 537a76739..0ded0b6fd 100755
--- a/hysop/operator/drag_and_lift.py
+++ b/hysop/operator/drag_and_lift.py
@@ -136,7 +136,7 @@ class MomentumForces(Forces):
# output setup
self._set_io('drag_and_lift', (1, 4))
- self.discrete_op.setWriter(self._writer)
+ self.discrete_op.set_writer(self._writer)
self._is_uptodate = True
def get_work_properties(self):
@@ -266,5 +266,5 @@ class NocaForces(Forces):
# output setup
self._set_io('drag_and_lift', (1, 4))
- self.discrete_op.setWriter(self._writer)
+ self.discrete_op.set_writer(self._writer)
self._is_uptodate = True
diff --git a/hysop/operator/energy_enstrophy.py b/hysop/operator/energy_enstrophy.py
index a59c1edeb..5251e6502 100644
--- a/hysop/operator/energy_enstrophy.py
+++ b/hysop/operator/energy_enstrophy.py
@@ -1,7 +1,8 @@
# -*- coding: utf-8 -*-
-"""
-@file energy_enstrophy.py
-Compute Energy and Enstrophy
+"""Compute Energy and Enstrophy from velocity and vorticity
+
+
+See :ref:`energy_enstrophy`.
"""
from hysop.operator.discrete.energy_enstrophy import EnergyEnstrophy as DEE
from hysop.operator.computational import Computational
@@ -10,78 +11,71 @@ from hysop.operator.continuous import opsetup
class EnergyEnstrophy(Computational):
"""
- Computes enstrophy and the kinetic energy
- \f{eqnarray*}
- enstrophy = \frac{1}{\Omega}\int_\Omega \omega^2 d\Omega
- \f} with \f$\Omega\f$ the volume or surface of the physical domain
- \f$\omega\f$ the vorticity and
- \f{eqnarray*}
- energy = \frac{1}{2\Omega}\int_\Omega v^2 d\Omega
- \f}
+ Computes enstrophy and kinetic energy
"""
def __init__(self, velocity, vorticity, is_normalized=True, **kwds):
"""
- Constructor.
- @param velocity field
- @param vorticity field
- @param isNormalized : boolean indicating whether the enstrophy
- and energy values have to be normalized by the domain lengths.
- Default file name = 'energy_enstrophy.dat'
- See hysop.tools.io_utils.Writer for details
+ Parameters
+ ----------
+ velocity, vorticity : :class:`~hysop.operator.field.continuous.Field`
+ velocity and vorticity fields
+ isNormalized : boolean
+ true if enstrophy and energy values have to be normalized
+ by the domain lengths.
+
+ Notes
+ -----
+ By default, values at each time step are appended
+ to energy_enstrophy file. See :ref:`operators_io`
+ for more details.
"""
assert 'variables' not in kwds, 'variables parameter is useless.'
super(EnergyEnstrophy, self).__init__(variables=[velocity, vorticity],
**kwds)
- ## velocity field
+ # velocity field
self.velocity = velocity
- ## vorticity field
+ # vorticity field
self.vorticity = vorticity
- ## are the energy end enstrophy values normalized by domain lengths ?
+ # are the energy end enstrophy values normalized by domain lengths ?
self.is_normalized = is_normalized
- ## self._buffer_1 = 0.
- ## self._buffer_2 = 0.
+ # self._buffer_1 = 0.
+ # self._buffer_2 = 0.
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(EnergyEnstrophy, self).get_work_properties()
vd = self.discreteFields[self.velocity]
wd = self.discreteFields[self.vorticity]
v_ind = vd.topology.mesh.compute_index
w_ind = wd.topology.mesh.compute_index
- shape_v = vd[0][v_ind].shape
- shape_w = wd[0][w_ind].shape
- if shape_v == shape_w:
- return {'rwork': [shape_v], 'iwork': None}
- else:
- return {'rwork': [shape_v, shape_w], 'iwork': None}
+ size_v = vd[0][v_ind].size
+ size_w = wd[0][w_ind].size
+ size_work = max(size_v, size_w)
+ # work array is a flat array, shared between energy and enstrophy
+ return {'rwork': [size_work], 'iwork': None}
@opsetup
def setup(self, rwork=None, iwork=None):
if not self._is_uptodate:
self.discrete_op = DEE(self.discreteFields[self.velocity],
- self.discreteFields[self.vorticity],
- self.is_normalized,
- rwork=rwork)
+ self.discreteFields[self.vorticity],
+ self.is_normalized,
+ rwork=rwork)
# Output setup
self._set_io('energy_enstrophy', (1, 3))
- self.discrete_op.setWriter(self._writer)
+ self.discrete_op.set_writer(self._writer)
self._is_uptodate = True
def energy(self):
- """
- Return last computed value of the energy
+ """Return last computed value of the energy
"""
return self.discrete_op.energy
def enstrophy(self):
- """
- Return last computed value of the enstrophy
+ """Return last computed value of the enstrophy
"""
return self.discrete_op.enstrophy
diff --git a/hysop/operator/monitoringPoints.py b/hysop/operator/monitoringPoints.py
index d67adc99f..4626d2289 100644
--- a/hysop/operator/monitoringPoints.py
+++ b/hysop/operator/monitoringPoints.py
@@ -63,6 +63,6 @@ class MonitoringPoints(Computational):
self.monitPt_coords, rwork=rwork)
# Output setup
self._set_io('monit', (1, 3))
- self.discrete_op.setWriter(self._writer)
+ self.discrete_op.set_writer(self._writer)
self._is_uptodate = True
diff --git a/hysop/operator/poisson.py b/hysop/operator/poisson.py
index 5151d2380..324ed1168 100644
--- a/hysop/operator/poisson.py
+++ b/hysop/operator/poisson.py
@@ -9,6 +9,7 @@ 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):
@@ -83,7 +84,7 @@ class Poisson(Computational):
def discretize(self):
# Poisson solver based on fftw
- if self.method[SpaceDiscretisation] is 'fftw':
+ if self.method[SpaceDiscretisation] is 'fftw' and __FFTW_ENABLED__:
super(Poisson, self)._fftw_discretize()
if self.withCorrection:
toporef = self.discreteFields[self.output_field].topology
@@ -102,7 +103,7 @@ class Poisson(Computational):
**self._config)
self.projection.discretize()
else:
- raise AttributeError("Method not yet implemented.")
+ raise AttributeError("Method not implemented.")
@debug
@opsetup
@@ -131,3 +132,6 @@ class Poisson(Computational):
formulation=self.formulation)
self._is_uptodate = True
+
+ def get_work_properties(self):
+ return {'rwork': None, 'iwork': None}
diff --git a/hysop/operator/profiles.py b/hysop/operator/profiles.py
index ac7dfbdfe..3f74bbb29 100644
--- a/hysop/operator/profiles.py
+++ b/hysop/operator/profiles.py
@@ -70,6 +70,6 @@ class Profiles(Computational):
rwork=rwork)
# Output setup
self._set_io('profile', (1, 9))
- self.discrete_op.setWriter(self._writer)
+ self.discrete_op.set_writer(self._writer)
self._is_uptodate = True
diff --git a/hysop/operator/reprojection.py b/hysop/operator/reprojection.py
index ded883620..d252314a2 100644
--- a/hysop/operator/reprojection.py
+++ b/hysop/operator/reprojection.py
@@ -43,7 +43,7 @@ class Reprojection(Computational):
self.frequency, rwork=rwork,
method=self.method)
self._set_io('reprojection', (1, 4))
- self.discrete_op.setWriter(self._writer)
+ self.discrete_op.set_writer(self._writer)
self._is_uptodate = True
def do_projection(self, ite):
diff --git a/hysop/operator/residual.py b/hysop/operator/residual.py
index 077f52632..3eeee3530 100644
--- a/hysop/operator/residual.py
+++ b/hysop/operator/residual.py
@@ -48,6 +48,6 @@ class Residual(Computational):
self.discrete_op.initialize_vortPrev()
# Output setup
self._set_io('residual', (1, 3))
- self.discrete_op.setWriter(self._writer)
+ self.discrete_op.set_writer(self._writer)
self._is_uptodate = True
diff --git a/hysop/operator/stretching.py b/hysop/operator/stretching.py
index 4b183ac67..7c58412dc 100755
--- a/hysop/operator/stretching.py
+++ b/hysop/operator/stretching.py
@@ -70,58 +70,26 @@ class Stretching(Computational):
self.output = [self.vorticity]
def get_work_properties(self):
- """
- Get properties of internal work arrays. Must be call after discretize
- but before setup.
-
- Returns
- -------
- dictionnary
- keys = 'rwork' and 'iwork', values = list of shapes of internal
- arrays required by this operator (real arrays for rwork, integer
- arrays for iwork).
-
- Example
- -------
-
- >> works_prop = op.get_work_properties()
- >> print works_prop
- {'rwork': [(12, 12), (45, 12, 33)], 'iwork': None}
-
- means that the operator requires two real arrays of
- shape (12,12) and (45, 12, 33) and no integer arrays
-
- """
- if not self._is_discretized:
- msg = 'The operator must be discretized '
- msg += 'before any call to this function.'
- raise RuntimeError(msg)
+ super(Stretching, self).get_work_properties()
# Get fields local shape.
vd = self.discreteFields[self.velocity]
- shape_v = vd[0][...].shape
# Collect info from numerical methods
# --> time-integrator required work space
ti = self.method[TimeIntegrator]
- rwork_length = ti.getWorkLengths(3)
+ topo = vd.topology
+ nbc = self.velocity.nb_components
+ res = ti.get_work_properties(nbc, topo)
# ---> differential operator work space
if self.formulation is GradUW:
- rwork_length += diff_op.GradVxW.get_work_length()
+ dop = diff_op.GradVxW
elif self.formulation is Conservative:
- rwork_length += diff_op.DivWV.get_work_length()
- res = {'rwork': [], 'iwork': None}
- # Set rwork shapes
- for _ in xrange(rwork_length):
- res['rwork'].append(shape_v)
+ dop = diff_op.DivWV
+ res['rwork'] += dop.get_work_properties(topo)['rwork']
return res
@debug
def discretize(self):
- if self.method[SpaceDiscretisation] is FDC4:
- nbghosts = 2
- else:
- raise ValueError("Unknown method for space discretization of the\
- stretching operator.")
-
+ nbghosts = self.method[SpaceDiscretisation].ghosts_layer_size
super(Stretching, self)._standard_discretize(nbghosts)
@debug
@@ -142,7 +110,6 @@ class StretchingLinearized(Stretching):
(\om_b \cdot \nabla)u
\f}
"""
-
@debug
def __init__(self, velocity_BF, vorticity_BF, **kwds):
"""
@@ -151,13 +118,12 @@ class StretchingLinearized(Stretching):
velocity_BF, vorticity_BF : base flow fields :
class:`~hysop.fields.continuous.Field`
**kwds : extra parameters for base class
-
+
Notes
-----
* The default formulation is the 'Conservative' one.
* The default solving method is finite differences, 4th order,
in space and Runge-Kutta 3 in time.
-
"""
super(StretchingLinearized, self).__init__(**kwds)
self.variables[velocity_BF] = None
@@ -167,7 +133,6 @@ class StretchingLinearized(Stretching):
self.velocity_BF = velocity_BF
# Base flow vorticity variable (vector)
self.vorticity_BF = vorticity_BF
-
# Usual stretching operator to compute the
# first term of the linearized rhs: (w'.grad)ub
self.usual_stretch = Stretching(self.velocity_BF, self.vorticity,
@@ -175,18 +140,13 @@ class StretchingLinearized(Stretching):
self.input.append(self.velocity_BF)
self.input.append(self.vorticity_BF)
-
+
@debug
def discretize(self):
- if self.method[SpaceDiscretisation] is FDC4:
- nbghosts = 2
- else:
- raise ValueError("Unknown method for space discretization of the\
- stretching operator.")
+ nbghosts = self.method[SpaceDiscretisation].ghosts_layer_size
self.usual_stretch.discretize()
super(StretchingLinearized, self)._standard_discretize(nbghosts)
-
@debug
@opsetup
def setup(self, rwork=None, iwork=None):
@@ -198,11 +158,9 @@ class StretchingLinearized(Stretching):
method_lin = self.method.copy()
method_lin[TimeIntegrator] = Euler
self.discrete_op =\
- StretchLinD(velocity=self.discreteFields[self.velocity],
- vorticity=self.discreteFields[self.vorticity],
- vorticity_BF=self.discreteFields[self.vorticity_BF],
- usual_op=self.usual_stretch.discrete_op,
- method=method_lin, rwork=rwork, iwork=iwork)
+ SLD(velocity=self.discreteFields[self.velocity],
+ vorticity=self.discreteFields[self.vorticity],
+ vorticity_BF=self.discreteFields[self.vorticity_BF],
+ usual_op=self.usual_stretch.discrete_op,
+ method=method_lin, rwork=rwork, iwork=iwork)
self._is_uptodate = True
-
-
diff --git a/hysop/operator/velocity_correction.py b/hysop/operator/velocity_correction.py
index 088eabd87..7cdd7e1c4 100755
--- a/hysop/operator/velocity_correction.py
+++ b/hysop/operator/velocity_correction.py
@@ -67,7 +67,7 @@ class VelocityCorrection(Computational):
iwork=iwork)
# Output setup
self._set_io('velocity_correction', (1, 2 + self.domain.dimension))
- self.discrete_op.setWriter(self._writer)
+ self.discrete_op.set_writer(self._writer)
self._is_uptodate = True
def computeCorrection(self):
diff --git a/hysop/operators.py b/hysop/operators.py
new file mode 100644
index 000000000..9dfb7a2eb
--- /dev/null
+++ b/hysop/operators.py
@@ -0,0 +1,29 @@
+"""Shortcuts for operators import
+
+Allows things like:
+
+from hysop.operators import Advection
+"""
+import hysop.operator.stretching
+import hysop.operator.differential
+import hysop.operator.analytic
+import hysop.operator.penalization
+import hysop.operator.drag_and_lift
+import hysop.operator.hdf_io
+import hysop.operator.poisson
+import hysop.operator.advection
+import hysop.operator.energy_enstrophy
+Stretching = hysop.operator.stretching.Stretching
+Curl = hysop.operator.differential.Curl
+Grad = hysop.operator.differential.Grad
+DivAdvection = hysop.operator.differential.DivAdvection
+Analytic = hysop.operator.analytic.Analytic
+Penalization = hysop.operator.penalization.Penalization
+PenalizeVorticity = hysop.operator.penalization.PenalizeVorticity
+MomentumForces = hysop.operator.drag_and_lift.MomentumForces
+NocaForces = hysop.operator.drag_and_lift.NocaForces
+HDF_Reader = hysop.operator.hdf_io.HDF_Reader
+HDF_Writer = hysop.operator.hdf_io.HDF_Writer
+Poisson = hysop.operator.poisson.Poisson
+Advection = hysop.operator.advection.Advection
+EnergyEnstrophy = hysop.operator.energy_enstrophy.EnergyEnstrophy
diff --git a/hysop/tools/misc.py b/hysop/tools/misc.py
index a065daba8..ce8c646df 100644
--- a/hysop/tools/misc.py
+++ b/hysop/tools/misc.py
@@ -120,6 +120,8 @@ class WorkSpaceTools(object):
assert len(np.where(
np.asarray(wk.shape) > 1)[0]) == 1, msg2
result.append(wk.ravel()[:subsize[i]].reshape(subshape[i]))
+ for i in xrange(len(result)):
+ assert npw.arrays_share_data(result[i], work[i])
except AttributeError:
# Work array has been replaced by an OpenCL Buffer
@@ -231,7 +233,9 @@ class WorkSpaceTools(object):
shapes = [(0,), ] * lwork
for prop in properties:
lp = len(prop)
- shapes[:lp] = np.maximum(shapes[:lp], prop)
+ for i in xrange(lp):
+ shapes[i] = tuple(np.maximum(shapes[i], np.prod(prop[i])))
+ print shapes
work = [npw.zeros(shape) for shape in shapes]
return work
diff --git a/hysop/tools/tests/test_work_arrays.py b/hysop/tools/tests/test_work_arrays.py
index dfb07f0e8..e6f79b488 100644
--- a/hysop/tools/tests/test_work_arrays.py
+++ b/hysop/tools/tests/test_work_arrays.py
@@ -8,6 +8,7 @@ from hysop import Field, Box
from hysop.operators import Stretching, EnergyEnstrophy, Curl
from hysop.tools.misc import WorkSpaceTools
import hysop.tools.numpywrappers as npw
+from hysop.constants import HYSOP_REAL, HYSOP_INTEGER
cos = np.cos
sin = np.sin
@@ -53,6 +54,70 @@ def init():
return op, simu
+def test_build_work_1():
+ """Build work array. Input = nb of required arrays
+ and a common shape for all of them.
+ """
+
+ # First case : create a work list
+ inshape = (12,)
+ work = WorkSpaceTools.check_work_array(3, inshape)
+ assert len(work) == 3
+ for wk in work:
+ assert wk.shape == inshape
+ assert wk.dtype == HYSOP_REAL
+ # Second case : 'reuse' memory from a predefined work array
+ inshape2 = (3, 4)
+ work2 = WorkSpaceTools.check_work_array(2, inshape2, work)
+ assert len(work2) == 2
+ for i in xrange(len(work2)):
+ assert work2[i].shape == inshape2
+ # ensure that work[i] and work2[i] use the same memory
+ assert npw.arrays_share_data(work2[i], work[i])
+
+
+def test_build_work_2():
+ """Build work array. Input = nb of required arrays
+ and list of shapes for those arrays.
+ """
+ inshapes = [(12,), (27,), (30,)]
+ work = WorkSpaceTools.check_work_array(3, inshapes)
+ assert len(work) == 3
+ for i in xrange(len(work)):
+ assert work[i].shape == inshapes[i]
+ assert work[i].dtype == HYSOP_REAL
+
+ # Second case : 'reuse' memory from a predefined work array
+ inshapes2 = [(3, 2, 2), (3, 8)]
+ work2 = WorkSpaceTools.check_work_array(2, inshapes2, work)
+ assert len(work2) == 2
+ for i in xrange(len(work2)):
+ assert work2[i].shape == inshapes2[i]
+ # ensure that work[i] and work2[i] use the same memory
+ assert npw.arrays_share_data(work2[i], work[i])
+
+
+def test_build_work_3():
+ """Same as above but for integers
+ """
+ inshapes = [(12,), (27,), (30,)]
+ work = WorkSpaceTools.check_work_array(3, inshapes,
+ data_type=HYSOP_INTEGER)
+ assert len(work) == 3
+ for i in xrange(len(work)):
+ assert work[i].shape == inshapes[i]
+ assert work[i].dtype == HYSOP_INTEGER
+
+ # Second case : 'reuse' memory from a predefined work array
+ inshapes2 = [(3, 2, 2), (3, 8)]
+ work2 = WorkSpaceTools.check_work_array(2, inshapes2, work)
+ assert len(work2) == 2
+ for i in xrange(len(work2)):
+ assert work2[i].shape == inshapes2[i]
+ # ensure that work[i] and work2[i] use the same memory
+ assert npw.arrays_share_data(work2[i], work[i])
+
+
def test_default_work():
"""Default case : each op has its own work (or no work at all)
"""
@@ -122,9 +187,3 @@ def test_common_work():
ll = min(lref, lw)
for i in xrange(ll):
assert npw.arrays_share_data(wkref[i], w[i])
-
-
-if __name__ == "__main__":
- test_default_work()
- test_no_common_work()
- test_common_work()
--
GitLab