diff --git a/Examples/dataTG.py b/Examples/dataTG.py
index 1e1a0031a3ca8ae912532c34b856792f630b595d..1b125b19060b233e31db6083507af82fe2db9bec 100644
--- a/Examples/dataTG.py
+++ b/Examples/dataTG.py
@@ -15,9 +15,9 @@ WITH_PROJ = True
# post treatment ...
OUTPUT_FREQ = 1
-FILENAME = './TG/energy2.dat'
+FILENAME = './TG/energyGUV.dat'
# simulation parameters
from parmepy.problem.simulation import Simulation
-simu = Simulation(tinit=0.0, tend=1., timeStep=0.002,
+simu = Simulation(tinit=0.0, tend=10., timeStep=0.002,
iterMax=1000000)
diff --git a/Examples/howto_integrators.py b/Examples/howto_integrators.py
index d15d5c092141479f581819ee21a4b1b6f2fa1557..ab4e27f9a42c3a14d22a979822e56e8913443848 100644
--- a/Examples/howto_integrators.py
+++ b/Examples/howto_integrators.py
@@ -7,7 +7,7 @@
# \todo : complete it with 2D and 3D cases.
#
-from parmepy.numerics.integrators import Euler, RK2, RK3, RK4
+from parmepy.methods import Euler, RK2, RK3, RK4
from parmepy.constants import WITH_GUESS
from parmepy.tools import numpywrappers as npw
import math
@@ -615,7 +615,7 @@ def test_2D_0(method, odemethod, optim):
errxodespy = npw.ones(nbdt)
erryodespy = npw.ones(nbdt)
dtodespy = npw.ones(nbdt)
- lwork, iwork = method.getWorkLengths(2)
+ lwork = method.getWorkLengths(2)
print lwork, nbdt
work = [npw.zeros((nb, nb)) for i in xrange(lwork)]
for i in xrange(nbdt):
@@ -652,4 +652,4 @@ print "\n====================== 1D ======================\n"
#test_1D_0()
print "\n====================== 2D ======================\n"
optim = WITH_GUESS
-test_2D_0(RK4, ode.RK2, optim)
+test_2D_0(RK4, ode.RK4, optim)
diff --git a/Examples/testPoisson.py b/Examples/testPoisson.py
index 728769663df9bc751e6f20386ec79b7133f83b5a..0fed7308ecf07e740969aec83ae748a275e9f3ed 100755
--- a/Examples/testPoisson.py
+++ b/Examples/testPoisson.py
@@ -17,7 +17,7 @@ cos = np.cos
dim = 3
LL = 2 * pi * np.ones((dim))
dom = pp.Box(dimension=dim, length=LL)
-resol = [17, 17, 17]
+resol = [65, 65, 65]
## Fields
velocity = pp.Field(domain=dom, isVector=True, name='Velocity')
@@ -99,6 +99,8 @@ assert np.allclose(ref.norm(), velocity.norm())
for i in range(dom.dimension):
assert np.allclose(vd[i], refD[i])
+print velocity.norm()
+
poisson.finalize()
#print poisson
diff --git a/HySoP/hysop/methods.py b/HySoP/hysop/methods.py
new file mode 100644
index 0000000000000000000000000000000000000000..6a9761308fd809ce9205fddedb4e14bc3f33c88d
--- /dev/null
+++ b/HySoP/hysop/methods.py
@@ -0,0 +1,44 @@
+"""
+@file methods.py
+Constant parameters names and values required for 'method' argument of
+the operators, numerical methods ...
+"""
+
+# TODO CHECK HOW THE LIST BELOW APPEARS IN DOXYGEN DOC
+
+# Dict Keys values:
+TimeIntegrator = 11111
+Remesh = 22222
+Interpolation = 33333
+Formulation = 44444
+SpaceDiscretisation = 55555
+
+# Dict authorized values:
+
+# Time integrators : method[TimeIntegrator]
+import parmepy.numerics.integrators.runge_kutta2 as runge_kutta2
+RK2 = runge_kutta2.RK2
+import parmepy.numerics.integrators.runge_kutta3 as runge_kutta3
+RK3 = runge_kutta3.RK3
+import parmepy.numerics.integrators.runge_kutta4 as runge_kutta4
+RK4 = runge_kutta4.RK4
+import parmepy.numerics.integrators.euler as euler
+Euler = euler.Euler
+
+# Remesh
+#import parmepy.numerics as numerics
+#L2_1 = numerics.remeshing.L2_1
+# A compléter ...
+
+# Interpolation
+#Linear = numerics.interpolation.Linear
+
+# Finite differences
+import parmepy.numerics.finite_differences as fd
+FD_C_4 = fd.FD_C_4
+FD_C_4_CAA = 1234
+
+# Stretching formulations
+import parmepy.operator.discrete.stretching as strd
+Conservative = strd.Conservative
+GradUW = strd.GradUW
diff --git a/HySoP/hysop/numerics/differential_operations.py b/HySoP/hysop/numerics/differential_operations.py
index 8849933c1789f2aee28eb9282c5bf1001a210a8e..3d30fa02ba9f144bb18423eb65ca843ecea12dfa 100755
--- a/HySoP/hysop/numerics/differential_operations.py
+++ b/HySoP/hysop/numerics/differential_operations.py
@@ -6,7 +6,7 @@ Library of functions used to perform classical vector calculus
"""
from parmepy.constants import debug, XDIR, YDIR, ZDIR
from abc import ABCMeta, abstractmethod
-from finite_differences import FD_C_4
+from parmepy.methods import FD_C_4_CAA, FD_C_4
import numpy as np
@@ -40,11 +40,11 @@ class Curl(DifferentialOperation):
"""
Computes \f$ \nabla \ times (f) \f$, f being a vector field.
"""
- def __init__(self, topo, work, method='FD_C_4'):
+ def __init__(self, topo, work, method=FD_C_4):
assert len(work) == self.getWorkLengths()
DifferentialOperation.__init__(self, topo, work)
- if method.find('FD_C_4') >= 0:
+ if method is FD_C_4:
# - 4th ordered FD,
# - work space provided by used during function call,
# i.e. memshape not required,
@@ -124,11 +124,11 @@ class DivV(DifferentialOperation):
- 1 seems to be faster than 2 for large systems (>200*3) and
less memory consumming. For smaller systems, well, it depends ...
"""
- def __init__(self, topo, work, method='FD_C_4'):
+ def __init__(self, topo, work, method=FD_C_4):
DifferentialOperation.__init__(self, topo, work)
- if method.find('FD_C_4_CAA') >= 0:
+ if method is FD_C_4_CAA:
# - 4th ordered FD,
# - fd.compute_and_add method
# declare and create fd scheme
@@ -140,7 +140,7 @@ class DivV(DifferentialOperation):
assert (topo.ghosts >= 2).all(),\
'you need a ghost layer for FD4 scheme.'
- elif method.find('FD_C_4') >= 0:
+ elif method is FD_C_4:
# - 4th ordered FD,
# - fd.compute method
# declare and create fd scheme
@@ -160,14 +160,14 @@ class DivV(DifferentialOperation):
self.fd_scheme.computeIndices(self._indices)
@staticmethod
- def getWorkLengths(nb_components=None, domain_dim=None, fd_method='C'):
+ def getWorkLengths(nb_components=None, domain_dim=None, fd_method=FD_C_4):
"""
- fd_method = 'C' --> call fd_compute
- fd_method = 'CAA' --> call fd.compute_and_add
+ fd_method = FD_C_4 --> call fd_compute
+ fd_method = FD_C_4_CAA --> call fd.compute_and_add
"""
if domain_dim == 1 or nb_components == 1:
return 1
- elif fd_method is 'CAA':
+ elif fd_method is FD_C_4_CAA:
return 1
else:
return 2
@@ -243,22 +243,14 @@ class DivT(DifferentialOperation):
\f$w\f$ and V some vector fields.
Methods :
- 1 - FD_C_4 (default) : 4th order, centered finite differences, with
- local workspace and based on fd.compute.
- init : func = DivT(topo, memshape=shape)
+ 1 - FD_C_4 (default) : 4th order, centered finite differences,
+ based on fd.compute.
+ init : func = DivT(topo, work)
call : result = func(var1, var2, result)
- 2 - FD_C_4_work : same as above but user must provide a work space
- during call (see DivT.FDCentral4_with_work)
- init : func = DivT(topo, method='FD_C_4_work')
- call : result = func(var1, var2, result, work)
- 3 - FD_C_4_CAA : 4th order, centered finite differences, with
- local workspace and based on fd.compute_and_add.
- init : func = DivT(topo, method='FD_C_4_CAA', memshape=shape)
+ 2 - FD_C_4_CAA : 4th order, centered finite differences,
+ based on fd.compute_and_add.
+ init : func = DivT(topo, work, method=FD_C_4_CAA)
call : result = func(var1, var2, result)
- 4 - FD_C_4_CAA_work : same as above but user must provide a work
- space during call (see DivT.FDCentral4_CAA_with_work)
- init : func = DivT(topo, method='FD_C_4_CAA_work')
- call : result = func(var1, var2, result, work)
var1 stands for W and var2 for V.
result must be a vector field (same number of components as var1 and var2),
@@ -267,9 +259,9 @@ class DivT(DifferentialOperation):
each component with the same size as var1[i].
"""
- def __init__(self, topo, work, method='FD_C_4'):
+ def __init__(self, topo, work, method=FD_C_4):
DifferentialOperation.__init__(self, topo, work)
- if method.find('FD_C_4') >= 0:
+ if method is FD_C_4 or method is FD_C_4_CAA:
# - 4th ordered FD,
# - fd.compute_and_add method
# connect call function
@@ -281,7 +273,7 @@ class DivT(DifferentialOperation):
raise ValueError("FD scheme Not yet implemented")
@staticmethod
- def getWorkLengths(nb_components, domain_dim, fd_method='C'):
+ def getWorkLengths(nb_components, domain_dim, fd_method):
return DivV.getWorkLengths(nb_components, domain_dim, fd_method)
def __call__(self, var1, var2, result):
@@ -305,8 +297,9 @@ class GradS(DifferentialOperation):
"""
Computes \f$ \nabla(\rho) \f$, \f$ \rho\f$ a scalar field
"""
- def __init__(self, topo, method='FD_C_4'):
- if method.find('FD_C_4') >= 0:
+ def __init__(self, topo, method=FD_C_4):
+
+ if method is FD_C_4:
self.fcall = self.FDCentral4
self.fd_scheme = FD_C_4(topo.mesh.space_step)
assert (topo.ghosts >= 2).all(),\
@@ -343,10 +336,10 @@ class GradV(DifferentialOperation):
Computes \f$ [\nabla(V)] \f$, with
\f$V\f$ a vector field.
"""
- def __init__(self, topo, method='FD_C_4'):
- if method.find('FD_C_4') >= 0:
+ def __init__(self, topo, method=FD_C_4):
+ if method is FD_C_4:
self.fcall = self.FDCentral4
- self.grad = GradS(topo, method='FD_C_4')
+ self.grad = GradS(topo, method)
else:
raise ValueError("FD scheme Not yet implemented")
@@ -373,10 +366,10 @@ class GradVxW(DifferentialOperation):
Computes \f$ [\nabla(V)][W] \f$, with
\f$V\f$ and \f$W\f$ some vector fields.
"""
- def __init__(self, topo, work, method='FD_C_4'):
+ def __init__(self, topo, work, method=FD_C_4):
DifferentialOperation.__init__(self, topo, work)
- if method.find('FD_C_4_diag') >= 0:
+ if method is FD_C_4:
# - 4th ordered FD, with diagnostics in output
# - fd.compute method
@@ -386,14 +379,6 @@ class GradVxW(DifferentialOperation):
assert (topo.ghosts >= 2).all(),\
'you need a ghost layer for FD4 scheme.'
self.fcall = self.FDCentral4_diag
-
- elif method.find('FD_C_4') >= 0:
- # - 4th ordered FD
- # - fd.compute method
- self.fd_scheme = FD_C_4(topo.mesh.space_step)
- self.fcall = self.FDCentral4
- assert (topo.ghosts >= 2).all(),\
- 'you need a ghost layer for FD4 scheme.'
else:
raise ValueError("FD scheme Not yet implemented")
@@ -441,21 +426,3 @@ class GradVxW(DifferentialOperation):
np.add(result[comp], self._work[0], result[comp])
diagnostics[1] = max(diagnostics[1], np.max(self._work[1]))
return result, diagnostics
-
- def FDCentral4(self, var1, var2, result, diagnostics=None):
- """
- """
- assert len(result) == self._dim
- nbc = len(var1)
- for comp in xrange(nbc):
- result[comp][...] = 0.0
- for cdir in xrange(self._dim):
- # self._work = d/dcdir (var1_comp)
- self.fd_scheme.compute(var1[comp], cdir, self._work[0])
- # compute self._work = self._work.var2[cdir]
- self._work[0][...] = self._work[0] * var2[cdir]
- # sum to obtain nabla(var_comp) . var2, saved into result[comp]
- np.add(result[comp], self._work[0], result[comp])
-
- return result
-
diff --git a/HySoP/hysop/numerics/integrators/__init__.py b/HySoP/hysop/numerics/integrators/__init__.py
index 183f62acb7d0f3a5b3147fea8531a6a3d2eb478c..d54697ff71812611d44f3637b7ecd3925db4d8da 100644
--- a/HySoP/hysop/numerics/integrators/__init__.py
+++ b/HySoP/hysop/numerics/integrators/__init__.py
@@ -21,13 +21,3 @@ Length of work list depends on integrator's type :
Work can be used to give extra parameters.
"""
-import runge_kutta2
-import runge_kutta3
-import runge_kutta4
-RK2 = runge_kutta2.RK2
-RK3 = runge_kutta3.RK3
-RK4 = runge_kutta4.RK4
-
-import euler
-Euler = euler.Euler
-
diff --git a/HySoP/hysop/operator/__init__.py b/HySoP/hysop/operator/__init__.py
index f4ebf2b9d73d370d5c0fa842e3b848d6edca3a65..09b8f2632344d26a2e52c155fcee2aab99c523e0 100644
--- a/HySoP/hysop/operator/__init__.py
+++ b/HySoP/hysop/operator/__init__.py
@@ -12,3 +12,15 @@
# ...
# advec.apply()
# \endcode
+#
+# \section aboutmeth About 'method' arg in operator init
+#
+# Method argument must be a dictionnary with some of the following keys/values:
+#
+# - method["TimeIntegrator"] = Euler, RK2, RK3, RK4
+# - method["Interpolation"] = Linear
+# - method["Remesh"] = L2_1, L2_2 ...
+# - method["SpaceDiscretization"] = FD_C_4
+# - method["Support"] = GPU, CPU
+# - method["Formulation"] = Conservative, GradUW (may depends on the operator)
+# - ...
diff --git a/HySoP/hysop/operator/continuous.py b/HySoP/hysop/operator/continuous.py
index a150fa698b5cfc7868375e49613b85b547700bd0..978e08fc4b31bceef359b350a7e50f27341d4739 100644
--- a/HySoP/hysop/operator/continuous.py
+++ b/HySoP/hysop/operator/continuous.py
@@ -31,13 +31,14 @@ class Operator(object):
@debug
@abstractmethod
- def __init__(self, variables, method=None, name=None,
+ def __init__(self, variables, method={},
topo=None, ghosts=None):
"""
Build the operator.
The only required parameter is a list of variables.
@param variables : list of fields on which this operator will apply.
- @param method : name of the method used to descretize
+ @param method : a dictionnary of methods.
+ See methods.py for authorized values.
the operator.
@param name : operator id.
@param topo : a predefined topology to discretize variables
@@ -60,8 +61,6 @@ class Operator(object):
must be defined on the same domain.'
## Object to describe the method of discretization of this operator.
self.discreteOperator = None
- ## Optional name to identify the operator.
- self.name = name
## bool to check if the setup function has been called for
## this operator
self._isUpToDate = False
@@ -90,11 +89,32 @@ class Operator(object):
else:
self.ghosts = topo.ghosts
+ @staticmethod
+ def getWorkLengths():
+ """
+ Return the length of working arrays lists required
+ for stretching discrete operator, depending on :
+ - the formulation (Conservative or GradUW)
+ - the time integrator (RK3, ...)
+ """
+ return 0, 0
+
+ def setWorks(self, rwork=[], iwork=[]):
+ self.discreteOperator.setWorks(rwork, iwork)
+
+ @abstractmethod
+ def discretize(self):
+ """
+ Build (or get) topologies required for this operator
+ and discretize all its variables.
+ """
@abstractmethod
def setUp(self):
"""
- Initialize the operator : compute or get topologies,
- discretize variables, create a discrete operator.
+ Last step of initializaton. After this, the operator must be
+ ready for apply call.
+
+ Main step : setup for discrete operators.
"""
@debug
diff --git a/HySoP/hysop/operator/discrete/analytic.py b/HySoP/hysop/operator/discrete/analytic.py
index e8dd941fdffa1c7ed8e20277f67005fae45aff59..66655aa6bfd2119406c36544bcc967fa4ced19cc 100644
--- a/HySoP/hysop/operator/discrete/analytic.py
+++ b/HySoP/hysop/operator/discrete/analytic.py
@@ -21,7 +21,7 @@ class Analytic_D(DiscreteOperator):
@param[in] formula : user defined formula to be applied.
@param method : Method used
"""
- DiscreteOperator.__init__(self, variables, method, name="Analytic")
+ DiscreteOperator.__init__(self, variables, method)
self.output = variables
self.vformula = np.vectorize(formula)
@@ -39,9 +39,3 @@ class Analytic_D(DiscreteOperator):
df.initialize(self.vformula, simulation.time)
for df in self.output:
df.contains_data = True
-
-
-if __name__ == "__main__":
- print __doc__
- print "- Provided class : Analytic_D"
- print Analytic_D.__doc__
diff --git a/HySoP/hysop/operator/discrete/discrete.py b/HySoP/hysop/operator/discrete/discrete.py
index 74f8e281ee477abd647baf8a17926f02281efd90..13aa3f9ac11fee84def208454c825a0e1d0d2dd2 100644
--- a/HySoP/hysop/operator/discrete/discrete.py
+++ b/HySoP/hysop/operator/discrete/discrete.py
@@ -23,7 +23,7 @@ class DiscreteOperator(object):
@debug
@abstractmethod
- def __init__(self, variables, method=None):
+ def __init__(self, variables, method={}):
"""
Create an empty discrete operator.
"""
@@ -38,8 +38,6 @@ class DiscreteOperator(object):
self.output = []
## Operator numerical method.
self.method = method
- ##
- self.numMethod = None
## Object to store computational times of lower level functions
self.timer = Timer(self)
## bool to check if the setup function has been called for
@@ -62,20 +60,15 @@ class DiscreteOperator(object):
def setWorks(self, rwork=[], iwork=[]):
- # Check inputs (list length)
- d_dim = self.variables[0].domain.dimension
- self._rwork_length, self._iwork_length = \
- self.getWorkLengths(domain_dim=d_dim)
- assert len(rwork) == self._rwork_length
- assert len(iwork) == self._iwork_length
# Set work arrays for real and int.
# Warning : no copy! We must have pointer
# links between work and self.work arrays.
self._rwork = rwork
self._iwork = iwork
self.hasExternalWork = True
+ # Warning : we do not check work lengths/size.
+ # This will be done during discreteOperator.setUp.
- @abstractmethod
def setUp(self):
"""
Update the operator --> after calling this function,
@@ -125,8 +118,3 @@ class DiscreteOperator(object):
for f in self.output:
s += str(f) + "\n"
return s
-
-if __name__ == "__main__":
- print __doc__
- print "- Provided class : DiscreteOperator"
- print DiscreteOperator.__doc__
diff --git a/HySoP/hysop/operator/discrete/stretching.py b/HySoP/hysop/operator/discrete/stretching.py
index 4db2f00a305dba353084ef0d2ed01eaa78ce7489..eeb26caec4b72ecf88cc2ff40c2f5302fd549eef 100755
--- a/HySoP/hysop/operator/discrete/stretching.py
+++ b/HySoP/hysop/operator/discrete/stretching.py
@@ -6,10 +6,11 @@ Discretisation of the stretching operator for two different formulations.
"""
from parmepy.constants import debug, WITH_GUESS
+from parmepy.methods import TimeIntegrator, SpaceDiscretisation,\
+ RK3, FD_C_4
from discrete import DiscreteOperator
from parmepy.numerics.integrators.euler import Euler
from parmepy.numerics.integrators.runge_kutta2 import RK2
-from parmepy.numerics.integrators.runge_kutta3 import RK3
from parmepy.numerics.integrators.runge_kutta4 import RK4
from parmepy.tools.timers import timed_function
from parmepy.numerics.differential_operations import DivT, GradVxW
@@ -24,14 +25,15 @@ class Stretching(DiscreteOperator):
"""
Abstract interface to stretching discrete operators.
Two formulations are available:
- - conservative, see ConservativeStretching
- - 'Grad(UxV), see GradStretching
+ - conservative, see operator.discrete.stretching.Conservative
+ - 'Grad(UxW), see operator.discrete.stretching.GradUW
"""
__metaclass__ = ABCMeta
@debug
- def __init__(self, velocity, vorticity, method='FD_C_4'):
+ def __init__(self, velocity, vorticity,
+ method={TimeIntegrator: RK3, SpaceDiscretisation: FD_C_4}):
"""
@param velocity : discrete field
@param vorticity : discrete field
@@ -57,18 +59,15 @@ class Stretching(DiscreteOperator):
## stability constant
self.cststretch = 0.
# Depends on time integration method
- if self.method.find('Euler') >= 0:
- self.num_method = Euler
+ timeint = self.method[TimeIntegrator]
+ if isinstance(timeint, Euler):
self.cststretch = 2.0
- elif self.method.find('RK2') >= 0:
- self.num_method = RK2
+ elif isinstance(timeint, RK2):
self.cststretch = 2.0
- elif self.method.find('RK4') >= 0:
- self.num_method = RK4
- self.cststretch = 2.7853
- else: # self.method.find('RK3') >= 0:
- self.num_method = RK3
+ elif isinstance(timeint, RK3):
self.cststretch = 2.5127
+ elif isinstance(timeint, RK4):
+ self.cststretch = 2.7853
# prepare ghost points synchro for velocity
self._synchronize = UpdateGhosts(self.velocity.topology,
@@ -96,19 +95,21 @@ class Stretching(DiscreteOperator):
self.str_work = self._rwork[self._work_length_ti:]
# A function to compute the gradient of a vector field
# Work vector is provided in input.
- self.strFunc = self.formulation(self.velocity.topology, self.str_work,
- method=self.str_method)
+ self.strFunc = \
+ self.formulation(self.velocity.topology,
+ self.str_work, method=
+ self.method[SpaceDiscretisation])
-class ConservativeStretching(Stretching):
+class Conservative(Stretching):
"""
Discretisation of the following problem :
\f{eqnarray*} \frac{\partial\omega}{\partial t} = \nabla.(\omega:v) \f}
"""
- def __init__(self, velocity, vorticity, method='FD_C_4'):
- Stretching.__init__(self, self.velocity, self.vorticity)
- self.formulation = DivT
- self.str_method = method
+ def __init__(self, velocity, vorticity,
+ method={TimeIntegrator: RK3, SpaceDiscretisation: FD_C_4}):
+ Stretching.__init__(self, velocity, vorticity, method)
+ self.formulation = Conservative
@staticmethod
def getWorkLengths(timeIntegrator):
@@ -126,11 +127,12 @@ class ConservativeStretching(Stretching):
return result
# Create the time integrator
- self.timeIntegrator = self.num_method(self.nbComponents,
- self.time_int_work,
- self.velocity.topology,
- f=rhs,
- optim=WITH_GUESS)
+ self.timeIntegrator = \
+ self.method[TimeIntegrator](self.nbComponents,
+ self.time_int_work,
+ self.velocity.topology,
+ f=rhs,
+ optim=WITH_GUESS)
self._isUpToDate = True
@timed_function
@@ -155,18 +157,16 @@ class ConservativeStretching(Stretching):
result=self.vorticity.data)
-class GradStretching(Stretching):
+class GradUW(Stretching):
"""
Discretisation of the following problem :
\f{eqnarray*} \frac{\partial\omega}{\partial t}=[\nabla(v)][\omega]\f}
"""
- def __init__(self, velocity, vorticity, method='FD_C_4'):
- Stretching.__init__(self, self.velocity, self.vorticity)
- self.formulation = GradVxW
- self.str_method = method + '_diag'
- ## a vector to save diagnostics computed from GradVxW (max div ...)
- self.diagnostics = npw.ones(2)
+ def __init__(self, velocity, vorticity,
+ method={TimeIntegrator: RK3, SpaceDiscretisation: FD_C_4}):
+ Stretching.__init__(self, velocity, vorticity, method)
+ self.formulation = GradUW
@staticmethod
def getWorkLengths(timeIntegrator):
@@ -177,19 +177,22 @@ class GradStretching(Stretching):
def setUp(self):
Stretching.setUp(self)
+ ## a vector to save diagnostics computed from GradVxW (max div ...)
+ self.diagnostics = npw.ones(2)
# Time integration scheme.
def rhs(t, y, result):
result, self.diagnostics =\
- self.graduv(self.velocity.data, y, result, self.diagnostics)
+ self.strFunc(self.velocity.data, y, result, self.diagnostics)
return result
# Create the time integrator
- self.timeIntegrator = self.num_method(self.nbComponents,
- self.time_int_work,
- self.velocity.topology,
- f=rhs,
- optim=WITH_GUESS)
+ self.timeIntegrator = \
+ self.method[TimeIntegrator](self.nbComponents,
+ self.time_int_work,
+ self.velocity.topology,
+ f=rhs,
+ optim=WITH_GUESS)
self._isUpToDate = True
@timed_function
diff --git a/HySoP/hysop/operator/stretching.py b/HySoP/hysop/operator/stretching.py
index 16624852b657cfc493aac11f1121e54eaef0e6e0..cbb7bbf7fce58546d4eb0eb20cdffb01c2e6c177 100755
--- a/HySoP/hysop/operator/stretching.py
+++ b/HySoP/hysop/operator/stretching.py
@@ -5,10 +5,10 @@
Computation of stretch. term in Navier-Stokes.
"""
-from parmepy.constants import debug, GRADUW, CONSERVATIVE
+from parmepy.constants import debug
+from parmepy.methods import TimeIntegrator, Formulation, SpaceDiscretisation,\
+ RK3, Conservative, FD_C_4
from parmepy.operator.continuous import Operator
-from parmepy.operator.discrete.stretching import ConservativeStretching
-from parmepy.operator.discrete.stretching import GradStretching
import numpy as np
@@ -19,7 +19,8 @@ class Stretching(Operator):
@debug
def __init__(self, velocity, vorticity, resolutions,
- formulation=None, method=None, topo=None, ghosts=None):
+ method={TimeIntegrator: RK3, Formulation: Conservative,
+ SpaceDiscretisation: FD_C_4}, topo=None, ghosts=None):
"""
Create a Stretching operator from given
velocity and vorticity variables.
@@ -44,32 +45,24 @@ class Stretching(Operator):
## Grid resolution for each variable (dictionnary)
self.resolutions = resolutions
## Numerical methods for time and space discretization
- if method is None:
- self.method = 'FD_C4_RK3'
- else:
- self.method = method
+ self.method = method
+ assert Formulation in self.method.keys()
+ assert SpaceDiscretisation in self.method.keys()
+ assert TimeIntegrator in self.method.keys()
+
## Formulation used for the stretching equation.
## Default = conservative form.
- if formulation is not None:
- self.formulation = formulation
- else:
- self.formulation = CONSERVATIVE
+ self.formulation = self.method[Formulation]
self.input = [self.velocity, self.vorticity]
self.output = [self.vorticity]
- @debug
- def setUp(self):
- print '================ Stretching scheme ====================='
- print ' Space discretization and integr. scheme:', self.method
- print '========================================================'
- nbGhosts = 0
- if self.method.find('FD_order2') >= 0:
- nbGhosts = 1
- elif self.method.find('FD_C4') >= 0:
+ def discretize(self):
+ if self.method[SpaceDiscretisation] is FD_C_4:
nbGhosts = 2
else:
- raise ValueError("Unknown method for stretching operator.")
+ raise ValueError("Unknown method for space discretization of the\
+ stretching operator.")
if self._predefinedTopo is not None:
assert (self._predefinedTopo.ghosts >= nbGhosts).all()
@@ -90,17 +83,27 @@ class Stretching(Operator):
ghosts=self.ghosts)
self.discreteFields[v] = v.discretize(topo)
- if self.formulation is CONSERVATIVE:
- self.discreteOperator =\
- ConservativeStretching(self.discreteFields[self.velocity],
- self.discreteFields[self.vorticity],
- method=self.method)
-
- elif self.formulation is GRADUW:
- self.discreteOperator =\
- GradStretching(self.discreteFields[self.velocity],
- self.discreteFields[self.vorticity],
- method=self.method)
+ @staticmethod
+ def getWorkLengths(formulation, timeIntegrator):
+ """
+ Return the length of working arrays lists required
+ for stretching discrete operator, depending on :
+ - the formulation (Conservative or GradUW)
+ - the time integrator (RK3, ...)
+ """
+ return formulation.getWorkLengths(timeIntegrator)
+
+ def setWorks(self, rwork=[], iwork=[]):
+ self.discreteOperator.setWorks(rwork, iwork)
+
+ @debug
+ def setUp(self):
+ submeth = self.method.copy()
+ del submeth[Formulation]
+ self.discreteOperator =\
+ self.formulation(self.discreteFields[self.velocity],
+ self.discreteFields[self.vorticity],
+ method=submeth)
self.discreteOperator.setUp()
self._isUpToDate = True