From 45be6b9b7b3339677c53322d68f1e01924c4b2f2 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Franck=20P=C3=A9rignon?= <franck.perignon@imag.fr>
Date: Thu, 4 Jul 2013 10:08:03 +0000
Subject: [PATCH] Update stretching (two formulations) and integrators
---
HySoP/CMakeLists.txt | 12 +-
HySoP/hysop/constants.py | 4 +
.../hysop/numerics/differential_operations.py | 153 ++++++++++-
HySoP/hysop/numerics/finite_differences.py | 7 +-
HySoP/hysop/numerics/integrators/euler.py | 30 +--
HySoP/hysop/numerics/integrators/odesolver.py | 26 +-
.../numerics/integrators/runge_kutta2.py | 156 +++++------
.../numerics/integrators/runge_kutta3.py | 57 ++--
.../numerics/integrators/runge_kutta4.py | 87 +++---
HySoP/hysop/operator/continuous.py | 1 -
HySoP/hysop/operator/curl.py | 21 +-
HySoP/hysop/operator/discrete/stretching.py | 251 ++++++++++--------
HySoP/hysop/operator/discrete/updateGhosts.py | 6 +-
.../hysop/operator/monitors/compute_forces.py | 2 +-
HySoP/hysop/operator/reprojection.py | 2 +-
HySoP/hysop/operator/stretching.py | 104 +++++---
HySoP/hysop/operator/tests/test_Stretching.py | 8 +-
HySoP/hysop/operator/updateGhosts.py | 2 +-
HySoP/hysop/problem/tests/test_problem.py | 0
19 files changed, 534 insertions(+), 395 deletions(-)
create mode 100644 HySoP/hysop/problem/tests/test_problem.py
diff --git a/HySoP/CMakeLists.txt b/HySoP/CMakeLists.txt
index cc669d2a9..6d3af4846 100644
--- a/HySoP/CMakeLists.txt
+++ b/HySoP/CMakeLists.txt
@@ -180,7 +180,7 @@ else()
endif()
# Target to remove parmepy from install-path.
-add_custom_target(python-cleaninstall COMMAND rm -rf ${${PROJECT_NAME}_INSTALL_DIR}/${PROJECT_NAME}*
+add_custom_target(python-cleaninstall COMMAND rm -rf ${CMAKE_INSTALL_PREFIX}/${PROJECT_NAME}*
COMMENT "remove parmepy package and its dependencies")
# Target to clean sources (remove .pyc files) and build dir.
@@ -246,9 +246,17 @@ if(VERBOSE_MODE)
message(STATUS " Project uses Scales : ${WITH_SCALES}")
message(STATUS " Project uses FFTW : ${WITH_FFTW}")
message(STATUS " Project uses GPU : ${WITH_GPU}")
- message(STATUS " Project will be built in ${${PROJECT_NAME}_PYTHON_BUILD_DIR}.")
+ message(STATUS " Project will be built in build/${${PROJECT_NAME}_PYTHON_BUILD_DIR}.")
message(STATUS " Python packages will be installed in ${CMAKE_INSTALL_PREFIX}")
message(STATUS " ${PROJECT_NAME} debug mode : ${DEBUG}")
message(STATUS "====================== ======= ======================")
+ message(STATUS " ")
+ message(STATUS "Try :")
+ message(STATUS " 'make -jN' to build the project, N being the number of available processes.")
+ message(STATUS " 'make python-install' to install python modules and their dependencies. ")
+ message(STATUS " 'make doc' to generate doxygen documentation for parmepy.")
+ message(STATUS " 'make test' to run some test (after the build!).")
+ message(STATUS " 'make clean' to clean build directory.")
+ message(STATUS " 'make python-cleaninstall' to clean install directory.")
endif()
diff --git a/HySoP/hysop/constants.py b/HySoP/hysop/constants.py
index ebfe50a4f..7bd6de176 100644
--- a/HySoP/hysop/constants.py
+++ b/HySoP/hysop/constants.py
@@ -31,6 +31,10 @@ ZDIR = 2
PERIODIC = 0
## Directions string
S_DIR = ["_X", "_Y", "_Z"]
+## Stretching formulation (div(w:u))
+CONSERVATIVE = 1
+## Stretching formulation ([grad(u)][w])
+GRADUW = 0
#define debug decorator:
diff --git a/HySoP/hysop/numerics/differential_operations.py b/HySoP/hysop/numerics/differential_operations.py
index cec05fb5e..d9df80826 100755
--- a/HySoP/hysop/numerics/differential_operations.py
+++ b/HySoP/hysop/numerics/differential_operations.py
@@ -1,4 +1,4 @@
-# -*- codingind utf-8 -*-
+# -*- coding: utf-8 -*-
"""
@file differential_operations.py
Library of functions used to perform classical vector calculus
@@ -7,6 +7,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
+import numpy as np
class DifferentialOperation(object):
@@ -18,15 +19,17 @@ class DifferentialOperation(object):
@debug
@abstractmethod
- def __init__(self):
+ def __init__(self, topo, method=None):
pass
-
+
class Curl(DifferentialOperation):
"""
Computes \f$ \nabla \ times (f) \f$, f being a vector field.
"""
def __init__(self, topo, method=None):
+
+ DifferentialOperation.__init__(self, topo)
if method is not None and method.find('FD_order2') >= 0:
raise ValueError("2nd order scheme Not yet implemented")
else:
@@ -36,19 +39,21 @@ class Curl(DifferentialOperation):
self.fd_scheme = FD_C_4(topo.mesh.space_step)
self.indices = topo.mesh.iCompute
self.fd_scheme.computeIndices(self.indices)
+ self.nbComponents = topo.domain.dimension
def __call__(self, variable, data=None):
return self.fcall(variable, data)
def FDCentral4(self, variable, data=None):
"""
+ @param variable : a list of numpy arrays
"""
# Get slices for subarrays ...
# init (if required) return data
if data is None:
data = []
- for i in xrange(variable.dimension):
- data.append(variable.data[i].copy())
+ for i in xrange(self.nbComponents):
+ data.append(variable[i].copy())
#data[i][...] = 0.0
#data[0][...] = 0.0
@@ -100,6 +105,7 @@ class DivV(DifferentialOperation):
self.fd_scheme.computeIndices(self.indices)
assert (topo.ghosts >= 2).all(),\
'you need a ghost layer for FD4 scheme.'
+ self.nbComponents = topo.domain.dimension
def __call__(self, var1, scal, data=None, work=None):
return self.fcall(var1, scal, data)
@@ -109,9 +115,9 @@ class DivV(DifferentialOperation):
Compute central finite difference scheme, order 4.
"""
if data is None:
- data = var1.data[0].copy()
+ data = np.zeros_like(var1[0])
if work is None:
- work = var1.data[0].copy()
+ work = np.zeros_like(var1[0])
#data[...] = 0.0
@@ -129,9 +135,9 @@ class DivV(DifferentialOperation):
"""
"""
if data is None:
- data = var1.data[0].copy()
-
- data[...] = 0.0
+ data = np.zeros_like(var1[0])
+ else:
+ data[...] = 0.0
# d/dx (scal * v1x)
self.fd_scheme.accumulate(var1[XDIR] * scal, XDIR, data)
@@ -154,6 +160,7 @@ class DivT(DifferentialOperation):
self.fcall = self.FDCentral4
self.div = DivV(topo, method='FD_C4')
self.divZ = DivV(topo, method='FD_C4_nowork')
+ self.nbComponents = topo.domain.dimension
def __call__(self, var1, var2, data=None):
return self.fcall(var1, var2, data)
@@ -164,8 +171,9 @@ class DivT(DifferentialOperation):
# init (if required) return data
if data is None:
data = []
- for i in xrange(var1.dimension):
- data.append(var1.data[i].copy())
+ for i in xrange(self.nbComponents):
+ data.append(np.zeros_like(var1[i]))
+
# we use data[Y]Â and data[Z] as work ...
data[XDIR] = self.div(var1, var2[XDIR],
data=data[XDIR], work=data[YDIR])
@@ -174,3 +182,124 @@ class DivT(DifferentialOperation):
# No work vector for zdir
data[ZDIR] = self.divZ(var1, var2[ZDIR], data=data[ZDIR])
return data
+
+
+class GradS(DifferentialOperation):
+ """
+ Computes \f$ \nabla(\rho) \f$, \f$ \rho\f$ a scalar field
+ """
+ def __init__(self, topo, method=None):
+ if method is not None and method.find('FD_order2') >= 0:
+ raise ValueError("2nd order scheme Not yet implemented")
+ elif method.find('FD_C4') >= 0:
+ self.fcall = self.FDCentral4
+ self.fd_scheme = FD_C_4((topo.mesh.space_step))
+ else:
+ raise ValueError("FD scheme Not yet implemented")
+
+ self.indices = topo.mesh.iCompute
+ self.fd_scheme.computeIndices(self.indices)
+ assert (topo.ghosts >= 2).all(),\
+ 'you need a ghost layer for FD4 scheme.'
+ ## number of components for vector fields
+ self.nbComponents = topo.domain.dimension
+
+ def __call__(self, scal, data=None, work=None):
+ return self.fcall(scal, data)
+
+ def FDCentral4(self, scal, data=None):
+ """
+ Compute central finite difference scheme, order 4.
+ @param scal a numpy array
+ """
+ if data is None:
+ data = []
+ for i in xrange(self.nbComponents):
+ data.append(np.zeros_like(scal))
+ #data[...] = 0.0
+
+ # d/dx (scal), saved in data
+ self.fd_scheme.inplace(scal, XDIR, data[XDIR])
+ # d/dy (scal), saved in data
+ self.fd_scheme.inplace(scal, YDIR, data[YDIR])
+ # d/dz(scal), saved in data
+ self.fd_scheme.inplace(scal, ZDIR, data[ZDIR])
+ return data
+
+
+class GradVxW(DifferentialOperation):
+ """
+ Computes \f$ [\nabla(V)][W] \f$, with
+ \f$V\f$ and \f$W\f$ some vector fields.
+ """
+ def __init__(self, topo, method=None):
+ if method is not None and method.find('FD_order2') >= 0:
+ raise ValueError("2nd order scheme Not yet implemented")
+ else:
+ self.fcall = self.FDCentral4
+ self.grad = GradS(topo, method='FD_C4')
+
+ self.divZ = DivV(topo, method='FD_C4_nowork')
+ self.nbComponents = topo.domain.dimension
+
+ if self.nbComponents == 1:
+ self.worklength = 1
+ elif self.nbComponents == 2:
+ self.worklength = 3
+ else:
+ self.worklength = 5
+ self.prod = BlockProd()
+
+ def __call__(self, var1, var2, diagnostics=None, data=None):
+ return self.fcall(var1, var2, diagnostics, data)
+
+ def FDCentral4(self, var1, var2, diagnostics=None, data=None):
+ """
+ """
+ # init (if required) return data
+ if data is None:
+ data = []
+ for i in xrange(self.worklength):
+ data.append(np.zeros_like(var1[0]))
+ else:
+ assert len(data) == self.worklength
+ if diagnostics is not None:
+ diagnostics[:] = 0.0
+
+ for direc in xrange(self.nbComponents):
+ data[direc:self.nbComponents + direc] = \
+ self.grad(var1[direc],
+ data=data[direc:self.nbComponents + direc])
+
+ if diagnostics is not None:
+ diagnostics[0] = max(diagnostics[0],
+ sum([abs(data[i])
+ for i in xrange(direc,
+ self.nbComponents
+ + direc)]))
+
+ # compute grad(Vx).var2
+ data[direc] = self.prod(data[direc:self.nbComponents + direc],
+ var2, data[direc])
+ return data
+
+
+class BlockProd(object):
+ """
+ Block product of two lists of numpy arrays
+ """
+ def __init__(self):
+ pass
+
+ def __call__(self, list1, list2, data=None):
+ assert len(list1) == len(list2)
+
+ if data is not None:
+ assert data.shape == list1[0].shape
+ data[...] = 0.0
+ else:
+ data = np.zeros_like(list1[0])
+
+ for i in xrange(len(list1)):
+ data += list1[i] * list2[i]
+ return data
diff --git a/HySoP/hysop/numerics/finite_differences.py b/HySoP/hysop/numerics/finite_differences.py
index 9d1dd160a..7f2ebf3e9 100644
--- a/HySoP/hysop/numerics/finite_differences.py
+++ b/HySoP/hysop/numerics/finite_differences.py
@@ -29,10 +29,11 @@ class FD_C_4(FiniteDifference):
def __init__(self, step):
"""
Constructor
- @param step : grid step size
+ @param step : a numpy array with grid step size in each dir.
"""
# dim of the field on which scheme will be applied
# (i.e dim of the domain)
+ step = np.asarray(step)
self._dim = step.size
# list of indices for index - 1
self._m1 = []
@@ -73,7 +74,7 @@ class FD_C_4(FiniteDifference):
def inplace(self, tab, cdir, result):
"""
- Compute \f$ \frac{\partial d tab()}{\partial cdir} \f$ using fd scheme,
+ Compute \f$ \frac{\partial tab()}{\partial cdir} \f$ using fd scheme,
The result is saved in input/ouput parameter result (in place!)
@param tab : a numpy array
@param cdir : direction of differentiation
@@ -88,7 +89,7 @@ class FD_C_4(FiniteDifference):
def accumulate(self, tab, cdir, result):
"""
- Compute \f$ \frac{\partial d tab()}{\partial cdir} \f$ using fd scheme,
+ Compute \f$ \frac{\partial tab()}{\partial cdir} \f$ using fd scheme,
The result is ADDED to input/ouput parameter result.
@param tab : a numpy array
@param cdir : direction of differentiation
diff --git a/HySoP/hysop/numerics/integrators/euler.py b/HySoP/hysop/numerics/integrators/euler.py
index 57bbb3f91..20102a9fb 100755
--- a/HySoP/hysop/numerics/integrators/euler.py
+++ b/HySoP/hysop/numerics/integrators/euler.py
@@ -4,7 +4,6 @@
Forward Euler method.
"""
-from parmepy.constants import np, PARMES_REAL, ORDER
from parmepy.numerics.integrators.odesolver import ODESolver
@@ -17,35 +16,18 @@ class Euler(ODESolver):
y(t_{n+1}) = y(t_n) + dt*f(y(t_n))
"""
- def __init__(self, f=None, dim=1):
+ def __init__(self, dim, f=None):
"""
- Constructor.
-
@param f function f(t,y) : Right hand side of the equation to solve.
@param dim : dimensions
- @param conditions : function to apply boundary conditions.
"""
- ODESolver.__init__(self, f)
+ ODESolver.__init__(self, dim, f)
self.name = 'euler'
- def __call__(self, f, t, dt, y):
- """
- Integration step for Euler Method.
- up = f[t, u(t,x,y,z)].
-
- @param u : function of position at time t.
- @param f : function of function of position at time t.
- @param t : current time.
- @param dt : time step.
-# """
-# result = y + dt * np.asarray(f(t, y), dtype=PARMES_REAL, order=ORDER)
-# return result
-
- RHS = f(t, y)
-
- return y[0] + dt * RHS[0], \
- y[1] + dt * RHS[1], \
- y[2] + dt * RHS[2]
+ def __call__(self, y, dt, t):
+
+ result = self.f(t, y)
+ return [y[i] + dt * result[i] for i in xrange(self._dim)]
def __str__(self):
"""ToString method"""
diff --git a/HySoP/hysop/numerics/integrators/odesolver.py b/HySoP/hysop/numerics/integrators/odesolver.py
index 7d91f24b4..df9538d53 100644
--- a/HySoP/hysop/numerics/integrators/odesolver.py
+++ b/HySoP/hysop/numerics/integrators/odesolver.py
@@ -12,46 +12,40 @@ class ODESolver(NumMethod):
"""
Abstract description for ODE solvers.
Solve the system :
- y'(t) = f(t,y)
+ \f[ \frac{dy(t)}{dt} = f(t,y) \f]
- By default, f = 0.
"""
__metaclass__ = ABCMeta
@abstractmethod
- def __init__(self, f=None):
+ def __init__(self, dim, f=None):
"""
- Constructor.
-
- @param f function f(t,y) : Right hand side of the equation to solve.
+ @param dim, vector field (f) dimension (number of component)
+ @param f : function f(t,y), right hand side of the equation to solve.
+ WARNING : y arg must be a list of numpy arrays and f must return
+ a list of numpy arrays.
"""
super(ODESolver, self).__init__()
## RHS.
if f is not None:
self.f = f
else:
- #[0. for x in y]
- self.f = lambda t, y: np.vectorize(f)(t, y)
+ self.f = lambda t, y: [np.zeros_like(y[i]) for i in xrange(dim)]
+
+ self._dim = dim
@abstractmethod
- def __call__(self, y, t, dt):
+ def __call__(self, y, dt, t):
"""
- Abstract method, apply operaton on a variable.
- Must be implemented by sub-class.
-
@param y : position at time t.
@param t : current time.
@param dt : time step.
@return y at t+dt.
"""
- pass
if __name__ == "__main__":
print __doc__
print "- Provided class : ODESolver (abstract)"
print ODESolver.__doc__
-
-
-
diff --git a/HySoP/hysop/numerics/integrators/runge_kutta2.py b/HySoP/hysop/numerics/integrators/runge_kutta2.py
index 3bad50427..df6793e02 100755
--- a/HySoP/hysop/numerics/integrators/runge_kutta2.py
+++ b/HySoP/hysop/numerics/integrators/runge_kutta2.py
@@ -3,7 +3,6 @@
RK2 method interface.
"""
-from parmepy.constants import np
from parmepy.numerics.integrators.odesolver import ODESolver
@@ -15,98 +14,79 @@ class RK2(ODESolver):
y(t_n+1)= y(t_n) + dt*y'[y(t_n)+dt/2*y'(y(t_n))].
"""
- def __init__(self, f=None, dim=3):
+ def __init__(self, dim, f=None):
"""
- Constructor.
-
+ @param dim, vector field (f) dimension (number of component)
@param f function f(t,y) : Right hand side of the equation to solve.
- @param dim : dimensions
"""
- ODESolver.__init__(self, f)
+ ODESolver.__init__(self, dim, f)
self.name = 'RK2'
- def __call__(self, f, t, dt, y):
-
- """
- Integration step for RK2 Method.
- up = f[t, Y(t,space)].
-
- @param f : right-hand side.
- @param t : current time.
- @param dt : time step.
- @param y : function Y(t,space)
- """
-
-# # K1 = dt *f[t,y(t,x,y,z)]
-# K1 = dt * np.asarray(f(t, y), dtype=PARMES_REAL, order=ORDER)
-# # result = y(t,x,y,z)] + dt *f[t + dt/2 , y(t,x,y,z) + k1/2.]
-
-# result = y + dt * np.asarray(f(t + dt / 2., y + K1 / 2.),
-# dtype=PARMES_REAL, order=ORDER)
-# return result
-
- RHS1 = f(t, y)
- K1 = np.array([dt * RHS1[0], dt * RHS1[1], dt * RHS1[2]])
-
- RHS2 = f(t + dt / 2., y + K1 / 2.)
-
- return y[0] + dt * RHS2[0], \
- y[1] + dt * RHS2[1], \
- y[2] + dt * RHS2[2]
-
-
-class RK2_scalar(ODESolver):
- """
- ODESolver implementation for solving an equation system with RK2 method
- for scalar quantities.
- y'(t) = f(t,y)
-
- y(t_n+1)= y(t_n) + dt*y'[y(t_n)+dt/2*y'(y(t_n))].
-
- """
- def __init__(self, f=None):
- """
- Constructor.
-
- @param f function f(t,y) : Right hand side of the equation to solve.
- @param dim : dimensions
- """
- ODESolver.__init__(self, f)
- self.name = 'RK2'
-
- def __call__(self, t, dt, y, res=None,
- work=None, y_temp=None,
- yprime=None, **f_kwargs):
-
- """
- Integration step for RK2 Method.
- up = f[t, Y(t,space)].
-
- @param t : current time.
- @param dt : time step.
- @param y : function Y(t,space).
- @param res : result (may be allocated when not None).
- @param work : allocated working array.
- @param yprime : y'(t) already computed (avoid 1st call to f)
- @param f_kwargs : Parameters for right-hand side computing.
- """
- if res is None:
- res = np.empty_like(y)
- if work is None:
- work = np.empty_like(y)
- if yprime is None:
- yprime = np.empty_like(y)
- yprime[...] = self.f(t, y, res=yprime, **f_kwargs)
-
- y_temp[...] = y + 0.5 * dt * yprime
- work[...] = self.f(t + 0.5 * dt, y_temp, res=work, **f_kwargs)
-
- res[...] = y
- res[...] += dt * work
-
- def __str__(self):
- """ToString method"""
- return "RK2 Method"
+ def __call__(self, y, dt, t):
+ # k1 = f(t,y)
+ # k2 = f(t + dt/2, y + dt/2 * k1)
+ # result = y + dt * k2
+
+ assert len(y) == self._dim
+ # k1 ... in k2
+ k2 = self.f(t, y)
+ # k2
+ k2 = [y[i] + dt / 2 * k2[i] for i in xrange(self._dim)]
+ k2 = self.f(t + dt / 2., k2)
+
+ return [y[i] + dt * k2[i] for i in xrange(self._dim)]
+
+
+## class RK2_scalar(ODESolver):
+## """
+## ODESolver implementation for solving an equation system with RK2 method
+## for scalar quantities.
+## y'(t) = f(t,y)
+
+## y(t_n+1)= y(t_n) + dt*y'[y(t_n)+dt/2*y'(y(t_n))].
+
+## """
+## def __init__(self, f=None):
+## """
+## Constructor.
+
+## @param f function f(t,y) : Right hand side of the equation to solve.
+## @param dim : dimensions
+## """
+## ODESolver.__init__(self, f)
+## self.name = 'RK2'
+
+## def __call__(self, f, t, dt, y, res=None, work=None, yprime=None,
+## **f_kwargs):
+
+## """
+## Integration step for RK2 Method.
+## up = f[t, Y(t,space)].
+
+## @param f : right-hand side.
+## @param t : current time.
+## @param dt : time step.
+## @param y : function Y(t,space).
+## @param res : result (may be allocated when not None).
+## @param work : allocated working array.
+## @param yprime : y'(t) already computed (avoid 1st call to f)
+## @param f_kwargs : Parameters for right-hand side computing.
+## """
+## if res is None:
+## res = np.empty_like(y)
+## if work is None:
+## work = np.empty_like(y)
+## if yprime is None:
+## work[...] = f(t, y, res=work, **f_kwargs)
+## else:
+## work[...] = yprime
+## work[...] *= 0.5 * dt
+## res[...] = y
+## res[...] += dt * f(t + 0.5 * dt, y + work, **f_kwargs)
+
+## def __str__(self):
+## """ToString method"""
+## return "RK2 Method"
if __name__ == "__main__":
diff --git a/HySoP/hysop/numerics/integrators/runge_kutta3.py b/HySoP/hysop/numerics/integrators/runge_kutta3.py
index d417ab7c2..905710465 100755
--- a/HySoP/hysop/numerics/integrators/runge_kutta3.py
+++ b/HySoP/hysop/numerics/integrators/runge_kutta3.py
@@ -3,64 +3,41 @@
RK3 method interface.
"""
-from parmepy.constants import np, PARMES_REAL, ORDER
from parmepy.numerics.integrators.odesolver import ODESolver
class RK3(ODESolver):
"""
ODESolver implementation for solving an equation system with RK3 method.
- y'(t) = f(t,y)
-
- y(t_n+1)= y(t_n) + dt*y'[y(t_n)+dt/2*y'(y(t_n))].
-
"""
- def __init__(self, f=None, dim=3):
+ def __init__(self, dim, f=None):
"""
- Constructor.
-
@param f function f(t,y) : Right hand side of the equation to solve.
@param dim : dimensions
- @param conditions : function to apply boundary conditions.
"""
- ODESolver.__init__(self, f)
+ ODESolver.__init__(self, dim, f)
self.name = 'RK3'
- def __call__(self, f, t, dt, y):
- """
- Integration step for RK2 Method.
- up = f[t, Y(t,space)].
-
- @param f : right-hand side.
- @param t : current time.
- @param dt : time step.
- @param y : function Y(t,space)
- """
-# # K1 = dt *f[t,y(t,x,y,z)]
-# K1 = dt * np.asarray(f(t, y), dtype=PARMES_REAL, order=ORDER)
-# # result = y(t,x,y,z)] + dt *f[t + dt/2 , y(t,x,y,z) + k1/2.]
-# K2 = dt * np.asarray(f(t + dt / 3., y + K1[:, :, :, :] / 3.),
-# dtype=PARMES_REAL, order=ORDER)
-# K3 = dt * np.asarray(
-# f(t + dt * 2. / 3., y + K2[:, :, :, :] * 2. / 3.),
-# dtype=PARMES_REAL, order=ORDER)
-# # result = y(t,x,y,z)] + 1/4 [ k1 + 3 k3]
-# result = y + 1. / 4. * (K1 + 3. * K3)
-# return result
+ def __call__(self, y, dt, t):
- RHS1 = f(t, y)
- K1 = np.array([dt * RHS1[0], dt * RHS1[1], dt * RHS1[2]])
+ # k1 = f(t,y)
+ # k2 = f(t+dt/3, y + dt/3 *k1))
+ # k3 = f(t + 2/3 *dt , y + 2/3 dt * k2))
+ # result = y + 0.25 * dt * (k1 + 3 * k3)
- RHS2 = f(t + dt / 3., y + K1 / 3.)
- K2 = np.array([dt * RHS2[0], dt * RHS2[1], dt * RHS2[2]])
+ assert len(y) == self._dim
- RHS3 = f(t + dt * 2. / 3., y + K2 * 2. / 3.)
- K3 = np.array([dt * RHS3[0], dt * RHS3[1], dt * RHS3[2]])
+ k1 = self.f(t, y)
+ # k2 ... in k3
+ k3 = [y[i] + dt / 3 * k1[i] for i in xrange(self._dim)]
+ k3 = self.f(t + dt / 3., k3)
+ # k3 ...
+ k3 = [y[i] + 2 * dt / 3 * k3[i] for i in xrange(self._dim)]
+ k3 = self.f(t + dt * 2. / 3., k3)
# result = y(t,x,y,z)] + 1/4 [ k1 + 3 k3]
- return y[0] + 1. / 4. * (K1[0] + 3. * K3[0]), \
- y[1] + 1. / 4. * (K1[1] + 3. * K3[1]), \
- y[2] + 1. / 4. * (K1[2] + 3. * K3[2])
+ return [y[i] + 0.25 * dt * (k1[i] + 3 * k3[i])
+ for i in xrange(self._dim)]
def __str__(self):
"""ToString method"""
diff --git a/HySoP/hysop/numerics/integrators/runge_kutta4.py b/HySoP/hysop/numerics/integrators/runge_kutta4.py
index c75c5545f..c5d33dece 100755
--- a/HySoP/hysop/numerics/integrators/runge_kutta4.py
+++ b/HySoP/hysop/numerics/integrators/runge_kutta4.py
@@ -4,75 +4,59 @@
RK4 method interface.
"""
-from parmepy.constants import np, PARMES_REAL, ORDER
from parmepy.numerics.integrators.odesolver import ODESolver
class RK4(ODESolver):
"""
ODESolver implementation for solving an equation system with RK4 method.
- y'(t) = f(t,y)
-
- y(t_n+1)= y(t_n) + dt*y'[y(t_n)+dt/2*y'(y(t_n))].
-
"""
- def __init__(self, f=None, dim=3):
+ def __init__(self, dim, f=None):
"""
- Constructor.
-
@param f function f(t,y) : Right hand side of the equation to solve.
@param dim : dimensions
- @param conditions : function to apply boundary conditions.
"""
- ODESolver.__init__(self, f)
+ ODESolver.__init__(self, dim, f)
self.name = 'RK4'
- def __call__(self, f, t, dt, y):
- """
- Integration step for RK4 Method.
- up = f[t, Y(t,space)].
-
- @param f : right-hand side.
- @param t : current time.
- @param dt : time step.
- @param y : function Y(t,space)
- """
-
-# # K1 = dt *f[t, y(t, x, y, z)]
-# K1 = dt * np.asarray(f(t, y), dtype=PARMES_REAL, order=ORDER)
-# K2 = dt * np.asarray(f(t + dt / 2., y + K1[:, :, :, :] / 2.),
-# dtype=PARMES_REAL, order=ORDER)
-# K3 = dt * np.asarray(f(t + dt / 2., y + K2[:, :, :, :] / 2.),
-# dtype=PARMES_REAL, order=ORDER)
-# K4 = dt * np.asarray(f(t + dt, y + K3[:, :, :, :]),
-# dtype=PARMES_REAL, order=ORDER)
-# # result = y(t, x, y, z)] + dt *f[t + dt/2 , y(t, x, y, z) + k1/2.]
-# result = y + 1. / 6. * (K1[:, :, :, :] + 2. * K2[:, :, :, :] +
-# 2. * K3[:, :, :, :] + K4[:, :, :, :])
-# return result
-
- RHS1 = f(t, y)
- K1 = np.array([dt * RHS1[0], dt * RHS1[1], dt * RHS1[2]])
-
- RHS2 = f(t + dt / 2., y + K1 / 1.)
- K2 = np.array([dt * RHS2[0], dt * RHS2[1], dt * RHS2[2]])
-
- RHS3 = f(t + dt * 2., y + K2 * 2.)
- K3 = np.array([dt * RHS3[0], dt * RHS3[1], dt * RHS3[2]])
-
- RHS4 = f(t + dt, y + K3)
- K4 = np.array([dt * RHS4[0], dt * RHS4[1], dt * RHS4[2]])
-
- # result = y(t, x, y, z)] + dt *f[t + dt/2 , y(t, x, y, z) + k1/2.]
- return y[0] + 1. / 6. * (K1[0] + 2. * K2[0] + 2. * K3[0] + K4[0]), \
- y[1] + 1. / 6. * (K1[1] + 2. * K2[1] + 2. * K3[1] + K4[1]), \
- y[2] + 1. / 6. * (K1[2] + 2. * K2[2] + 2. * K3[2] + K4[2])
+ def __call__(self, y, dt, t):
+
+ # k1 = f(t,y)
+ # k2 = f(t + dt/2, y + dt/2 * k1)
+ # k3 = f(t + dt/2, y + dt/2 * k2)
+ # k4 = f(t + dt, y + dt * k3)
+ # result = y + dt/6( k1 + 2 * k2 + 2 * k3 + k4)
+
+ # k1 in result
+ result = self.f(t, y)
+ # k2 in work
+ work = [y[i] + dt / 2 * result[i] for i in xrange(self._dim)]
+ work = self.f(t + dt / 2, work)
+
+ # k1 + 2 * k2
+ for i in xrange(self._dim):
+ result[i] += 2 * work[i]
+ # k3 in work
+ work = [y[i] + dt / 2 * work[i] for i in xrange(self._dim)]
+ work = self.f(t + dt / 2, work)
+ # k1 + 2 * k2 + 2 * k3
+ for i in xrange(self._dim):
+ result[i] += 2 * work[i]
+ # k4 in work
+ work = [y[i] + dt * work[i] for i in xrange(self._dim)]
+ work = self.f(t + dt, work)
+ # k1 + 2 * k2 + 2 * k3 + k4
+ for i in xrange(self._dim):
+ result[i] += work[i]
+ result[i] *= dt / 6
+ result[i] += y[i]
+
+ return result
def __str__(self):
"""ToString method"""
return "RK4 Method"
-
class RK4_scalar(ODESolver):
"""
ODESolver implementation for solving an equation system with RK4 method
@@ -132,7 +116,6 @@ class RK4_scalar(ODESolver):
"""ToString method"""
return "RK4 Method"
-
if __name__ == "__main__":
print __doc__
print "- Provided class : RK4"
diff --git a/HySoP/hysop/operator/continuous.py b/HySoP/hysop/operator/continuous.py
index 2e827ca14..7b36c2cbc 100644
--- a/HySoP/hysop/operator/continuous.py
+++ b/HySoP/hysop/operator/continuous.py
@@ -75,7 +75,6 @@ class Operator(object):
## Number of points in the ghost layer
if ghosts is not None:
self.ghosts = np.asarray(ghosts)
- print "ooioii", self.ghosts, self.ghosts.size, self.domain.dimension
assert self.ghosts.size == self.domain.dimension
else:
self.ghosts = None
diff --git a/HySoP/hysop/operator/curl.py b/HySoP/hysop/operator/curl.py
index 7a320a2b0..11a353503 100644
--- a/HySoP/hysop/operator/curl.py
+++ b/HySoP/hysop/operator/curl.py
@@ -6,6 +6,7 @@ Curl of a field.
from parmepy.constants import debug, np
from parmepy.operator.continuous import Operator
from parmepy.operator.discrete.curl import CurlFD
+from parmepy.operator.updateGhosts import UpdateGhosts
class Curl(Operator):
@@ -21,7 +22,7 @@ class Curl(Operator):
@param invar : the input vector field
@param outvar : curl of input vector field
@param resolution : resolution (global) of the fields
- @param method : Method used (default = finite difference, 4th order)
+ @param method : method used (default = finite difference, 4th order)
@param topo : a predefined topology to discretize invar/outvar
"""
@@ -35,15 +36,15 @@ class Curl(Operator):
self.invar = invar
self.outvar = outvar
if self.method is None:
- self.method = 'FD_order4'
+ self.method = 'FD_C4'
+ ## An operator to update ghost points of input var.
+ self.ghostsOperator = None
@debug
def setUp(self):
-
- print "set up for curl ..."
if self.method.find('FD_order2') >= 0:
nbGhosts = 1
- elif self.method.find('FD_order4') >= 0:
+ elif self.method.find('FD_C4') >= 0:
nbGhosts = 2
else:
nbGhosts = 0
@@ -75,8 +76,18 @@ class Curl(Operator):
self.discreteFields[self.outvar],
method=self.method)
self.discreteOperator.setUp()
+ # ghost op
+ self.ghostsOperator = UpdateGhosts(self.input, topo=topo)
+
self._isUpToDate = True
+ def apply(self, simulation=None):
+ # update ghost points, if required
+ if self.ghostsOperator is not None:
+ self.ghostsOperator.apply()
+ # computation ...
+ self.discreteOperator.apply(simulation)
+
if __name__ == "__main__":
print __doc__
print "- Provided class : Curl"
diff --git a/HySoP/hysop/operator/discrete/stretching.py b/HySoP/hysop/operator/discrete/stretching.py
index 48c6fd9b8..b46cb2c61 100755
--- a/HySoP/hysop/operator/discrete/stretching.py
+++ b/HySoP/hysop/operator/discrete/stretching.py
@@ -2,50 +2,64 @@
"""
@file operator/discrete/stretching.py
-Discrete stretching representation
+Discretisation of the stretching operator for two different formulations.
"""
from parmepy.constants import debug
from discrete import DiscreteOperator
-from parmepy.numerics.differentialOperator import DifferentialOperator
-from parmepy.numerics.fct2op import Fct2Op
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.operator.discrete.synchronizeGhosts import SynchronizeGhosts_d
+from parmepy.numerics.differential_operations import DivT, GradVxW
+import numpy as np
+import math
+ceil = math.ceil
-class Stretching_d(DiscreteOperator):
+class ConservativeStretching(DiscreteOperator):
"""
- Stretching operator representation.
+ Discretisation of the following problem :
+ \f{eqnarray*} \frac{\partial\omega}{\partial t} = \nabla.(\omega:v) \f}
"""
@debug
- def __init__(self, velocity, vorticity,
- method='FD_order4 RK3', propertyOp='divConservation'):
+ def __init__(self, velocity, vorticity, method=None):
"""
- Constructor.
- Create a Stretching operator on a given continuous domain.
- Work on a given velocity and vorticity to compute the stretching term.
-
- @param stretch : Stretching operator.
- @param method : Method to use
+ @param velocity : discrete field
+ @param vorticity : discrete field
+ @param method : numerical method for space/time discretizations
"""
- DiscreteOperator.__init__(self, [velocity, vorticity], method,
- name="Stretching")
- ## Velocity.
+ ## velocity discrete field
self.velocity = velocity
- ## Vorticity.
+ ## vorticity discrete field
self.vorticity = vorticity
- # ## input fields
- # self.input = [self.velocity, self.vorticity]
- self.propertyOp = propertyOp
+
+ DiscreteOperator.__init__(self, [self.velocity, self.vorticity],
+ method,
+ name="ConservativeStretching")
+
self.input = [self.velocity, self.vorticity]
self.output = [self.vorticity]
+ # \todo multiresolution case
+ assert self.velocity.topology == self.vorticity.topology,\
+ 'Multiresolution case not yet implemented.'
+
+ # We must have a topo for this operator ...
+ # A function to compute the gradient of a vector field
+ self.divFunc = DivT(self.velocity.topology,
+ method=method)
+
+ # Time integration scheme.
+ def rhs(t, y):
+ return self.divFunc(self.velocity.data, y)
+
+ # TOdo : check fd call with numpy and/or parmes fields.
+ # stability constant
self.cststretch = 0.
+ # set time integration method
if self.method.find('Euler') >= 0:
self.num_method = Euler
self.cststretch = 2.0
@@ -62,99 +76,124 @@ class Stretching_d(DiscreteOperator):
self.num_method = RK3
self.cststretch = 2.5127
+ self.timeIntegrator = self.num_method(self.vorticity.nbComponents, rhs)
+
@timed_function
def apply(self, simulation):
+ if simulation is None:
+ raise ValueError("Missing simulation value for computation.")
+
+ # time step
+ dt = simulation.timeStep
+ # current time
+ t = simulation.time
+ # Call time integrator
+ self.vorticity.data[0][...], self.vorticity.data[1][...], \
+ self.vorticity.data[2][...] = \
+ self.timeIntegrator(self.vorticity.data, dt, t)
+
+
+class GradStretching(DiscreteOperator):
+ """
+ Discretisation of the following problem :
+ \f{eqnarray*} \frac{\partial\omega}{\partial t}=[\nabla(v)][\omega]\f}
+ """
+
+ @debug
+ def __init__(self, velocity, vorticity, method=None):
"""
- Apply Stretching operator.
-
- Stretching algorithm:
- @li 1. discretization of the divergence operator (div wu)
- or w grad(u) : \n
- @li 2. Time integration. Performs a method choose to resolove
- the equation dw/dt = div(wu) or dw/dt = w grad(u).
- -define the function of ODEsolver as the result of the first step
- -integrate with the chosen method (RK4. RK2. euler)
+ @param velocity : discrete field
+ @param vorticity : discrete field
+ @param method : numerical method for space/time discretizations
"""
+ ## velocity discrete field
+ self.velocity = velocity
+ ## vorticity discrete field
+ self.vorticity = vorticity
+
+ DiscreteOperator.__init__(self, [self.velocity, self.vorticity],
+ method,
+ name="ConservativeStretching")
+
+ self.input = [self.velocity, self.vorticity]
+ self.output = [self.vorticity]
+
+ # \todo multiresolution case
+ assert self.velocity.topology == self.vorticity.topology,\
+ 'Multiresolution case not yet implemented.'
+
+ # We must have a topo for this operator ...
+ # A function to compute the gradient of a vector field
+ self.graduv = GradVxW(self.velocity.topology, method=method)
+
+ # Time integration scheme.
+ def rhs(t, y):
+ return self.graduv(self.velocity.data, y)
+
+ # TOdo : check fd call with numpy and/or parmes fields.
+ # stability constant
+ self.cststretch = 0.
+ # set time integration method
+ if self.method.find('Euler') >= 0:
+ self.num_method = Euler
+ self.cststretch = 2.0
+ elif self.method.find('RK2') >= 0:
+ self.num_method = RK2
+ self.cststretch = 2.0
+ elif self.method.find('RK3') >= 0:
+ self.num_method = RK3
+ self.cststretch = 2.5127
+ elif self.method.find('RK4') >= 0:
+ self.num_method = RK4
+ self.cststretch = 2.7853
+ else:
+ self.num_method = RK3
+ self.cststretch = 2.5127
+
+ ## time integrator
+ self.timeIntegrator = self.num_method(self.vorticity.nbComponents, rhs)
+ ## a vector of float with diagnostics computed
+ ## from GradVxW (max div ...)
+ self.diagnostics = np.ones((5))
+
+ @timed_function
+ def apply(self, simulation):
if simulation is None:
raise ValueError("Missing simulation value for computation.")
+ # time step
dt = simulation.timeStep
- if self.num_method is not None:
- condA = self.propertyOp.find('divConservation') < 0
- if condA and self.propertyOp.find('massConservation') < 0:
- print 'No or wrong stretching operator property specified, \
- use default divConservation'
- self.propertyOp = 'divConservation'
-
- if self.propertyOp.find('massConservation') >= 0:
- ## Determination of Stretching time step (stability condition)
- self.ndt = 1
- self.dtstretch = dt
-
- ## fct2Op
- self.function = Fct2Op(self.vorticity,
- self.velocity,
- choice='divWU',
- topology=self.velocity.topology)
-
- else: # self.propertyOp == 'divConservation'
-
- ## Space discretization of grad(u)
- self.stretchOp = DifferentialOperator(
- self.vorticity, self.velocity,
- method='FD_order4',
- choice='gradV')
- #une synchro
- self.synchro = SynchronizeGhosts_d(
- self.vorticity.data,
- self.vorticity.topology, 'greaterSend')
- self.synchro.setUp()
- gradResult, maxArray, divergence = self.stretchOp(self.synchro)
-
- print "maxadv, maxgersh", maxArray[1], maxArray[0]
- print "dtAdvec", 0.5 / maxArray[1]
-
- ## Determination of Stretching time step (stability condition)
- self.ndt = self.stabilityTest(dt, maxArray[0],
- self.num_method, self.cststretch)
- self.dtstretch = dt / float(self.ndt)
-
- ## fct2Op
- self.function = Fct2Op(self.vorticity,
- gradResult,
- method='FD_order4',
- choice='gradV',
- topology=self.vorticity.topology)
- self.synchro.apply()
-
- ## Time integration
- t = simulation.time
- for i in range(self.ndt):
- methodInt = self.num_method(self.function)
-# self.vorticity.data[...] = methodInt(self.function.apply,
-# t, self.dtstretch,
-# self.vorticity.data)
- self.vorticity.data[0][...], self.vorticity.data[1][...], \
- self.vorticity.data[2][...] = methodInt(
- self.function.apply,
- t, self.dtstretch,
- self.vorticity.data)
- self.synchro.apply()
-# print 'norm omega', LA.norm(self.vorticity.data)
-
- def stabilityTest(self, dt, maxi, num_method, cststretch):
- dtstretch = dt
- dtstretch = min(dtstretch, cststretch / maxi)
- if dt == dtstretch:
- print "dt, ndt , dtstretch, upperbound", dt, int(dt / dtstretch),\
- dt / float(int(dt / dtstretch)), cststretch / maxi
- return int(dt / dtstretch)
- else:
- print "dt, ndt , dtstretch, upperbound", dt, int(dt / dtstretch)\
- + 1, dt / float(int(dt / dtstretch) + 1), cststretch / maxi
- return int(dt / dtstretch) + 1
+ # current time
+ t = simulation.time
+ self.graduv(self.velocity, self.vorticity, self.diagnostics)
+ # Compute the number of required subcycles
+ ndt, subdt = self.checkStability(dt)
+ # Call time integrator
+ print "time steps ...", dt, subdt
+ assert sum(subdt) == dt
+ for i in xrange(ndt):
+ self.vorticity.data[0], self.vorticity.data[1],\
+ self.vorticity.data[2] =\
+ self.timeIntegrator(self.vorticity.data, subdt[i], t)
+
+ def checkStability(self, dt):
+ """
+ Computes a stability condition depending on some
+ diagnostics (from GradVxW)
+ @param current time step
+ @return int,float : the number of required subcycles and
+ the subcycle time-step.
+ """
+ dt_stab = min(dt, self.cststretch / self.diagnostics[0])
+ nb_cycles = int(ceil(dt / dt_stab))
+ subdt = np.zeros((nb_cycles))
+ subdt[:] = dt_stab
+ subdt[-1] = dt - (nb_cycles - 1) * dt_stab
+ return nb_cycles, subdt
if __name__ == "__main__":
print __doc__
- print "- Provided class : Stretching_d"
- print Stretching_d.__doc__
+ print "- Provided class : ConservativeStretching, GradStretching"
+ print ConservativeStretching.__doc__
+
+
diff --git a/HySoP/hysop/operator/discrete/updateGhosts.py b/HySoP/hysop/operator/discrete/updateGhosts.py
index 888e5cf0f..7c7b40306 100644
--- a/HySoP/hysop/operator/discrete/updateGhosts.py
+++ b/HySoP/hysop/operator/discrete/updateGhosts.py
@@ -117,8 +117,10 @@ class UpdateGhosts_d(DiscreteOperator):
comm = self.topology.comm
rank = self.topology.rank
neighbours = self.topology.neighbours
- exchange_dir = [d for d in xrange(self._dim)
- if self.topology.cutdir[d]]
+ exchange_dir = []
+ if self.topology.size > 1:
+ exchange_dir = [d for d in xrange(self._dim)
+ if self.topology.cutdir[d]]
i = 0
for d in exchange_dir:
# 1 - Fill in buffers
diff --git a/HySoP/hysop/operator/monitors/compute_forces.py b/HySoP/hysop/operator/monitors/compute_forces.py
index 5e99259b0..de36f634e 100644
--- a/HySoP/hysop/operator/monitors/compute_forces.py
+++ b/HySoP/hysop/operator/monitors/compute_forces.py
@@ -64,7 +64,7 @@ class Compute_forces(Monitoring):
nbGhosts = 0
if self.method.find('FD_order2') >= 0:
nbGhosts = 1
- elif self.method.find('FD_order4') >= 0:
+ elif self.method.find('FD_C4') >= 0:
nbGhosts = 2
else:
nbGhosts = 2
diff --git a/HySoP/hysop/operator/reprojection.py b/HySoP/hysop/operator/reprojection.py
index 87e3a2e99..ead25a59f 100644
--- a/HySoP/hysop/operator/reprojection.py
+++ b/HySoP/hysop/operator/reprojection.py
@@ -48,7 +48,7 @@ class Reprojection_criterion(Monitoring):
nbGhosts = 0
if self.method.find('FD_order2') >= 0:
nbGhosts = 1
- elif self.method.find('FD_order4') >= 0:
+ elif self.method.find('FD_C4') >= 0:
nbGhosts = 2
else:
nbGhosts = 2
diff --git a/HySoP/hysop/operator/stretching.py b/HySoP/hysop/operator/stretching.py
index 4b381ab67..fd888b29d 100755
--- a/HySoP/hysop/operator/stretching.py
+++ b/HySoP/hysop/operator/stretching.py
@@ -2,33 +2,40 @@
"""
@file operator/stretching.py
-Stretching operator representation
-
+Computation of stretch. term in Navier-Stokes.
"""
-from parmepy.constants import debug
+from parmepy.constants import debug, GRADUW, CONSERVATIVE
from parmepy.operator.continuous import Operator
-from parmepy.operator.discrete.stretching import Stretching_d
+from parmepy.operator.discrete.stretching import ConservativeStretching
+from parmepy.operator.discrete.stretching import GradStretching
+from parmepy.operator.updateGhosts import UpdateGhosts
import numpy as np
class Stretching(Operator):
"""
- Stretching operator representation
\todo write latex formulas
"""
@debug
- def __init__(self, velocity, vorticity,
- resolutions=None, method='FD_order4 RK3',
- topo=None, ghosts=None, **other_config):
+ def __init__(self, velocity, vorticity, resolutions,
+ formulation=None, method=None, topo=None, ghosts=None,
+ ghostUpdate=None):
"""
- Constructor.
Create a Stretching operator from given
velocity and vorticity variables.
+ 'FD_C4 RK3',
- @param velocity ContinuousVectorField : velocity variable.
- @param vorticity ContinuousVectorField : vorticity variable.
+ @param velocity field
+ @param vorticity field
+ @param resolutions : grid resolution of velocity and vorticity
+ @param method : solving method
+ (default = finite differences, 4th order, in space
+ and Runge-Kutta 3 in time.)
+ @param topo : a predefined topology to discretize velocity/vorticity
+ @param ghosts : number of ghosts points. Default depends on the method.
+ Autom. computed if not set.
"""
Operator.__init__(self, [velocity, vorticity], method, topo=topo,
ghosts=ghosts, name="Stretching")
@@ -36,64 +43,91 @@ class Stretching(Operator):
self.velocity = velocity
## vorticity variable (vector)
self.vorticity = vorticity
+ ## Grid resolution for each variable (dictionnary)
self.resolutions = resolutions
- self.config = other_config
+ ## Numerical methods for time and space discretization
+ if method is None:
+ self.method = 'FD_C4_RK3'
+ else:
+ self.method = method
+ ## Formulation used for the stretching equation.
+ ## Default = conservative form.
+ if formulation is not None:
+ self.formulation = formulation
+ else:
+ self.formulation = CONSERVATIVE
+ ## Do this operator deals with ghost points update?
+ ## default = yes
+ if ghostUpdate is None:
+ self.ghostUpdate = True
+ else:
+ self.ghostUpdate = ghostUpdate
self.input = [self.velocity, self.vorticity]
self.output = [self.vorticity]
+ ## operator to update ghost points of the input variables
+ self.synchronize = None
@debug
def setUp(self):
- """
- Stretching operator discretization method.
- Create a discrete Stretching operator from given specifications.
-
- @param idVelocityD : Index of velocity discretisation to use.
- @param idVorticityD : Index of vorticity discretisation to use.
- @param propertyOp : property garantied by the stretching Op
- (divConservation, massConservation) default : divConservation
- @param method : the method to use :
- space : (FD_order2, FD_order4, ..) default : FD_order4
- and time : (Euler, RK, ..) default : RK3
- """
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_order4') >= 0:
+ elif self.method.find('FD_C4') >= 0:
nbGhosts = 2
else:
- nbGhosts = 2
+ raise ValueError("Unknown method for stretching operator.")
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)
+
else:
# Variables discretization
if self.ghosts is not None:
self.ghosts[self.ghosts < nbGhosts] = nbGhosts
else:
- self.ghosts = np.zeros((self.domain.dimension)) * nbGhosts
+ self.ghosts = np.ones((self.domain.dimension)) * nbGhosts
+
for v in self.variables:
- # the topology for v ...
topo = self.domain.getOrCreateTopology(self.domain.dimension,
self.resolutions[v],
ghosts=self.ghosts)
-
- # ... and the corresponding discrete field
self.discreteFields[v] = v.discretize(topo)
- self.discreteOperator = \
- Stretching_d(self.discreteFields[self.velocity],
- self.discreteFields[self.vorticity],
- method=self.method, **self.config)
-
+ if self.formulation is CONSERVATIVE:
+ self.discreteOperator =\
+ ConservativeStretching(self.discreteFields[self.velocity],
+ self.discreteFields[self.vorticity],
+ method=self.method)
+ if self.ghostUpdate:
+ topov = self.discreteFields[self.velocity].topology
+ self.synchronize = UpdateGhosts(self.input, topo=topov)
+ self.synchronize.setUp()
+ elif self.formulation is GRADUW:
+ self.discreteOperator =\
+ GradStretching(self.discreteFields[self.velocity],
+ self.discreteFields[self.vorticity],
+ method=self.method)
+ if self.ghostUpdate:
+ topov = self.discreteFields[self.velocity].topology
+ self.synchronize = UpdateGhosts(self.velocity, topo=topov)
+ self.synchronize.setUp()
self.discreteOperator.setUp()
self._isUpToDate = True
+ def apply(self, simulation=None):
+ # update ghost points, if required
+ if self.synchronize is not None:
+ self.synchronize.apply()
+ # computation ...
+ self.discreteOperator.apply(simulation)
+
+
if __name__ == "__main__":
print __doc__
print "- Provided class : Stretching"
diff --git a/HySoP/hysop/operator/tests/test_Stretching.py b/HySoP/hysop/operator/tests/test_Stretching.py
index db2a4786d..a6d044def 100755
--- a/HySoP/hysop/operator/tests/test_Stretching.py
+++ b/HySoP/hysop/operator/tests/test_Stretching.py
@@ -1,13 +1,9 @@
# -*- coding: utf-8 -*-
-import time
-
import parmepy as pp
import numpy as np
from parmepy.fields.continuous import Field
from parmepy.operator.stretching import Stretching
-from parmepy.problem.navier_stokes import NSProblem
from parmepy.problem.simulation import Simulation
-from parmepy.mpi.main_var import main_comm, main_rank
def computeVel(x, y, z):
@@ -77,7 +73,7 @@ def test_Stretching():
stretch = Stretching(velo, vorti,
resolutions={velo: nbElem,
vorti: nbElem},
- method='FD_order4 RK3',
- propertyOp='divConservation')
+ method='FD_C4_RK3')
stretch.setUp()
simulation = Simulation(0, timeStep, 20)
+ stretch.apply(simulation)
diff --git a/HySoP/hysop/operator/updateGhosts.py b/HySoP/hysop/operator/updateGhosts.py
index ea04b96f3..bbab6918e 100644
--- a/HySoP/hysop/operator/updateGhosts.py
+++ b/HySoP/hysop/operator/updateGhosts.py
@@ -23,7 +23,7 @@ class UpdateGhosts(Operator):
"""
Operator.__init__(self, variables, method, topo=topo)
- self.topo = topo
+ self.topo = self._predefinedTopo
@debug
def setUp(self):
diff --git a/HySoP/hysop/problem/tests/test_problem.py b/HySoP/hysop/problem/tests/test_problem.py
new file mode 100644
index 000000000..e69de29bb
--
GitLab