diff --git a/Examples/NSDebug.py b/Examples/NSDebug.py
index aeb96c8ca3d965af4c5cae2eab9d99e97f10b80c..be6101baa26a8348346a9d964144d96b5c13241a 100755
--- a/Examples/NSDebug.py
+++ b/Examples/NSDebug.py
@@ -46,7 +46,7 @@ boxorigin = npw.realarray([0., 0., 0.])
box = pp.Box(dim, length=boxlength, origin=boxorigin)
## Global resolutionsimu
-nbElem = [33,]*3
+nbElem = [129]*3
# Upstream flow velocity
uinf = 1.0
@@ -72,17 +72,9 @@ velo = Field(domain=box, formula=computeVel,
vorti = Field(domain=box, formula=computeVort,
name='Vorticity', isVector=True)
-## Usual Cartesian topology definition
-# At the moment we use two (or three?) topologies :
-# - "topo" for Stretching and all operators based on finite differences.
-# --> ghost layer = 2
-# - topo from Advection operator for all the other operators.
-# --> no ghost layer
-# - topo from fftw for Poisson and Diffusion.
-# Todo : check compat between scales and fft operators topologies.
+# Topo for operators that required ghosts points (stretching)
ghosts = np.ones((box.dimension)) * NBGHOSTS
-topo = Cartesian(box, box.dimension, nbElem,
- ghosts=ghosts)
+topostr = Cartesian(box, box.dimension, nbElem, ghosts=ghosts)
## === Obstacles (sphere + up and down plates) ===
sphere_pos = (boxlength - boxorigin) * 0.5
@@ -98,7 +90,7 @@ op['advection'] = Advection(velo, vorti,
op['stretching'] = Stretching(velo, vorti,
resolutions={velo: nbElem, vorti: nbElem},
- topo=topo)
+ topo=topostr)
op['diffusion'] = Diffusion(vorti, resolution=nbElem, viscosity=VISCOSITY)
@@ -134,18 +126,8 @@ distr = {}
distr['fft2curl'] = Redistribute([velo], op['penalization'], op['curl'])
distr['curl2fft'] = Redistribute([vorti], op['curl'], op['poisson'])
-
-# 1 - Curl to stretching (velocity only)
-#distr['curl2str'] = Redistribute([velo], op['curl'], op['stretching'])
-# 2 - Advection to stretching
-# velocity only
-#distr['adv2str_v'] = Redistribute([velo], op['advection'], op['stretching'])
-# vorticity only
-#distr['adv2str_w'] = Redistribute([vorti], op['advection'], op['stretching'])
-# Both
distr['adv2str'] = Redistribute([velo, vorti],
op['advection'], op['stretching'])
-# 3 - Stretching to Diffusion
distr['str2diff'] = Redistribute([vorti], op['stretching'], op['diffusion'])
distr['curl2adv'] = Redistribute([velo, vorti], op['curl'], op['advection'])
distr['curl2str'] = Redistribute([velo, vorti], op['curl'], op['stretching'])
@@ -160,16 +142,19 @@ simu = Simulation(tinit=0.0, tend=10., nbIter=2000, iterMax=1000000)
monitors = {}
# Save velo and vorti values on topofft
monitors['printerFFT'] = Printer(variables=[velo, vorti], topo=topofft,
- formattype=HDF5)
+ formattype=HDF5, prefix='fft_io')
monitors['printerCurl'] = Printer(variables=[velo, vorti], topo=topocurl,
- formattype=HDF5, prefix='curl_io_vorti')
-monitors['printerAdvec'] = Printer(variables=[velo, vorti], topo=topofft,
+ formattype=HDF5, prefix='curl_io')
+monitors['printerAdvec'] = Printer(variables=[velo, vorti], topo=topoadvec,
formattype=HDF5, prefix='advec_io',
xmfalways=True)
+monitors['printerStr'] = Printer(variables=[velo, vorti], topo=topostr,
+ formattype=HDF5, prefix='str_io',
+ xmfalways=True)
# Compute and save enstrophy/Energy, on topofft
io_default = {}
-monitors['energy'] = Energy_enstrophy(velo, vorti, topo=topofft,
+monitors['energy'] = Energy_enstrophy(velo, vorti, topo=topoadvec,
io_params={},
viscosity=VISCOSITY, isNormalized=False)
# Setup for reprojection
@@ -199,7 +184,7 @@ print coeffForce
## topo, cb,
## obstacles=[sphere], io_params={})
monitors['forces'] = DragAndLift(velo, vorti, VISCOSITY, coeffForce,
- topo, cb, io_params={})
+ topostr, cb, io_params={})
# === Setup for all declared operators ===
for ope in op.values():
@@ -220,16 +205,21 @@ allops = dict(op.items() + distr.items() + monitors.items())
# - compute vorticity with curl
def initFields_mode1():
# The velocity is initialized on the fftw topology
- # penalized and distributed on curl topo
+ # penalized and distributed to curl topo
velo.initialize(topo=topofft)
op['penalization'].apply(simu)
distr['fft2curl'].apply(simu)
# Vorticity is initialized after penalization of velocity
op['curl'].apply(simu)
- # Distribute vorticity on topofft
+ # Distribute vorticity and velocity on topostr
distr['curl2str'].apply(simu)
+ # Warning : distribute operators do not update ghost points. This is done
+ # by computational operators. So a synchronize call may be required
+ # to monitor proper values.
+ #op['stretching'].discreteOperator._synchronize(op['stretching'].discreteOperator.velocity.data + op['stretching'].discreteOperator.vorticity.data)
+ op['stretching'].updateGhosts()
# From this point both velocity and vorticity are initialized
- # on topocurl and topoadvec. velocity is also init. on topofft.
+ # on topostr and topocurl. velocity is also init. on topofft.
## Call initialization
initFields_mode1()
@@ -237,14 +227,21 @@ initFields_mode1()
##### Check initialize process results #####
norm_vel_fft = velo.norm(topofft)
norm_vel_curl = velo.norm(topocurl)
-norm_vort_fft = vorti.norm(topofft)
+norm_vel_str = velo.norm(topostr)
+
norm_vort_curl = vorti.norm(topocurl)
+norm_vort_str = vorti.norm(topostr)
-print norm_vort_curl, norm_vort_fft
+print norm_vort_curl, norm_vort_str
print norm_vel_curl, norm_vel_fft
-assert np.allclose(norm_vort_curl, norm_vort_fft)
+assert np.allclose(norm_vort_curl, norm_vort_str)
assert np.allclose(norm_vel_curl, norm_vel_fft)
+assert np.allclose(norm_vel_curl, norm_vel_str)
+wref = np.asarray([0., 4354.98713193, 612.8061093], dtype=np.float64)
+vref = np.asarray([1418.04337028, 0., 0.], dtype=np.float64)
+assert np.allclose(norm_vel_curl, vref)
+assert np.allclose(norm_vort_str, wref)
# Print initial state
for mon in monitors.values():
mon.apply(simu)
@@ -252,39 +249,48 @@ for mon in monitors.values():
## Note Franck : tests init ok, same values on both topologies, for v and w.
# === After init ===
#
-# v init on topocurl (which may be equal to topofft)
+# v,w init on topostr
#
# === Definition of the 'one-step' sequence of operators ===
#
-# 1 - w = curl(v), topocurl
-# 2 - w = advection(v,w), topo-advec
-# 3 - w = stretching(v,w), topostr
-# 4 - w = diffusion(w), topofft
-# 5 - v = poisson(w), topofft
-# 6 - v = correction(v, w), topofft
-# 7 - v = penalization(v), topofft
+# A - w = stretching(v,w), topostr
+# B
+# 1 - w = diffusion(w), topofft
+# 2 - w = reprojection(w), topofft
+# 3 - v = poisson(w), topofft
+# 4 - v = correction(v, w), topofft
+# 5 - v = penalization(v), topofft
+# C - w = curl(v), topocurl
+# D - w = advection(v,w), topoadvec
#
-# Required bridges :
-# 1 --> 2 : (v,w) on topo-advec
-# 2 --> 3 : (v,w) on topostr
-# 3 --> 4 : w on topofft
-# 7 --> 1 : v on topocurl
-
-## fullseq = ['advection', # in: v,w out: w, on topoadvec
-## 'adv2str', 'stretching', # in: v,w out: w on topostr
-## # in: w, out: v, w on topofft
-## 'str2diff', 'diffusion', 'reprojection', 'poisson', 'correction',
-## 'penalization',
-## 'fft2curl', 'curl',
-## 'curl2adv', 'energy', 'forces', 'printerAdvec',
-## ]
-fullseq = ['stretching', # in: v,w out: w on topostr
- # in: w, out: v, w on topofft
+
+# Required bridges and in/out var:
+# A, IN : v, w on topostr
+# OUT : w on topostr
+# A - B bridge topostr-topofft for w
+# B, IN : w
+# OUT : w, v on topofft
+# B - C : bridge topofft - topocurl for v
+# C, IN : v on topocurl
+# OUT : w on topocurl
+# C - D : bridge topocurl - topoadvec for v and w
+# D, IN : v,w on topoadvec
+# OUT : w on topoadvec
+# D - A : bridge topoadvec - topostr for v and w
+
+
+# Monitors :
+# - energy(v,w), works for any topo
+# - reprojection(w), needs ghost points
+# - printer(v,w), works for any topo
+# - forces(v,w), needs ghost points
+
+fullseq = ['stretching',
'str2diff', 'diffusion', 'reprojection', 'poisson', 'correction',
'penalization',
'fft2curl', 'curl',
'curl2adv',
- 'advection', # in: v,w out: w, on topoadvec
+ 'advection',
'energy', 'printerAdvec',
'adv2str', 'forces']
@@ -296,12 +302,12 @@ def run(sequence):
allops[name].apply(simu)
seq = fullseq
-
+seq = []
simu.initialize()
while not simu.isOver:
- if topocurl.rank == 0:
- simu.printState()
+ #if topocurl.rank == 0:
+ # simu.printState()
run(seq)
simu.advance()
@@ -337,3 +343,21 @@ for ope in distr.values():
ope.finalize()
for monit in monitors.values():
monit.finalize()
+
+from abc import ABCMeta, abstractmethod
+
+
+class A(object):
+
+ def toto(self):
+ print "toto A "
+
+ def titi(self):
+ self.toto()
+
+
+class B(A):
+
+ def toto(self):
+ print "toto B"
+
diff --git a/Examples/dataNS_bb.py b/Examples/dataNS_bb.py
index 3652e6358dd1668bb82ab5d9150bb8f21e948669..da4f93cb60d734b77b9b2f9a449d6a92ead83dd7 100644
--- a/Examples/dataNS_bb.py
+++ b/Examples/dataNS_bb.py
@@ -36,7 +36,7 @@ ADVECTION_METHOD = {Scales: 'p_M6', Splitting: 'classic'}
# Splitting: 'o2_FullHalf'}
# Curl method
#CURL_METHOD = {SpaceDiscretisation: FD_C_4, GhostUpdate: True}
-CURL_METHOD = {SpaceDiscretisation: fftw2py}
+CURL_METHOD = {SpaceDiscretisation: 'fftw'}
# Viscosity of the fluid
VISCOSITY = 1. / 300.
# Vorticity projection (list of 2 elements):
diff --git a/HySoP/CMakeLists.txt b/HySoP/CMakeLists.txt
index 4682220ffe121e308f29eff3d9739b4342367845..7373a839ebbaa45d961a2a5eab246c141f43ed2f 100644
--- a/HySoP/CMakeLists.txt
+++ b/HySoP/CMakeLists.txt
@@ -31,7 +31,7 @@ option(USE_MPI "compile and link parmes with mpi when this mode is enable. Defau
option(WITH_PPM "link Parmes with PPM library (core component) - Deprecated. Default = off." OFF)
option(WITH_PPMNumerics "link Parmes with PPM-numerics - Deprecated. Default = off" OFF)
option(WITH_TESTS "Enable testing. Default = off" ON)
-option(BUILD_SHARED_LIBS "Enable dynamic library build, default = OFF." OFF)
+option(BUILD_SHARED_LIBS "Enable dynamic library build, default = OFF." ON)
option(WITH_LIB_FORTRAN "Generate libparmes from fortran files in src, wrapped into parmepy.f2py module. Default = ON." ON)
option(WITH_SCALES "compile/create parmesscales lib and link it with Parmes. Default = ON." ON)
option(WITH_FFTW "Link with fftw library (required for some Parmes solvers), default = ON" ON)
@@ -175,17 +175,21 @@ add_custom_target(python-build ALL COMMAND ${PYTHON_EXECUTABLE} ${CMAKE_CURRENT_
# To install python package AND parmes library and modules
if(WITH_LIB_FORTRAN)
add_custom_target(python-install COMMAND ${PYTHON_EXECUTABLE} ${CMAKE_CURRENT_BINARY_DIR}/setup.py install ${install-opt} config_fc --f90exec=${CMAKE_Fortran_COMPILER}
- COMMAND make install
+ #COMMAND make install
WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR} COMMENT "build/install parmepy package")
else()
add_custom_target(python-install COMMAND ${PYTHON_EXECUTABLE} ${CMAKE_CURRENT_BINARY_DIR}/setup.py install ${install-opt} config_fc --f90exec=${CMAKE_Fortran_COMPILER}
- COMMAND make install
+ #COMMAND make install
WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR} COMMENT "build/install parmepy package")
endif()
+install(CODE "execute_process(COMMAND ${CMAKE_BUILD_TOOL} python-install WORKING_DIRECTORY \"${CMAKE_CURRENT_BINARY_DIR}\")")
+
# Target to remove parmepy from install-path.
-add_custom_target(python-cleaninstall COMMAND rm -rf ${CMAKE_INSTALL_PREFIX}*
- COMMENT "remove parmepy package and its dependencies")
+# Remark : this does not work properly if a PREFIX is given by user.
+# Todo : fix this.
+add_custom_target(python-cleaninstall COMMAND rm -rf ${CMAKE_INSTALL_PREFIX}
+ COMMENT "remove parmepy package and its dependencies in ${CMAKE_INSTALL_PREFIX}.")
# Target to clean sources (remove .pyc files) and build dir.
file(GLOB_RECURSE PYCFILES "${CMAKE_SOURCE_DIR}/*.pyc")
@@ -209,9 +213,10 @@ else()
endif()
-add_dependencies(python-build ${PARMES_LIBRARY_NAME})
-add_dependencies(python-install ${PARMES_LIBRARY_NAME})
-
+if(WITH_LIB_FORTRAN)
+ add_dependencies(python-build ${PARMES_LIBRARY_NAME})
+ add_dependencies(python-install ${PARMES_LIBRARY_NAME})
+endif()
file(GLOB _DIR_FILES_EXT ${_DIR}/${_EXT})
diff --git a/HySoP/hysop/default_methods.py b/HySoP/hysop/default_methods.py
new file mode 100644
index 0000000000000000000000000000000000000000..8acaa3aedc4e8d41ff808b6a322065bdf431afad
--- /dev/null
+++ b/HySoP/hysop/default_methods.py
@@ -0,0 +1,33 @@
+"""
+@file default_methods.py
+Default parameter values for methods in operators.
+"""
+from parmepy.methods_keys import TimeIntegrator, Interpolation, GhostUpdate,\
+ Remesh, Support, Splitting, MultiScale, Formulation, SpaceDiscretisation, \
+ dtAdvecCrit
+from parmepy.numerics.integrators.runge_kutta2 import RK2
+from parmepy.numerics.integrators.runge_kutta3 import RK3
+from parmepy.numerics.interpolation import Linear
+from parmepy.numerics.remeshing import L2_1
+
+
+#class DefaultMethods(object):
+
+ADVECTION = {TimeIntegrator: RK2, Interpolation: Linear,
+ Remesh: L2_1, Support: '', Splitting: 'o2', MultiScale: L2_1}
+
+from parmepy.numerics.finite_differences import FD_C_4
+
+DIFFERENTIAL = {SpaceDiscretisation: FD_C_4, GhostUpdate: True}
+
+ADAPT_TIME_STEP = {TimeIntegrator: RK3, SpaceDiscretisation: FD_C_4,
+ dtAdvecCrit: 'vort'}
+
+BAROCLINIC = {SpaceDiscretisation: FD_C_4}
+
+DIFFUSION = {SpaceDiscretisation: 'fftw', GhostUpdate: True}
+
+POISSON = {SpaceDiscretisation: 'fftw', GhostUpdate: True}
+
+STRETCHING = {TimeIntegrator: RK3, Formulation: "Conservative",
+ SpaceDiscretisation: FD_C_4}
diff --git a/HySoP/hysop/f2py/__init__.py b/HySoP/hysop/fakef2py/__init__.py
similarity index 100%
rename from HySoP/hysop/f2py/__init__.py
rename to HySoP/hysop/fakef2py/__init__.py
diff --git a/HySoP/hysop/f2py/fftw2py/__init__.py b/HySoP/hysop/fakef2py/fftw2py/__init__.py
similarity index 100%
rename from HySoP/hysop/f2py/fftw2py/__init__.py
rename to HySoP/hysop/fakef2py/fftw2py/__init__.py
diff --git a/HySoP/hysop/f2py/scales2py/__init__.py b/HySoP/hysop/fakef2py/scales2py/__init__.py
similarity index 100%
rename from HySoP/hysop/f2py/scales2py/__init__.py
rename to HySoP/hysop/fakef2py/scales2py/__init__.py
diff --git a/HySoP/hysop/fields/continuous.py b/HySoP/hysop/fields/continuous.py
index 1e5007e8a6fd06a35c36794ef18c44a348ff55fd..cbafa8fcae7705895084eda56822287fb861a5ed 100644
--- a/HySoP/hysop/fields/continuous.py
+++ b/HySoP/hysop/fields/continuous.py
@@ -6,8 +6,6 @@ Continuous variable description.
from parmepy.constants import debug
from parmepy.fields.discrete import DiscreteField
from parmepy.mpi import main_rank
-import parmepy.tools.numpywrappers as npw
-import numpy as np
class Field(object):
@@ -74,7 +72,7 @@ class Field(object):
## Is this field a vector field?
self.isVector = isVector
## Analytic formula.
- self._formula = formula
+ self.formula = formula
## A list of extra parameters, used in _formula
self.extraParameters = tuple([])
## Number of components of the field
@@ -120,7 +118,7 @@ class Field(object):
(i.e. is of type 'user_func_2', see parmepy.fields for details)
Warning: if set, this formula will overwrite any previous setting.
"""
- self._formula = formula
+ self.formula = formula
self.doVectorize = doVectorize
@debug
@@ -138,7 +136,7 @@ class Field(object):
else:
self.topoInit = topo
df = self.discretization(topo)
- df.initialize(self._formula, self.doVectorize, currentTime,
+ df.initialize(self.formula, self.doVectorize, currentTime,
*self.extraParameters)
def value(self, *pos):
@@ -150,8 +148,8 @@ class Field(object):
--> depends on formula's definition ...
@return formula(pos): value of the formula at the given position.
"""
- if self._formula is not None:
- return self._formula(*pos)
+ if self.formula is not None:
+ return self.formula(*pos)
else:
return 0.0
@@ -160,9 +158,9 @@ class Field(object):
s = '[' + str(main_rank) + '] " ' + self.name + ' ", '
s += str(self.dimension) + 'D '
- if self._formula is not None:
+ if self.formula is not None:
s += 'an analytic formula is associated to this field:'
- s += str(self._formula.func_name)
+ s += str(self.formula.func_name)
if self.isVector:
s += 'vector field '
else:
diff --git a/HySoP/hysop/numerics/differential_operations.py b/HySoP/hysop/numerics/differential_operations.py
index 1575827da586daae58298878362d80a76d17871c..ed9b566fd11a6b19ab3a264a3817632e06421a4c 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 parmepy.numerics.finite_differences import FD_C_4, FD2_C_2
+from parmepy.numerics.finite_differences import FD_C_4, FD_C_2, FD2_C_2
import numpy as np
@@ -19,7 +19,7 @@ class DifferentialOperation(object):
@debug
@abstractmethod
- def __init__(self, topo, work):
+ def __init__(self, topo, work=None):
if work is None:
work = []
self._work = work
@@ -36,6 +36,9 @@ class DifferentialOperation(object):
fields on which this method operates.
@return length of list of work arrays of reals.
"""
+ assert nb_components is None
+ assert domain_dim is None
+ assert fd_method is None
return 0
@@ -62,7 +65,7 @@ class Curl(DifferentialOperation):
assert (topo.ghosts >= 2).all(),\
'you need a ghost layer for FD4 scheme.'
elif method is FD_C_2:
- # - 4th ordered FD,
+ # - 2nd ordered FD,
# - work space provided by used during function call,
# i.e. memshape not required,
# - fd.compute method
@@ -376,20 +379,17 @@ class GradV(DifferentialOperation):
\f$V\f$ a vector field.
"""
def __init__(self, topo, method=FD_C_4):
+ DifferentialOperation.__init__(self, topo)
if method is FD_C_4:
self.fcall = self.FDCentral4
self.grad = GradS(topo, method)
else:
raise ValueError("FD scheme Not yet implemented")
- self._dim = topo.domain.dimension
-
def __call__(self, var1, result):
return self.fcall(var1, result)
def FDCentral4(self, var1, result):
- """
- """
nbc = len(var1)
assert len(result) == nbc * self._dim
for cdir in xrange(self._dim):
diff --git a/HySoP/hysop/operator/__init__.py b/HySoP/hysop/operator/__init__.py
index 8cebe57a80a8ba046ec546043d544d3a651adec8..3fdf73ccb7942c7f8cf48d15c2192d80b3973bf9 100644
--- a/HySoP/hysop/operator/__init__.py
+++ b/HySoP/hysop/operator/__init__.py
@@ -17,11 +17,7 @@
#
# 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)
-# - ...
-
+# Methods have default values, taken from parmepy.default_methods
+# preset dictionnaries.
+#
+# Keys in methods dict are given in parmepy.method_keys.
diff --git a/HySoP/hysop/operator/adapt_timestep.py b/HySoP/hysop/operator/adapt_timestep.py
index 70662d439d2d4f9830c97a33ff96ac5dd463d4ac..e1239bcb0651d490f358b2d1d84c545da71ed528 100755
--- a/HySoP/hysop/operator/adapt_timestep.py
+++ b/HySoP/hysop/operator/adapt_timestep.py
@@ -8,12 +8,12 @@ Definition of the adaptative time step according to the flow fields.
from parmepy.constants import debug
from parmepy.methods_keys import TimeIntegrator, SpaceDiscretisation,\
dtAdvecCrit
-from parmepy.numerics.integrators.runge_kutta3 import RK3
from parmepy.numerics.finite_differences import FD_C_4
from parmepy.operator.discrete.adapt_timestep import AdaptTimeStep_D
from parmepy.operator.continuous import Operator
from parmepy.variable_parameter import VariableParameter
import numpy as np
+import parmepy.default_methods as default
class AdaptTimeStep(Operator):
@@ -54,9 +54,7 @@ class AdaptTimeStep(Operator):
"""
if method is None:
- method = {TimeIntegrator: RK3,
- SpaceDiscretisation: FD_C_4,
- dtAdvecCrit: 'vort'}
+ method = default.ADAPT_TIME_STEP
Operator.__init__(self, [velocity, vorticity], method, topo=topo,
ghosts=ghosts, task_id=task_id, comm=comm)
## velocity variable (vector)
diff --git a/HySoP/hysop/operator/advection.py b/HySoP/hysop/operator/advection.py
index 2ed83e4643cda77619f1e267416322bb05fdf021..21e4cdae2e1a1cd61a770de6cf48564f14bcf742 100644
--- a/HySoP/hysop/operator/advection.py
+++ b/HySoP/hysop/operator/advection.py
@@ -7,12 +7,11 @@ from parmepy.constants import debug, np, PARMES_INDEX, S_DIR
from parmepy.operator.continuous import Operator
from parmepy.methods_keys import Scales, TimeIntegrator, Interpolation,\
Remesh, Support, Splitting, MultiScale
-from parmepy.numerics.integrators.runge_kutta2 import RK2
-from parmepy.numerics.interpolation import Linear
from parmepy.numerics.remeshing import L2_1
from parmepy.fields.continuous import Field
from parmepy.operator.redistribute import Redistribute
from parmepy.operator.advectionDir import AdvectionDir
+import parmepy.default_methods as default
class Advection(Operator):
@@ -65,12 +64,7 @@ class Advection(Operator):
@param topo : a predefined topology to discretize variables
"""
if method is None:
- method = {TimeIntegrator: RK2,
- Interpolation: Linear,
- Remesh: L2_1,
- Support: '',
- Splitting: 'o2',
- MultiScale: L2_1}
+ method = default.ADVECTION
v = [velocity]
if isinstance(advectedFields, list):
self.advectedFields = advectedFields
@@ -145,6 +139,8 @@ class Advection(Operator):
task_id=task_id, name_suffix=vars_str[0:-1] + ')',
**self.config)
+ # function to switch between CPU or GPU setup.
+ self._my_setUp = None
@debug
def discretize(self):
"""
@@ -288,20 +284,34 @@ class Advection(Operator):
self._my_setUp = self.setUp_Python
@staticmethod
- def getWorkLengths(time_integrator, interpolation, remesh, domain_dim):
+ def getWorkLengths(method=None, domain_dim=None):
"""
Return the length of working arrays lists required
for advction discrete operator, depending on :
- the time integrator (RK2, ...)
- the interpolation (which depends on domain dimension)
- the remeshing (which depends on domain dimension)
+ @param method : the dict of parameters for the operator.
+ Default = parmepy.default_methods.ADVECTION
"""
- tw = time_integrator.getWorkLengths(1)
- iw, iiw = interpolation.getWorkLengths(domain_dim=domain_dim)
- rw, riw = remesh.getWorkLengths(domain_dim=domain_dim)
+ if method is None:
+ method = default.ADVECTION
+ assert Interpolation in method,\
+ 'An interpolation is required for the advection method.'
+ assert TimeIntegrator in method,\
+ 'A time integrator is required for the advection method.'
+ assert Remesh in method,\
+ 'A remesh is required for the advection method.'
+ tw = method[TimeIntegrator].getWorkLengths(1)
+ iw, iiw = method[Interpolation].getWorkLengths(domain_dim=domain_dim)
+ rw, riw = method[Remesh].getWorkLengths(domain_dim=domain_dim)
return max(tw + iw, rw), max(iiw, riw)
- def setWorks(self, rwork=[], iwork=[]):
+ def setWorks(self, rwork=None, iwork=None):
+ if rwork is None:
+ rwork = []
+ if iwork is None:
+ iwork = []
if not self._isSCALES:
for i in xrange(self._dim):
self.advecDir[i].setWorks(rwork, iwork)
diff --git a/HySoP/hysop/operator/advection_dir.py b/HySoP/hysop/operator/advection_dir.py
index 61bdd2ae9ba726c17010ea236b0e94a654dfc045..97b76f9c7c91d9f315582e8d84727fa6562cf02b 100644
--- a/HySoP/hysop/operator/advection_dir.py
+++ b/HySoP/hysop/operator/advection_dir.py
@@ -5,12 +5,13 @@ Advection of a field in a single direction.
"""
from parmepy.constants import debug, np, S_DIR, PARMES_INDEX
from parmepy.methods_keys import TimeIntegrator, Interpolation, Remesh, \
- Support, Splitting, MultiScale
-from parmepy.numerics.integrators.runge_kutta2 import RK2
+ Support, MultiScale
from parmepy.numerics.interpolation import Linear
from parmepy.numerics.remeshing import L2_1, L4_2, L4_4
from parmepy.operator.continuous import Operator
from parmepy.tools.timers import Timer
+# To get default method values for advection:
+import parmepy.default_methods as default
class AdvectionDir(Operator):
@@ -35,12 +36,7 @@ class AdvectionDir(Operator):
task_id=None, name_suffix='', comm=None,
**other_config):
if method is None:
- method = {TimeIntegrator: RK2,
- Interpolation: Linear,
- Remesh: L2_1,
- Support: '',
- Splitting: 'o2',
- MultiScale: L2_1}
+ method = default.ADVECTION
v = [velocity]
[v.append(f) for f in advectedFields]
Operator.__init__(self, v, method, topo=topo,
@@ -65,7 +61,7 @@ class AdvectionDir(Operator):
raise ValueError("Multiscale advection is not supported in "
"Python yet, user should use Scales or GPU.")
if not MultiScale in method.keys():
- method[MultiScale] == L2_1
+ method[MultiScale] = L2_1
if method[MultiScale] == Linear:
self._v_ghosts = np.array([1, ] * self.domain.dimension,
dtype=PARMES_INDEX)
@@ -130,7 +126,7 @@ class AdvectionDir(Operator):
shape[self.dir] > 1:
topodims[-1] = size
else:
- topodims[...] = shape
+ topodims[...] = shape
else:
# MPI topology depends on direction
if self.dir == self.domain.dimension - 1:
@@ -154,20 +150,34 @@ class AdvectionDir(Operator):
self.discreteFields[v] = v.discretize(topo)
@staticmethod
- def getWorkLengths(time_integrator, interpolation, remesh, domain_dim):
+ def getWorkLengths(method=None, domain_dim=None):
"""
Return the length of working arrays lists required
for advction discrete operator, depending on :
- the time integrator (RK2, ...)
- the interpolation (which depends on domain dimension)
- the remeshing (which depends on domain dimension)
+ @param method : the dict of parameters for the operator.
+ Default = parmepy.default_methods.ADVECTION
"""
- tw = time_integrator.getWorkLengths(1)
- iw, iiw = interpolation.getWorkLengths(domain_dim=domain_dim)
- rw, riw = remesh.getWorkLengths(domain_dim=domain_dim)
+ if method is None:
+ method = default.ADVECTION
+ assert Interpolation in method,\
+ 'An interpolation is required for the advection method.'
+ assert TimeIntegrator in method,\
+ 'A time integrator is required for the advection method.'
+ assert Remesh in method,\
+ 'A remesh is required for the advection method.'
+ tw = method[TimeIntegrator].getWorkLengths(1)
+ iw, iiw = method[Interpolation].getWorkLengths(domain_dim=domain_dim)
+ rw, riw = method[Remesh].getWorkLengths(domain_dim=domain_dim)
return max(tw + iw, rw), max(iiw, riw)
- def setWorks(self, rwork=[], iwork=[]):
+ def setWorks(self, rwork=None, iwork=None):
+ if rwork is None:
+ rwork = []
+ if iwork is None:
+ iwork = []
self.discreteOperator.setWorks(rwork, iwork)
def setUp(self):
diff --git a/HySoP/hysop/operator/analytic.py b/HySoP/hysop/operator/analytic.py
index 78a176c42ea9eae4db6cfe8b2c958b1df0aaf8fb..cb50737fb27b9b010b467a9183f7e54ba6c600e5 100644
--- a/HySoP/hysop/operator/analytic.py
+++ b/HySoP/hysop/operator/analytic.py
@@ -36,11 +36,11 @@ class Analytic(Operator):
v.setFormula(formula, doVectorize)
self.doVectorize = doVectorize
else:
- self.formula = self.variables[0]._formula
+ self.formula = self.variables[0].formula
self.doVectorize = self.variables[0].doVectorize
# Only one formula allowed per operator
for v in self.variables:
- assert v._formula is self.formula
+ assert v.formula is self.formula
## Dict of resolutions (one per variable)
self.resolutions = resolutions
diff --git a/HySoP/hysop/operator/baroclinic.py b/HySoP/hysop/operator/baroclinic.py
index 824050d6074a4b957a2b4138905fd6159802fb2e..535f644d345c95a74b2bb0f46d6f0358db870631 100644
--- a/HySoP/hysop/operator/baroclinic.py
+++ b/HySoP/hysop/operator/baroclinic.py
@@ -9,6 +9,7 @@ from parmepy.operator.discrete.baroclinic import Baroclinic_d
from parmepy.methods_keys import SpaceDiscretisation
from parmepy.numerics.finite_differences import FD_C_4
from parmepy.constants import debug, np
+import parmepy.default_methods as default
class Baroclinic(Operator):
@@ -19,7 +20,7 @@ class Baroclinic(Operator):
@debug
def __init__(self, velocity, vorticity, density, viscosity,
resolutions=None, topo=None, ghosts=None,
- method={SpaceDiscretisation: FD_C_4},
+ method=None,
task_id=None, comm=None, **other_config):
"""
Constructor.
@@ -32,10 +33,13 @@ class Baroclinic(Operator):
@param resolutions : grid resolution of velocity, vorticity, density
@param method : solving method
(default = finite differences, 4th order, in space)
- @param topo : a predefined topology to discretize velocity/vorticity/density
+ @param topo : a predefined topology to discretize
+ velocity/vorticity/density
@param ghosts : number of ghosts points. Default depends on the method.
Autom. computed if not set.
"""
+ if method is None:
+ method = default.BAROCLINIC
Operator.__init__(self, [velocity, vorticity, density], method,
topo=topo, ghosts=ghosts, task_id=task_id, comm=comm)
diff --git a/HySoP/hysop/operator/continuous.py b/HySoP/hysop/operator/continuous.py
index 7a8c0b5467c3196806704620259d0d0544d1e7f9..5739b941e9ac21c241bd0e97a8a34c80ed72cdbb 100644
--- a/HySoP/hysop/operator/continuous.py
+++ b/HySoP/hysop/operator/continuous.py
@@ -103,13 +103,13 @@ class Operator(object):
self.time_info = None
@staticmethod
- def getWorkLengths():
+ def getWorkLengths(method=None, domain_dim=None):
"""
Return the length of working arrays lists required
- for stretching discrete operator, depending on :
- - the formulation (Conservative or GradUW)
- - the time integrator (RK3, ...)
+ for the discrete operator.
"""
+ assert method is None
+ assert domain_dim is None
return 0, 0
def setWorks(self, rwork=None, iwork=None):
@@ -192,6 +192,12 @@ class Operator(object):
s = shortName + " operator. Not discretised."
return s + "\n"
+ def updateGhosts(self):
+ """
+ Update ghost points values, if any.
+ """
+ self.discreteOperator.updateGhosts()
+
class EmptyOperator(object):
"""Empty operator"""
diff --git a/HySoP/hysop/operator/curlAndDiffusion.py b/HySoP/hysop/operator/curlAndDiffusion.py
index 55ac3e977875ae5466a3de18e1e7dfe388f499f4..808d27b695d047251eef15917ad04e78e13829ea 100644
--- a/HySoP/hysop/operator/curlAndDiffusion.py
+++ b/HySoP/hysop/operator/curlAndDiffusion.py
@@ -6,7 +6,10 @@ Operator for diffusion problem.
"""
from parmepy.operator.continuous import Operator
-from parmepy.f2py import fftw2py
+try:
+ from parmepy.f2py import fftw2py
+except ImportError:
+ from parmepy.fakef2py import fftw2py
from parmepy.operator.discrete.diffusion_fft import DiffusionFFT
from parmepy.constants import debug
@@ -37,10 +40,13 @@ class Diffusion(Operator):
Operator.__init__(self, [velocity, vorticity], method,
name="Curl and Diffusion",
task_id=task_id, comm=comm)
+
self.velocity = velocity
self.vorticity = vorticity
self.resolutions = resolutions
self.config = other_config
+ raise ValueError("This operator is obsolete and must be reviewed.\
+ Do not use it.")
@debug
def setUp(self):
diff --git a/HySoP/hysop/operator/differential.py b/HySoP/hysop/operator/differential.py
index f9627250d4118be45edf69cd26b2e82f043d3031..2c55b1a487860bd1ae943b7c0fb8cf121b3b9899 100644
--- a/HySoP/hysop/operator/differential.py
+++ b/HySoP/hysop/operator/differential.py
@@ -5,11 +5,15 @@ Differential operators
"""
from parmepy.constants import debug, np
from parmepy.operator.continuous import Operator
-from parmepy.operator.discrete.differential import Curl_d, GradFD
-from parmepy.methods_keys import SpaceDiscretisation, GhostUpdate
+from parmepy.operator.discrete.differential import CurlFFT, CurlFD, GradFD
+from parmepy.methods_keys import SpaceDiscretisation
from parmepy.numerics.finite_differences import FD_C_4
-from parmepy.f2py import fftw2py
-from abc import ABCMeta
+try:
+ from parmepy.f2py import fftw2py
+except ImportError:
+ from parmepy.fakef2py import fftw2py
+import parmepy.default_methods as default
+from abc import ABCMeta, abstractmethod
class Differential(Operator):
@@ -24,7 +28,7 @@ class Differential(Operator):
method=None, topo=None, ghosts=None, task_id=None, comm=None):
if method is None:
- method = {SpaceDiscretisation: FD_C_4, GhostUpdate: True}
+ method = default.DIFFERENTIAL
Operator.__init__(self, [invar, outvar], method, topo=topo,
ghosts=ghosts, task_id=task_id, comm=comm)
## input variable
@@ -33,8 +37,8 @@ class Differential(Operator):
self.outvar = outvar
## Grid resolution for each variable (dictionnary)
self.resolutions = resolutions
- # The only available method at the time is fftw
- if self.method[SpaceDiscretisation] is fftw2py:
+
+ if self.method[SpaceDiscretisation] is 'fftw':
self.resolution = self.resolutions[self.outvar]
assert self.resolution == self.resolutions[self.invar],\
'for fftw method, all variables must have\
@@ -43,61 +47,91 @@ class Differential(Operator):
self.input = [invar]
def discretize(self):
- if self.method[SpaceDiscretisation] is FD_C_4:
- nbGhosts = 2
- elif self.method[SpaceDiscretisation] is fftw2py:
- nbGhosts = 0
- else:
- raise ValueError("Unknown method for space discretization of the\
- stretching operator.")
- # get (or create) the topology
+ domdim = self.domain.dimension
+ # Get main communicator
if self._predefinedTopo is not None:
- assert (self._predefinedTopo.ghosts >= nbGhosts).all()
- topo = self._predefinedTopo
- for v in self.variables:
- self.discreteFields[v] = v.discretize(topo)
-
+ comm = self._predefinedTopo.comm
+ commsize = self._predefinedTopo.size
+ elif self._comm is not None:
+ comm = self._comm
+ commsize = self._comm.Get_size()
else:
- if self.method[SpaceDiscretisation] is fftw2py:
- if self._comm is not None:
- main_size = self._comm.Get_size()
- comm = self._comm
- else:
- from parmepy.mpi.main_var import main_size
- from parmepy.mpi.main_var import main_comm as comm
-
- # Compute local resolution/distribution of data
- # according to fftw requirements.
- localres, localoffset = fftw2py.init_fftw_solver(
- self.resolution, self.domain.length, comm=comm.py2f())
+ from parmepy.mpi.main_var import main_comm as comm
+ from parmepy.mpi.main_var import main_size as commsize
+
+ if self.method[SpaceDiscretisation] is 'fftw':
+ # FFTW case : init fftw (fortran) solver
+ localres, localoffset = fftw2py.init_fftw_solver(
+ self.resolution, self.domain.length, comm=comm.py2f())
+ topodims = np.ones(domdim)
+ topodims[-1] = commsize
+
+ # Case 1 : topology provided by user at init:
+ if self._predefinedTopo is not None:
+ # Check if input topo is complient with fftw topo
+ topo = self._predefinedTopo
+ assert (topo.shape == topodims).all(), 'input topology is\
+ not compliant with fftw.'
+ for v in self.variables:
+ self.discreteFields[v] = v.discretize(topo)
- topodims = np.ones((self.domain.dimension))
- topodims[-1] = main_size
- #variables discretization
+ else:
for v in self.variables:
topo = self.domain.getOrCreateTopology(
self.domain.dimension,
- self.resolution, topodims, precomputed=True,
- offset=localoffset, localres=localres,
- ghosts=self.ghosts, comm=self._comm)
+ self.resolution, topodims,
+ precomputed=True,
+ offset=localoffset,
+ localres=localres,
+ ghosts=self.ghosts,
+ comm=self._comm)
+ self.discreteFields[v] = v.discretize(topo)
+
+ elif self.method[SpaceDiscretisation] is FD_C_4:
+ # Finite differences method
+ # Minimal number of ghost points
+ nbGhosts = 2
+
+ # Case 1 : topology provided by user at init:
+ if self._predefinedTopo is not None:
+ # Check if input topo is complient with fftw topo
+ topo = self._predefinedTopo
+ # check if topo ghost nb is complient with numerical method
+ assert (topo.ghosts >= nbGhosts).all()
+ for v in self.variables:
self.discreteFields[v] = v.discretize(topo)
+
else:
- # default topology constructor, with global resolution
- # and domain
- dim = self.domain.dimension
if self.ghosts is None:
- self.ghosts = np.asarray([nbGhosts] * dim)
+ self.ghosts = np.asarray([nbGhosts] * domdim)
+ else:
+ assert (self.ghosts >= nbGhosts).all()
+
for v in self.variables:
- topo = self.domain.getOrCreateTopology(
- dim, self.resolutions[v],
- ghosts=self.ghosts, comm=self._comm)
+ topo = self.domain.getOrCreateTopology(domdim,
+ self.resolutions[v],
+ ghosts=self.ghosts,
+ comm=self._comm)
self.discreteFields[v] = v.discretize(topo)
+ else:
+ raise ValueError("Unknown method for space discretization of the\
+ differential operator.")
+
assert self.discreteFields[self.invar].topology == \
self.discreteFields[self.outvar].topology, \
'Operator not yet implemented for multiple resolutions.'
+ @abstractmethod
+ def setUp(self):
+ """
+ Last step of initialization. After this, the operator must be
+ ready for apply call.
+
+ Main step : setup for discrete operators.
+ """
+
class Curl(Differential):
"""
@@ -107,19 +141,35 @@ class Curl(Differential):
@debug
def setUp(self):
# Create a discrete operator, according to the chosen method
- self.discreteOperator = Curl_d(self.discreteFields[self.invar],
- self.discreteFields[self.outvar],
- method=self.method)
+ assert self.domain.dimension == 3, "Not yet implemented for dim < 3"
+ # todo : implement a 'false' 2D curl.
+ if self.method[SpaceDiscretisation] is 'fftw':
+ self.discreteOperator = CurlFFT(self.discreteFields[self.invar],
+ self.discreteFields[self.outvar],
+ method=self.method)
+ elif self.method[SpaceDiscretisation] is FD_C_4:
+ self.discreteOperator = CurlFD(self.discreteFields[self.invar],
+ self.discreteFields[self.outvar],
+ method=self.method)
+ else:
+ raise ValueError("The required Space Discretisation is\
+ not available for Curl.")
self.discreteOperator.setUp()
self._isUpToDate = True
@staticmethod
- def getWorkLengths():
+ def getWorkLengths(method=None, domain_dim=None):
"""
Return the length of working arrays lists required
for curl discrete operator
+ @param method: useless, just there to fit with base class interface.
+ @param domain_dim : useless, just there
+ to fit with base class interface.
"""
- return Curl_d.getWorkLengths()
+ if method[SpaceDiscretisation] is FD_C_4:
+ return CurlFD.getWorkLengths()
+ else:
+ return 0
def setWorks(self, rwork=None, iwork=None):
self.discreteOperator.setWorks(rwork, iwork)
@@ -132,6 +182,10 @@ class Grad(Differential):
@debug
def setUp(self):
+ if self.method[SpaceDiscretisation] is not FD_C_4:
+ raise ValueError("Grad operator only\
+ implemented with finite differences. Please change\
+ method[SpaceDiscretisation] value.")
# Create a discrete operator, according to the chosen method
# (only finite differences at the time).
self.discreteOperator = GradFD(self.discreteFields[self.invar],
diff --git a/HySoP/hysop/operator/diffusion.py b/HySoP/hysop/operator/diffusion.py
index c1eb3ab8ec9a8f2119342ae1486fba1fe7b0eb23..55ff1a7f73804acaa91309583d6fe4c528f6e4eb 100644
--- a/HySoP/hysop/operator/diffusion.py
+++ b/HySoP/hysop/operator/diffusion.py
@@ -6,11 +6,15 @@ Operator for diffusion problem.
"""
from parmepy.operator.continuous import Operator
-from parmepy.f2py import fftw2py
+try:
+ from parmepy.f2py import fftw2py
+except ImportError:
+ from parmepy.fakef2py import fftw2py
from parmepy.operator.discrete.diffusion_fft import DiffusionFFT
from parmepy.constants import debug
-from parmepy.methods_keys import SpaceDiscretisation, GhostUpdate
+from parmepy.methods_keys import SpaceDiscretisation
import numpy as np
+import parmepy.default_methods as default
class Diffusion(Operator):
@@ -26,8 +30,7 @@ class Diffusion(Operator):
@debug
def __init__(self, vorticity, resolution, viscosity, method=None,
- topo=None, ghosts=None, task_id=None, comm=None,
- **other_config):
+ topo=None, ghosts=None, task_id=None, comm=None):
"""
Constructor for the diffusion operator.
@param[in,out] vorticity : field \f$ \omega \f$
@@ -36,7 +39,7 @@ class Diffusion(Operator):
"""
# The only available method at the time is fftw
if method is None:
- method = {SpaceDiscretisation: 'fftw', GhostUpdate: True}
+ method = default.DIFFUSION
## input/output field, solution of the problem
self.vorticity = vorticity
@@ -46,7 +49,7 @@ class Diffusion(Operator):
self.viscosity = viscosity
Operator.__init__(self, [vorticity], method=method, ghosts=ghosts,
topo=topo, task_id=task_id, comm=comm)
- self.config = other_config
+
self.input = [self.vorticity]
self.output = [self.vorticity]
@@ -104,7 +107,7 @@ class Diffusion(Operator):
"""
self.discreteOperator = DiffusionFFT(
self.discreteFields[self.vorticity], self.viscosity,
- method=self.method, **self.config)
+ method=self.method)
self.discreteOperator.setUp()
self._isUpToDate = True
diff --git a/HySoP/hysop/operator/discrete/curlAndDiffusion_fft.py b/HySoP/hysop/operator/discrete/curlAndDiffusion_fft.py
index 51350f1d57348143924411251a8ae14e7ac12162..591a524b7019db64211087b00fc26e046b6bb4cf 100644
--- a/HySoP/hysop/operator/discrete/curlAndDiffusion_fft.py
+++ b/HySoP/hysop/operator/discrete/curlAndDiffusion_fft.py
@@ -3,11 +3,13 @@
@file diffusion_fft.py
Discrete Diffusion operator using FFTW (fortran)
"""
-from parmepy.f2py import fftw2py
-from discrete import DiscreteOperator
+try:
+ from parmepy.f2py import fftw2py
+except ImportError:
+ from parmepy.fakef2py import fftw2py
+from parmepy.operator.discrete.discrete import DiscreteOperator
from parmepy.constants import debug
from parmepy.tools.timers import timed_function
-from parmepy.mpi import MPI, main_rank
class DiffusionFFT(DiscreteOperator):
diff --git a/HySoP/hysop/operator/discrete/differential.py b/HySoP/hysop/operator/discrete/differential.py
index 52a16ddee97f6b42d60f3f2d1441627f2d5f5582..0ca0def6704752bbbe69540c953704f0b4318921 100644
--- a/HySoP/hysop/operator/discrete/differential.py
+++ b/HySoP/hysop/operator/discrete/differential.py
@@ -4,16 +4,18 @@
Discretisation of the differential operators (curl, grad ...)
"""
-from parmepy.constants import debug, np
-from discrete import DiscreteOperator
-from parmepy.tools.timers import timed_function
+from parmepy.constants import debug
+from parmepy.operator.discrete.discrete import DiscreteOperator
from parmepy.numerics.differential_operations import Curl, GradV
-from parmepy.numerics.finite_differences import FD_C_4
import parmepy.tools.numpywrappers as npw
-from abc import ABCMeta
+from abc import ABCMeta, abstractmethod
from parmepy.numerics.updateGhosts import UpdateGhosts
-from parmepy.methods_keys import SpaceDiscretisation, GhostUpdate
-from parmepy.f2py import fftw2py
+from parmepy.methods_keys import SpaceDiscretisation
+try:
+ from parmepy.f2py import fftw2py
+except ImportError:
+ from parmepy.fakef2py import fftw2py
+import parmepy.default_methods as default
class Differential(DiscreteOperator):
@@ -36,97 +38,80 @@ class Differential(DiscreteOperator):
self.invar = invar
self.outvar = outvar
if method is None:
- method = {SpaceDiscretisation: FD_C_4, GhostUpdate: True}
+ method = default.DIFFERENTIAL
DiscreteOperator.__init__(self, [self.invar, self.outvar],
method=method)
self.input = [self.invar]
self.output = [self.outvar]
self._dim = self.invar.topology.domain.dimension
- # prepare ghost points synchro for velocity
- if self.method[GhostUpdate]:
- self._synchronize = UpdateGhosts(self.invar.topology,
- self.invar.nbComponents)
- else:
- def noghost(x):
- pass
+ self._synchronize = None
- self._synchronize = noghost
+ @staticmethod
+ def getWorkLengths(method=None, domain_dim=None):
+ assert method is None
+ assert domain_dim is None
+ return 0.
- @debug
- @timed_function
+ @abstractmethod
def apply(self, simulation=None):
- # Calling for requirements completion
- DiscreteOperator.apply(self, simulation)
-
- if isinstance(self, Curl_d) and \
- self.method[SpaceDiscretisation] is fftw2py:
-
- self.outvar.data[0], self.outvar.data[1],\
- self.outvar.data[2] = \
- fftw2py.solve_curl_3d(self.invar.data[0],
- self.invar.data[1],
- self.invar.data[2],
- self.outvar.data[0],
- self.outvar.data[1],
- self.outvar.data[2],
- self.ghosts_in,
- self.ghosts_out)
-
-# self.outvar.data[0], self.outvar.data[1], \
-# self.outvar.data[2] = \
-# fftw2py.projection_om_3d(self.outvar.data[0],
-# self.outvar.data[1],
-# self.outvar.data[2],
-# self.ghosts_out)
+ """
+ """
- else:
- self._synchronize(self.invar.data)
-# vd_init=[]
-# for i in xrange(3):
-# vd_init.append(np.empty_like(self.invar.data))
-# vd_init[i][...] = self.invar.data[i][...]
+class CurlFFT(Differential):
+ """
+ Compute the curl of a discrete field, using Fourier fftw
+ """
- self.outvar.data = self._function(self.invar.data, self.outvar.data)
+ def apply(self, simulation=None):
+ DiscreteOperator.apply(self, simulation)
+ ghosts_in = self.invar.topology.ghosts
+ ghosts_out = self.outvar.topology.ghosts
+ self.outvar.data[0], self.outvar.data[1], self.outvar.data[2] = \
+ fftw2py.solve_curl_3d(self.invar.data[0], self.invar.data[1],
+ self.invar.data[2], self.outvar.data[0],
+ self.outvar.data[1], self.outvar.data[2],
+ ghosts_in, ghosts_out)
-# for i in xrange(3):
-# print 'erreur', np.max(np.abs(vd_init[i][...] - self.invar.data[i][...]))
+ def finalize(self):
+ """
+ Clean memory (fftw plans and so on)
+ """
+ fftw2py.clean_fftw_solver(self.outvar.dimension)
-class Curl_d(Differential):
+class CurlFD(Differential):
"""
Compute the curl of a discrete field, using finite differences.
"""
def setUp(self):
- if self.method[SpaceDiscretisation] is fftw2py:
- self.ghosts_in = self.invar.topology.ghosts
- self.ghosts_out = self.outvar.topology.ghosts
+ # prepare ghost points synchro for velocity
+ self._synchronize = UpdateGhosts(self.invar.topology,
+ self.invar.nbComponents)
+
+ worklength = self.getWorkLengths()
+ memshape = self.invar.data[0].shape
+ if not self.hasExternalWork:
+ self._rwork = [npw.zeros(memshape) for i in xrange(worklength)]
else:
- self._worklength = self.getWorkLengths()
- memshape = self.invar.data[0].shape
- if not self.hasExternalWork:
- self._rwork = [npw.zeros(memshape)
- for i in xrange(self._worklength)]
- else:
- assert len(self._rwork) == self._work_length
- for wk in self._rwork:
- assert wk.shape == memshape
- assert self.outvar.nbComponents == self.invar.nbComponents
-
- self._function = Curl(self.invar.topology,
- self._rwork,
- self.method[SpaceDiscretisation])
+ assert len(self._rwork) == worklength
+ for wk in self._rwork:
+ assert wk.shape == memshape
+ assert self.outvar.nbComponents == self.invar.nbComponents
+
+ self._function = Curl(self.invar.topology, self._rwork,
+ self.method[SpaceDiscretisation])
@staticmethod
- def getWorkLengths():
+ def getWorkLengths(method=None, domain_dim=None):
return Curl.getWorkLengths()
- def finalize(self):
- """
- Clean memory (fftw plans and so on)
- """
- fftw2py.clean_fftw_solver(self.vorticity.dimension)
+ def apply(self, simulation=None):
+ # Calling for requirements completion
+ DiscreteOperator.apply(self, simulation)
+ self._synchronize(self.invar.data)
+ self.outvar.data = self._function(self.invar.data, self.outvar.data)
class GradFD(Differential):
@@ -135,6 +120,21 @@ class GradFD(Differential):
"""
def setUp(self):
+ # prepare ghost points synchro for velocity
+ self._synchronize = UpdateGhosts(self.invar.topology,
+ self.invar.nbComponents)
assert self.outvar.nbComponents == self._dim * self.invar.nbComponents
self._function = GradV(self.invar.topology,
self.method[SpaceDiscretisation])
+
+ def curlFD(self, simulation=None):
+ # Calling for requirements completion
+ DiscreteOperator.apply(self, simulation)
+ self._synchronize(self.invar.data)
+ self.outvar.data = self._function(self.invar.data, self.outvar.data)
+
+ def apply(self, simulation=None):
+ # Calling for requirements completion
+ DiscreteOperator.apply(self, simulation)
+ self._synchronize(self.invar.data)
+ self.outvar.data = self._function(self.invar.data, self.outvar.data)
diff --git a/HySoP/hysop/operator/discrete/diffusion_fft.py b/HySoP/hysop/operator/discrete/diffusion_fft.py
index 51709c0687e74c235749d6f84cb5cbf5eff16990..99c0cbb1bc88936960d7ca3eb27a4f53a886c604 100644
--- a/HySoP/hysop/operator/discrete/diffusion_fft.py
+++ b/HySoP/hysop/operator/discrete/diffusion_fft.py
@@ -3,7 +3,10 @@
@file diffusion_fft.py
Discrete Diffusion operator using FFTW (fortran)
"""
-from parmepy.f2py import fftw2py
+try:
+ from parmepy.f2py import fftw2py
+except ImportError:
+ from parmepy.fakef2py import fftw2py
from parmepy.operator.discrete.discrete import DiscreteOperator
from parmepy.constants import debug
from parmepy.mpi import MPI
@@ -36,7 +39,7 @@ class DiffusionFFT(DiscreteOperator):
self.output = [self.vorticity]
@debug
- def apply(self, simulation):
+ def apply(self, simulation=None):
if simulation is None:
raise ValueError("Missing dt value for diffusion computation.")
# Calling for requirements completion
@@ -68,4 +71,3 @@ class DiffusionFFT(DiscreteOperator):
Clean memory (fftw plans and so on)
"""
fftw2py.clean_fftw_solver(self.vorticity.dimension)
-
diff --git a/HySoP/hysop/operator/discrete/discrete.py b/HySoP/hysop/operator/discrete/discrete.py
index 96fe01654db125aba97bc0cee186d5b105d6c581..a3d8fa7d56f1c940f63a6fa5011f3b531d4bfe9a 100644
--- a/HySoP/hysop/operator/discrete/discrete.py
+++ b/HySoP/hysop/operator/discrete/discrete.py
@@ -67,6 +67,8 @@ class DiscreteOperator(object):
@return length of list of work arrays of reals.
@return length of list of work arrays of int.
"""
+ assert nb_components is None
+ assert domain_dim is None
return 0, 0
def setWorks(self, rwork=None, iwork=None):
@@ -134,3 +136,14 @@ class DiscreteOperator(object):
for f in self.output:
s += str(f) + "\n"
return s
+
+ def updateGhosts(self):
+ """
+ Update ghost points values, if any.
+ This function must be implemented in the discrete
+ operator if it may be useful to ask for ghost
+ points update without a call to apply.
+ For example for monitoring purpose before
+ operator apply.
+ """
+ pass
diff --git a/HySoP/hysop/operator/discrete/particle_advection.py b/HySoP/hysop/operator/discrete/particle_advection.py
index 1e2b44fe27132da8de5bc090faa9e9dd5078a639..050f16fb07d6040d543e9cd95cf024d62ecf5680 100644
--- a/HySoP/hysop/operator/discrete/particle_advection.py
+++ b/HySoP/hysop/operator/discrete/particle_advection.py
@@ -5,16 +5,12 @@ Advection solver, particular method, pure-python version.
"""
from parmepy.constants import debug, WITH_GUESS, PARMES_INDEX
-from parmepy.methods_keys import TimeIntegrator, Interpolation, Remesh, \
- Support, Splitting
-from parmepy.numerics.integrators.runge_kutta2 import RK2
-from parmepy.numerics.interpolation import Linear
-from parmepy.numerics.remeshing import L2_1
+from parmepy.methods_keys import TimeIntegrator, Interpolation, Remesh
from parmepy.operator.discrete.discrete import DiscreteOperator
from parmepy.fields.continuous import Field
-from parmepy.tools.timers import timed_function
import parmepy.tools.numpywrappers as npw
from parmepy.mpi import MPI
+import parmepy.default_methods as default
class ParticleAdvection(DiscreteOperator):
@@ -25,11 +21,7 @@ class ParticleAdvection(DiscreteOperator):
@debug
def __init__(self, velocity, advectedFields, d,
part_position=None, part_advectedFields=None,
- method={TimeIntegrator: RK2,
- Interpolation: Linear,
- Remesh: L2_1,
- Support: '',
- Splitting: 'o2'},
+ method=None,
isMultiScale=False):
"""
Constructor.
@@ -64,7 +56,8 @@ class ParticleAdvection(DiscreteOperator):
self.advectedFields = advectedFields
# update the list of discrete variables for this operator
[variables.append(advecF) for advecF in self.advectedFields]
-
+ if method is None:
+ method = default.ADVECTION
DiscreteOperator.__init__(self, variables, method)
self.input = self.variables
@@ -77,10 +70,27 @@ class ParticleAdvection(DiscreteOperator):
self._isMultiScale = isMultiScale
@staticmethod
- def getWorkLengths(time_integrator, interpolation, remesh, domain_dim):
- tw = time_integrator.getWorkLengths(1)
- iw, iiw = interpolation.getWorkLengths(domain_dim=domain_dim)
- rw, riw = remesh.getWorkLengths(domain_dim=domain_dim)
+ def getWorkLengths(method=None, domain_dim=None):
+ """
+ Return the length of working arrays lists required
+ for advction discrete operator, depending on :
+ - the time integrator (RK2, ...)
+ - the interpolation (which depends on domain dimension)
+ - the remeshing (which depends on domain dimension)
+ @param method : the dict of parameters for the operator.
+ Default = parmepy.default_methods.ADVECTION
+ """
+ if method is None:
+ method = default.ADVECTION
+ assert Interpolation in method,\
+ 'An interpolation is required for the advection method.'
+ assert TimeIntegrator in method,\
+ 'A time integrator is required for the advection method.'
+ assert Remesh in method,\
+ 'A remesh is required for the advection method.'
+ tw = method[TimeIntegrator].getWorkLengths(1)
+ iw, iiw = method[Interpolation].getWorkLengths(domain_dim=domain_dim)
+ rw, riw = method[Remesh].getWorkLengths(domain_dim=domain_dim)
return max(tw + iw, rw), max(iiw, riw)
@debug
@@ -104,7 +114,7 @@ class ParticleAdvection(DiscreteOperator):
Rmsh = self.method[Remesh]
memshape = self.part_position.data[0].shape
# Get work lengths
- w, iw = self.getWorkLengths(RKn, Interpol, Rmsh, self._dim)
+ w, iw = self.getWorkLengths(self.method, self._dim)
w_rk = RKn.getWorkLengths(nb_components=1)
w_interp, iw_interp = Interpol.getWorkLengths(domain_dim=self._dim)
w_remesh, iw_remesh = Rmsh.getWorkLengths(domain_dim=self._dim)
diff --git a/HySoP/hysop/operator/discrete/penalization.py b/HySoP/hysop/operator/discrete/penalization.py
index 5d3db9b7cb3ce63e37b859427ba344c22986de5c..dbc5e9f5d0e88e6ebf6180d4711ffe83b5d6160d 100644
--- a/HySoP/hysop/operator/discrete/penalization.py
+++ b/HySoP/hysop/operator/discrete/penalization.py
@@ -16,7 +16,7 @@ class Penalization_d(DiscreteOperator):
"""
@debug
- def __init__(self, variables, obstacles, factor, method=None):
+ def __init__(self, variables, obstacles, factor):
"""
Constructor.
@param[in,out] variables : list of discrete fields to be penalized
@@ -24,7 +24,7 @@ class Penalization_d(DiscreteOperator):
is applied.
@param[in] factor : penalization factoricient
"""
- DiscreteOperator.__init__(self, variables, method)
+ DiscreteOperator.__init__(self, variables)
## Penalization parameter
self.factor = np.asarray(factor)
@@ -51,7 +51,7 @@ class Penalization_d(DiscreteOperator):
@debug
@timed_function
- def apply(self, simulation):
+ def apply(self, simulation=None):
if simulation is None:
raise ValueError("Missing dt value for penalization computation.")
diff --git a/HySoP/hysop/operator/discrete/poisson_fft.py b/HySoP/hysop/operator/discrete/poisson_fft.py
index 9e24ff4f38e846060aa2f2ad2fe38388e50b26c8..28911c2cae934f93fd2398e700da73b89a134ccc 100644
--- a/HySoP/hysop/operator/discrete/poisson_fft.py
+++ b/HySoP/hysop/operator/discrete/poisson_fft.py
@@ -4,7 +4,11 @@
Discrete operator for Poisson problem (fftw based)
"""
import parmepy.tools.numpywrappers as npw
-from parmepy.f2py import fftw2py
+try:
+ from parmepy.f2py import fftw2py
+except ImportError:
+ from parmepy.fakef2py import fftw2py
+
from parmepy.operator.discrete.discrete import DiscreteOperator
from parmepy.constants import debug, prof
from parmepy.mpi import MPI
diff --git a/HySoP/hysop/operator/discrete/scales_advection.py b/HySoP/hysop/operator/discrete/scales_advection.py
index 08bca68a441d3b7cb90022fcd8358aa4cf444a68..1564a954fe6dfebd10c8820d0d52cd2e770548d9 100644
--- a/HySoP/hysop/operator/discrete/scales_advection.py
+++ b/HySoP/hysop/operator/discrete/scales_advection.py
@@ -4,13 +4,13 @@
Discrete Advection operator based on scales library (Jean-Baptiste)
"""
-from parmepy.f2py import scales2py
-from discrete import DiscreteOperator
-from parmepy.constants import debug, np
-from parmepy.numerics.differential_operations import GradVxW
+try:
+ from parmepy.f2py import scales2py
+except ImportError:
+ from parmepy.fakef2py import scales2py
+from parmepy.operator.discrete.discrete import DiscreteOperator
+from parmepy.constants import debug
from parmepy.mpi import MPI
-from parmepy.numerics.updateGhosts import UpdateGhosts
-import parmepy.tools.numpywrappers as npw
import math
ceil = math.ceil
diff --git a/HySoP/hysop/operator/discrete/stretching.py b/HySoP/hysop/operator/discrete/stretching.py
index be7845162f7d95ca5af99936b8f9553922463605..85c213a3b05f0c99e4645904b89be36e6f5f3e1b 100755
--- a/HySoP/hysop/operator/discrete/stretching.py
+++ b/HySoP/hysop/operator/discrete/stretching.py
@@ -7,7 +7,6 @@ Discretisation of the stretching operator for two different formulations.
from parmepy.constants import debug, WITH_GUESS
from parmepy.methods_keys import TimeIntegrator, SpaceDiscretisation
-from parmepy.numerics.finite_differences import FD_C_4
from parmepy.operator.discrete.discrete import DiscreteOperator
from parmepy.numerics.integrators.euler import Euler
from parmepy.numerics.integrators.runge_kutta2 import RK2
@@ -17,7 +16,8 @@ import parmepy.numerics.differential_operations as diff_op
import parmepy.tools.numpywrappers as npw
from parmepy.numerics.updateGhosts import UpdateGhosts
from parmepy.mpi import MPI
-from abc import ABCMeta
+from abc import ABCMeta, abstractmethod
+import parmepy.default_methods as default
import math
ceil = math.ceil
@@ -45,7 +45,7 @@ class Stretching(DiscreteOperator):
## vorticity discrete field
self.vorticity = vorticity
if method is None:
- method = {TimeIntegrator: RK3, SpaceDiscretisation: FD_C_4}
+ method = default.STRETCHING
DiscreteOperator.__init__(self, [self.velocity, self.vorticity],
method=method)
@@ -104,6 +104,17 @@ class Stretching(DiscreteOperator):
self.str_work, method=
self.method[SpaceDiscretisation])
+ def updateGhosts(self):
+ """
+ Update ghost points values
+ """
+ self._synchronize(self.velocity.data + self.vorticity.data)
+
+ @abstractmethod
+ def apply(self, simulation=None):
+ """
+ """
+
class Conservative(Stretching):
"""
@@ -115,9 +126,11 @@ class Conservative(Stretching):
self.formulation = diff_op.DivWV
@staticmethod
- def getWorkLengths(timeIntegrator):
- # Stretching only in 3D
- rwork_length = timeIntegrator.getWorkLengths(3)
+ def getWorkLengths(method=None, domain_dim=3):
+ assert TimeIntegrator in method,\
+ 'A time integrator is required for the stretching.'
+ # Stretching occurs only in 3D
+ rwork_length = method[TimeIntegrator].getWorkLengths(domain_dim)
rwork_length += diff_op.DivWV.getWorkLengths()
return rwork_length, 0
@@ -166,7 +179,6 @@ class Conservative(Stretching):
class GradUW(Stretching):
""" Discretisation of the following problem:
-
\f{eqnarray*} \frac{\partial\omega}{\partial t}=[\nabla(v)][\omega]\f}
"""
@@ -175,9 +187,20 @@ class GradUW(Stretching):
self.formulation = diff_op.GradVxW
@staticmethod
- def getWorkLengths(timeIntegrator):
- # Stretching only in 3D
- rwork_length = timeIntegrator.getWorkLengths(3)
+ def getWorkLengths(method=None, domain_dim=3):
+ """
+ Return the length of working arrays lists required
+ for stretching discrete operator, depending on :
+ - the formulation (Conservative or GradUW)
+ - the time integrator (RK3, ...)
+ @param method : the dict of parameters for the operator.
+ """
+ if method is None:
+ method = default.STRETCHING
+ assert TimeIntegrator in method,\
+ 'A time integrator is required for the stretching.'
+ # Stretching occurs only in 3D
+ rwork_length = method[TimeIntegrator].getWorkLengths(domain_dim)
rwork_length += diff_op.GradVxW.getWorkLengths()
return rwork_length, 0
diff --git a/HySoP/hysop/operator/monitors/compute_forces.py b/HySoP/hysop/operator/monitors/compute_forces.py
index 64e6db21e78d6ca6d3b9cc4e2519551fe7a6b28f..561362e48d14f939661027ab472d0eda3e877393 100644
--- a/HySoP/hysop/operator/monitors/compute_forces.py
+++ b/HySoP/hysop/operator/monitors/compute_forces.py
@@ -99,6 +99,7 @@ class DragAndLift(Monitoring):
self._step = np.asarray(self.topo.mesh.space_step)
self._dvol = np.prod(self._step)
self._work = npw.zeros(self.topo.mesh.resolution)
+ assert (self.topo.ghosts >= 1).all()
# function to compute the laplacian of a
# scalar field. Default fd scheme. (See Laplacian)
self._laplacian = Laplacian(self.topo)
diff --git a/HySoP/hysop/operator/penalization.py b/HySoP/hysop/operator/penalization.py
index 6023001768ec0bd866b90520be351f7187c63349..d74e68f9465914c99fe5c79aa98d165f7dbed236 100644
--- a/HySoP/hysop/operator/penalization.py
+++ b/HySoP/hysop/operator/penalization.py
@@ -23,8 +23,8 @@ class Penalization(Operator):
@debug
def __init__(self, variables, obstacles, coeff,
- resolutions=None, method=None, topo=None,
- ghosts=None, task_id=None, comm=None, **other_config):
+ resolutions=None, topo=None,
+ ghosts=None, task_id=None, comm=None):
"""
Constructor.
@param[in,out] variables : list of fields to be penalized
@@ -33,10 +33,7 @@ class Penalization(Operator):
@param[in] resolutions : list of resolutions (one per variable)
"""
- if method is None:
- method = {}
- Operator.__init__(self, variables, method=method,
- topo=topo, ghosts=ghosts,
+ Operator.__init__(self, variables, topo=topo, ghosts=ghosts,
task_id=task_id, comm=comm)
## domain where penalization must be applied.
@@ -52,7 +49,6 @@ class Penalization(Operator):
# all variables must have the same resolution
for resol in resolutions.values():
assert resol == self.resolution
- self.config = other_config
self.input = self.output = self.variables
def discretize(self):
@@ -76,10 +72,7 @@ class Penalization(Operator):
"""
"""
self.discreteOperator = Penalization_d(self.discreteFields.values(),
- self.obstacles, self.coeff,
- method=self.method,
- **self.config)
+ self.obstacles, self.coeff)
self.discreteOperator.setUp()
self._isUpToDate = True
-
diff --git a/HySoP/hysop/operator/poisson.py b/HySoP/hysop/operator/poisson.py
index a9c7a454f8eb8709a0ee6b28dbf8c8ca433118de..085d20fc76d196de86468e49b8123fd79057d9f7 100644
--- a/HySoP/hysop/operator/poisson.py
+++ b/HySoP/hysop/operator/poisson.py
@@ -5,12 +5,18 @@ Poisson problem.
"""
from parmepy.operator.continuous import Operator
-from parmepy.f2py import fftw2py
+try:
+ from parmepy.f2py import fftw2py
+except ImportError:
+ from parmepy.fakef2py import fftw2py
+
from parmepy.operator.discrete.poisson_fft import PoissonFFT
from parmepy.constants import debug
import numpy as np
from parmepy.mpi.main_var import main_size
from parmepy.operator.velocity_correction import VelocityCorrection
+from parmepy.methods_keys import SpaceDiscretisation
+import parmepy.default_methods as default
class Poisson(Operator):
@@ -52,14 +58,16 @@ class Poisson(Operator):
Operator.__init__(self, [velocity, vorticity], method, ghosts=ghosts,
topo=topo, task_id=task_id, comm=comm)
if method is None:
- self.method = 'fftw'
- if self.method is 'fftw':
+ self.method = default.POISSON
+
+ if self.method[SpaceDiscretisation] is 'fftw':
self.resolution = self.resolutions[self.velocity]
assert self.resolution == self.resolutions[self.vorticity],\
'Poisson error : for fftw, all variables must have\
the same global resolution.'
else:
raise AttributeError("Method not yet implemented.")
+
self.input = [self.vorticity]
self.output = [self.velocity]
if flowrate is not None:
diff --git a/HySoP/hysop/operator/redistribute.py b/HySoP/hysop/operator/redistribute.py
index d805e01eaae9915a081a0d6ac59ca2a91fbf35d6..b78094507e1d710b7c8e5884f1516cb01edf28f4 100644
--- a/HySoP/hysop/operator/redistribute.py
+++ b/HySoP/hysop/operator/redistribute.py
@@ -23,7 +23,6 @@ from parmepy import __VERBOSE__
from parmepy.constants import debug, PARMES_MPI_REAL, ORDERMPI, np, S_DIR
from parmepy.operator.continuous import Operator
from parmepy.mpi.topology import Bridge
-from parmepy.mpi import main_rank, main_comm
from parmepy.methods_keys import Support
@@ -89,7 +88,6 @@ class Redistribute(Operator):
# function of this redistribute. This operator have to know self.
self.opTo.addRedistributeRequirement(self)
-
def discretize(self):
pass
@@ -298,12 +296,13 @@ class Redistribute(Operator):
for rk in br.recvFrom.keys():
self.r_request[v_name + str(rk)] = \
self._main_comm.Irecv([vTo, 1, self._r_types[v][rk]],
- source=rk, tag=rk)
+ source=rk, tag=rk)
self._hasRequests = True
for rk in br.sendTo.keys():
self.s_request[v_name + str(rk)] = \
- self._main_comm.Issend([vFrom, 1, self._s_types[v][rk]],
- dest=rk, tag=self._main_rank)
+ self._main_comm.Issend([vFrom, 1,
+ self._s_types[v][rk]],
+ dest=rk, tag=self._main_rank)
self._hasRequests = True
def _wait_host(self):
diff --git a/HySoP/hysop/operator/redistribute_intercomm.py b/HySoP/hysop/operator/redistribute_intercomm.py
index 959e4d1506b54427ec257275290691e356984892..030d68f0d11fa8ca015b16de6d343a7ce6c76749 100644
--- a/HySoP/hysop/operator/redistribute_intercomm.py
+++ b/HySoP/hysop/operator/redistribute_intercomm.py
@@ -11,7 +11,6 @@ from parmepy.constants import debug, PARMES_MPI_REAL, ORDERMPI, S_DIR, np
from parmepy import __VERBOSE__
from parmepy.operator.continuous import Operator
from parmepy.mpi.topology import Bridge_intercomm
-from parmepy.operator.redistribute import Redistribute
from parmepy.methods_keys import Support
@@ -83,9 +82,10 @@ class RedistributeIntercomm(Operator):
self._s_types[v] = {}
self._toHost_fields = []
self._toDevice_fields = []
+ self._parent_rank = self.parent_comm.Get_rank()
+ self._my_rank = None
def discretize(self):
- self._parent_rank = self.parent_comm.Get_rank()
for v in self.variables:
if self.topo is None:
@@ -345,4 +345,3 @@ class RedistributeIntercomm(Operator):
if not res:
return res
return res
-
diff --git a/HySoP/hysop/operator/stretching.py b/HySoP/hysop/operator/stretching.py
index 0725416b7dd0fd75b2509eefdc403768a809053f..ce310da71e34b4bca421e84a68b877911122a8db 100755
--- a/HySoP/hysop/operator/stretching.py
+++ b/HySoP/hysop/operator/stretching.py
@@ -8,10 +8,10 @@ Computation of stretch. term in Navier-Stokes.
from parmepy.constants import debug
from parmepy.methods_keys import TimeIntegrator, Formulation, \
SpaceDiscretisation
-from parmepy.numerics.integrators.runge_kutta3 import RK3
from parmepy.numerics.finite_differences import FD_C_4
-from parmepy.operator.discrete.stretching import Conservative
from parmepy.operator.continuous import Operator
+from parmepy.operator.discrete.stretching import Conservative, GradUW
+import parmepy.default_methods as default
import numpy as np
@@ -49,8 +49,7 @@ class Stretching(Operator):
self.resolutions = resolutions
## Numerical methods for time and space discretization
if method is None:
- method = {TimeIntegrator: RK3, Formulation: Conservative,
- SpaceDiscretisation: FD_C_4}
+ method = default.STRETCHING
self.method = method
assert Formulation in self.method.keys()
assert SpaceDiscretisation in self.method.keys()
@@ -58,7 +57,10 @@ class Stretching(Operator):
## Formulation used for the stretching equation.
## Default = conservative form.
- self.formulation = self.method[Formulation]
+ if self.method[Formulation] == "GradUW":
+ self.formulation = GradUW
+ else:
+ self.formulation = Conservative
self.input = [self.velocity, self.vorticity]
self.output = [self.vorticity]
@@ -91,14 +93,27 @@ class Stretching(Operator):
self.discreteFields[v] = v.discretize(topo)
@staticmethod
- def getWorkLengths(formulation, timeIntegrator):
+ def getWorkLengths(method=None, domain_dim=None):
"""
Return the length of working arrays lists required
for stretching discrete operator, depending on :
- the formulation (Conservative or GradUW)
- the time integrator (RK3, ...)
+ @param method : the dict of parameters for the operator.
+ Default = {TimeIntegrator: RK3, Formulation: Conservative,
+ SpaceDiscretisation: FD_C_4}
"""
- return formulation.getWorkLengths(timeIntegrator)
+ if method is None:
+ method = default.STRETCHING
+ assert Formulation in method,\
+ 'A formulation is required for the stretching.'
+ assert TimeIntegrator in method,\
+ 'A time integrator is required for the stretching.'
+ if method[Formulation] == "GradUW":
+ formulation = GradUW
+ else:
+ formulation = Conservative
+ return formulation.getWorkLengths(method[TimeIntegrator])
def setWorks(self, rwork=None, iwork=None):
if rwork is None:
@@ -109,12 +124,10 @@ class Stretching(Operator):
@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)
+ method=self.method)
self.discreteOperator.setUp()
self._isUpToDate = True
diff --git a/HySoP/hysop/problem/tests/test_simulation.py b/HySoP/hysop/problem/tests/test_simulation.py
new file mode 100644
index 0000000000000000000000000000000000000000..38efa286c5c0689e406d979944ebbb7b1aafce9e
--- /dev/null
+++ b/HySoP/hysop/problem/tests/test_simulation.py
@@ -0,0 +1,52 @@
+"""
+@file test_simulation.py
+tests simulation incr and io_utils writer
+"""
+from parmepy.problem.simulation import Simulation
+from parmepy.tools.io_utils import Writer
+
+simu = Simulation(tinit=0.0, tend=1.0, nbIter=10)
+
+
+def test_simu_incr():
+ wr = Writer({'frequency': 2})
+ assert wr.doWrite(simu.currentIteration)
+
+ simu.initialize()
+
+ assert not wr.doWrite(simu.currentIteration)
+
+ count = 1
+ while not simu.isOver:
+ if count % 2 == 0:
+ assert wr.doWrite(simu.currentIteration)
+ else:
+ assert not wr.doWrite(simu.currentIteration)
+ simu.printState()
+ simu.advance()
+ count += 1
+ assert simu.currentIteration == 10
+ simu.finalize()
+ assert wr.doWrite(simu.currentIteration)
+
+
+def test_simu_incr2():
+ wr = Writer({'frequency': 3})
+ assert wr.doWrite(simu.currentIteration)
+ simu.timeStep = 0.10000000001
+ simu.initialize()
+
+ assert not wr.doWrite(simu.currentIteration)
+
+ count = 1
+ while not simu.isOver:
+ if count % 3 == 0:
+ assert wr.doWrite(simu.currentIteration)
+ else:
+ assert not wr.doWrite(simu.currentIteration)
+ simu.printState()
+ simu.advance()
+ count += 1
+ assert simu.currentIteration == 10
+ simu.finalize()
+ assert wr.doWrite(simu.currentIteration)
diff --git a/HySoP/setup.py.in b/HySoP/setup.py.in
index 3c75c97226a0edb1984fd3a8d8261a6be509fcda..f771350fcb802f916433b5e6029c257c9ebc57ea 100644
--- a/HySoP/setup.py.in
+++ b/HySoP/setup.py.in
@@ -74,12 +74,17 @@ if(enable_fortran is "ON"):
parmeslib.append('fftw3')
parmeslib.append('fftw3_mpi')
parmes_libdir.append(fftwdir)
-
+ else:
+ packages.append('parmepy.fakef2py')
+ packages.append('parmepy.fakef2py.fftw2py')
withscales = '@WITH_SCALES@'
if(withscales is "ON"):
if(withfftw is "OFF"):
fortran_src.append(fortran_dir+'parameters.f90')
fortran_src.append(fortran_dir+'scales2py.f90')
+ else:
+ packages.append('parmepy.fakef2py')
+ packages.append('parmepy.fakef2py.scales2py')
options = [('F2PY_REPORT_ON_ARRAY_COPY', '1')]
if(os.uname()[0] == 'Linux'):
@@ -94,9 +99,9 @@ if(enable_fortran is "ON"):
define_macros=options)
ext_modules.append(parpyModule)
else:
- packages.append('parmepy.f2py')
- packages.append('parmepy.f2py.scales2py')
- packages.append('parmepy.f2py.fftw2py')
+ packages.append('parmepy.fakef2py')
+ packages.append('parmepy.fakef2py.scales2py')
+ packages.append('parmepy.fakef2py.fftw2py')
data_files = []
diff --git a/HySoP/src/CMakeLists.txt b/HySoP/src/CMakeLists.txt
index c03b48f1d21157d920cf0b3356bb1ad10bee29a4..b758e96894258ee8475adf2a90aa4f8a812389e8 100644
--- a/HySoP/src/CMakeLists.txt
+++ b/HySoP/src/CMakeLists.txt
@@ -34,7 +34,7 @@ if(WITH_SCALES)
${SCALES_DIR}/layout
${SCALES_DIR}/particles
${SCALES_DIR}/particles/advec_line
- ${SCALES_DIR}/output
+ # ${SCALES_DIR}/output
)
endif()
diff --git a/HySoP/src/fftw/Diffusion.f90 b/HySoP/src/fftw/Diffusion.f90
index fe5cc9316c030749c777030ec33cf3d0847341ef..e9178a16f1d4b7f00d8282cd650e936004e3f3ea 100755
--- a/HySoP/src/fftw/Diffusion.f90
+++ b/HySoP/src/fftw/Diffusion.f90
@@ -12,21 +12,21 @@ module diffusion
contains
- subroutine solveDiffusion2D(nudt,velocity_x,velocity_y,omega)
-
- real(mk), intent(in) :: nudt
- real(mk),dimension(:,:),intent(in) :: velocity_x,velocity_y
- real(mk),dimension(:,:),intent(inout) :: omega
-
- !! Compute fftw forward transform
- !! Omega is used to initialize the fftw buffer for input field.
- call r2c_diffusion_2d(velocity_x,velocity_y)
-
- call filter_diffusion_2d(nudt)
-
- call c2r_diffusion_2d(omega)
-
- end subroutine solveDiffusion2D
+!!$ subroutine solveDiffusion2D(nudt,velocity_x,velocity_y,omega)
+!!$
+!!$ real(mk), intent(in) :: nudt
+!!$ real(mk),dimension(:,:),intent(in) :: velocity_x,velocity_y
+!!$ real(mk),dimension(:,:),intent(inout) :: omega
+!!$
+!!$ !! Compute fftw forward transform
+!!$ !! Omega is used to initialize the fftw buffer for input field.
+!!$ call r2c_diffusion_2d(velocity_x,velocity_y)
+!!$
+!!$ call filter_diffusion_2d(nudt)
+!!$
+!!$ call c2r_diffusion_2d(omega)
+!!$
+!!$ end subroutine solveDiffusion2D
subroutine solveDiffusion3D(nudt,velocity_x,velocity_y,velocity_z,omega_x,omega_y,omega_z, ghosts_v, ghosts_w)
real(mk), intent(in) :: nudt
@@ -59,7 +59,8 @@ contains
call mpi_comm_dup(parent_comm, main_comm, ierr)
if(pbdim == 2) then
- call init_c2c_diffusion_2d(resolution,lengths)
+ !!call init_c2c_diffusion_2d(resolution,lengths)
+ print * , "WARNING, diffusion not yet implemented in 2D"
elseif(pbdim == 3) then
call init_c2c_3d(resolution,lengths)
endif