From 8ac30eed4313dff25926ffce735f30cbf0ccb377 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Franck=20P=C3=A9rignon?= <franck.perignon@imag.fr>
Date: Mon, 30 May 2016 16:28:08 +0200
Subject: [PATCH] Fix stretching + rewiew working arrays
---
hysop/default_methods.py | 5 +-
hysop/operator/adapt_timestep.py | 174 +++++---
hysop/operator/advection.py | 5 +-
hysop/operator/advection_dir.py | 55 +--
hysop/operator/analytic.py | 3 +
hysop/operator/computational.py | 31 +-
hysop/operator/density.py | 4 +
hysop/operator/diffusion.py | 6 +-
hysop/operator/discrete/adapt_timestep.py | 386 +++++++++---------
hysop/operator/discrete/differential.py | 8 +-
hysop/operator/discrete/poisson_fft.py | 37 +-
hysop/operator/discrete/stretching.py | 18 +-
hysop/operator/hdf_io.py | 3 +
hysop/operator/multiphase_gradp.py | 3 +
hysop/operator/multiresolution_filter.py | 46 +--
hysop/operator/penalization.py | 3 +
hysop/operator/reprojection.py | 3 +
hysop/operator/stretching.py | 4 +-
hysop/operator/tests/test_Stretching.py | 231 ++++++-----
.../operator/tests/test_adaptive_time_step.py | 337 +++++++++++----
hysop/operator/velocity_correction.py | 3 +
hysop/testsenv.py | 2 +-
hysop/tools/misc.py | 1 -
23 files changed, 838 insertions(+), 530 deletions(-)
diff --git a/hysop/default_methods.py b/hysop/default_methods.py
index d1e20b5c0..012ee19c3 100644
--- a/hysop/default_methods.py
+++ b/hysop/default_methods.py
@@ -2,7 +2,7 @@
"""
from hysop.methods_keys import TimeIntegrator, Interpolation, GhostUpdate,\
Remesh, Support, Splitting, MultiScale, Formulation, SpaceDiscretisation, \
- dtCrit, Precision
+ Precision
from hysop.constants import HYSOP_REAL
from hysop.numerics.odesolvers import RK2, RK3
from hysop.numerics.interpolation import Linear
@@ -21,8 +21,7 @@ DIFFERENTIAL = {SpaceDiscretisation: FDC4, GhostUpdate: True}
"""Differential operators (finite differencies based) default setup
"""
-ADAPT_TIME_STEP = {TimeIntegrator: RK3, SpaceDiscretisation: FDC4,
- dtCrit: 'vort'}
+ADAPT_TIME_STEP = {TimeIntegrator: RK3, SpaceDiscretisation: FDC4}
"""Adaptative time step default setup
"""
diff --git a/hysop/operator/adapt_timestep.py b/hysop/operator/adapt_timestep.py
index 1e0750355..4d76f3516 100755
--- a/hysop/operator/adapt_timestep.py
+++ b/hysop/operator/adapt_timestep.py
@@ -1,67 +1,105 @@
# -*- coding: utf-8 -*-
-"""
-@file operator/adapt_timestep.py
+"""Update time-step, depending on the flow field.
-Definition of the adaptative time step according to the flow fields.
+See :ref:`adaptive_time_step` for details.
"""
from hysop.constants import debug
-from hysop.methods_keys import TimeIntegrator, SpaceDiscretisation,\
- dtCrit
-from hysop.numerics.finite_differences import FDC4
-from hysop.operator.discrete.adapt_timestep import AdaptTimeStep_D
+from hysop.methods_keys import TimeIntegrator, SpaceDiscretisation
+from hysop.operator.discrete.adapt_timestep import AdaptiveTimeStepD
from hysop.operator.continuous import opsetup
from hysop.operator.computational import Computational
import hysop.default_methods as default
from hysop.mpi import main_comm, MPI
+from hysop.numerics.differential_operations import MaxDiagGradV, \
+ DiagAndStretch, StrainCriteria, StretchLike, StrainAndStretch
-class AdaptTimeStep(Computational):
- """
- The adaptative Time Step is computed according
+class AdaptiveTimeStep(Computational):
+ """The adaptative Time Step is computed according
to the following expression :
+
dt_adapt = min (dt_advection, dt_stretching, dt_cfl)
+
"""
+ authorized_criteria = ['vort', 'gradU', 'stretch', 'cfl', 'strain']
+
@debug
def __init__(self, velocity, vorticity, simulation,
- time_range=None, lcfl=0.125, cfl=0.5,
- maxdt=9999., **kwds):
+ lcfl=0.125, cfl=0.5, time_range=None,
+ maxdt=9999., criteria=None, **kwds):
"""
- Create a time_step-evaluation operator from given
- velocity and vorticity variables.
-
- Note : If cfl is None, the computation of the adaptative time step
- is only based on dt_advection and dt_stretching, taking the minimum
- value of the two.
-
- @param velocity field
- @param vorticity field
- @param dt_adapt : adaptative timestep variable
- @param time_range : [start, end] use to define a 'window' in which
- the current operator is applied. Outside start-end, this operator
- has no effect. Start/end are iteration numbers.
- Default = [2, endofsimu]
+ Parameters
+ ----------
+ velocity, vorticity : :class:`~hysop.fields.continuous.DiscreteFields`
+ fields used to update the time steps. Read only.
+ simulation : :class:`hysop.fields.problem.simulation.Simulation`
+ time loop parameters description
+ lcfl : double, optional
+ the lagrangian CFL coefficient used for advection stability
+ cfl : double, optional
+ CFL coefficient
+ criteria : list of strings, optional
+ name of the criteria used to compute the new time step.
+ See notes below.
+ time_range : list of two integers, optional
+ use to define a 'window' in which the current operator is applied.
+ time_range =Â [start, end], outside this range, the operator
+ has no effect. Start/end are iteration numbers.
+ Default = [2, end of simu]
+ maxdt : double, optional
+ maximum value allowed for time step. Default = 9999.
+ kwds : arguments passed to base class.
+
+ Notes
+ -----
+ * This operator has no effect on input variables.
+ * Authorized criteria are:
+ * criteria for the advection scheme, one of :
+ 1. 'vort' (default)
+ 2. 'gradU'
+ 3. 'strain'
+ * criteria for the stretching scheme :
+ 4. 'stretch'
+ * cfl-like :
+ 5. 'cfl'
+ * Computes a 'diagnostics' vector :
+ diagnostics = (time, time_step, c1, c2, ...)
+ ci = time step computed from the criterium number i in the list
+ above.
+ * dt_adapt = min(diagnostics, dtmax)
"""
assert 'variables' not in kwds, 'variables parameter is useless.'
- super(AdaptTimeStep, self).__init__(variables=[velocity, vorticity],
- **kwds)
+ super(AdaptiveTimeStep, self).__init__(variables=[velocity, vorticity],
+ **kwds)
if self.method is None:
self.method = default.ADAPT_TIME_STEP
assert SpaceDiscretisation in self.method.keys()
assert TimeIntegrator in self.method.keys()
- if dtCrit not in self.method.keys():
- self.method[dtCrit] = 'vort'
-
- ## velocity variable (vector)
+ # Definition of criterion for dt_advec computation
+ if criteria is None:
+ criteria = ['vort']
+ msg = 'criteria arg must be a list.'
+ assert isinstance(criteria, list), msg
+ msg = 'Unknow criteria : '
+ for cr in criteria:
+ assert cr in self.authorized_criteria, msg + cr
+ msg = 'criteria list is empty!'
+ assert len(criteria) >= 1, msg
+ unknowns = [k for k in criteria if k not in self.authorized_criteria]
+ msg = 'Unknown criteria ' + str(unknowns)
+ assert len(unknowns) == 0, msg
+
+ # list of criteria used to compute dt
+ self.criteria = criteria
+
+ # velocity variable (vector)
self.velocity = velocity
- ## vorticity variable (vector)
+ # vorticity variable (vector)
self.vorticity = vorticity
- ## adaptative time step variable ("variable" object)
+ # simulation
self.simulation = simulation
- #assert isinstance(self.dt_adapt, VariableParameter)
- # Check if 'dt' key is present in dt_adapt dict
- #assert 'dt' in self.dt_adapt.data
self.input = self.variables
self.output = []
@@ -71,8 +109,7 @@ class AdaptTimeStep(Computational):
self._set_inter_comm()
def _set_inter_comm(self):
- """
- Create intercommunicators, if required (i.e. if there are several
+ """Create intercommunicators, if required (i.e. if there are several
tasks defined in the domain).
"""
task_is_source = self._mpis.task_id == self.domain.current_task()
@@ -88,39 +125,41 @@ class AdaptTimeStep(Computational):
0, main_comm, rk)
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)
-
- vd = self.discreteFields[self.velocity]
- shape_v = vd[0].shape
- rwork_length = self.velocity.nb_components ** 2
- res = {'rwork': [], 'iwork': None}
- for _ in xrange(rwork_length):
- res['rwork'].append(shape_v)
- return res
+ super(AdaptiveTimeStep, self).get_work_properties()
+ diffop = None
+ topo = self.discreteFields[self.velocity].topology
+ if 'gradU' in self.criteria:
+ diffop = MaxDiagGradV
+ if 'stretch' in self.criteria:
+ diffop = StretchLike
+ if 'strain' in self.criteria:
+ diffop = StrainCriteria
+ if 'stretch' in self.criteria and 'gradU' in self.criteria:
+ diffop = DiagAndStretch
+ if 'stretch' in self.criteria and 'strain' in self.criteria:
+ diffop = StrainAndStretch
+ if diffop is not None:
+ return diffop.get_work_properties(topo)
+ else:
+ return {'rwork': None, 'iwork': None}
def discretize(self):
- if self.method[SpaceDiscretisation] is FDC4:
- nbGhosts = 2
- else:
- raise ValueError("Unknown method for space discretization of the\
- time step-evaluation operator.")
- super(AdaptTimeStep, self)._standard_discretize(nbGhosts)
+ nb_ghosts = self.method[SpaceDiscretisation].ghosts_layer_size
+ super(AdaptiveTimeStep, self)._standard_discretize(nb_ghosts)
@debug
@opsetup
def setup(self, rwork=None, iwork=None):
self.discrete_op =\
- AdaptTimeStep_D(self.discreteFields[self.velocity],
- self.discreteFields[self.vorticity],
- self.simulation, method=self.method,
- time_range=self.time_range,
- lcfl=self.lcfl,
- cfl=self.cfl,
- maxdt=self.maxdt,
- rwork=rwork, iwork=iwork)
+ AdaptiveTimeStepD(self.discreteFields[self.velocity],
+ self.discreteFields[self.vorticity],
+ self.simulation, method=self.method,
+ time_range=self.time_range,
+ criteria=self.criteria,
+ lcfl=self.lcfl,
+ cfl=self.cfl,
+ maxdt=self.maxdt,
+ rwork=rwork, iwork=iwork)
# Output setup
self._set_io('dt_adapt', (1, 7))
self.discrete_op.set_writer(self._writer)
@@ -140,3 +179,8 @@ class AdaptTimeStep(Computational):
else:
dt = self._intercomms[rk].bcast(dt, root=0)
self.simulation.update_time_step(dt)
+
+ def diagnostics(self):
+ """Get the list of computed dt (for each criteria)
+ """
+ return self.discrete_op.diagnostics
diff --git a/hysop/operator/advection.py b/hysop/operator/advection.py
index fd5abc420..85883c51a 100644
--- a/hysop/operator/advection.py
+++ b/hysop/operator/advection.py
@@ -1,7 +1,6 @@
-"""
-@file advection.py
+"""Advection of a field.
-Advection of a field.
+See :ref:`advection` in user guide.
"""
from hysop.constants import debug, S_DIR, ZDIR
from hysop.operator.computational import Computational
diff --git a/hysop/operator/advection_dir.py b/hysop/operator/advection_dir.py
index 4f84e1cfd..513949dd1 100644
--- a/hysop/operator/advection_dir.py
+++ b/hysop/operator/advection_dir.py
@@ -193,36 +193,45 @@ class AdvectionDir(Computational):
[ interp , part_positions, fields_on_particles]
interp part is used also for remesh and time-integrator.
"""
- if not self._is_discretized:
- msg = 'The operator must be discretized '
- msg += 'before any call to this function.'
- raise RuntimeError(msg)
- dimension = self.domain.dimension
- res = {'rwork': [], 'iwork': []}
- if self.method[Support].find('gpu') < 0:
- tiw = self.method[TimeIntegrator].getWorkLengths(1)
- iw, iiw = \
- self.method[Interpolation].getWorkLengths(domain_dim=dimension)
- rw, riw = Remeshing.getWorkLengths(domain_dim=dimension)
- iwl = max(iiw, riw)
- rw = max(tiw + iw, rw)
- else:
- # For GPU version, no need of numerics works
- iwl, rw = 0, 0
+ assert self._is_discretized
+
# Shape of reference comes from fields, not from velocity
advected_discrete_fields = [self.discreteFields[v]
for v in self.variables
if v is not self.velocity]
+ topo = advected_discrete_fields[0].topology
+ # number of components of the rhs (time integration)
+ rhs_size = 1
+
+ if self.method[Support].find('gpu') < 0:
+ # -- pure python advection --
+ # work array shape depends on the time integrator
+ # interpolation scheme and remeshing scheme
+ ti_work = self.method[TimeIntegrator].get_work_properties(
+ rhs_size, topo)
+ ti_rwork_length = len(ti_work['rwork'])
+ iw_prop = self.method[Interpolation].get_work_properties(topo)
+ rw_prop = Remeshing.get_work_properties(topo)
+ interp_iwork_length = len(iw_prop['iwork'])
+ interp_rwork_length = len(iw_prop['rwork'])
+ remesh_iwork_length = len(rw_prop['iwork'])
+ remesh_rwork_length = len(rw_prop['rwork'])
+ iwork_length = max(interp_iwork_length, remesh_iwork_length)
+ rwork_length = max(ti_rwork_length + interp_rwork_length,
+ remesh_rwork_length)
+ else:
+ # -- GPU advection --
+ # no work array
+ iwork_length, rwork_length = 0, 0
+
memshape = advected_discrete_fields[0].topology.mesh.resolution
- rw += np.sum([f.nb_components for f in self.advected_fields])
+ rwork_length += np.sum([f.nb_components for f in self.advected_fields])
if self.method[Support].find('gpu') < 0 or \
self.method[Support].find('gpu_2k') >= 0:
- rw += 1 # positions
- for i in xrange(rw):
- res['rwork'].append(memshape)
- for i in xrange(iwl):
- res['iwork'].append(memshape)
- return res
+ rwork_length += 1 # work array for position
+ memsize = np.prod(topo.mesh.resolution)
+ return {'rwork': [(memsize,)] * rwork_length,
+ 'iwork': [(memsize,)] * iwork_length}
@debug
@opapply
diff --git a/hysop/operator/analytic.py b/hysop/operator/analytic.py
index 4e605264c..9e3cf0c17 100644
--- a/hysop/operator/analytic.py
+++ b/hysop/operator/analytic.py
@@ -58,3 +58,6 @@ class Analytic(Computational):
def get_profiling_info(self):
pass
+
+ def get_work_properties(self):
+ return {'rwork': None, 'iwork': None}
diff --git a/hysop/operator/computational.py b/hysop/operator/computational.py
index ee6db6c9f..186e27489 100755
--- a/hysop/operator/computational.py
+++ b/hysop/operator/computational.py
@@ -87,6 +87,7 @@ class Computational(Operator):
# Turn to true after self._discretize_vars call.
self._is_discretized = False
+ @abstractmethod
def get_work_properties(self):
"""Get properties of internal work arrays. Must be call after discretize
but before setup.
@@ -109,7 +110,9 @@ class Computational(Operator):
shape (12,12) and (45, 12, 33) and no integer arrays
"""
- return {'rwork': None, 'iwork': None}
+ msg = 'The operator must be discretized '
+ msg += 'before any call to this function.'
+ assert self._is_discretized, msg
def discretize(self):
"""
@@ -137,7 +140,10 @@ class Computational(Operator):
Check variables and discretization parameters
Set single_topo: if true all fields are discretized with the
same topo
- @return build_topos : a dict (key == field), if build_topos[v] is true,
+
+ Returns
+ -------
+ build_topos : a dict (key == field), if build_topos[v] is true,
a topology must be built for v. In that case, the discretization has
been saved in self.variables[v] during init. In the other case
self.variables[v] is the required topology.
@@ -210,8 +216,9 @@ class Computational(Operator):
self.variables[v] = topo
else:
# Topo is already built, just check its ghosts
- topo = self.variables.values()[0].mesh.discretization.ghosts
- assert (topo >= min_ghosts).all()
+ msg = 'The proposed ghost layer is not large enough.'
+ ghosts = self.variables.values()[0].mesh.discretization.ghosts
+ assert (ghosts >= min_ghosts).all(), msg
else:
# ... or one topo for each field.
@@ -219,13 +226,18 @@ class Computational(Operator):
if build_topos[v]:
self.variables[v] = self._build_topo(self.variables[v],
min_ghosts)
+ build_topos[v] = False
else:
- assert (self.variables[v].ghosts >= min_ghosts).all()
+ assert (self.variables[v].ghosts() >= min_ghosts).all()
# All topos are built, we can discretize fields.
self._discretize_vars()
def _build_topo(self, discretization, min_ghosts):
+ """Build an mpi topology and its mesh from a given
+ discretization.
+
+ """
# Reset ghosts if necessary
ghosts = discretization.ghosts
ghosts[ghosts < min_ghosts] = min_ghosts
@@ -328,8 +340,8 @@ class Computational(Operator):
self.time_info = self.discrete_op.time_info
else:
from hysop.mpi import main_rank
- shortName = str(self.__class__).rpartition('.')[-1][0:-2]
- s = '[' + str(main_rank) + '] ' + shortName
+ short_name = str(self.__class__).rpartition('.')[-1][0:-2]
+ s = '[' + str(main_rank) + '] ' + short_name
s += " : operator not discretized --> no computation, time = 0."
print s
@@ -344,13 +356,14 @@ class Computational(Operator):
"""
Common printings for operators
"""
- shortName = str(self.__class__).rpartition('.')[-1][0:-2]
+ short_name = str(self.__class__).rpartition('.')[-1][0:-2]
if self.discrete_op is not None:
s = str(self.discrete_op)
else:
- s = shortName + " operator. Not discretised."
+ s = short_name + " operator. Not discretised."
return s + "\n"
def get_profiling_info(self):
+ """Update profiler"""
if self.discrete_op is not None:
self.profiler += self.discrete_op.profiler
diff --git a/hysop/operator/density.py b/hysop/operator/density.py
index f30d56e5d..bca083361 100644
--- a/hysop/operator/density.py
+++ b/hysop/operator/density.py
@@ -43,3 +43,7 @@ class DensityVisco(Computational):
viscosity=self.discreteFields[self.viscosity],
method=self.method)
self._is_uptodate = True
+
+ def get_work_properties(self):
+ return {'rwork': None, 'iwork': None}
+
diff --git a/hysop/operator/diffusion.py b/hysop/operator/diffusion.py
index cbf96d9ab..c72f69564 100644
--- a/hysop/operator/diffusion.py
+++ b/hysop/operator/diffusion.py
@@ -11,6 +11,7 @@ from hysop.constants import debug
from hysop.operator.continuous import opsetup
from hysop.methods_keys import SpaceDiscretisation
import hysop.default_methods as default
+from hysop import __FFTW_ENABLED__
class Diffusion(Computational):
@@ -54,7 +55,7 @@ class Diffusion(Computational):
self.output = [self.vorticity]
def discretize(self):
- if self.method[SpaceDiscretisation] is 'fftw':
+ if self.method[SpaceDiscretisation] is 'fftw' and __FFTW_ENABLED__:
super(Diffusion, self)._fftw_discretize()
elif self.method[SpaceDiscretisation] is 'fd':
super(Diffusion, self)._standard_discretize()
@@ -78,3 +79,6 @@ class Diffusion(Computational):
viscosity=self.viscosity,
**kw)
self._is_uptodate = True
+
+ def get_work_properties(self):
+ return {'rwork': None, 'iwork': None}
diff --git a/hysop/operator/discrete/adapt_timestep.py b/hysop/operator/discrete/adapt_timestep.py
index 85e1ade95..ac32343fa 100755
--- a/hysop/operator/discrete/adapt_timestep.py
+++ b/hysop/operator/discrete/adapt_timestep.py
@@ -1,69 +1,103 @@
# -*- coding: utf-8 -*-
-"""
-Evaluation of the adaptative time step according to the flow fields.
+"""Discrete operator to Compute the time step according to the flow fields.
"""
from hysop.constants import debug
-from hysop.methods_keys import TimeIntegrator, SpaceDiscretisation,\
- dtCrit
+from hysop.methods_keys import TimeIntegrator, SpaceDiscretisation
from hysop.operator.discrete.discrete import DiscreteOperator
-from hysop.numerics.differential_operations import GradV
+from hysop.numerics.differential_operations import MaxDiagGradV, \
+ DiagAndStretch, StrainCriteria, StretchLike, StrainAndStretch
import hysop.tools.numpywrappers as npw
from hysop.numerics.update_ghosts import UpdateGhosts
from hysop.mpi import MPI
from hysop.constants import np, HYSOP_MPI_REAL
from hysop.tools.profiler import profile
+from hysop.problem.simulation import Simulation
+from hysop.tools.misc import WorkSpaceTools
-class AdaptTimeStep_D(DiscreteOperator):
+class AdaptiveTimeStepD(DiscreteOperator):
"""
- The adaptative Time Step is computed according
- to the following expression :
+ The adaptative Time Step is computed as:
+
dt_adapt = min (dt_advection, dt_stretching, dt_cfl)
"""
+ authorized_criteria = ['vort', 'gradU', 'stretch', 'cfl', 'deform']
+
@debug
def __init__(self, velocity, vorticity, simulation,
- lcfl=0.125, cfl=0.5, time_range=None, maxdt=9999., **kwds):
+ lcfl=0.125, cfl=0.5, criteria=None,
+ time_range=None, maxdt=9999., **kwds):
"""
- @param velocity : discretization of the velocity field
- @param vorticity : discretization of the vorticity field
- @param dt_adapt : adaptative timestep
- (a hysop.variable_parameter.VariableParameter)
- @param lcfl : the lagrangian CFL coefficient used
- for advection stability
- @param cfl : the CFL coefficient.
- @param time_range : [start, end] use to define a 'window' in which
- the current operator is applied. Outside start-end, this operator
- has no effect. Start/end are iteration numbers.
- Default = [2, endofsimu]
+ Parameters
+ ----------
+ velocity, vorticity : :class:`~hysop.fields.discrete.DiscreteFields`
+ discrete fields used to update the time steps. Read only.
+ simulation : :class:`~hysop.problem.simulation.Simulation`
+ lcfl : double, optional
+ the lagrangian CFL coefficient used for advection stability
+ cfl : double, optional
+ CFL coefficient
+ criteria : list of strings, optional
+ name of the criteria used to compute the new time step.
+ See notes below.
+ time_range : list of two integers, optional
+ use to define a 'window' in which the current operator is applied.
+ time_range =Â [start, end], outside this range, the operator
+ has no effect. Start/end are iteration numbers.
+ Default = [2, end of simu]
+ maxdt : double, optional
+ maximum value allowed for time step. Default = 9999.
+ kwds : arguments passed to base class.
+
+ Notes
+ -----
+ * This operator has no effect on input variables.
+ * Authorized criteria are:
+ * criteria for the advection scheme, one of :
+ 1. 'vort' (default)
+ 2. 'gradU'
+ 3. 'strain'
+ * criteria for the stretching scheme :
+ 4. 'stretch'
+ * cfl-like :
+ 5. 'cfl'
+ * Computes a 'diagnostics' vector :
+ diagnostics = (time, time_step, c1, c2, ...)
+ ci = time step computed from the criterium number i in the list
+ above.
"""
- ## velocity discrete field
+ # velocity discrete field
self.velocity = velocity
- ## vorticity discrete field
+ # vorticity discrete field
self.vorticity = vorticity
- ## adaptative time step variable
- from hysop.problem.simulation import Simulation
+ # adaptative time step variable
assert isinstance(simulation, Simulation)
self.simulation = simulation
+ # Build the user required function list
+ self._used_functions = {}
+ # local buffer used to compute max values.
+ self._result = None
+ # position in result (only needed when stretch criterium is used)
+ self._pos = 0
+ # name of (optional) differential operator required to compute
+ # the criteria.
+ self._diffop_name = self._init_differential_operator(criteria)
assert 'variables' not in kwds, 'variables parameter is useless.'
- super(AdaptTimeStep_D, self).__init__(variables=[velocity, vorticity],
- **kwds)
+ super(AdaptiveTimeStepD, self).__init__(
+ variables=[velocity, vorticity], **kwds)
self.input = self.variables
self.output = []
- ## Courant Fredrich Levy coefficient
+ # Courant Fredrich Levy coefficient
self.cfl = cfl
- ## Lagrangian CFL coefficient
+ # Lagrangian CFL coefficient
self.lcfl = lcfl
- ## Max. timestep
+ # Max. timestep
self.maxdt = maxdt
- # Definition of criterion for dt_advec computation
- self.dtCrit = self.method[dtCrit]
- if not isinstance(self.dtCrit, list):
- self.dtCrit = [self.dtCrit]
- ## Time range
+ # Time range
if time_range is None:
time_range = [2, np.infty]
self.time_range = time_range
@@ -72,188 +106,168 @@ class AdaptTimeStep_D(DiscreteOperator):
# [time, dt, d1, d2, d3, d4, d5]
# for d1...d5 see computation details in apply.
self.diagnostics = npw.zeros((7))
- self._t_diagnostics = npw.zeros_like(self.diagnostics)
+ # list of indices of diagnostics of interest.
+ self._diag_ind = []
+ # True if ghost points synchro is required
+ self._needs_synchro = False
+ # The (optional) differential operation used to
+ # compute the criterium
+ self._diff_op = None
+ # Coefficient used to compute streching stability criterium
+ self.coeff_stretch = self.method[TimeIntegrator].stability_coeff()
+ # Init functions ...
+ if self._diffop_name is not None:
+ self._diff_op = self._diffop_name(
+ topo=self.velocity.topology,
+ method=self.method[SpaceDiscretisation],
+ work=self._rwork)
+ self._needs_synchro = True
+ # Ghost synchronisation operator.
+ if 'vort' in criteria and len(criteria) == 1:
+ self._synchronize = UpdateGhosts(self.vorticity.topology,
+ self.vorticity.nb_components)
+ self._synchr_buff = self.vorticity.data
+ elif 'vort' in criteria:
+ self._synchronize = UpdateGhosts(
+ self.velocity.topology,
+ self.vorticity.nb_components + self.velocity.nb_components)
+ self._synchr_buff = self.velocity.data + self.vorticity.data
- # All diagnostcs function definition:
- # (Index in self.diagnostics, function, is gradU needed)
- self._all_functions = {
- 'gradU': (2, self._compute_gradU, True),
- 'stretch': (3, self._compute_stretch, True),
- 'cfl': (4, self._compute_cfl, False),
- 'vort': (5, self._compute_vort, False),
- 'deform': (6, self._compute_deform, True),
- }
-
- # Build the user required function list
- self._used_functions = []
- self._is_gradU_needed = False
- for crit in self.dtCrit:
- self._is_gradU_needed = \
- self._is_gradU_needed or self._all_functions[crit][-1]
- self._used_functions.append(self._all_functions[crit])
- assert len(self._used_functions) >= 1, "You must specify at least " + \
- "one criterion among: " + str(self._all_functions.keys())
-
- # Definition of dt:
- self.get_all_dt = []
- self._prepare_dt_list()
-
- # prepare ghost points synchro for velocity
- self._synchronize = UpdateGhosts(self.velocity.topology,
- self.velocity.nb_components)
- # gradU function
- self._function = GradV(topo=self.velocity.topology,
- method=self.method[SpaceDiscretisation])
+ else:
+ self._synchronize = UpdateGhosts(
+ self.velocity.topology,
+ self.velocity.nb_components)
+ self._synchr_buff = self.velocity.data
def _set_work_arrays(self, rwork=None, iwork=None):
- memshape = self.velocity.data[0].shape
- worklength = self.velocity.nb_components ** 2
- # rwork is used to save gradU
- if rwork is None:
- self._rwork = [npw.zeros(memshape) for _ in xrange(worklength)]
-
- else:
- assert isinstance(rwork, list), 'rwork must be a list.'
- self._rwork = rwork
- assert len(self._rwork) == worklength
- for wk in self._rwork:
- assert wk.shape == memshape
+ """Allocate local work spaces (if needed)
+ """
+ if self._diffop_name is not None:
+ topo = self.velocity.topology
+ wkp = self._diffop_name.get_work_properties(topo)
+ lwork = len(wkp['rwork'])
+ subshape = wkp['rwork'][0]
+ self._rwork = WorkSpaceTools.check_work_array(lwork,
+ subshape, rwork)
- @staticmethod
- def _compute_stability_coeff(timeint):
- 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
- # Definition of stability coefficient for stretching operator
- coef_stretch = 0.0
- classtype = timeint.mro()[0]
- if classtype is Euler:
- coef_stretch = 2.0
- elif classtype is RK2:
- coef_stretch = 2.0
- elif classtype is RK3:
- coef_stretch = 2.5127
- elif classtype is RK4:
- coef_stretch = 2.7853
- return coef_stretch
+ def _init_differential_operator(self, criteria):
+ """Select the required differential operator (if any)
+ depending on the chosen criteria
- def _prepare_dt_list(self):
- # definition of dt_advection
- if 'gradU' in self.dtCrit:
- # => based on gradU
- self.get_all_dt.append(
- lambda diagnostics: self.lcfl / diagnostics[2])
- if 'deform' in self.dtCrit:
- # => based on the deformations
- self.get_all_dt.append(
- lambda diagnostics: self.lcfl / diagnostics[6])
- if 'vort' in self.dtCrit:
- # => based on the vorticity
- self.get_all_dt.append(
- lambda diagnostics: self.lcfl / diagnostics[5])
- if 'stretch' in self.dtCrit:
- coeff_stretch = self._compute_stability_coeff(
- self.method[TimeIntegrator])
- self.get_all_dt.append(
- lambda diagnostics: coeff_stretch / diagnostics[3])
- if 'cfl' in self.dtCrit:
- h = self.velocity.topology.mesh.space_step[0]
- self.get_all_dt.append(
- lambda diagnostics: (self.cfl * h) / diagnostics[4])
+ Notes:
+ * required to check/allocate internal work arrays
+ * required to connect the functions needed to compute the criteria
+ * MUST be called before parent class init
+ """
+ # Dictionnary of the available diagnostics functions
+ # all_functions[name] =
+ # (Index in self.diagnostics, function)
+ all_functions = {
+ 'vort': (2, self._compute_vort),
+ 'gradU': (3, self._compute_grad_u),
+ 'strain': (4, self._compute_strain),
+ 'stretch': (5, self._compute_stretch),
+ 'cfl': (6, self._compute_cfl)
+ }
- def _gradU(self):
- # Synchronize ghost points of velocity
- self._synchronize(self.velocity.data)
- # gradU computation
- self._rwork = self._function(self.velocity.data, self._rwork)
+ # select functions of interest, from input criteria list.
+ self._used_functions = {all_functions[k][0]: all_functions[k][1]
+ for k in criteria if k in all_functions}
+ diffop = None
+ dimension = self.velocity.domain.dimension
+ if 'gradU' in criteria:
+ diffop = MaxDiagGradV
+ self._result = npw.zeros(dimension)
+ if 'stretch' in criteria:
+ diffop = StretchLike
+ self._result = npw.zeros(dimension)
+ self._pos = 0
+ if 'strain' in criteria:
+ diffop = StrainCriteria
+ self._result = npw.zeros(dimension)
+ if 'stretch' in criteria and 'gradU' in criteria:
+ diffop = DiagAndStretch
+ self._result = npw.zeros(2 * dimension)
+ self._pos = dimension
+ if 'stretch' in criteria and 'strain' in criteria:
+ diffop = StrainAndStretch
+ self._result = npw.zeros(2 * dimension)
+ self._pos = dimension
+ return diffop
- def _compute_gradU(self):
- res = 0.
- nb_components = self.velocity.nb_components
- for direc in xrange(nb_components):
- # maxima of partial derivatives of velocity :
- # needed for advection stability condition (1st option)
- res = max(res, np.max(abs(self._rwork[(nb_components + 1)
- * direc])))
- return res
+ def _compute_grad_u(self):
+ """Compute criterium ...
+ """
+ self._result = self._diff_op(self.velocity.data, self._result)
+ return self.lcfl / max(self._result[:self._dim])
def _compute_stretch(self):
- res = 0.
- nb_components = self.velocity.nb_components
- for direc in xrange(nb_components):
- # maxima of partial derivatives of velocity:
- # needed for stretching stability condition
- tmp = np.max(sum([abs(self._rwork[i])
- for i in xrange(nb_components * direc,
- nb_components * (direc + 1))]))
- res = max(res, tmp)
- return res
+ """Compute criterium ...
+ """
+ self._result = self._diff_op(self.velocity.data, self._result)
+ return self.coeff_stretch / max(
+ self._result[self._pos:self._pos + self._dim])
def _compute_cfl(self):
- # maxima of velocity : needed for CFL based time step
- return np.max([np.max(np.abs(v_c)) for v_c in self.velocity.data])
+ """maxima of velocity : needed for CFL based time step
+ """
+ coef = self.cfl * self.velocity.topology.mesh.space_step[0]
+ return coef / np.max([np.max(np.abs(v_c))
+ for v_c in self.velocity.data])
def _compute_vort(self):
- # maxima of vorticity :
- # needed for advection stability condition (2nd option)
- return np.max([np.max(np.abs(w_c)) for w_c in self.vorticity.data])
+ """maxima of vorticity :
+ needed for advection stability condition
+ """
+ return self.lcfl / np.max([np.max(np.abs(w_c))
+ for w_c in self.vorticity.data])
- def _compute_deform(self):
- # 1/2(gradU + gradU^T) computation
- self._rwork[1] += self._rwork[3]
- self._rwork[2] += self._rwork[6]
- self._rwork[5] += self._rwork[7]
- self._rwork[1] *= 0.5
- self._rwork[2] *= 0.5
- self._rwork[5] *= 0.5
- self._rwork[3][...] = self._rwork[1][...]
- self._rwork[6][...] = self._rwork[2][...]
- self._rwork[7][...] = self._rwork[5][...]
- # maxima of deformation tensor:
- # needed for advection stability condition (3rd option)
- res = 0.
- nb_components = self.velocity.nb_components
- for direc in xrange(nb_components):
- tmp = np.max(sum([abs(self._rwork[i])
- for i in xrange(nb_components * direc,
- nb_components * (direc + 1))
- ]))
- res = max(res, tmp)
- return res
+ def _compute_strain(self):
+ """maxima of strain tensor:
+ needed for advection stability condition
+ """
+ self._result = self._diff_op(self.velocity.data, self._result)
+ return self.lcfl / np.max(self._result[:self._dim])
@debug
@profile
def apply(self, simulation=None):
-# if simulation is not None:
-# assert self.simulation is simulation
-
- # current time
+ if simulation is not None:
+ assert self.simulation is simulation
+ ind = self.velocity.topology.mesh.compute_index
+ # current time and iteration number
time = self.simulation.time
iteration = self.simulation.current_iteration
- Nmax = min(self.simulation.max_iter, self.time_range[1])
- self.diagnostics[0] = time
- if iteration >= self.time_range[0] and iteration <= Nmax:
- if self._is_gradU_needed:
- self._gradU()
+ # upper value allowed for iterations
+ nmax = min(self.simulation.max_iter, self.time_range[1])
+ buff = npw.zeros(7)
+ buff[0] = time
+ if iteration >= self.time_range[0] and iteration <= nmax:
+ # Apply each function listed in used_functions
+ if self._needs_synchro:
+ # Synchronize ghost points of velocity
+ self._synchronize(self._synchr_buff)
+
for func in self._used_functions:
# func: (Index in self.diagnostics, function, is gradU needed)
- self.diagnostics[func[0]] = func[1]()
+ buff[func] = self._used_functions[func]()
+ # mpi reduction
self.velocity.topology.comm.Allreduce(
- sendbuf=[self.diagnostics, 7, HYSOP_MPI_REAL],
- recvbuf=[self._t_diagnostics, 7, HYSOP_MPI_REAL],
+ sendbuf=[buff, 7, HYSOP_MPI_REAL],
+ recvbuf=[self.diagnostics, 7, HYSOP_MPI_REAL],
op=MPI.MAX)
- self.diagnostics[...] = self._t_diagnostics
+ ind = self._used_functions.keys()
+ time_step = np.min(list(self.diagnostics[ind]) +
+ [self.maxdt])
+ self.diagnostics[0] = time
+ self.diagnostics[1] = time_step
- dt = np.min([dt(self.diagnostics) for dt in self.get_all_dt] +
- [self.maxdt])
- self.diagnostics[1] = dt
if self._writer is not None and self._writer.do_write(iteration):
self._writer.buffer[0, :] = self.diagnostics
self._writer.write()
# Update simulation time step with the new dt
- self.simulation.update_time_step(dt)
+ self.simulation.update_time_step(time_step)
# Warning this update is done only for the current MPI task!
# See wait function in base class.
diff --git a/hysop/operator/discrete/differential.py b/hysop/operator/discrete/differential.py
index 424a8491d..1046ac0c5 100644
--- a/hysop/operator/discrete/differential.py
+++ b/hysop/operator/discrete/differential.py
@@ -16,13 +16,12 @@ from hysop.numerics.differential_operations import Curl, GradV,\
from abc import ABCMeta, abstractmethod
from hysop.numerics.update_ghosts import UpdateGhosts
from hysop.methods_keys import SpaceDiscretisation
-try:
- from hysop.f2hysop import fftw2py
-except ImportError:
- from hysop.fakef2py import fftw2py
import hysop.default_methods as default
from hysop.tools.profiler import profile
from hysop.tools.misc import WorkSpaceTools
+from hysop import __FFTW_ENABLED__
+if __FFTW_ENABLED__:
+ from hysop.f2hysop import fftw2py
class Differential(DiscreteOperator):
@@ -70,6 +69,7 @@ class CurlFFT(Differential):
"""
def __init__(self, **kwds):
+ assert __FFTW_ENABLED__, "Run hysop with fftw to use this class."
super(CurlFFT, self).__init__(**kwds)
if self.domain.dimension == 3:
self._apply = self._apply_3d
diff --git a/hysop/operator/discrete/poisson_fft.py b/hysop/operator/discrete/poisson_fft.py
index 1cd215f14..3027058fd 100644
--- a/hysop/operator/discrete/poisson_fft.py
+++ b/hysop/operator/discrete/poisson_fft.py
@@ -1,18 +1,17 @@
# -*- coding: utf-8 -*-
-"""
-@file poisson_fft.py
-Discrete operator for Poisson problem (fftw based)
+"""Discrete operator for Poisson problem (fftw based)
"""
import hysop.tools.numpywrappers as npw
-try:
- from hysop.f2hysop import fftw2py
-except ImportError:
- from hysop.fakef2py import fftw2py
from hysop.operator.discrete.discrete import DiscreteOperator
from hysop.operator.discrete.reprojection import Reprojection
from hysop.constants import debug
from hysop.tools.profiler import profile
+from hysop import __FFTW_ENABLED__
+
+if __FFTW_ENABLED__:
+ from hysop.f2hysop import fftw2py
+
class PoissonFFT(DiscreteOperator):
@@ -43,13 +42,13 @@ class PoissonFFT(DiscreteOperator):
assert 'variables' not in kwds, 'variables parameter is useless.'
super(PoissonFFT, self).__init__(variables=[output_field, input_field],
**kwds)
- ## Solution field
+ # Solution field
self.output_field = output_field
- ## RHS field
+ # RHS field
self.input_field = input_field
- ## Solenoidal projection of input_field ?
+ # Solenoidal projection of input_field ?
self.projection = projection
- ## Filter size array = domainLength/(CoarseRes-1)
+ # Filter size array = domainLength/(CoarseRes-1)
self.filterSize = filterSize
# If 2D problem, input_field must be a scalar
self.dim = self.output_field.domain.dimension
@@ -60,7 +59,7 @@ class PoissonFFT(DiscreteOperator):
self.input = [self.input_field]
self.output = [self.output_field]
- ## The function called during apply
+ # The function called during apply
self.solve = None
# a sub function ...
self._solve = None
@@ -74,7 +73,7 @@ class PoissonFFT(DiscreteOperator):
f(output.nbComponents) + pb type 'pressure-poisson'
"""
- ## Multiresolution ?
+ # Multiresolution ?
multires = self.output_field.topology.mesh != self.input_field.topology.mesh
# connexion to the required apply function
@@ -109,6 +108,16 @@ class PoissonFFT(DiscreteOperator):
self.solve = self._solve_and_correct
else:
self.solve = self._solve
+
+ if not __FFTW_ENABLED__:
+ self.solve = self._solve_error
+
+ def _solve_error():
+ """Error message when fftw fortran interface is not available.
+ """
+ msg = "Warning: fftw fortran interface is not available,"
+ msg += "the fft solver does nothing."
+ raise NotImplementedError('msg')
def do_projection_with_op(self, simu):
self.projection.apply(simu)
@@ -197,7 +206,7 @@ class PoissonFFT(DiscreteOperator):
self.output_field.data[2], ghosts_w, ghosts_v)
def _solve_and_correct(self, simu):
- self._solve(simu)
+ self._solve(simu.current_iteration)
self.correction.apply(simu)
diff --git a/hysop/operator/discrete/stretching.py b/hysop/operator/discrete/stretching.py
index 8e3d2e16d..7e96a901d 100755
--- a/hysop/operator/discrete/stretching.py
+++ b/hysop/operator/discrete/stretching.py
@@ -155,7 +155,6 @@ class Conservative(Stretching):
@profile
def __init__(self, **kwds):
- # Right-hand side for time integration
def rhs(t, y, result):
"""rhs used in time integrator
"""
@@ -201,7 +200,7 @@ class GradUW(Stretching):
assert sum(subdt) == dt
for i in xrange(ndt):
- self._integrate(t, subdt[i])
+ self._integrate(subdt[i], t)
def _check_stability(self, dt):
"""Computes a stability condition depending on some
@@ -224,7 +223,7 @@ class GradUW(Stretching):
return nb_cycles, subdt
-class StretchingLinearized(Conservative):
+class StretchingLinearized(Stretching):
"""Conservative formulation of the linearized stretching.
"""
@@ -237,24 +236,26 @@ class StretchingLinearized(Conservative):
self.vorticity_BF.nb_components)
self.usual_op = usual_op
- # Right-hand side for time integration
def rhs(t, y, result, form):
+ """rhs used in time integrator
+ """
if form == "div(w:u)":
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)
-
+
def _integrate(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._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],
+ "div(w:u)")
# perform integration and save result in-place
self.vorticity.data = self.timeIntegrator(t, self.vorticity.data, dt,
result=self.vorticity.data)
@@ -263,7 +264,8 @@ class StretchingLinearized(Conservative):
# rhs of the integrator
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],
+ "div(u:w)")
# perform integration and save result in-place
self.vorticity.data = self.timeIntegrator(t, self.vorticity.data, dt,
result=self.vorticity.data)
diff --git a/hysop/operator/hdf_io.py b/hysop/operator/hdf_io.py
index 6936a3370..5ccd13a73 100755
--- a/hysop/operator/hdf_io.py
+++ b/hysop/operator/hdf_io.py
@@ -169,6 +169,9 @@ class HDF_IO(Computational):
Abstract interface to read/write process
"""
+ def get_work_properties(self):
+ return {'rwork': None, 'iwork': None}
+
class HDF_Writer(HDF_IO):
"""
diff --git a/hysop/operator/multiphase_gradp.py b/hysop/operator/multiphase_gradp.py
index a9411d346..d1415f0fb 100644
--- a/hysop/operator/multiphase_gradp.py
+++ b/hysop/operator/multiphase_gradp.py
@@ -75,3 +75,6 @@ class MultiphaseGradP(Computational):
def initialize_velocity(self):
self.discrete_op.initialize_velocity()
+
+ def get_work_properties(self):
+ return {'rwork': None, 'iwork': None}
diff --git a/hysop/operator/multiresolution_filter.py b/hysop/operator/multiresolution_filter.py
index 09127899c..01cc7d0a6 100644
--- a/hysop/operator/multiresolution_filter.py
+++ b/hysop/operator/multiresolution_filter.py
@@ -1,52 +1,45 @@
-"""
-@file operator/multiresolution_filter.py
-Filter values between grids of different resolution.
+"""Filter values between grids of different resolution.
+
"""
from hysop.constants import debug
from hysop.operator.continuous import opsetup
from hysop.operator.computational import Computational
import hysop.default_methods as default
from hysop.methods_keys import Support
-from hysop.tools.parameters import Discretization
class MultiresolutionFilter(Computational):
- """
- Apply a filter from a fine grid to a coarser grid.
+ """Interpolation from fine grid to coarse grid
+
"""
@debug
def __init__(self, d_in, d_out, **kwds):
"""
- Operator to apply Apply a filter from a fine grid to a coarser grid.
+ Parameters
+ ----------
+ d_in, d_out : :class:`hysop.mpi.topology.Cartesian`
+ or :class:`tools.parameters.Discretization`
+ data distribution for source (in) and target (out) grids.
+ kwds : base class parameters
+
"""
if 'method' not in kwds:
kwds['method'] = default.MULTIRESOLUTION_FILER
super(MultiresolutionFilter, self).__init__(**kwds)
self.d_in, self.d_out = d_in, d_out
+ self.input = self.variables
self.output = self.variables
+ self._df_in = []
+ self._df_out = []
def discretize(self):
super(MultiresolutionFilter, self)._standard_discretize()
- if isinstance(self.d_in, Discretization):
- topo_in = self._build_topo(self.d_in, 0)
- else:
- topo_in = self.d_in
- if isinstance(self.d_out, Discretization):
- topo_out = self._build_topo(self.d_out, 0)
- else:
- topo_out = self.d_out
- self._df_in = []
- self._df_out = []
+ topo_in = self._build_topo(self.d_in, 0)
+ topo_out = self._build_topo(self.d_out, 0)
for v in self.variables:
- if not topo_in in v.discreteFields.keys():
- self._df_in.append(v.discretize(topo_in))
- else:
- self._df_in.append(v.discreteFields[topo_in])
- if not topo_out in v.discreteFields.keys():
- self._df_out.append(v.discretize(topo_out))
- else:
- self._df_out.append(v.discreteFields[topo_out])
+ self._df_in.append(v.discretize(topo_in))
+ self._df_out.append(v.discretize(topo_out))
@opsetup
def setup(self, rwork=None, iwork=None):
@@ -61,3 +54,6 @@ class MultiresolutionFilter(Computational):
field_in=self._df_in, field_out=self._df_out,
method=self.method, rwork=rwork, iwork=iwork)
self._is_uptodate = True
+
+ def get_work_properties(self):
+ return {'rwork': None, 'iwork': None}
diff --git a/hysop/operator/penalization.py b/hysop/operator/penalization.py
index bb5e73c66..e9bf4653a 100644
--- a/hysop/operator/penalization.py
+++ b/hysop/operator/penalization.py
@@ -99,6 +99,9 @@ class Penalization(Computational):
self._is_uptodate = True
+ def get_work_properties(self):
+ return {'rwork': None, 'iwork': None}
+
class PenalizeVorticity(Penalization):
"""
diff --git a/hysop/operator/reprojection.py b/hysop/operator/reprojection.py
index d252314a2..eb52b3b64 100644
--- a/hysop/operator/reprojection.py
+++ b/hysop/operator/reprojection.py
@@ -51,3 +51,6 @@ class Reprojection(Computational):
True if projection must be done
"""
return self.discrete_op.do_projection(ite)
+
+ def get_work_properties(self):
+ return {'rwork': None, 'iwork': None}
diff --git a/hysop/operator/stretching.py b/hysop/operator/stretching.py
index 7c58412dc..c0e242cbd 100755
--- a/hysop/operator/stretching.py
+++ b/hysop/operator/stretching.py
@@ -103,8 +103,8 @@ class Stretching(Computational):
class StretchingLinearized(Stretching):
- """
- Solve the linearized stretching equation, i.e:
+ """Solve the linearized stretching equation, i.e:
+
\f{eqnarray*}
\frac{\partial \omega}{\partial t} &=& (\om \cdot \nabla)u_b +
(\om_b \cdot \nabla)u
diff --git a/hysop/operator/tests/test_Stretching.py b/hysop/operator/tests/test_Stretching.py
index ad76dca56..3a0f89b2f 100755
--- a/hysop/operator/tests/test_Stretching.py
+++ b/hysop/operator/tests/test_Stretching.py
@@ -1,21 +1,21 @@
# -*- coding: utf-8 -*-
-import hysop as pp
+"""Tests for stretching component resolution
+"""
+from hysop import Field, Box
import numpy as np
-from hysop.fields.continuous import Field
-from hysop.operator.stretching import Stretching, \
- StretchingLinearized
+from hysop.operator.stretching import Stretching, StretchingLinearized
from hysop.problem.simulation import Simulation
from hysop.methods_keys import TimeIntegrator, Formulation,\
SpaceDiscretisation
-from hysop.methods import Euler, RK3, FDC4, Conservative
+from hysop.methods import RK3, FDC4
from hysop.tools.parameters import Discretization
-import hysop.tools.numpywrappers as npw
+from hysop.tools.misc import WorkSpaceTools
pi = np.pi
cos = np.cos
sin = np.sin
-def computeVel(res, x, y, z, t):
+def compute_vel(res, x, y, z, t):
amodul = cos(pi * 1. / 3.)
pix = pi * x
piy = pi * y
@@ -32,7 +32,7 @@ def computeVel(res, x, y, z, t):
return res
-def computeVort(res, x, y, z, t):
+def compute_vort(res, x, y, z, t):
amodul = cos(pi * 1. / 3.)
pix = pi * x
piy = pi * y
@@ -41,16 +41,16 @@ def computeVort(res, x, y, z, t):
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))
+ (- 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))
+ (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))
+ (cos(pi2x) * sin(piy) * sin(piy) +
+ sin(pix) * sin(pix) * cos(pi2y))
return res
@@ -80,138 +80,143 @@ def computeVortBF(res, x, y, z, t):
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))
-
+ (- 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))
-
+ (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))
-
- return res
+ (cos(pi2x) * sin(piy) * sin(piy) +
+ sin(pix) * sin(pix) * cos(pi2y))
+ return res
-def test_stretching():
- # Parameters
+def init(method=None, work=False):
+ """Build, init, setup for operator
+ """
nb = 33
- boxLength = [1., 1., 1.]
- boxMin = [0., 0., 0.]
- nbElem = Discretization([nb, nb, nb], [2, 2, 2])
- time_step = 0.05
+ box_length = [1., 1., 1.]
+ box_min = [0., 0., 0.]
+ nb_elem = Discretization([nb, nb, nb], [2, 2, 2])
# Domain
- box = pp.Box(length=boxLength, origin=boxMin)
+ box = Box(length=box_length, origin=box_min)
# Fields
velo = Field(
- domain=box, formula=computeVel,
+ domain=box, formula=compute_vel,
name='Velocity', is_vector=True)
vorti = Field(
- domain=box, formula=computeVort,
+ domain=box, formula=compute_vort,
name='Vorticity', is_vector=True)
+ op = Stretching(velo, vorti, discretization=nb_elem, method=method)
+ op.discretize()
+ rwork = None
+ if work:
+ wkp = op.get_work_properties()['rwork']
+ 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)
+ return op, simulation
+
- # Operators
+def test_stretching():
+ """Default case
+ """
#method = {TimeIntegrator: RK3, Formulation: Conservative,
# SpaceDiscretisation: FDC4}
- stretch = Stretching(velo, vorti, discretization=nbElem)
- stretch.discretize()
- topo = stretch.discreteFields[velo].topology
- velo.initialize(topo=topo)
- vorti.initialize(topo=topo)
- stretch.setup()
- simulation = Simulation(start=0, end=1., time_step=time_step)
- stretch.apply(simulation)
+ op, simu = init()
+ op.apply(simu)
+
+
+def test_stretching_external_work():
+ """User-defined work arrays
+ """
+ op, simu = init(work=True)
+ op.apply(simu)
+
-def test_stretchingLinearized():
-
- # Parameters
+def test_stretching_graduw():
+ """GradUW formulation
+ """
+ method = {TimeIntegrator: RK3, Formulation: "GradUW",
+ SpaceDiscretisation: FDC4}
+ op, simu = init(method=method)
+ op.apply(simu)
+
+
+def init_linearized(method=None, work=False):
+ """Build, init, setup for operator
+ """
nb = 33
- boxLength = [1., 1., 1.]
- boxMin = [0., 0., 0.]
- nbElem = Discretization([nb, nb, nb], [2, 2, 2])
- time_step = 0.05
-
+ box_length = [1., 1., 1.]
+ box_min = [0., 0., 0.]
+ nb_elem = Discretization([nb, nb, nb], [2, 2, 2])
+
# Domain
- box = pp.Box(length=boxLength, origin=boxMin)
-
+ box = Box(length=box_length, origin=box_min)
+
# Fields
velo = Field(
- domain=box, formula=computeVel,
+ domain=box, formula=compute_vel,
name='Velocity', is_vector=True)
vorti = Field(
- domain=box, formula=computeVort,
+ domain=box, formula=compute_vort,
name='Vorticity', is_vector=True)
- veloBF = Field(
+ velo_bf = Field(
domain=box, formula=computeVelBF,
name='VelocityBF', is_vector=True)
- vortiBF = Field(
+ vorti_bf = Field(
domain=box, formula=computeVortBF,
name='VorticityBF', is_vector=True)
-
- # Operators
+
# Usual stretching operator
- stretch1 = Stretching(velo, vorti, discretization=nbElem)
- # Linearized stretching operator
- stretch2 = StretchingLinearized(velocity=velo,
- vorticity=vorti,
- velocity_BF=veloBF,
- vorticity_BF=vortiBF,
- discretization=nbElem)
+ stretch1 = Stretching(velo, vorti, discretization=nb_elem)
+ stretch2 = StretchingLinearized(velocity=velo, vorticity=vorti,
+ velocity_BF=velo_bf,
+ vorticity_BF=vorti_bf,
+ discretization=nb_elem, method=method)
stretch1.discretize()
stretch2.discretize()
+ rwork = None
+ if work:
+ wkp = stretch2.get_work_properties()['rwork']
+ rwork = WorkSpaceTools.check_work_array(len(wkp), wkp[0])
topo = stretch1.discreteFields[velo].topology
velo.initialize(topo=topo)
vorti.initialize(topo=topo)
- veloBF.initialize(topo=topo)
- vortiBF.initialize(topo=topo)
+ velo_bf.initialize(topo=topo)
+ vorti_bf.initialize(topo=topo)
stretch1.setup()
- stretch2.setup()
- simulation = Simulation(start=0, end=1., time_step=time_step)
- stretch1.apply(simulation)
- print 'norm vorti (usual):', vorti.norm(topo)
-
- stretch2.apply(simulation)
- print 'norm vorti (lin):', vorti.norm(topo)
-
-
-def test_stretching_external_work():
- # Parameters
- nb = 33
- boxLength = [1., 1., 1.]
- boxMin = [0., 0., 0.]
- nbElem = Discretization([nb, nb, nb], [2, 2, 2])
- time_step = 0.05
-
- # Domain
- box = pp.Box(length=boxLength, origin=boxMin)
-
- # Fields
- velo = Field(
- domain=box, formula=computeVel,
- name='Velocity', is_vector=True)
- vorti = Field(
- domain=box, formula=computeVort,
- name='Vorticity', is_vector=True)
-
- # Operators
- #method = {TimeIntegrator: RK3, Formulation: Conservative,
- # SpaceDiscretisation: FDC4}
- stretch = Stretching(velo, vorti, discretization=nbElem)
- stretch.discretize()
- wk_p = stretch.get_work_properties()
- rwork = []
- wk_length = len(wk_p['rwork'])
- for i in xrange(wk_length):
- memshape = wk_p['rwork'][i]
- rwork.append(npw.zeros(memshape))
-
- topo = stretch.discreteFields[velo].topology
- velo.initialize(topo=topo)
- vorti.initialize(topo=topo)
- stretch.setup(rwork=rwork)
- simulation = Simulation(start=0, end=1., time_step=time_step)
- stretch.apply(simulation)
+ stretch2.setup(rwork=rwork)
+ simulation = Simulation(start=0, end=1., time_step=0.05)
+ return stretch1, stretch2, simulation
+
+
+def test_stretching_linearized():
+ """Linearized formulation for stretching
+ """
+ str1, str2, simu = init_linearized()
+ str1.apply(simu)
+ str2.apply(simu)
+
+
+def test_stretching_external_work_graduv():
+ """User-defined work arrays for GradUW formulation
+ """
+ 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/hysop/operator/tests/test_adaptive_time_step.py b/hysop/operator/tests/test_adaptive_time_step.py
index 54284a67b..c3b7c3709 100755
--- a/hysop/operator/tests/test_adaptive_time_step.py
+++ b/hysop/operator/tests/test_adaptive_time_step.py
@@ -1,98 +1,291 @@
# -*- coding: utf-8 -*-
-import hysop as pp
-from hysop.operator.adapt_timestep import AdaptTimeStep
+from hysop import Field, Box
+from hysop.operator.adapt_timestep import AdaptiveTimeStep
from hysop.problem.simulation import Simulation
from hysop.tools.parameters import Discretization
from hysop.mpi import main_comm, main_size
import numpy as np
import hysop.tools.numpywrappers as npw
-import os
+from hysop.tools.misc import WorkSpaceTools
+from hysop.numerics.tests.test_differential_operations import compute_vel,\
+ compute_vort, diag_grad_func_3d, compute_max, grad_v_func_3d, strain_f_3d
+from math import pi
+from hysop.methods_keys import TimeIntegrator
+
sin = np.sin
cos = np.cos
+ldom = npw.asrealarray([pi * 2., 4. * pi, 2. * pi])
+xdom = [0., 0., 0.]
+
+d3d = Discretization([65, 123, 33], [2, 2, 2])
+
-d3d = Discretization([33, 33, 33], [2, 2, 2])
+def init(criteria=None, work=False):
+ """build fields and simu"""
+ box = Box(length=ldom, origin=xdom)
+ velo = Field(domain=box, formula=compute_vel,
+ name='Velocity', is_vector=True)
+ vorti = Field(domain=box, formula=compute_vort,
+ name='Vorticity', is_vector=True)
+ simu = Simulation(nb_iter=4)
+ simu.initialize()
+ assert simu.current_iteration == 0
+ op = AdaptiveTimeStep(velo, vorti, simulation=simu,
+ criteria=criteria,
+ discretization=d3d, lcfl=0.125, cfl=0.5,
+ time_range=[0, 4])
+ op.discretize()
+ rwork = None
+ if work is True:
+ wkp = op.get_work_properties()
+ lwork = len(wkp['rwork'])
+ subshape = wkp['rwork'][0]
+ rwork = WorkSpaceTools.check_work_array(lwork, subshape)
+ op.setup(rwork=rwork)
+ vorti.initialize(time=simu.time, topo=op.discrete_op.vorticity.topology)
+ velo.initialize(time=simu.time, topo=op.discrete_op.vorticity.topology)
+ op.apply()
+ op.wait()
+ return simu, op
-def computeVel(res, x, y, z, t):
- res[0][...] = sin(x * t) * cos(y) * cos(z)
- res[1][...] = - cos(x) * sin(y) * cos(z)
- res[2][...] = 0.
- return res
+def compute_dt_ref(topo, simu, op):
+ """Compute some reference values for dt for
+ each criteria
+ """
+ ref = Field(domain=topo.domain, formula=grad_v_func_3d,
+ name='Analytical', nb_components=topo.domain.dimension ** 2)
+ rd = ref.discretize(topo)
+ ref.initialize(time=simu.time, topo=topo)
+ ref2 = Field(domain=topo.domain, formula=strain_f_3d,
+ name='Analytical',
+ nb_components=topo.domain.dimension *
+ (topo.domain.dimension + 1) / 2)
+ rd2 = ref2.discretize(topo)
+ ref2.initialize(time=simu.time, topo=topo)
+ dd = topo.domain.dimension
+ ll = [rd[i * (dd + 1)] for i in xrange(dd)]
+ ic = topo.mesh.compute_index
+ refmax_adv = refmax_str = refmax_strain = 0.
+ ll2 = [None, ] * dd
+ ll2[0] = [rd2[0], rd2[3], rd2[4]]
+ ll2[1] = [rd2[1], rd2[3], rd2[5]]
+ ll2[2] = [rd2[2], rd2[4], rd2[5]]
+ for d in xrange(dd):
+ refmax_adv = max(compute_max([ll[d]], ic), refmax_adv)
+ refmax_str = max(compute_max(rd[d * dd: (d + 1) * dd], ic), refmax_str)
+ refmax_strain = max(compute_max(ll2[d], ic), refmax_strain)
+ dt = {}
+ dt['gradU'] = op.lcfl / refmax_adv
+ dt['stretch'] = op.method[TimeIntegrator].stability_coeff() / refmax_str
+ refmax = max([np.max(np.abs(v[ic]))
+ for v in op.velocity.discretize(topo).data])
+ dt['cfl'] = op.cfl * topo.mesh.space_step[0] / refmax
+ dt['strain'] = op.lcfl / refmax_strain
+ refmax = max([np.max(np.abs(w[ic]))
+ for w in op.vorticity.discretize(topo).data])
+ dt['vort'] = op.lcfl / refmax
+ return dt
-def computeVort(res, x, y, z, t):
- res[0][...] = 0.
- res[1][...] = 0.
- res[2][...] = sin(x) * cos(y * t) * cos(z)
- return res
+def test_adapt_vort():
+ """
+ Test 'vort' criteria (default config)
+ """
+ simu, op = init()
+ assert 'vort' in op.criteria
+ dt0 = simu.time_step
+ topo = op.discrete_op.vorticity.topology
+ wd = op.vorticity.discretize(topo)
+ ic = topo.mesh.compute_index
+ ref = max([np.max(np.abs(w[ic])) for w in wd.data])
+ assert np.allclose(op.diagnostics()[2], op.lcfl / ref)
+ assert np.allclose(simu.time_step, min(dt0, op.diagnostics()[2]))
-def init():
- box = pp.Box(length=[2.] * 3, origin=[0.0] * 3)
- velo = pp.Field(domain=box, formula=computeVel,
- name='Velocity', is_vector=True)
- vorti = pp.Field(domain=box, formula=computeVort,
- name='Vorticity', is_vector=True)
- return velo, vorti
+def test_adapt_fail():
+ """ Check wrong crit. value
+ """
+ fail = False
+ try:
+ init(['vorti', 'stretch'])
+ except:
+ fail = True
+ assert fail
-def test_adapt():
+def test_adapt_gradu():
"""
- Todo : write proper tests.
- Here we just check if discr/setup/apply process goes well.
+ Test 'gradu' criteria
"""
- velo, vorti = init()
- simu = Simulation(nb_iter=2)
- op = AdaptTimeStep(velo, vorti, simulation=simu,
- discretization=d3d, lcfl=0.125, cfl=0.5)
- op.discretize()
- op.setup()
- op.apply()
- op.wait()
+ simu, op = init(['gradU'])
+ dt0 = simu.time_step
+ topo = op.discrete_op.velocity.topology
+ ref = Field(domain=topo.domain, formula=diag_grad_func_3d,
+ name='Analytical', is_vector=True)
+ rd = ref.discretize(topo)
+ ref.initialize(time=simu.time, topo=topo)
+ dd = topo.domain.dimension
+ refmax = max([compute_max([rd[d]], topo.mesh.compute_index)
+ for d in xrange(dd)])
+ assert np.allclose(op.diagnostics()[3], op.lcfl / refmax,
+ rtol=topo.mesh.space_step.min())
+ assert np.allclose(simu.time_step, min(dt0, op.diagnostics()[3]))
-def test_adapt_2():
+def test_adapt_strain():
"""
- The same but with file output
+ Test 'strain' criteria
"""
- velo, vorti = init()
- simu = Simulation(nb_iter=2)
- op = AdaptTimeStep(velo, vorti, simulation=simu, io_params=True,
- discretization=d3d, lcfl=0.125, cfl=0.5)
- op.discretize()
- op.setup()
- op.apply()
- op.wait()
- filename = op.io_params.filename
- assert os.path.exists(filename)
+ simu, op = init(['strain'])
+ dt0 = simu.time_step
+ topo = op.discrete_op.velocity.topology
+ resdim = topo.domain.dimension * (topo.domain.dimension + 1) / 2
+ ref = Field(domain=topo.domain, formula=strain_f_3d,
+ name='Analytical', nb_components=resdim)
+ rd = ref.discretize(topo)
+ ref.initialize(time=simu.time, topo=topo)
+ dd = topo.domain.dimension
+ ll = [None, ] * dd
+ ll[0] = [rd[0], rd[3], rd[4]]
+ ll[1] = [rd[1], rd[3], rd[5]]
+ ll[2] = [rd[2], rd[4], rd[5]]
+ ic = topo.mesh.compute_index
+ refmax = 0.
+ for d in xrange(dd):
+ refmax = max(compute_max(ll[d], ic), refmax)
+ dt = op.lcfl / refmax
+ assert np.allclose(op.diagnostics()[4], dt,
+ rtol=topo.mesh.space_step.min())
+ assert np.allclose(simu.time_step, min(dt0, op.diagnostics()[4]))
-def test_adapt_3():
+def test_adapt_stretch():
"""
- The same but with external work vector
+ Test 'stretch' criteria
"""
- velo, vorti = init()
- simu = Simulation(nb_iter=2)
- op = AdaptTimeStep(velo, vorti, simulation=simu, io_params=True,
- discretization=d3d, lcfl=0.125, cfl=0.5)
- op.discretize()
- wk_p = op.get_work_properties()
- rwork = []
- wk_length = len(wk_p['rwork'])
- for i in xrange(wk_length):
- memshape = wk_p['rwork'][i]
- rwork.append(npw.zeros(memshape))
+ simu, op = init(['stretch'])
+ dt0 = simu.time_step
+ topo = op.discrete_op.velocity.topology
+ ref = Field(domain=topo.domain, formula=grad_v_func_3d,
+ name='Analytical', nb_components=topo.domain.dimension ** 2)
+ rd = ref.discretize(topo)
+ ref.initialize(time=simu.time, topo=topo)
+ dd = topo.domain.dimension
+ refmax = 0.
+ ic = topo.mesh.compute_index
+ for d in xrange(dd):
+ pos = d * dd
+ refmax = max(compute_max(rd[pos: pos + dd], ic), refmax)
+ cs = op.method[TimeIntegrator].stability_coeff()
+ dt = cs / refmax
+ assert np.allclose(op.diagnostics()[5], dt,
+ rtol=topo.mesh.space_step.min())
+ assert np.allclose(simu.time_step, min(dt0, op.diagnostics()[5]))
- op.setup(rwork=rwork)
- op.apply(simu)
- op.wait()
- filename = op.io_params.filename
- assert os.path.exists(filename)
+
+def test_adapt_cfl():
+ """
+ Test 'stretch' criteria
+ """
+ simu, op = init(['cfl'])
+ dt0 = simu.time_step
+ topo = op.discrete_op.velocity.topology
+ dt0 = simu.time_step
+ vd = op.velocity.discretize(topo)
+ ic = topo.mesh.compute_index
+ ref = max([np.max(np.abs(v[ic])) for v in vd.data])
+ dt = op.cfl * topo.mesh.space_step[0] / ref
+ assert np.allclose(op.diagnostics()[6], dt)
+ assert np.allclose(simu.time_step, min(dt0, op.diagnostics()[6]))
+
+
+def test_adapt_multi_crit():
+ """
+ Test combination of 3 criteria ('grad_u', 'stretch', 'cfl')
+ """
+ criteria = ['gradU', 'stretch', 'cfl']
+ simu, op = init(criteria)
+ topo = op.discrete_op.velocity.topology
+ dtref = compute_dt_ref(topo, simu, op)
+ dt0 = simu.time_step
+
+ assert np.allclose(op.diagnostics()[3], dtref['gradU'],
+ rtol=topo.mesh.space_step.min())
+ assert np.allclose(op.diagnostics()[5], dtref['stretch'],
+ rtol=topo.mesh.space_step.min())
+ assert np.allclose(op.diagnostics()[6], dtref['cfl'],
+ rtol=topo.mesh.space_step.min())
+ dtref_min = min([dtref[c] for c in criteria])
+ assert np.allclose(op.diagnostics()[0], simu.time)
+ assert np.allclose(op.diagnostics()[1], simu.time_step)
+ assert simu.time_step <= op.maxdt
+ assert np.allclose(simu.time_step, min(dt0, dtref_min))
+
+
+def test_adapt_all_crit():
+ """
+ Test combination of all criteria
+ """
+ criteria = ['gradU', 'stretch', 'cfl', 'strain', 'vort']
+ simu, op = init(criteria)
+ topo = op.discrete_op.velocity.topology
+ dtref = compute_dt_ref(topo, simu, op)
+ dt0 = simu.time_step
+
+ assert np.allclose(op.diagnostics()[0], simu.time)
+ assert np.allclose(op.diagnostics()[1], simu.time_step)
+ assert simu.time_step <= op.maxdt
+ assert np.allclose(op.diagnostics()[2], dtref['vort'],
+ rtol=topo.mesh.space_step.min())
+ assert np.allclose(op.diagnostics()[3], dtref['gradU'],
+ rtol=topo.mesh.space_step.min())
+ assert np.allclose(op.diagnostics()[4], dtref['strain'],
+ rtol=topo.mesh.space_step.min())
+ assert np.allclose(op.diagnostics()[5], dtref['stretch'],
+ rtol=topo.mesh.space_step.min())
+ assert np.allclose(op.diagnostics()[6], dtref['cfl'],
+ rtol=topo.mesh.space_step.min())
+ dtref_min = min([dtref[c] for c in criteria])
+ assert np.allclose(op.diagnostics()[0], simu.time)
+ assert np.allclose(op.diagnostics()[1], simu.time_step)
+ assert simu.time_step <= op.maxdt
+ assert np.allclose(simu.time_step, min(dt0, dtref_min))
+
+
+def test_adapt_all_crit_with_work():
+ """
+ Test combination of all criteria
+ """
+ criteria = ['gradU', 'stretch', 'cfl', 'strain', 'vort']
+ simu, op = init(criteria, work=True)
+ topo = op.discrete_op.velocity.topology
+ dtref = compute_dt_ref(topo, simu, op)
+ dt0 = simu.time_step
+
+ assert np.allclose(op.diagnostics()[0], simu.time)
+ assert np.allclose(op.diagnostics()[1], simu.time_step)
+ assert simu.time_step <= op.maxdt
+ assert np.allclose(op.diagnostics()[2], dtref['vort'],
+ rtol=topo.mesh.space_step.min())
+ assert np.allclose(op.diagnostics()[3], dtref['gradU'],
+ rtol=topo.mesh.space_step.min())
+ assert np.allclose(op.diagnostics()[4], dtref['strain'],
+ rtol=topo.mesh.space_step.min())
+ assert np.allclose(op.diagnostics()[5], dtref['stretch'],
+ rtol=topo.mesh.space_step.min())
+ assert np.allclose(op.diagnostics()[6], dtref['cfl'],
+ rtol=topo.mesh.space_step.min())
+ dtref_min = min([dtref[c] for c in criteria])
+ assert np.allclose(op.diagnostics()[0], simu.time)
+ assert np.allclose(op.diagnostics()[1], simu.time_step)
+ assert simu.time_step <= op.maxdt
+ assert np.allclose(simu.time_step, min(dt0, dtref_min))
-def test_adapt_4():
+def test_adapt_multi_tasks():
"""
- The same but with external work vector
+ Multi mpi tasks version
"""
# MPI procs are distributed among two tasks
GPU = 4
@@ -105,18 +298,18 @@ def test_adapt_4():
proc_tasks[2] = GPU
proc_tasks[1] = VISU
- dom = pp.Box(dimension=3, proc_tasks=proc_tasks)
- velo = pp.Field(domain=dom, formula=computeVel,
- name='Velocity', is_vector=True)
- vorti = pp.Field(domain=dom, formula=computeVort,
- name='Vorticity', is_vector=True)
+ dom = Box(dimension=3, proc_tasks=proc_tasks)
+ velo = Field(domain=dom, formula=compute_vel,
+ name='Velocity', is_vector=True)
+ vorti = Field(domain=dom, formula=compute_vort,
+ name='Vorticity', is_vector=True)
from hysop.tools.parameters import MPIParams
cpu_task = MPIParams(comm=dom.comm_task, task_id=CPU)
simu = Simulation(nb_iter=4)
- op = AdaptTimeStep(velo, vorti, simulation=simu, io_params=True,
- discretization=d3d, lcfl=0.125, cfl=0.5,
- mpi_params=cpu_task)
+ op = AdaptiveTimeStep(velo, vorti, simulation=simu, io_params=True,
+ discretization=d3d, lcfl=0.125, cfl=0.5,
+ mpi_params=cpu_task)
simu.initialize()
if dom.is_on_task(CPU):
op.discretize()
diff --git a/hysop/operator/velocity_correction.py b/hysop/operator/velocity_correction.py
index 7cdd7e1c4..9c863f19a 100755
--- a/hysop/operator/velocity_correction.py
+++ b/hysop/operator/velocity_correction.py
@@ -76,3 +76,6 @@ class VelocityCorrection(Computational):
but do not apply it onto velocity.
"""
self.discrete_op.computeCorrection()
+
+ def get_work_properties(self):
+ return {'rwork': None, 'iwork': None}
diff --git a/hysop/testsenv.py b/hysop/testsenv.py
index b3b39590e..d677d10b6 100755
--- a/hysop/testsenv.py
+++ b/hysop/testsenv.py
@@ -1,6 +1,6 @@
"""Set some functions and variables useful to run tests.
"""
-from hysop import __FFTW_ENABLED__, __GPU_ENABLED__, __SCALES_ENABLED__
+from hysop import __FFTW_ENABLED__, __SCALES_ENABLED__
import pytest
import shutil
from hysop.tools.io_utils import IO
diff --git a/hysop/tools/misc.py b/hysop/tools/misc.py
index ce8c646df..66f6fcadb 100644
--- a/hysop/tools/misc.py
+++ b/hysop/tools/misc.py
@@ -235,7 +235,6 @@ class WorkSpaceTools(object):
lp = len(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
--
GitLab