From 5227f5edf5d8d9628dbd045cc86f720ceaf0978a Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Franck=20P=C3=A9rignon?= <franck.perignon@imag.fr>
Date: Wed, 22 May 2013 15:59:35 +0000
Subject: [PATCH] Operators review (Jean-Matthieu approved ...)
---
HySoP/hysop/domain/obstacle/cylinder2d.py | 4 +-
HySoP/hysop/domain/obstacle/sphere.py | 8 +-
HySoP/hysop/fields/continuous.py | 11 +-
HySoP/hysop/operator/advection.py | 5 +-
HySoP/hysop/operator/analytic.py | 1 +
HySoP/hysop/operator/continuous.py | 13 +--
HySoP/hysop/operator/diffusion.py | 2 +
HySoP/hysop/operator/discrete/analytic.py | 18 ++-
.../operator/discrete/curlAndDiffusion_fft.py | 4 -
.../hysop/operator/discrete/diffusion_fft.py | 18 +--
HySoP/hysop/operator/discrete/discrete.py | 15 ++-
HySoP/hysop/operator/discrete/multiphase.py | 23 ++--
.../operator/discrete/particle_advection.py | 7 +-
HySoP/hysop/operator/discrete/penalization.py | 12 +-
HySoP/hysop/operator/discrete/poisson_fft.py | 9 +-
.../operator/discrete/scales_advection.py | 104 ++----------------
HySoP/hysop/operator/discrete/stretching.py | 32 +++---
HySoP/hysop/operator/monitors/monitoring.py | 6 +-
HySoP/hysop/operator/monitors/printer.py | 37 ++++---
HySoP/hysop/operator/multiphase.py | 2 +
HySoP/hysop/operator/penalization.py | 1 +
HySoP/hysop/operator/redistribute.py | 3 +-
HySoP/hysop/operator/stretching.py | 4 +-
HySoP/hysop/problem/problem.py | 88 ++++-----------
HySoP/hysop/problem/tests/test_transport.py | 6 +-
HySoP/hysop/problem/transport.py | 5 +-
26 files changed, 157 insertions(+), 281 deletions(-)
diff --git a/HySoP/hysop/domain/obstacle/cylinder2d.py b/HySoP/hysop/domain/obstacle/cylinder2d.py
index 5d7f8b978..dbded5d02 100644
--- a/HySoP/hysop/domain/obstacle/cylinder2d.py
+++ b/HySoP/hysop/domain/obstacle/cylinder2d.py
@@ -20,8 +20,8 @@ class Cylinder2D(Sphere):
@param radius : sphere radius, default = 1
@param porousLayers : a list of thicknesses
for successive porous layers
- radius is the outside sphere radius and thicknesses are given from
- outside layer to inside one.
+ radius is the inside sphere radius and thicknesses are given from
+ inside layer to outside one.
@param vd : velocity of the disk (considered as a rigid body),
default = 0.
"""
diff --git a/HySoP/hysop/domain/obstacle/sphere.py b/HySoP/hysop/domain/obstacle/sphere.py
index 99fa9cb31..ae2420223 100644
--- a/HySoP/hysop/domain/obstacle/sphere.py
+++ b/HySoP/hysop/domain/obstacle/sphere.py
@@ -21,8 +21,8 @@ class Sphere(Obstacle):
@param radius : sphere radius, default = 1
@param porousLayers : a list of thicknesses
for successive porous layers
- radius is the outside sphere radius and thicknesses are given from
- outside layer to inside one.
+ radius is the inside sphere radius and thicknesses are given from
+ inside layer to outside one.
@param vd : velocity of the sphere (considered as a rigid body),
default = 0.
"""
@@ -89,8 +89,8 @@ class HemiSphere(Sphere):
@param radius : sphere radius, default = 1
@param porousLayers : a list of thicknesses
for successive porous layers
- radius is the outside sphere radius and thicknesses are given from
- outside layer to inside one.
+ radius is the inside sphere radius and thicknesses are given from
+ inside layer to outside one.
@param vd : velocity of the sphere (considered as a rigid body),
default = 0.
"""
diff --git a/HySoP/hysop/fields/continuous.py b/HySoP/hysop/fields/continuous.py
index a5ad06085..521f79235 100644
--- a/HySoP/hysop/fields/continuous.py
+++ b/HySoP/hysop/fields/continuous.py
@@ -67,6 +67,8 @@ class Field(object):
self.isVector = isVector
## Analytic formula.
self._formula = formula
+ ## A list of extra parameters, used in _formula
+ self.extraParameters = tuple([])
@debug
def discretize(self, topo):
@@ -114,7 +116,7 @@ class Field(object):
@param *args : see parmepy.fields
"""
for dF in self.discreteFields.values():
- dF.initialize(self._formula, *args)
+ dF.initialize(self._formula, *(args + self.extraParameters))
def value(self, *pos):
"""
@@ -160,6 +162,13 @@ class Field(object):
else:
return self.discreteFields.values()[0].norm()
+ def setExtraParameters(self, *args):
+ """
+ Set values for (optional) list of extra parameters
+ that may be required in formula.
+ """
+ self.extraParameters = args
+
if __name__ == "__main__":
print __doc__
print "- Provided class : ContinuousField"
diff --git a/HySoP/hysop/operator/advection.py b/HySoP/hysop/operator/advection.py
index 03ff69959..c0dd4da43 100644
--- a/HySoP/hysop/operator/advection.py
+++ b/HySoP/hysop/operator/advection.py
@@ -1,7 +1,7 @@
"""
@file advection.py
-Advection of a scalar or vector field.
+Advection of a field.
"""
from parmepy.constants import debug, np, PARMES_INDEX
from parmepy.operator.continuous import Operator
@@ -31,8 +31,10 @@ class Advection(Operator):
v = [velocity]
if isinstance(advectedFields, list):
[v.append(f) for f in advectedFields]
+ self.output = advectedFields
else:
v.append(advectedFields)
+ self.output = [advectedFields]
Operator.__init__(self, v, method)
## Transport velocity
self.velocity = velocity
@@ -45,6 +47,7 @@ class Advection(Operator):
self.resolutions = resolutions
## Extra parameters (depend on the method)
self.config = other_config
+ self.input = [self.velocity]
@debug
def setUp(self):
diff --git a/HySoP/hysop/operator/analytic.py b/HySoP/hysop/operator/analytic.py
index e7c23ea0b..202695b57 100644
--- a/HySoP/hysop/operator/analytic.py
+++ b/HySoP/hysop/operator/analytic.py
@@ -43,6 +43,7 @@ class Analytic(Operator):
## Dict of resolutions (one per variable)
self.resolutions = resolutions
self.config = other_config
+ self.output = self.variables
@debug
def setUp(self):
diff --git a/HySoP/hysop/operator/continuous.py b/HySoP/hysop/operator/continuous.py
index dd2ca52c2..9ed546a13 100644
--- a/HySoP/hysop/operator/continuous.py
+++ b/HySoP/hysop/operator/continuous.py
@@ -85,15 +85,14 @@ class Operator(object):
self.discreteOperator.finalize()
@debug
- def apply(self, t=None, dt=None, ite=None):
+ def apply(self, simulation=None):
"""
- apply the operator on its variables.
- @param t : Current simulation time.
- @param dt : Time step
- @param ite : Iteration number
- @return duration of computation.
+ Apply this operator to its variables.
+ @param simulation : object that describes the simulation
+ parameters (time, time step, iteration number ...), see
+ parmepy.problem.simulation.Simulation for details.
"""
- return self.discreteOperator.apply(t, dt, ite)
+ return self.discreteOperator.apply(simulation)
def printComputeTime(self):
""" Time monitoring."""
diff --git a/HySoP/hysop/operator/diffusion.py b/HySoP/hysop/operator/diffusion.py
index e334b65ca..2c7bffaef 100644
--- a/HySoP/hysop/operator/diffusion.py
+++ b/HySoP/hysop/operator/diffusion.py
@@ -44,6 +44,8 @@ class Diffusion(Operator):
# The only available method at the time is fftw
if method is None:
self.method = 'fftw'
+ self.input = [self.vorticity]
+ self.output = [self.vorticity]
@debug
def setUp(self):
diff --git a/HySoP/hysop/operator/discrete/analytic.py b/HySoP/hysop/operator/discrete/analytic.py
index 575c9fd09..80fcc6d7a 100644
--- a/HySoP/hysop/operator/discrete/analytic.py
+++ b/HySoP/hysop/operator/discrete/analytic.py
@@ -1,5 +1,5 @@
"""
-@file operator/discrete/analytic.py
+@file analytic.py
Apply user-defined formula to a set of discrete variables.
"""
@@ -25,21 +25,19 @@ class Analytic_D(DiscreteOperator):
else:
DiscreteOperator.__init__(self, [variables], method)
- self.input = variables
self.output = variables
self.name = "analytic"
- self.formula = formula
+ self.vformula = np.vectorize(formula)
@debug
- def setUp(self):
- # numpy formula vectorization
- self.vformula = np.vectorize(self.formula)
-
- @debug
- def apply(self, t, dt, ite):
+ def apply(self, simulation):
+ """
+ Initialize is always called with coords + current time.
+ Any extra parameters must be set using field.setExtraParameters.
+ """
self.compute_time = MPI.Wtime()
for df in self.variables:
- df.initialize(self.vformula, t, dt, ite)
+ df.initialize(self.vformula, simulation.time)
for df in self.output:
df.contains_data = True
diff --git a/HySoP/hysop/operator/discrete/curlAndDiffusion_fft.py b/HySoP/hysop/operator/discrete/curlAndDiffusion_fft.py
index e546f2f49..6bede01e0 100644
--- a/HySoP/hysop/operator/discrete/curlAndDiffusion_fft.py
+++ b/HySoP/hysop/operator/discrete/curlAndDiffusion_fft.py
@@ -34,10 +34,6 @@ class DiffusionFFT(DiscreteOperator):
self.with_curl = with_curl
self.compute_time = 0.
- @debug
- def setUp(self):
- pass
-
@debug
def apply(self, t=None, dt=None, ite=None):
"""
diff --git a/HySoP/hysop/operator/discrete/diffusion_fft.py b/HySoP/hysop/operator/discrete/diffusion_fft.py
index d43d54174..0da5487b4 100644
--- a/HySoP/hysop/operator/discrete/diffusion_fft.py
+++ b/HySoP/hysop/operator/discrete/diffusion_fft.py
@@ -31,22 +31,16 @@ class DiffusionFFT(DiscreteOperator):
raise AttributeError("Wrong problem dimension: only 2D \
and 3D cases are implemented.")
# Base class initialisation
- DiscreteOperator.__init__(self, [vorticity], method)
+ DiscreteOperator.__init__(self, [self.vorticity], method)
+ self.input = [self.vorticity]
+ self.output = [self.vorticity]
@debug
- def setUp(self):
- pass
-
- @debug
- def apply(self, t, dt, ite):
- """
- Apply operator.
- @param dt : current time step
- """
+ def apply(self, simulation):
self.compute_time = MPI.Wtime()
- if dt is None:
+ if simulation is None:
raise ValueError("Missing dt value for diffusion computation.")
-
+ dt = simulation.timeStep
if (self.vorticity.dimension == 2):
self.vorticity.data = fftw2py.solve_diffusion_2d(
self.viscosity * dt, self.vorticity.data)
diff --git a/HySoP/hysop/operator/discrete/discrete.py b/HySoP/hysop/operator/discrete/discrete.py
index 0fb1d4d81..ff6ab109d 100644
--- a/HySoP/hysop/operator/discrete/discrete.py
+++ b/HySoP/hysop/operator/discrete/discrete.py
@@ -57,12 +57,21 @@ class DiscreteOperator(object):
s += str(self.compute_time)
print s
+ def setUp(self):
+ """
+ Update the operator --> after calling this function,
+ the operator is ready to use.
+ """
+ pass
+
@debug
@abstractmethod
- def apply(self, t, dt, ite):
+ def apply(self, simulation=None):
"""
- Abstract method, apply operaton on the of discrete variables of
- this operator.
+ Apply this operator to its variables.
+ @param simulation : object that describes the simulation
+ parameters (time, time step, iteration number ...), see
+ parmepy.problem.simulation.Simulation for details.
"""
@debug
diff --git a/HySoP/hysop/operator/discrete/multiphase.py b/HySoP/hysop/operator/discrete/multiphase.py
index 14fbcb283..172aaf677 100644
--- a/HySoP/hysop/operator/discrete/multiphase.py
+++ b/HySoP/hysop/operator/discrete/multiphase.py
@@ -14,7 +14,7 @@ class Pressure_d(DiscreteOperator):
"""
@debug
- def __init__(self, velocity, vorticity, density, method=''):
+ def __init__(self, velocity, vorticity, density, method=None):
"""
Constructor.
Create the old velocity field to compute the -Dp/rho
@@ -28,25 +28,16 @@ class Pressure_d(DiscreteOperator):
self.vorticity = vorticity
self.density = density
self.compute_time = 0.
+ self.input = [self.velocity_new, self.vorticity, self.density]
+ self.output = [self.vorticity]
@debug
- def setUp(self):
- pass
- ## La topologie se recupere depuis les variables discretes
- ## self.topology = self.velocity_old.topology
- ## self.topology = self.velocity_new.topology
- ## self.topology = self.vorticity.topology
- ## self.topology = self.vorticity.density
+ def apply(self, simulation=None):
+ if simulation is None:
+ raise ValueError("Missing dt value for diffusion computation.")
- ## self.resolution = self.topology.mesh.resolution
- ## self.meshSize = self.topology.mesh.size
-
- @debug
- def apply(self, t=None, dt=None, ite=None):
- """
- Apply operator.
- """
self.compute_time = time.time()
+ dt = simulation.timeStep
result = (self.velocity_new - self.velocity_old) / dt
gradU = DifferentialOperator(self.velocity_new, self.velocity_new,
choice='gradV',
diff --git a/HySoP/hysop/operator/discrete/particle_advection.py b/HySoP/hysop/operator/discrete/particle_advection.py
index ed760bfcd..aa8b05b86 100644
--- a/HySoP/hysop/operator/discrete/particle_advection.py
+++ b/HySoP/hysop/operator/discrete/particle_advection.py
@@ -147,8 +147,11 @@ class ParticleAdvection(DiscreteOperator):
p_adF.data[...], d)
@debug
- def apply(self, t=None, dt=None, ite=None):
- self.numMethod(t, dt)
+ def apply(self, simulation=None):
+ if simulation is None:
+ raise ValueError("Missing dt value for diffusion computation.")
+
+ self.numMethod(simulation.time, simulation.timeStep)
if __name__ == "__main__":
diff --git a/HySoP/hysop/operator/discrete/penalization.py b/HySoP/hysop/operator/discrete/penalization.py
index eeca52502..2b8b94702 100644
--- a/HySoP/hysop/operator/discrete/penalization.py
+++ b/HySoP/hysop/operator/discrete/penalization.py
@@ -50,22 +50,20 @@ class Penalization_d(DiscreteOperator):
for c in xrange(lsize + 1, len(cond)):
self.cond[c] = cond[c]
assert self.factor.size == len(self.cond)
+ self.input = self.output = self.variables
@debug
- def setUp(self):
- pass
-
- @debug
- def apply(self, t=None, dt=0.0, ite=None):
+ def apply(self, simulation):
+ if simulation is None:
+ raise ValueError("Missing dt value for penalization computation.")
self.compute_time = MPI.Wtime()
-
+ dt = simulation.timeStep
coef = 1.0 / (1.0 + dt * self.factor)
for v in self.variables:
i = 0
for cond in self.cond:
- #i = self.cond.index(cond)
v[cond] *= coef[i]
i += 1
diff --git a/HySoP/hysop/operator/discrete/poisson_fft.py b/HySoP/hysop/operator/discrete/poisson_fft.py
index e3b878d93..ac679e22d 100644
--- a/HySoP/hysop/operator/discrete/poisson_fft.py
+++ b/HySoP/hysop/operator/discrete/poisson_fft.py
@@ -46,16 +46,9 @@ class PoissonFFT(DiscreteOperator):
self.input = [self.vorticity]
self.output = [self.velocity]
- @debug
- def setUp(self):
- pass
-
@debug
@prof
- def apply(self, t, dt, ite):
- """
- Apply operator.
- """
+ def apply(self, simulation=None):
self.compute_time = MPI.Wtime()
if self.case is '2D':
diff --git a/HySoP/hysop/operator/discrete/scales_advection.py b/HySoP/hysop/operator/discrete/scales_advection.py
index 26ba0e2f7..3c4aaa9b6 100644
--- a/HySoP/hysop/operator/discrete/scales_advection.py
+++ b/HySoP/hysop/operator/discrete/scales_advection.py
@@ -40,20 +40,12 @@ class ScalesAdvection(DiscreteOperator):
self.compute_time = 0.
@debug
- def setUp(self):
- pass
- #self.stab_coeff = stab_coeff
-# self.topology = self.velocity.topology
-# self.ghosts = self.topology.ghosts
-# self.resolution = self.topology.mesh.resolution
-# self.dim = self.topology.dim
+ def apply(self, simulation=None):
+ if simulation is None:
+ raise ValueError("Missing dt value for diffusion computation.")
- @debug
- def apply(self, t=None, dt=None, ite=None):
- """
- Apply operator.
- """
self.compute_time = time.time()
+ dt = simulation.timeStep
for adF in self.advectedFields:
if isinstance(adF, VectorField):
for d in xrange(adF.dimension):
@@ -68,93 +60,15 @@ class ScalesAdvection(DiscreteOperator):
self.velocity.data[1],
self.velocity.data[2],
adF.data)
- # ind0a = self.ghosts[0]
- # ind0b = self.resolution[0] - self.ghosts[0]
- # ind1a = self.ghosts[1]
- # ind1b = self.resolution[1] - self.ghosts[1]
- # ind2a = self.ghosts[2]
- # ind2b = self.resolution[2] - self.ghosts[2]
-
- # velocityNoG = [np.zeros((self.resolution - 2 * self.ghosts),
- # dtype=PARMES_REAL, order=ORDER)
- # for d in xrange(self.dim)]
- # for i in xrange(self.dim):
- # velocityNoG[i][...] = self.velocity[i][ind0a:ind0b,
- # ind1a:ind1b,
- # ind2a:ind2b]
- # if self.advectedField.vector:
-
- # advecFieldNoG = [np.zeros((self.resolution - 2 * self.ghosts),
- # dtype=PARMES_REAL, order=ORDER)
- # for d in xrange(self.dim)]
- # for i in xrange(self.dim):
- # advecFieldNoG[i][...] = self.advectedField[i][ind0a:ind0b,
- # ind1a:ind1b,
- # ind2a:ind2b]
-
- # # Vector field advection
- # advecFieldNoG[0][...] = \
- # scales2py.solve_advection(dt, velocityNoG[0][...],
- # velocityNoG[1][...],
- # velocityNoG[2][...],
- # advecFieldNoG[0][...])
- # advecFieldNoG[1][...] = \
- # scales2py.solve_advection(dt, velocityNoG[0][...],
- # velocityNoG[1][...],
- # velocityNoG[2][...],
- # advecFieldNoG[1][...])
- # advecFieldNoG[2][...] = \
- # scales2py.solve_advection(dt, velocityNoG[0][...],
- # velocityNoG[1][...],
- # velocityNoG[2][...],
- # advecFieldNoG[2][...])
-
- # for i in xrange(self.dim):
- # self.advectedField[i][ind0a:ind0b,
- # ind1a:ind1b, ind2a:ind2b] = \
- # advecFieldNoG[i][...]
-
- # # Computation of advected vector field norm --> addField :
- # self.addField[ind0a:ind0b,
- # ind1a:ind1b, ind2a:ind2b] = \
- # np.sqrt(advecFieldNoG[0][...] ** 2 + \
- # advecFieldNoG[1][...] ** 2 + \
- # advecFieldNoG[2][...] ** 2)
-
- # addFieldNoG = np.zeros((self.resolution - 2 * self.ghosts),
- # dtype=PARMES_REAL, order=ORDER)
- # addFieldNoG[...] = self.addField[ind0a:ind0b,
- # ind1a:ind1b,
- # ind2a:ind2b]
-
- # # Added Scalar field advection
- # addFieldNoG[...] = \
- # scales2py.solve_advection(dt, velocityNoG[0][...],
- # velocityNoG[1][...],
- # velocityNoG[2][...],
- # addFieldNoG[...])
-
- # self.addField[ind0a:ind0b, ind1a:ind1b, ind2a:ind2b] = \
- # addFieldNoG[...]
-
- # else:
- # advecFieldNoG = np.zeros((self.resolution - 2 * self.ghosts),
- # dtype=PARMES_REAL, order=ORDER)
- # advecFieldNoG[...] = self.advectedField[ind0a:ind0b,
- # ind1a:ind1b,
- # ind2a:ind2b]
- # # Scalar field advection
- # advecFieldNoG[...] = \
- # scales2py.solve_advection(dt, velocityNoG[0][...],
- # velocityNoG[1][...],
- # velocityNoG[2][...],
- # advecFieldNoG[...])
- # self.advectedField[ind0a:ind0b, ind1a:ind1b, ind2a:ind2b] = \
- # advecFieldNoG[...]
self.compute_time = time.time() - self.compute_time
self.total_time += self.compute_time
return self.compute_time
+ def finalize(self):
+ """
+ \todo check memory deallocation in scales???
+ """
+ pass
if __name__ == "__main__":
print __doc__
diff --git a/HySoP/hysop/operator/discrete/stretching.py b/HySoP/hysop/operator/discrete/stretching.py
index b7169629f..b26a9419a 100755
--- a/HySoP/hysop/operator/discrete/stretching.py
+++ b/HySoP/hysop/operator/discrete/stretching.py
@@ -40,10 +40,9 @@ class Stretching_d(DiscreteOperator):
# ## input fields
# self.input = [self.velocity, self.vorticity]
self.propertyOp = propertyOp
- self.input = [velocity, vorticity]
+ self.input = [self.velocity, self.vorticity]
+ self.output = [self.vorticity]
- @debug
- def setUp(self):
self.cststretch = 0.
if self.method.find('Euler') >= 0:
self.num_method = Euler
@@ -63,7 +62,7 @@ class Stretching_d(DiscreteOperator):
self.cststretch = self.cststretch / 2.
@debug
- def apply(self, t=None, dt=None, ite=None):
+ def apply(self, simulation):
"""
Apply Stretching operator.
@@ -79,13 +78,16 @@ class Stretching_d(DiscreteOperator):
-define the function of ODEsolver as the result of the first step
-integrate with the chosen method (RK4. RK2. euler)
"""
- self.compute_time = time.time()
+ if simulation is None:
+ raise ValueError("Missing dt value for diffusion computation.")
+ self.compute_time = time.time()
+ dt = simulation.timeStep
if self.num_method is not None:
-
- if (self.propertyOp.find('divConservation') < 0
- and self.propertyOp.find('massConservation') < 0):
- print 'No or wrong stretching operator property specified, use default divConservation'
+ condA = self.propertyOp.find('divConservation') < 0
+ if condA and self.propertyOp.find('massConservation') < 0:
+ print 'No or wrong stretching operator property specified, \
+ use default divConservation'
self.propertyOp = 'divConservation'
if self.propertyOp == 'massConservation':
@@ -99,7 +101,7 @@ class Stretching_d(DiscreteOperator):
choice='divWU',
topology=self.velocity.topology)
- else: # self.propertyOp == 'divConservation'
+ else: # self.propertyOp == 'divConservation'
## Space discretization of grad(u)
self.stretchOp = DifferentialOperator(
@@ -113,7 +115,8 @@ class Stretching_d(DiscreteOperator):
print "dtAdvec", 0.5 / maxgersh[1]
## Determination of Stretching time step (stability condition)
- self.ndt = self.stabilityTest(dt, maxgersh[0], self.num_method, self.cststretch)
+ self.ndt = self.stabilityTest(dt, maxgersh[0],
+ self.num_method, self.cststretch)
self.dtstretch = dt / float(self.ndt)
## fct2Op
@@ -124,11 +127,12 @@ class Stretching_d(DiscreteOperator):
topology=self.vorticity.topology)
## Time integration
+ t = simulation.time
for i in range(self.ndt):
methodInt = self.num_method(self.function)
- self.vorticity.data = methodInt(self.function.apply,
- t, self.dtstretch,
- self.vorticity.data)
+ self.vorticity.data[...] = methodInt(self.function.apply,
+ t, self.dtstretch,
+ self.vorticity.data)
# print 'norm omega', LA.norm(self.vorticity.data)
self.compute_time = time.time() - self.compute_time
diff --git a/HySoP/hysop/operator/monitors/monitoring.py b/HySoP/hysop/operator/monitors/monitoring.py
index 750913167..35d7a1870 100644
--- a/HySoP/hysop/operator/monitors/monitoring.py
+++ b/HySoP/hysop/operator/monitors/monitoring.py
@@ -1,19 +1,21 @@
"""
@file monitoring.py
+Global interface for monitors.
"""
from abc import ABCMeta, abstractmethod
from parmepy.operator.continuous import Operator
class Monitoring(Operator):
- """Abstract description a monitor. """
+ """Abstract interface to monitoring operators."""
__metaclass__ = ABCMeta
@abstractmethod
def __init__(self, variables, frequency):
""" Constructor
- @param dimension integer : domain dimension.
+ @param variables : list of fields to monitor
+ @param frequency : output rate
"""
Operator.__init__(self, variables)
## Monitor frequency
diff --git a/HySoP/hysop/operator/monitors/printer.py b/HySoP/hysop/operator/monitors/printer.py
index 524a908ce..89ccb20e8 100644
--- a/HySoP/hysop/operator/monitors/printer.py
+++ b/HySoP/hysop/operator/monitors/printer.py
@@ -7,7 +7,7 @@ from parmepy.constants import np, S_DIR, PARMES_REAL
from parmepy.operator.monitors.monitoring import Monitoring
import evtk.hl as evtk
from parmepy.mpi import main_rank
-import time
+from parmepy.mpi import MPI
class Printer(Monitoring):
@@ -17,18 +17,19 @@ class Printer(Monitoring):
Performs outputs in VTK images.
"""
- def __init__(self, frequency=0, fields=[], outputPrefix='./out_'):
+ def __init__(self, frequency=0, fields=[], prefix='./out_', ext='.vtk'):
"""
Create a results printer for given fields, filename
prefix (relative path) and an output frequency.
- @param fields : Fields to output.
- @param frequency : output file producing frequency.
- @param outputPrefix : file name prefix, contains relative path.
+ @param fields : list of variables to export.
+ @param frequency : output rate (output every freq iteration)
+ @param prefix : output file name prefix.
+ @param ext : output file extension, default=.vtk.
"""
Monitoring.__init__(self, fields, frequency)
- ## Filename prefix
- self.outputPrefix = outputPrefix + str(main_rank)
+ ## output file name prefix
+ self.prefix = prefix + str(main_rank)
## Method to collect data in case of distributed data
self.get_data_method = None
if self.freq != 0:
@@ -38,7 +39,8 @@ class Printer(Monitoring):
self.step = self._passStep
## Printing compute time
self.compute_time = 0.
- self.call_number = 0
+ ## Internal counter
+ self._call_number = 0
def apply(self, t=None, dt=None, ite=None):
self.step(ite)
@@ -78,21 +80,20 @@ class Printer(Monitoring):
def _printStep(self, ite):
"""
- Printer core method.
-
- Write fields data to file on necessary.\n
- VTK images are writed. If fails, classical ascii output is performed.
+ Try to write data into VTK files.
+ If fails, turn to classical ascii output.
"""
- t = time.time()
+ t = MPI.Wtime()
if (ite % self.freq) == 0:
- self.call_number += 1
- print "== IO"
+ self._call_number += 1
+ # Transfer from GPU to CPU if required
for f in self.variables:
for df in f.discreteFields.values():
try:
df.toHost()
except AttributeError:
pass
+ # Set output file name
filename = self.outputPrefix + "results_{0:05d}.dat".format(ite)
print "Print file : " + filename
## VTK output \todo: Need fix in 2D, getting an IOError.
@@ -149,10 +150,10 @@ class Printer(Monitoring):
def printComputeTime(self):
self.time_info[0] = "\"IO\" "
self.time_info[1] = str(self.compute_time) + " "
- print "IO total time : ", self.compute_time, self.call_number
+ print "IO total time : ", self.compute_time, self._call_number
if __name__ == "__main__":
print __doc__
- print "- Provided class : " + providedClass
- print eval(providedClass).__doc__
+ print "- Provided class : " + Printer
+ print Printer.__doc__
diff --git a/HySoP/hysop/operator/multiphase.py b/HySoP/hysop/operator/multiphase.py
index 3d896572d..82e38078d 100644
--- a/HySoP/hysop/operator/multiphase.py
+++ b/HySoP/hysop/operator/multiphase.py
@@ -32,6 +32,8 @@ class Pressure(Operator):
self.viscosity = viscosity
self.resolutions = resolutions
self.config = other_config
+ self.input = [self.velocity, self.vorticity, self.density]
+ self.output = [self.vorticity]
@debug
def setUp(self):
diff --git a/HySoP/hysop/operator/penalization.py b/HySoP/hysop/operator/penalization.py
index a25ab4fb6..a1da8d482 100644
--- a/HySoP/hysop/operator/penalization.py
+++ b/HySoP/hysop/operator/penalization.py
@@ -52,6 +52,7 @@ class Penalization(Operator):
for resol in resolutions.values():
assert resol == self.resolution
self.config = other_config
+ self.input = self.output = self.variables
@debug
def setUp(self):
diff --git a/HySoP/hysop/operator/redistribute.py b/HySoP/hysop/operator/redistribute.py
index 3cee0259a..33a415591 100644
--- a/HySoP/hysop/operator/redistribute.py
+++ b/HySoP/hysop/operator/redistribute.py
@@ -62,6 +62,7 @@ class Redistribute(Operator):
assert v in opFrom.variables and v in opTo.variables, \
'Redistribute error : one of the variable is not present\
in both source and target operator.'
+ self.input = self.output = self.variables
@debug
def setUp(self):
@@ -81,7 +82,7 @@ class Redistribute(Operator):
self._isUpToDate = True
- def apply(self):
+ def apply(self, simulation=None):
"""
Proceed with data redistribution from opFrom to opTo
"""
diff --git a/HySoP/hysop/operator/stretching.py b/HySoP/hysop/operator/stretching.py
index dd3cd807c..d7d060577 100755
--- a/HySoP/hysop/operator/stretching.py
+++ b/HySoP/hysop/operator/stretching.py
@@ -15,6 +15,7 @@ import numpy as np
class Stretching(Operator):
"""
Stretching operator representation
+ \todo write latex formulas
"""
@debug
@@ -36,6 +37,8 @@ class Stretching(Operator):
self.vorticity = vorticity
self.resolutions = resolutions
self.config = other_config
+ self.input = [self.velocity, self.vorticity]
+ self.output = [self.vorticity]
@debug
def setUp(self):
@@ -51,7 +54,6 @@ class Stretching(Operator):
space : (FD_order2, FD_order4, ..) default : FD_order4
and time : (Euler, RK, ..) default : RK3
"""
- Operator.setUp(self)
print '================ Stretching scheme ====================='
print ' Space discretization and integr. scheme:', self.method
print '========================================================'
diff --git a/HySoP/hysop/problem/problem.py b/HySoP/hysop/problem/problem.py
index ca3593b82..d5d961c63 100644
--- a/HySoP/hysop/problem/problem.py
+++ b/HySoP/hysop/problem/problem.py
@@ -1,9 +1,8 @@
"""
-@package parmepy.problem.problem
+@file problem.py
Complete problem description.
"""
-from abc import ABCMeta, abstractmethod
from parmepy.constants import debug
from parmepy import __VERBOSE__
from parmepy.operator.monitors.monitoring import Monitoring
@@ -25,19 +24,18 @@ class Problem(object):
OpenGLRendering object), The loop over time-steps is passed to Qt4
"""
- __metaclass__ = ABCMeta
-
@debug
def __new__(cls, *args, **kw):
return object.__new__(cls, *args, **kw)
@debug
- @abstractmethod
- def __init__(self, operators, monitors=[]):
+ def __init__(self, operators, simulation, monitors=[]):
"""
Create a transport problem instance.
@param operators : list of operators.
+ @param simulation : a parmepy.simulation.Simulation object
+ to describe simulation parameters.
@param monitors : list of monitors.
"""
## Problem operators
@@ -45,7 +43,7 @@ class Problem(object):
for m in monitors:
self.operators.append(m)
## Computes time step and manage iterations
- self.timer = None
+ self.simulation = simulation
## Timings informations
self.time_info = ["", ""]
self.domain = self.operators[0].variables[0].domain
@@ -54,13 +52,14 @@ class Problem(object):
if not self.domain is v.domain:
raise ValueError("Problem must have only one " +
"domain for variables.")
+ ## A list of variables that must be initialized before
+ ## any call to op.apply()
+ self.input = []
@debug
- def setUp(self, t_end, dt):
+ def setUp(self):
"""
Set solver.
- @param t_end : duration of the simulation.
- @param dt : simulation time step.
"""
for op in self.operators:
if not isinstance(op, Monitoring):
@@ -72,16 +71,13 @@ class Problem(object):
if isinstance(op, Redistribute):
op.setUp()
- self.timer = Timer(t_end, dt)
if __VERBOSE__:
print "==== Variables initialization ===="
- var_list = []
- for op in self.operators:
- for v in op.variables:
- var_list.append(v)
- # List of unique elements (initialize once each variable)
- for v in set(var_list):
+
+ # \todo find a way to set self.input properly.
+ for v in self.input:
v.initialize()
+
if __VERBOSE__:
print "===="
for op in self.operators:
@@ -105,19 +101,18 @@ class Problem(object):
print "\n\n Start solving ..."
for op in self.operators:
if isinstance(op, Monitoring):
- op.apply(self.timer.t, self.timer.dt, ite=0)
+ op.apply(self.simulation)
try:
## Pass the main loop to the OpenGL viewer
if not self.glRenderer is None:
self.glRenderer.startMainLoop()
except AttributeError:
- while not self.timer.end:
- if (self.domain.topologies.values()[0].rank == 0):
- print "==== Iteration : {0:3d} t={1:6.3f} ====".format(
- self.timer.ite, self.timer.t + self.timer.dt)
+ while not self.simulation.isOver:
+ self.simulation.printState()
for op in self.operators:
- op.apply(self.timer.t, self.timer.dt, self.timer.ite)
- self.timer.step()
+ op.apply(self.simulation)
+ self.simulation.advance()
+
print "\n\n End solving\n"
print "=== Timings ==="
if (self.domain.topologies.values()[0].rank == 0):
@@ -160,50 +155,7 @@ class Problem(object):
return s
-
-class Timer:
- """
- Manage time steps and simulation end.
-
- """
-
- def __init__(self, t_end, dt, t_init=0.):
- """
- Creates a Timer.
-
- @param t_end : Simulation final time.
- @param dt : Time step.
- @param t_init : Simulation starting time.
- """
- ## Simulation final time
- self.t_end = t_end
- ## Simulation time step
- self.dt = dt
- ## Simulation current time
- self.t = t_init
- ## Is simulation is terminated
- self.end = False
- ## Iteration counter
- self.ite = 1
- ## Max iteration
- self.ite_max = int(round(self.t_end / self.dt))
-
- def step(self):
- """
- Update current time.
- """
- if self.ite < self.ite_max:
- self.ite += 1
- self.t += self.dt
- else:
- self.end = True
-
- def __str__(self):
- s = "Timer (DiscreteOperator). "
- return s
-
if __name__ == "__main__":
print __doc__
- print "- Provided class : problem, Timer"
+ print "- Provided class : problem"
print Problem.__doc__
- print Timer.__doc__
diff --git a/HySoP/hysop/problem/tests/test_transport.py b/HySoP/hysop/problem/tests/test_transport.py
index b0df7fe42..cb4885251 100644
--- a/HySoP/hysop/problem/tests/test_transport.py
+++ b/HySoP/hysop/problem/tests/test_transport.py
@@ -7,6 +7,7 @@ from parmepy.domain.box import Box
from parmepy.fields.continuous import Field
from parmepy.operator.advection import Advection
from parmepy.problem.transport import TransportProblem
+from parmepy.problem.simulation import Simulation
def cosinus_product_2D(x, y):
@@ -54,8 +55,9 @@ def assert_nullVelocity(dim, boxLength, boxMin, nbElem, finalTime, timeStep,
resolutions={velo: nbElem,
scal: nbElem},
)
- pb = TransportProblem([advec])
- pb.setUp(finalTime, timeStep)
+ simu = Simulation(tend=finalTime, nbiter=1)
+ pb = TransportProblem([advec], simu)
+ pb.setUp()
initial_scalar = scal.discreteFields.values()[0].data.copy()
pb.solve()
return np.allclose(initial_scalar, scal.discreteFields.values()[0].data,
diff --git a/HySoP/hysop/problem/transport.py b/HySoP/hysop/problem/transport.py
index 1eda4cbd7..3a1f8cdff 100644
--- a/HySoP/hysop/problem/transport.py
+++ b/HySoP/hysop/problem/transport.py
@@ -10,8 +10,8 @@ class TransportProblem(Problem):
"""
Transport problem description.
"""
- def __init__(self, operators, monitors=[]):
- Problem.__init__(self, operators, monitors)
+ def __init__(self, operators, simulation, monitors=[]):
+ Problem.__init__(self, operators, simulation, monitors)
self.advection, self.velocity = None, None
for op in self.operators:
if isinstance(op, Advection):
@@ -20,4 +20,3 @@ class TransportProblem(Problem):
self.velocity = op
if self.advection is None:
raise ValueError("Transport problem with no Advection operator")
-
--
GitLab