diff --git a/CMake/FindFFTW.cmake b/CMake/FindFFTW.cmake
index d553d45365ccf65b2ef0c5b0ee843a60b8a6aeb1..ffcb3e4ab80501af1cb0ae4d3de4c3aec73476d2 100644
--- a/CMake/FindFFTW.cmake
+++ b/CMake/FindFFTW.cmake
@@ -38,9 +38,22 @@ find_library(
)
find_library(FFTW_LIBRARY
- NAMES fftw3
+ NAMES fftw3
)
+# Finally the library itself, mpi version
+find_library(
+ FFTW_MPI_LIBRARY
+ NAMES fftw3_mpi
+ PATHS ${fftw_DIR} ${FFTW_INCLUDE_DIR}/.. ENV LIBRARY_PATH ENV LD_LIBRARY_PATH ${FFTW_PKGCONF_LIBRARY_DIRS}
+ PATH_SUFFIXES lib
+ NO_DEFAULT_PATH
+)
+
+find_library(FFTW_MPI_LIBRARY
+ NAMES fftw3_mpi
+)
+
find_library(
FFTWFloat_LIBRARY
NAMES fftw3f
@@ -54,5 +67,5 @@ find_library(FFTWFloat_LIBRARY
set(FFTW_PROCESS_INCLUDES FFTW_INCLUDE_DIR)
-set(FFTW_PROCESS_LIBS FFTW_LIBRARY FFTWFloat_LIBRARY)
+set(FFTW_PROCESS_LIBS FFTW_LIBRARY FFTWFloat_LIBRARY FFTW_MPI_LIBRARY)
libfind_process(FFTW)
diff --git a/HySoP/hysop/box.py b/HySoP/hysop/box.py
deleted file mode 100644
index a4a7fc07c75e7e50310a61d7fdd5534b741276c5..0000000000000000000000000000000000000000
--- a/HySoP/hysop/box.py
+++ /dev/null
@@ -1,58 +0,0 @@
-# -*- coding: utf-8 -*-
-"""
-@package Domain
-Physical domain representation.
-"""
-from ..Param import *
-from ContinuousDomain import ContinuousDomain
-from CartesianGrid import CartesianGrid
-
-
-class Box(ContinuousDomain):
- """
- Periodic box representation.
- ContinuousDomain.ContinuousDomain specialization.
- """
-
- def __init__(self, dimension=1, minPosition=[0.], length=[1.]):
- """
- Constructor.
- Create a Box from a dimension, length and minimum position.
- Parameters dimensions must coincide. Raise ValueError in case of inconsistent parameters dimensions.
-
- @param dimension integer : Box dimension.
- @param length : Box length.
- @param minPosition : Box minimum position.
- """
- if not (dimension == len(length) and dimension == len(minPosition)):
- raise ValueError("Box parameters inconsistents dimensions")
- ContinuousDomain.__init__(self, dimension)
- ## Box length. length = max - min.
- self.length = np.array(length)
- ## Minimum Box position.
- self.min = np.array(minPosition)
- ## Maximum Box position.
- self.max = self.min + self.length
-
- def discretize(self, spec=None):
- """
- Box discretization method.
- Create a Grig.Grid from the Box according to discretization specifications.
-
- @param spec : discretization specifications, element number for each directions.
- """
- if self.discreteDomain is None:
- self.discreteDomain = CartesianGrid(self, spec)
-
- def __str__(self):
- """
- ToString method.
- """
- s = " Box (ContinuousDomain) \n"
- s += " dimension={0}, minPosition={1}, maxPosition={2}, length={3}".format(self.dimension, self.min, self.max, self.length)
- return s
-
-if __name__ == "__main__":
- print __doc__
- print "- Provided class : " + providedClass
- print eval(providedClass).__doc__
diff --git a/HySoP/hysop/domain/__init__.py b/HySoP/hysop/domain/__init__.py
index 53ad3f95bedeafddc49219886cee57692694e86b..aae2d3b688f41c120ab2c4611d911042633aff48 100644
--- a/HySoP/hysop/domain/__init__.py
+++ b/HySoP/hysop/domain/__init__.py
@@ -1,5 +1,6 @@
-from ContinuousDomain import ContinuousDomain
-from Box import Box
-from DiscreteDomain import DiscreteDomain
-from CartesianGrid import CartesianGrid
-from ParticleField import ParticleField
+"""
+@package ParMePy.domain
+
+Everything concerning physical domains, their discretizations and MPI topologies.
+
+"""
diff --git a/HySoP/hysop/domain/box.py b/HySoP/hysop/domain/box.py
new file mode 100644
index 0000000000000000000000000000000000000000..2905edf4ba03929082e1cf5221893e502d969387
--- /dev/null
+++ b/HySoP/hysop/domain/box.py
@@ -0,0 +1,51 @@
+# -*- coding: utf-8 -*-
+from ..constants import PERIODIC
+from .domain import Domain
+from .grid import CartesianGrid
+import numpy as np
+
+class Box(Domain):
+ """
+ Periodic box representation.
+ Note FP : BC as parameter may be better?
+
+ """
+
+ def __init__(self,dimension=3,length=[1.0,1.0,1.0], minPosition=[0.,0.,0.]):
+ """
+ Constructor.
+ Create a Box from a dimension, length and minimum position.
+ Parameters dimensions must coincide. Raise ValueError in case of inconsistent parameters dimensions.
+
+ @param length : Box length.
+ @param minPosition : Box minimum position.
+ """
+ if not (dimension == len(length) and dimension == len(minPosition)):
+ raise ValueError("Box parameters inconsistents dimensions")
+ Domain.__init__(self, dimension)
+ ## Box length. length = max - min.
+ self.length = np.array(length)
+ ## Minimum Box position.
+ self.min = np.array(minPosition)
+ ## Maximum Box position.
+ self.max = self.min + self.length
+ # Boundary conditions type :
+ self.boundaries = np.zeros((self.dimension))
+ self.boundaries[:] = PERIODIC
+
+ def discretize(self,resolution):
+ """Box discretization method."""
+
+ self.discreteDomain = CartesianGrid(resolution,self)
+ return self.discreteDomain
+
+ def __str__(self):
+ """ Doc display. """
+ s = str(self.dimension)+"D box (parallelepipedic or rectangular) domain : \n"
+ s += " minPosition : "+str(self.min)+", maxPosition :"+str(self.max)+", lengthes :"+str(self.length)+"."
+ return s
+
+if __name__ == "__main__":
+ print __doc__
+ print "- Provided class : Box."
+ print Box.__doc__
diff --git a/HySoP/hysop/domain/discrete.py b/HySoP/hysop/domain/discrete.py
index cb88b1d837dda4520a85d49a0e94171bc1826485..9ddaba5f52138dc13fcc15a0b60b8499c0a07d7e 100644
--- a/HySoP/hysop/domain/discrete.py
+++ b/HySoP/hysop/domain/discrete.py
@@ -1,28 +1,32 @@
# -*- coding: utf-8 -*-
"""
-@package Domain
-Physical domain representation.
+@package ParMePy.domain.discrete
+
+Physical domain discretization
+
"""
-from ..Param import *
+from abc import ABCMeta,abstractmethod
+
+class DiscreteDomain(object):
+ """ Abstract description of a discretized physical domain. """
-class DiscreteDomain:
- """
- Abstract description of a discretized physical domain.
- """
+ __metaclass__=ABCMeta
- def __init__(self, dimension):
+ @abstractmethod
+ def __init__(self,domain=None):
"""
Constructor.
Create a DiscreteDomain with a given discretization.
"""
+ self.domain = domain
+ if(domain is not None):
+ self.dimension = domain.dimension
+
## Discretization specifications. Discrete elements number.
- self.elementNumber = None
+ #self.elementNumber = None
## Discretization specifications. Discrete elements size.
- self.elementSize = None
- ## Function to create all Points depending on Box.Box dimension
- self.dimension = dimension
-
+ #self.elementSize = None
if __name__ == "__main__":
print __doc__
diff --git a/HySoP/hysop/domain/domain.py b/HySoP/hysop/domain/domain.py
index ab64a02f4f51d11a5ad92f5e9633ffcf178f904e..126107e90672b891591588280761b7b338890616 100644
--- a/HySoP/hysop/domain/domain.py
+++ b/HySoP/hysop/domain/domain.py
@@ -1,38 +1,34 @@
# -*- coding: utf-8 -*-
-"""
-@package Domain
-Physical domain representation.
-"""
-from ..Param import *
+"""@package ParMePy.domain.domain
+Classes for physical domains description. """
-class ContinuousDomain:
- """
- Abstract description of physical domain.
- """
+from abc import ABCMeta,abstractmethod
+class Domain:
+ """Abstract description of physical domain. """
+
+ __metaclass__=ABCMeta
+
+ @abstractmethod
def __init__(self, dimension):
"""
- Constructor.
- Create a ContinuousDomain with given dimension.
-
+
@param dimension integer : domain dimension.
+
"""
## Domain dimension.
self.dimension = dimension
## Domain discretization.
self.discreteDomain = None
- def discretize(self, spec=None):
+ @abstractmethod
+ def discretize(self,resolution=None):
"""
- Abstract method.
- Must be implemented by sub-class.
-
- @param spec : discretization specifications.
+ @param resolution : discretization specifications.
"""
- raise NotImplementedError("Need to override method in a subclass of " + providedClass)
if __name__ == "__main__":
print __doc__
- print "- Provided class : ContinuousDomain"
- print ContinuousDomain.__doc__
+ print "- Provided class : Domain (abstract)."
+ print Domain.__doc__
diff --git a/HySoP/hysop/domain/grid.py b/HySoP/hysop/domain/grid.py
index 3797c4538e632b19e1b1a2d46c70e73fa2a9b482..db8b12bf502c8d7a1bf639f2e66df6a9262d3dbe 100644
--- a/HySoP/hysop/domain/grid.py
+++ b/HySoP/hysop/domain/grid.py
@@ -1,38 +1,36 @@
# -*- coding: utf-8 -*-
"""
-@package Domain
-Physical domain representation.
-"""
-from ..Param import *
-from DiscreteDomain import DiscreteDomain
+@package ParMePy.domain.grid
+Cartesian grid definition
+"""
+from .discrete import DiscreteDomain
+import numpy as np
class CartesianGrid(DiscreteDomain):
- """
- Discrete periodic box representation as a cartesian grid.
- DiscreteDomain.DiscreteDomain specialization.
- """
+ """ Discrete box representation as a cartesian grid. """
- def __init__(self, box, spec):
- """
- Constructor.
- Create a CartesianGrid from a Box.Box and a discretization specification.
- Select functions regarding dimension.
+ def __init__(self,resolution,domain=None):
+ if(domain is not None):
+ DiscreteDomain.__init__(self,domain)
+ assert(self.dimension==len(resolution))
+ else:
+ self.dimension = len(resolution)
+ DiscreteDomain.__init__(self)
+ self.resolution = resolution
+ ## self.elementNumber = spec
+ ## self.shape = np.append(np.asarray(self.elementNumber), box.dimension)
+ ## self.elementSize = box.length / self.elementNumber
+ ## ## Minimum Grid position.
+ ## self.min = box.min
+ ## ## Maximum Grid position.
+ ## self.max = box.max
+ ## ## Grid length.
+ ## self.length = box.length
+ ## self.axes = [np.asarray(np.arange(self.min[i], self.max[i], self.elementSize[i]), dtype=dtype_real, order=order) for i in xrange(self.dimension)]
- @param box Box.Box : ContinuousDomain.ContinuousDomain discretized.
- @param spec : element number for each directions.
- """
- DiscreteDomain.__init__(self, box.dimension)
- self.elementNumber = spec
- self.shape = np.append(np.asarray(self.elementNumber), box.dimension)
- self.elementSize = box.length / self.elementNumber
- ## Minimum Grid position.
- self.min = box.min
- ## Maximum Grid position.
- self.max = box.max
- ## Grid length.
- self.length = box.length
- self.axes = [np.asarray(np.arange(self.min[i], self.max[i], self.elementSize[i]), dtype=dtype_real, order=order) for i in xrange(self.dimension)]
+ def update(start):
+ self.start = start
def applyConditions(self, pos, dir):
"""
@@ -53,14 +51,14 @@ class CartesianGrid(DiscreteDomain):
@param index list : index to get
@return point with given index
"""
- return self.points[index]
+ return self.__coords.__getitem__(i)
def __str__(self):
"""
ToString method.
"""
- s = "# CartesianGrid of {0} Points".format(np.prod(self.elementNumber))
- return s + "\n"
+ s = str(self.dimension)+"D cartesian grid with resolution equal to : "+str(self.resolution)+".\n"
+ return s
if __name__ == "__main__":
print __doc__
diff --git a/HySoP/hysop/fields/__init__.py b/HySoP/hysop/fields/__init__.py
new file mode 100644
index 0000000000000000000000000000000000000000..e69de29bb2d1d6434b8b29ae775ad8c2e48c5391
diff --git a/HySoP/hysop/fields/analytical.py b/HySoP/hysop/fields/analytical.py
new file mode 100644
index 0000000000000000000000000000000000000000..b38fd01b77acdec6ba15c284d858b1aaffa273a3
--- /dev/null
+++ b/HySoP/hysop/fields/analytical.py
@@ -0,0 +1,59 @@
+# -*- coding: utf-8 -*-
+"""
+@package ParMePy.fields.analytical
+
+"""
+from .continuous import ContinuousField
+from .discrete import DiscreteField
+
+class AnalyticalField(ContinuousField):
+ """
+ A variable (continuous) defined with an analytic formula.
+
+ """
+
+ def __init__(self, dimension, domain, formula=None, scalar=False, name=""):
+ """
+ Constructor.
+ Create an AnalyticalField from a formula.
+
+ @param domain : application domain of the variable.
+ @param dimension : variable dimension.
+ @param formula : function defining the variable.
+ """
+ ContinuousField.__init__(self, dimension, domain)
+ ## Analytic formula.
+ self.formula = formula
+ ## Is scalar variable
+ self.scalar = scalar
+ self.name = name
+
+ def discretize(self, d_spec=None):
+ """
+ AnalyticalField discretization method.
+ Create an DiscreteField.DiscreteField from given specifications.
+
+ @param d_spec : discretization specifications, not used.
+ """
+ if self.discreteField is None:
+ self.discreteField = DiscreteField(domain=self.domain.discreteDomain, dimension=self.dimension, initialValue=self.formula, scalar=self.scalar, spec=d_spec, name=self.name)
+
+ def value(self, pos):
+ """
+ Evaluation of the variable at a given position.
+
+ @param pos : Position to evaluate.
+ @return : value of the formula at the given position.
+ """
+ return self.formula(pos)
+
+ def __str__(self):
+ """ToString method"""
+ s = " Analytical variable (ContinuousField)\n"
+ s += " dimension={0}\n".format(self.dimension)
+ return s + "\n"
+
+if __name__ == "__main__":
+ print __doc__
+ print "- Provided class : AnalyticalField"
+ print AnalyticalField.__doc__
diff --git a/HySoP/hysop/fields/continuous.py b/HySoP/hysop/fields/continuous.py
new file mode 100644
index 0000000000000000000000000000000000000000..bad14fe29943c5239c56d071405d7b596e475a74
--- /dev/null
+++ b/HySoP/hysop/fields/continuous.py
@@ -0,0 +1,41 @@
+# -*- coding: utf-8 -*-
+"""
+@package ParMePy.fields.continuous
+Continuous variable description
+"""
+
+class ContinuousField(object):
+ """
+ A variable defined on a physical domain
+
+ """
+
+ def __init__(self, dimension, domain):
+ """
+ Constructor.
+ Create a ContinuousField.
+
+ @param dimension integer : variable dimension.
+ @param domain : variable application domain.
+ """
+ ## Application domain of the variable.
+ self.domain = domain
+ ## Field dimension.
+ self.dimension = dimension
+ ## Discretisation of this field
+ self.discreteField = None
+ self.name = ""
+
+ def discretize(self, spec=None):
+ """
+ Abstract method.
+ Must be implemented by sub-class.
+
+ @param spec : discretization specifications
+ """
+ raise NotImplementedError("Need to override method in a subclass of " + providedClass)
+
+if __name__ == "__main__":
+ print __doc__
+ print "- Provided class : ContinuousField"
+ print ContinuousField.__doc__
diff --git a/HySoP/hysop/fields/discrete.py b/HySoP/hysop/fields/discrete.py
new file mode 100644
index 0000000000000000000000000000000000000000..3b2641b3169f5d94fc751b6df2b4d359ed7b015e
--- /dev/null
+++ b/HySoP/hysop/fields/discrete.py
@@ -0,0 +1,78 @@
+# -*- coding: utf-8 -*-
+"""
+@package ParMePy.fields.discrete
+
+"""
+
+import numpy as np
+from ..constants import *
+
+class ScalarField(object):
+ """
+ A scalar field, defined on each grid of each subdomain of the given topology.
+
+ """
+
+ def __init__(self,topology,name=""):
+ self.topology = topology
+ self.name = name
+ self.dimension = topology.domain.dimension
+ resolution = self.topology.mesh.resolution
+ self.data = np.zeros((resolution),dtype=PARMES_REAL)
+
+ def __str__(self):
+ """ Display class information """
+ s = self.name+": "+str(self.dimension)+"D scalar field of shape " + str(self.data.shape)
+ return s + "\n"
+
+ def __getitem__(self,i):
+ """ Access to the content of the field """
+ return self.data.__getitem__(i)
+
+ def __setitem__(self,i,value):
+ """
+ Set the content of the field component at position i
+ Usage :
+ A = ScalarField(topo)
+ A[2,1,1] = 12.0 # Calls A.data[2,1,1] = 12.0
+ """
+ self.data.__setitem__(i,value)
+
+class VectorField(object):
+ """
+ A vector field, defined on each grid of each subdomain of the given topology.
+
+ """
+
+ def __init__(self,topology,name=""):
+ self.topology = topology
+ self.name = name
+ self.dimension = topology.domain.dimension
+ resolution = self.topology.mesh.resolution
+ self.data = np.zeros((resolution),dtype=PARMES_REAL)
+
+ def __str__(self):
+ """ Display class information """
+ s = self.name+": "+str(self.dimension)+"D scalar field of shape " + str(self.data.shape)
+ return s + "\n"
+
+ def __getitem__(self,i):
+ """ Access to the content of the field """
+ return self.data.__getitem__(i)
+
+ def __setitem__(self,i,value):
+ """
+ Set the content of the field component at position i
+ Usage :
+ A = ScalarField(topo)
+ A[2,1,1] = 12.0 # Calls A.data[2,1,1] = 12.0
+ """
+ self.data.__setitem__(i,value)
+
+
+
+
+if __name__ == "__main__":
+ print __doc__
+ print "- Provided class : ScalarField"
+ print ScalarField.__doc__
diff --git a/HySoP/hysop/fields/topology.py b/HySoP/hysop/fields/topology.py
new file mode 100644
index 0000000000000000000000000000000000000000..14ad0c79cc3737fa1ece293960f435bd01e8267c
--- /dev/null
+++ b/HySoP/hysop/fields/topology.py
@@ -0,0 +1,133 @@
+"""
+@package ParMePy.domain.topology
+MPI Topologies
+"""
+import numpy as np
+import mpi4py.MPI as MPI
+from ..constants import *
+
+class CartesianTopology(object):
+ """
+
+ Define an MPI cartesian topology with an associated mesh for each sub-domain.
+ Keyword arguments :
+ - domain : the geometry; it must be a box.
+ - resolution : global resolution
+ - dim : dimension of the topology
+ - comm : MPI communicator used to create this topology
+ - periods : periodicity of the topology (in each direction)
+ - ghosts : number of poinst in the ghost layer
+
+ """
+ def __init__(self,domain,resolution,dim=3,comm=MPI.COMM_WORLD,periods=3*(True,),ghosts=np.zeros((3))):
+
+ # Topology dimension
+ self.dim=dim
+ ## Associated domain
+ self.domain=domain
+ # Compute the "optimal" processus distribution for each direction of the grid topology
+ nbprocs = comm.Get_size()
+ self.dims=np.asarray(MPI.Compute_dims(nbprocs,dim))
+ # Define the topology
+ self.topo = comm.Create_cart(self.dims,periods=periods[0:dim])
+ # Size of the topology
+ self.size = self.topo.Get_size()
+ # Rank of the current process, in the new topology
+ self.rank = self.topo.Get_rank()
+ # Coordinates of the current process
+ self.coords = self.topo.Get_coords(self.rank)
+ # Neighbours
+ # self.neighbours = np.zeros((dim,2))
+ self.down,self.up=self.topo.Shift(XDIR,1)
+ if(dim>1) :
+ self.west,self.east=self.topo.Shift(YDIR,1)
+ if(dim>2) :
+ self.south,self.north=self.topo.Shift(ZDIR,1)
+ # ghost layer
+ self.ghosts = ghosts[0:dim]
+ # Global resolution
+ self.resolution = resolution
+ nbpoints = self.resolution[:self.dim]//self.dims
+ remainingPoints = self.resolution[:self.dim]%self.dims
+ localResolution = self.resolution.copy()
+ localResolution[:self.dim] = nbpoints
+ for i in range(dim):
+ if(self.coords[i] == self.dims[i]-1):
+ localResolution[i] = nbpoints[i]+remainingPoints[i]
+ start = np.zeros((domain.dimension))
+ start[:self.dim]=self.coords*nbpoints
+ self.mesh = LocalMesh(self.rank,resolution=localResolution,start=start,ghosts=self.ghosts)
+
+ def __str__(self):
+ """ Topology info display """
+ s ='=======================================================\n'
+ if(self.rank==0):
+ s += str(self.dim)+'D cartesian Topology of size '+str(self.dims)+'\n'
+ s += str(self.mesh)
+ s +='=======================================================\n'
+ return s
+
+class FFTTopology(object):
+ """ A 2D or 3D MPI cartesian topology
+
+ """
+
+ def __init__(self, dim=2,comm=MPI.COMM_WORLD,nbprocs=MPI.COMM_WORLD.Get_size(),periods=3*(True,),dims=(1,1,1)):
+
+ # Topology dimension
+ self.dim=dim
+ # Compute the "optimal" processus distribution for each direction of the grid topology
+ dims=MPI.Compute_dims(nbprocs,dim)
+ self.dims=np.ones((dim+1))
+ self.dims[0:dim] = np.asarray(dims)
+ # Define the topology
+ self.topo = comm.Create_cart(self.dims,periods=periods[0:dim+1])
+ # Size of the topology
+ sizeT = self.topo.Get_size()
+ # Rank of the current process
+ self.rank = self.topo.Get_rank()
+ # Coordinates of the current process
+ self.coords = self.topo.Get_coords(self.rank)
+ # Neighbours
+ # self.neighbours = np.zeros((dim,2))
+ self.down,self.up=self.topo.Shift(XDIR,1)
+ if(dim>1) :
+ self.west,self.east=self.topo.Shift(YDIR,1)
+ if(dim>2) :
+ self.south,self.north=self.topo.Shift(ZDIR,1)
+
+ def __str__(self):
+ """ Topology info display """
+ s ='======================================================='
+ s += str(self.dim)+'D FFT Topology of size '+str(self.dims)
+ s+=self.mesh.__str__
+ s+='======================================================='
+ return s
+
+class LocalMesh(object):
+ """
+ """
+
+ def __init__(self,rank,resolution,start,ghosts=np.zeros((3))) :
+ # Current process rank in the topology that owns this mesh
+ self.rank = rank
+ # Local resolution
+ self.resolution = resolution.copy()+2*ghosts
+ # index of the lower point (in each dir) in the global mesh
+ self.start = start
+ # index of the upper point (in each dir), global mesh
+ self.end = self.start + self.resolution - 1
+ # Mesh step size
+
+ def __str__(self):
+ """ Local mesh display """
+ s = '['+str(self.rank)+'] Local mesh resolution :'+str(self.resolution)+'\n'
+ s+= 'Global indices :'+str(self.start)+'\n'
+ return s
+
+
+if __name__ == "__main__" :
+ print "This module defines the following classes:"
+ print "- CartesianTopology: ", CartesianTopology.__doc__
+ print "- FFTTopology: ", FFTTopology.__doc__
+ print "- LocalMesh: ", LocalMesh.__doc__