Skip to content
Snippets Groups Projects
Commit 30e77327 authored by Jean-Matthieu Etancelin's avatar Jean-Matthieu Etancelin
Browse files

Cleanning parmepy and add some documentation.

parent 88d1147c
No related branches found
No related tags found
No related merge requests found
Showing
with 1196 additions and 968 deletions
This diff is collapsed.
""" """
@package ParMePy @package parmepy
Python package dedicated to flow simulation using particular methods on hybrid architectures (MPI-GPU) Python package dedicated to flow simulation using particular methods on hybrid architectures (MPI-GPU)
...@@ -25,8 +25,8 @@ Box = domain.box.Box ...@@ -25,8 +25,8 @@ Box = domain.box.Box
## Cartesian grid ## Cartesian grid
## MPI topologies and associated meshes ## MPI topologies and associated meshes
import fields.topology import domain.topology
CartesianTopology = fields.topology.CartesianTopology CartesianTopology = domain.topology.CartesianTopology
## Fields ## Fields
import fields.discrete import fields.discrete
...@@ -37,10 +37,10 @@ ContinuousField = fields.continuous.ContinuousField ...@@ -37,10 +37,10 @@ ContinuousField = fields.continuous.ContinuousField
AnalyticalField = fields.analytical.AnalyticalField AnalyticalField = fields.analytical.AnalyticalField
## Operators ## Operators
import operator.Transport import operator.transport
import operator.Velocity import operator.velocity
Transport = operator.Transport.Transport Transport = operator.transport.Transport
Velocity = operator.Velocity.Velocity Velocity = operator.velocity.Velocity
## Problem ## Problem
import problem.problem import problem.problem
......
# -*- coding: utf-8 -*-
""" """
@package parmepy.constants @package parmepy.constants
Constant parameters required for the parmepy package (internal use). Constant parameters required for the parmepy package (internal use).
""" """
...@@ -10,7 +10,6 @@ import math ...@@ -10,7 +10,6 @@ import math
from parmepy.particular_solvers import __path__ as solver_path from parmepy.particular_solvers import __path__ as solver_path
import os import os
#
PI = math.pi PI = math.pi
# Set default type for real and integer numbers # Set default type for real and integer numbers
PARMES_REAL = np.float64 PARMES_REAL = np.float64
......
""" """
@package parmepy.domain @package parmepy.domain
@file parmepy/domain/__init__.py
Everything concerning physical domains, their discretizations and MPI topologies. Everything concerning physical domains, their discretizations and MPI topologies.
......
# -*- coding: utf-8 -*- """@package parmepy.domain.box
Classes for box domains description.
"""
from ..constants import * from ..constants import *
from .domain import Domain from .domain import Domain
from .grid import CartesianGrid from .grid import CartesianGrid
import numpy as np
class Box(Domain): class Box(Domain):
""" """
Periodic box representation. Periodic box representation.
Note FP : BC as parameter may be better? @note FP : BC as parameter may be better?
@todo Have different BC
""" """
def __init__(self, dimension=3, length=[1.0, 1.0, 1.0], origin=[0., 0., 0.]): def __init__(self, dimension=3, length=[1.0, 1.0, 1.0], origin=[0., 0., 0.]):
""" """
Constructor. Create a Periodic Box from a dimension, length and origin.
Create a Box from a dimension, length and minimum position.
Parameters dimensions must coincide. Raise ValueError in case of inconsistent parameters dimensions. Parameters dimensions must coincide. Raise ValueError in case of inconsistent parameters dimensions.
@param length : Box length. @par
@param origin : Box minimum position. By defaults, it creates a \f$[0;1]^3\f$ Box.
@param dimension : Box dimension. Default: 3
@param length : Box length. Default [1.0, 1.0, 1.0]
@param origin : Box minimum position. Default [0., 0., 0.]
""" """
if not (dimension == len(length) and dimension == len(origin)): if not (dimension == len(length) and dimension == len(origin)):
raise ValueError("Box parameters inconsistents dimensions") raise ValueError("Box parameters inconsistents dimensions")
Domain.__init__(self, dimension) Domain.__init__(self, dimension)
## Box length. length = max - min. ## Box length.
self.length = np.asarray(length, dtype=PARMES_REAL) self.length = np.asarray(length, dtype=PARMES_REAL)
## Minimum Box position. ## Box origin
self.origin = np.asarray(origin, dtype=PARMES_REAL) self.origin = np.asarray(origin, dtype=PARMES_REAL)
## Maximum Box position. ## Maximum Box position. max = origin + length
self.max = self.origin + self.length self.max = self.origin + self.length
# Boundary conditions type : ## Boundary conditions type
self.boundaries = np.zeros((self.dimension), dtype=PARMES_INTEGER) self.boundaries = np.zeros((self.dimension), dtype=PARMES_INTEGER)
self.boundaries[:] = PERIODIC self.boundaries[:] = PERIODIC
self.discreteDomain = None
def discretize(self, resolution): def discretize(self, resolution):
"""Box discretization method. """
Box discretization method.
Creates a CartesianGrid from the Box and a resolution.
@param resolution : resolution of the discretization.
@return CartesianGrid as discretized Box.
Use an array for discreteDomain??? i.e to associate several discretisations for the same domain. Use an array for discreteDomain??? i.e to associate several discretisations for the same domain.
Or maybe this should be a member of the fields?? Or maybe this should be a member of the fields??
""" """
self.discreteDomain = CartesianGrid(resolution, self) self.discreteDomain = CartesianGrid(resolution, self)
return self.discreteDomain return self.discreteDomain
def __str__(self): def __str__(self):
""" Doc display. """ """
Informations display.
@return Informations
"""
s = str(self.dimension) + "D box (parallelepipedic or rectangular) domain : \n" s = str(self.dimension) + "D box (parallelepipedic or rectangular) domain : \n"
s += " origin : " + str(self.origin) + ", maxPosition :" + str(self.max) + ", lengthes :" + str(self.length) + "." s += " origin : " + str(self.origin) + ", maxPosition :" + str(self.max) + ", lengthes :" + str(self.length) + "."
if self.discreteDomain is not None: if self.discreteDomain is not None:
......
# -*- coding: utf-8 -*- """@package parmepy.domain.discrete
"""
@package parmepy.domain.discrete
Physical domain discretization
Classes for physical domains discretization description.
""" """
from abc import ABCMeta, abstractmethod from abc import ABCMeta, abstractmethod
class DiscreteDomain: class DiscreteDomain:
""" Abstract description of a discretized physical domain. """ """ Abstract description of a discretized physical domain."""
__metaclass__ = ABCMeta __metaclass__ = ABCMeta
@abstractmethod @abstractmethod
def __init__(self, domain=None): def __init__(self, domain=None):
""" """
Constructor.
Create a DiscreteDomain with a given discretization. Create a DiscreteDomain with a given discretization.
@param domain : Domain to discretize
""" """
## Domain to discretize
self.domain = domain self.domain = domain
if(domain is not None): if(domain is not None):
## Domain dimension
self.dimension = domain.dimension self.dimension = domain.dimension
## Discretization specifications. Discrete elements number.
#self.elementNumber = None
## Discretization specifications. Discrete elements size.
#self.elementSize = None
if __name__ == "__main__": if __name__ == "__main__":
print __doc__ print __doc__
......
# -*- coding: utf-8 -*-
"""@package parmepy.domain.domain """@package parmepy.domain.domain
Classes for physical domains description. Classes for physical domains description.
...@@ -13,10 +12,8 @@ class Domain: ...@@ -13,10 +12,8 @@ class Domain:
@abstractmethod @abstractmethod
def __init__(self, dimension): def __init__(self, dimension):
""" """ Constructor
@param dimension integer : domain dimension. @param dimension integer : domain dimension.
""" """
## Domain dimension. ## Domain dimension.
self.dimension = dimension self.dimension = dimension
......
# -*- coding: utf-8 -*- """@package parmepy.domain.grid
"""
@package parmepy.domain.grid
Cartesian grid definition Cartesian grid definition.
""" """
from ..constants import *
from .discrete import DiscreteDomain from .discrete import DiscreteDomain
import numpy as np
class CartesianGrid(DiscreteDomain): class CartesianGrid(DiscreteDomain):
""" Discrete box representation as a cartesian grid. """ """ Discrete box representation as a cartesian grid."""
def __init__(self, resolution, domain): def __init__(self, resolution, domain):
"""
Creates a CartesianGrid.
@param resolution : resolution used
@param domain : Domain to discretize
"""
DiscreteDomain.__init__(self, domain) DiscreteDomain.__init__(self, domain)
assert(self.dimension == len(resolution)) assert(self.dimension == len(resolution))
## lowest point of the grid ## lowest point of the grid
...@@ -20,24 +25,11 @@ class CartesianGrid(DiscreteDomain): ...@@ -20,24 +25,11 @@ class CartesianGrid(DiscreteDomain):
## number of points in each direction ## number of points in each direction
self.resolution = resolution self.resolution = resolution
## space step size ## space step size
self.step = self.length/(self.resolution-1) self.step = self.length / (self.resolution - 1)
## 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): def update(start):
self.start = start self.start = start
def applyConditions(self, pos, dir):
"""
Apply periodic boundary conditions of the periodic domain.
@param pos position : position to check.
@return position inside the box boundary.
"""
## ------------------- NumPy Optimisation
pos[..., dir] = np.where(pos[..., dir] >= self.max[dir], pos[..., dir] - self.max[dir] + self.min[dir], pos[..., dir])
pos[..., dir] = np.where(pos[..., dir] < self.min[dir], pos[..., dir] + self.max[dir] - self.min[dir], pos[..., dir])
## -------------------
def __getitem__(self, index): def __getitem__(self, index):
""" """
Get iter overriding Get iter overriding
......
...@@ -2,7 +2,6 @@ ...@@ -2,7 +2,6 @@
@package parmepy.domain.topology @package parmepy.domain.topology
MPI Topologies MPI Topologies
""" """
import numpy as np
import mpi4py.MPI as MPI import mpi4py.MPI as MPI
from ..constants import * from ..constants import *
...@@ -65,7 +64,12 @@ class CartesianTopology(object): ...@@ -65,7 +64,12 @@ class CartesianTopology(object):
localResolution[i] = nbpoints[i] + remainingPoints[i] localResolution[i] = nbpoints[i] + remainingPoints[i]
start = np.zeros((domain.dimension), dtype=PARMES_INTEGER) start = np.zeros((domain.dimension), dtype=PARMES_INTEGER)
start[:self.dim] = self.coords * nbpoints start[:self.dim] = self.coords * nbpoints
self.mesh = LocalMesh(self.rank, resolution=localResolution, start=start, dom_size=self.domain.length / self.resolution,dom_origin=self.domain.origin, ghosts=self.ghosts) self.mesh = LocalMesh(self.rank,
resolution=localResolution,
start=start,
dom_size=self.domain.length / self.resolution,
dom_origin=self.domain.origin,
ghosts=self.ghosts)
def __str__(self): def __str__(self):
""" Topology info display """ """ Topology info display """
......
...@@ -2,5 +2,4 @@ ...@@ -2,5 +2,4 @@
@package parmepy.fields @package parmepy.fields
Everything concerning Fields. Everything concerning Fields.
""" """
# -*- coding: utf-8 -*-
""" """
@package parmepy.fields.analytical @package parmepy.fields.analytical
Continuous variable description defined with an analytic formula.
""" """
from .continuous import ContinuousField from .continuous import ContinuousField
from .discrete import ScalarField from .discrete import ScalarField
...@@ -16,12 +16,14 @@ class AnalyticalField(ContinuousField): ...@@ -16,12 +16,14 @@ class AnalyticalField(ContinuousField):
def __init__(self, domain, formula, name="", vector=False): def __init__(self, domain, formula, name="", vector=False):
""" """
Constructor.
Create an AnalyticalField from a formula. Create an AnalyticalField from a formula.
@param domain : application domain of the variable. @param domain : variable application domain.
@param dimension : variable dimension. @param formula : python function
@param formula : function defining the variable. @param name : name of the variable (used for vtk output).
@param vector : is variable is a vector.
@note formula is used in ScalarField or VectorField as vectorized function by numpy.
""" """
ContinuousField.__init__(self, domain, name, vector) ContinuousField.__init__(self, domain, name, vector)
## Analytic formula. ## Analytic formula.
...@@ -32,7 +34,7 @@ class AnalyticalField(ContinuousField): ...@@ -32,7 +34,7 @@ class AnalyticalField(ContinuousField):
Evaluation of the variable at a given position. Evaluation of the variable at a given position.
@param pos : Position to evaluate. @param pos : Position to evaluate.
@return : value of the formula at the given position. @return formula(pos): value of the formula at the given position.
""" """
return self.formula(*pos) return self.formula(*pos)
......
# -*- coding: utf-8 -*-
""" """
@package parmepy.fields.continuous @package parmepy.fields.continuous
Continuous variable description Continuous variable description
""" """
from .discrete import ScalarField from .discrete import ScalarField
...@@ -12,11 +12,11 @@ class ContinuousField(object): ...@@ -12,11 +12,11 @@ class ContinuousField(object):
def __init__(self, domain, name="?", vector=False): def __init__(self, domain, name="?", vector=False):
""" """
Constructor.
Create a ContinuousField. Create a ContinuousField.
@param domain : variable application domain. @param domain : variable application domain.
@param name : name of the variable (used for vtk output). @param name : name of the variable (used for vtk output).
@param vector : is variable is a vector.
""" """
## Application domain of the variable. ## Application domain of the variable.
self.domain = domain self.domain = domain
...@@ -26,17 +26,18 @@ class ContinuousField(object): ...@@ -26,17 +26,18 @@ class ContinuousField(object):
self.discreteField = None self.discreteField = None
## Name of this field ## Name of this field
self.name = name self.name = name
# number of different discretizations # Number of different discretizations
self._fieldId = -1 self._fieldId = -1
## list of the various discretizations of this field ## List of the various discretizations of this field
self.discreteField = [] self.discreteField = []
## is field is a vector field ## Is field is a vector field
self.vector = vector self.vector = vector
def discretize(self, topology=None): def discretize(self, topology=None):
""" """
discretization of the field on a topology Discretization of the field on a topology
@param topology : Topology definition to use
@return discrete field and its index
""" """
self._fieldId += 1 self._fieldId += 1
if self.vector: if self.vector:
...@@ -50,12 +51,17 @@ class ContinuousField(object): ...@@ -50,12 +51,17 @@ class ContinuousField(object):
return self.discreteField[self._fieldId], self._fieldId return self.discreteField[self._fieldId], self._fieldId
def initialize(self): def initialize(self):
"""
Initialize a field.
@see VectorField.initialize
@see ScalarField.initialize
"""
if self._fieldId == -1: if self._fieldId == -1:
raise ValueError("Cannot initialise analytical field non discretized.") raise ValueError("Cannot initialise analytical field non discretized.")
print self.name, print self.name,
for dF in self.discreteField: for dF in self.discreteField:
dF.initialize() dF.initialize()
print "Done" print " .Done"
def __str__(self): def __str__(self):
""" Field info display """ """ Field info display """
......
# -*- coding: utf-8 -*-
""" """
@package parmepy.fields.discrete @package parmepy.fields.discrete
Discretes variables (scalar and vectors) descriptions.
""" """
import numpy as np import numpy as np
from ..constants import * from ..constants import *
...@@ -15,25 +15,42 @@ class ScalarField(object): ...@@ -15,25 +15,42 @@ class ScalarField(object):
""" """
def __init__(self, parent, topology=None, name="?", idFromParent=None): def __init__(self, parent, topology=None, name="?", idFromParent=None):
"""
Creates a ScalaField.
@param parent : Continuous field.
@param topology : Topology informations
@param name : Field name
@param idFromParent : Index in the parent's discrete fields
"""
if(topology is not None): if(topology is not None):
## Topology informations.
self.topology = topology self.topology = topology
else: else:
raise NotImplementedError() raise NotImplementedError()
## Field name.
self.name = name self.name = name
## Field dimension.
self.dimension = topology.domain.dimension self.dimension = topology.domain.dimension
## Field resolution.
self.resolution = self.topology.mesh.resolution self.resolution = self.topology.mesh.resolution
## \todo is numpy array zeros or empty??? ## Field data numpy array. \todo is numpy array zeros or empty???
self.data = np.zeros((self.resolution), dtype=PARMES_REAL, order=ORDER) self.data = np.zeros((self.resolution), dtype=PARMES_REAL, order=ORDER)
## Pointer to gpu data array if needed.
self.gpu_data = None self.gpu_data = None
## @private Continuous field
self.__parentField = parent self.__parentField = parent
## @private Index in parent's field discrete fields
self.__idFromParent = idFromParent self.__idFromParent = idFromParent
## Application domain of the variable.
self.domain = self.__parentField.domain.discreteDomain self.domain = self.__parentField.domain.discreteDomain
## Is self.data contains data
self.contains_data = False self.contains_data = False
## Scalar is not vector
self.vector = False self.vector = False
def __str__(self): def __str__(self):
""" Display class information """ """ Display class information. """
s = '[' + str(self.topology.rank) + '] ' s = '[' + str(self.topology.rank) + '] '
s += " id " + str(self.__idFromParent) + ", " s += " id " + str(self.__idFromParent) + ", "
s += self.name + " " + str(self.dimension) + 'D discrete field with resolution ' s += self.name + " " + str(self.dimension) + 'D discrete field with resolution '
...@@ -41,21 +58,36 @@ class ScalarField(object): ...@@ -41,21 +58,36 @@ class ScalarField(object):
return s + "\n" return s + "\n"
def __getitem__(self, i): def __getitem__(self, i):
""" Access to the content of the field """ """
Access to the content of the field.
@param i : requested index.
@return data[i].
"""
return self.data.__getitem__(i) return self.data.__getitem__(i)
def __setitem__(self, i, value): def __setitem__(self, i, value):
""" """
Set the content of the field component at position i Set the content of the field component at position i.
Usage :
A = ScalarField(topo) @param i : requested index.
@param value : value to set.
Usage :\n
A = ScalarField(topo)\n
A[2,1,1] = 12.0 # Calls A.data[2,1,1] = 12.0 A[2,1,1] = 12.0 # Calls A.data[2,1,1] = 12.0
""" """
self.data.__setitem__(i, value) self.data.__setitem__(i, value)
def initialize(self, formula=None): def initialize(self, formula=None):
""" """
Initialize values with given formula Initialize values with given formula.
@param formula : formula to apply on coordinates.
Formula need to be vectorized by numpy. Therefore formula must have the following signature:
@li 3D : float formula(float x, float y, float z)
@li 2D : float formula(float x, float y)
where x,y,z are point coordinates.
""" """
if formula is not None: if formula is not None:
print "...", print "...",
...@@ -73,9 +105,6 @@ class ScalarField(object): ...@@ -73,9 +105,6 @@ class ScalarField(object):
else: else:
print "No formula", print "No formula",
def tofile(self, filename):
evtk.imageToVTK(filename, pointData={self.name: self.data})
class VectorField(object): class VectorField(object):
""" """
...@@ -84,21 +113,38 @@ class VectorField(object): ...@@ -84,21 +113,38 @@ class VectorField(object):
""" """
def __init__(self, parent, topology=None, name="?", idFromParent=None): def __init__(self, parent, topology=None, name="?", idFromParent=None):
"""
Creates a VectorField.
@param parent : Continuous field.
@param topology : Topology informations
@param name : Field name
@param idFromParent : Index in the parent's discrete fields
"""
if(topology is not None): if(topology is not None):
## Topology informations.
self.topology = topology self.topology = topology
else: else:
raise NotImplementedError() raise NotImplementedError()
## Field name.
self.name = name self.name = name
## Field dimension.
self.dimension = topology.domain.dimension self.dimension = topology.domain.dimension
## Field resolution.
self.resolution = self.topology.mesh.resolution self.resolution = self.topology.mesh.resolution
## \todo is numpy array zeros or empty??? ## Field data numpy array. \todo is numpy array zeros or empty???
self.data = [np.zeros((self.resolution), dtype=PARMES_REAL, order=ORDER) for d in xrange(self.dimension)] self.data = [np.zeros((self.resolution), dtype=PARMES_REAL, order=ORDER) for d in xrange(self.dimension)]
## Pointer to gpu data arrays if needed.
self.gpu_data = [None for d in xrange(self.dimension)] self.gpu_data = [None for d in xrange(self.dimension)]
## @private Continuous field
self.__parentField = parent self.__parentField = parent
## @private Index in parent's field discrete fields
self.__idFromParent = idFromParent self.__idFromParent = idFromParent
## Application domain of the variable.
self.domain = self.__parentField.domain.discreteDomain self.domain = self.__parentField.domain.discreteDomain
## Is self.data contains data
self.contains_data = False self.contains_data = False
## Scalar is not vector
self.vector = True self.vector = True
def __str__(self): def __str__(self):
...@@ -110,12 +156,17 @@ class VectorField(object): ...@@ -110,12 +156,17 @@ class VectorField(object):
return s + "\n" return s + "\n"
def __getitem__(self, i): def __getitem__(self, i):
""" Access to the content of the field """ Access to the content of the field.
Usage (3D): A = VectorField(...), We have 3 components len(a.data) == 3.
@param i : requested index.
@return data[i] depending on i.
Usage (3D): \n
A = VectorField(...), We have 3 components len(a.data) == 3. \n
Following instructions access to index 2,1,1 of y component: Following instructions access to index 2,1,1 of y component:
A[2,1,1,1] @li A[2,1,1,1]
A[1][2,1,1] @li A[1][2,1,1]
Access to whole vector of index 2,1,1: A[2,1,1] @li Access to whole vector of index 2,1,1: A[2,1,1]
""" """
try: try:
if len(i) == len(self.data): if len(i) == len(self.data):
...@@ -127,11 +178,15 @@ class VectorField(object): ...@@ -127,11 +178,15 @@ class VectorField(object):
def __setitem__(self, i, value): def __setitem__(self, i, value):
""" """
Set the content of the vector field component at position i Set the content of the vector field component at position i.
Usage :
A = VectorField(topo) @param i : requested index.
A[2,1,1] = 12.0 # Calls A.data[d][2,1,1] = 12.0 for all components d. @param value : value to set.
A[2,1,1] = [12.0, 13.0, 14.0] # Calls A.data[0][2,1,1] = 12.0, A.data[1][2,1,1] = 13.0 and A.data[2][2,1,1] = 14.0
Usage :\n
A = VectorField(topo)\n
A[2,1,1] = 12.0 # Calls A.data[d][2,1,1] = 12.0 for all components d.\n
A[2,1,1] = [12.0, 13.0, 14.0] # Calls A.data[0][2,1,1] = 12.0, A.data[1][2,1,1] = 13.0 and A.data[2][2,1,1] = 14.0\n
A[1][2,1,1] = 13.0 # Calls A.data[1][2,1,1] = 12.0 A[1][2,1,1] = 13.0 # Calls A.data[1][2,1,1] = 12.0
""" """
if len(i) == len(self.data): if len(i) == len(self.data):
...@@ -144,9 +199,15 @@ class VectorField(object): ...@@ -144,9 +199,15 @@ class VectorField(object):
def initialize(self, formula=None): def initialize(self, formula=None):
""" """
Initialize values with given formula Initialize values with given formula.
@param formula : formula to apply on coordinates.
Formula need to be vectorized by numpy. Therefore formula must have the following signature:
@li 3D : float, float, float formula(float x, float y, float z)
@li 2D : float, float formula(float x, float y)
where x,y,z are point coordinates.
""" """
#print self.topology.mesh.coords
if formula is not None: if formula is not None:
print "...", print "...",
v_formula = np.vectorize(formula) v_formula = np.vectorize(formula)
...@@ -163,5 +224,5 @@ class VectorField(object): ...@@ -163,5 +224,5 @@ class VectorField(object):
if __name__ == "__main__": if __name__ == "__main__":
print __doc__ print __doc__
print "- Provided class : ScalarField" print "- Provided class : ScalarField, VectorField"
print ScalarField.__doc__ print ScalarField.__doc__
...@@ -2,5 +2,4 @@ ...@@ -2,5 +2,4 @@
@package parmepy.operator @package parmepy.operator
Everything concerning operators. Everything concerning operators.
""" """
# -*- coding: utf-8 -*-
""" """
@package parmepy.operator.continuous @package parmepy.operator.continuous
Operator representation Operator representation.
""" """
from abc import ABCMeta, abstractmethod from abc import ABCMeta, abstractmethod
...@@ -17,47 +16,58 @@ class ContinuousOperator: ...@@ -17,47 +16,58 @@ class ContinuousOperator:
@abstractmethod @abstractmethod
def __init__(self): def __init__(self):
""" """
Constructor.
Create a ContinuousOperator. Create a ContinuousOperator.
""" """
## Variables ## Variables
self.variables = [] self.variables = []
## Operator discretization. ## Operator discretization.
self.discreteOperator = None self.discreteOperator = None
## Is need to split operator
self.needSplitting = False
def addVariable(self, cVariable): def addVariable(self, cVariable):
""" """
Add an continuous variables to the operator. Add an continuous variables to the operator.
Also add variables' domains to the operator.
@param cVariable ContinuousVariable.ContinuousVariable : variables to add. @param cVariable : list of variables to add.
""" """
for v in cVariable: for v in cVariable:
if self.variables.count(v) == 0: if self.variables.count(v) == 0:
self.variables.append(v) self.variables.append(v)
def apply(self, *args): def apply(self, *args):
"""
Apply operator.
@param *args : Operator arguments.
@return computing time.
"""
return self.discreteOperator.apply(*args) return self.discreteOperator.apply(*args)
def setMethod(self, method): def setMethod(self, method):
"""
Sets method to use in apply method.
@param method : the method to use.
"""
if self.discreteOperator is not None: if self.discreteOperator is not None:
self.discreteOperator.setMethod(method) self.discreteOperator.setMethod(method)
else: else:
raise ValueError("Cannot set mumerical method to non discretized operator") raise ValueError("Cannot set mumerical method to non discretized operator")
def printComputeTime(self): def printComputeTime(self):
"""Displays compute time for operator."""
if self.discreteOperator is not None: if self.discreteOperator is not None:
self.discreteOperator.printComputeTime() self.discreteOperator.printComputeTime()
else: else:
raise ValueError("Cannot print compute time of a non discretized operator") raise ValueError("Cannot print compute time of a non discretized operator")
@abstractmethod @abstractmethod
def discretize(self, spec=None): def discretize(self, *spec):
""" """
Abstract method. Abstract method.
Must be implemented by sub-class. Must be implemented by sub-class.
@param spec : discretization specifications. @param *spec : discretization specifications.
""" """
raise NotImplementedError("Need to override method in a subclass of " + providedClass) raise NotImplementedError("Need to override method in a subclass of " + providedClass)
......
# -*- coding: utf-8 -*-
""" """
@package operator @package parmepy.operator.discrete
Discrete operator representation
Discrete operator representation.
""" """
from abc import ABCMeta, abstractmethod from abc import ABCMeta, abstractmethod
...@@ -16,19 +16,22 @@ class DiscreteOperator: ...@@ -16,19 +16,22 @@ class DiscreteOperator:
@abstractmethod @abstractmethod
def __init__(self): def __init__(self):
""" """
Constructor.
Create an empty discrete operator. Create an empty discrete operator.
""" """
self.input, self.output = None, None ## Input variables
self.input = None
## Output variables
self.output = None
## variables ## variables
self.variables = [] self.variables = []
## Operator numerical method. ## Operator numerical method.
self.numMethod = [] self.numMethod = None
self.needSplitting = False ## DiscreteOperator is a discrete operator
self.gpu_kernel = None
self.gpu_kernel_name = ""
self.discreteOperator = self self.discreteOperator = self
## Total compute time
self.total_time = 0. self.total_time = 0.
## Operator name
self.name = "?"
def setMethod(self, method): def setMethod(self, method):
""" """
...@@ -41,9 +44,8 @@ class DiscreteOperator: ...@@ -41,9 +44,8 @@ class DiscreteOperator:
def addVariable(self, cVariable): def addVariable(self, cVariable):
""" """
Add an continuous variables to the operator. Add an continuous variables to the operator.
Also add variables' domains to the operator.
@param cVariable ContinuousVariable.ContinuousVariable : variables to add. @param cVariable : list of variables to add.
""" """
for v in cVariable: for v in cVariable:
if self.variables.count(v) == 0: if self.variables.count(v) == 0:
...@@ -51,9 +53,7 @@ class DiscreteOperator: ...@@ -51,9 +53,7 @@ class DiscreteOperator:
@abstractmethod @abstractmethod
def printComputeTime(self): def printComputeTime(self):
""" """Print total computing time."""
Print total computing time
"""
raise NotImplementedError("Need to override method in a subclass of " + providedClass) raise NotImplementedError("Need to override method in a subclass of " + providedClass)
@abstractmethod @abstractmethod
...@@ -61,8 +61,6 @@ class DiscreteOperator: ...@@ -61,8 +61,6 @@ class DiscreteOperator:
""" """
Abstract method, apply operaton on a variable. Abstract method, apply operaton on a variable.
Must be implemented by sub-class. Must be implemented by sub-class.
@param point DiscreteVariable : apply operator on this point.
""" """
raise NotImplementedError("Need to override method in a subclass of " + providedClass) raise NotImplementedError("Need to override method in a subclass of " + providedClass)
......
# -*- coding: utf-8 -*-
""" """
@package Operator @package parmepy.operator.remeshing
Operator representation
Discrete remeshing representation
""" """
from .discrete import DiscreteOperator from .discrete import DiscreteOperator
from ..constants import * from ..constants import *
...@@ -12,20 +12,16 @@ import time ...@@ -12,20 +12,16 @@ import time
class Remeshing(DiscreteOperator): class Remeshing(DiscreteOperator):
""" """
Remeshing operator representation. Remeshing operator representation.
DiscreteOperator.DiscreteOperator specialization.
""" """
def __init__(self, partPositions, partScalar, resscal, method): def __init__(self, partPositions, partScalar, resscal, method=None):
""" """
Constructor. Create a Remeshing operator.
Create a Remeshing operator on a given discrete domain.
Work with a particle field (positions, scalar, type and tag) to set scalar on a grid.
@param partPositions : particles positions. @param partPositions : particles positions.
@param partScalar : particles scalar. @param partScalar : particles scalar.
@param partBloc : particles bloc number. @param resscal : result grid scalar values.
@param partTag : particle tag number. @param method : the method to use.
@param resscal DiscreteVariable.DiscreteVariable : new grid scalar values.
""" """
DiscreteOperator.__init__(self) DiscreteOperator.__init__(self)
## Particles positions. ## Particles positions.
...@@ -35,18 +31,25 @@ class Remeshing(DiscreteOperator): ...@@ -35,18 +31,25 @@ class Remeshing(DiscreteOperator):
## Result scalar ## Result scalar
self.res_scalar = resscal self.res_scalar = resscal
self.addVariable([self.ppos, self.pscal, self.res_scalar]) self.addVariable([self.ppos, self.pscal, self.res_scalar])
self.numMethod = method
## input fields
self.input = [self.ppos, self.pscal] self.input = [self.ppos, self.pscal]
## output fields
self.output = [self.res_scalar] self.output = [self.res_scalar]
self.needSplitting = True self.method = method
## Compute time detailed per directions
self.compute_time = [0., 0., 0.] self.compute_time = [0., 0., 0.]
self.name = "remeshing" self.name = "remeshing"
def apply(self, t, dt, splittingDirection): def apply(self, t, dt, splittingDirection):
""" """
Apply Remeshing operator. Apply Remeshing operator.
@param t : current time.
@param dt : time step.
@param splittingDirection : Direction of splitting.
Remeshing algorithm:
@li 1. Remesh particles on grid
- Use a M'6 formula
@li 2. Profile timings of OpenCL kernels.
""" """
c_time = 0. c_time = 0.
if self.numMethod is not None: if self.numMethod is not None:
...@@ -58,6 +61,7 @@ class Remeshing(DiscreteOperator): ...@@ -58,6 +61,7 @@ class Remeshing(DiscreteOperator):
PARMES_REAL_GPU(self.res_scalar.domain.step[splittingDirection])) PARMES_REAL_GPU(self.res_scalar.domain.step[splittingDirection]))
for df in self.output: for df in self.output:
df.contains_data = False df.contains_data = False
# Get timpings from OpenCL events
self.numMethod.queue.finish() self.numMethod.queue.finish()
c_time += (evt.profile.end - evt.profile.start) * 1e-9 c_time += (evt.profile.end - evt.profile.start) * 1e-9
self.compute_time[splittingDirection] += c_time self.compute_time[splittingDirection] += c_time
......
# -*- coding: utf-8 -*-
""" """
@package Operator @package parmepy.operator.splitting
Operator representation
Splitting operator representation
""" """
from .discrete import DiscreteOperator from .discrete import DiscreteOperator
import time import time
...@@ -10,18 +10,33 @@ import time ...@@ -10,18 +10,33 @@ import time
class Splitting(DiscreteOperator): class Splitting(DiscreteOperator):
""" """
Splitting operator representation. Splitting operator representation.
DiscreteOperator.DiscreteOperator specialization.
Operator of operators applying in a Strang splitting.
Implements a 2nd order splitting in 3D:
@li X-dir, half time step
@li Y-dir, half time step
@li Z-dir, full time step
@li Y-dir, half time step
@li X-dir, half time step
\n
Implements a 2nd order splitting in 2D:
@li X-dir, half time step
@li Y-dir, full time step
@li X-dir, half time step
""" """
def __init__(self, operators=[], dim=3): def __init__(self, operators=[], dim=3):
""" """
Constructor. Create a Splitting operator on a given list of operators and a dimension.
Create a Splittinf operator on a given discrete domain.
Work with a particle field (positions, scalar, type and tag) to set scalar on a grid.
@param operators : list of operators to split.
@param dim : problem dimension.
""" """
DiscreteOperator.__init__(self) DiscreteOperator.__init__(self)
## Operators to split
self.operators = operators self.operators = operators
## Splitting at 2nd order.
self.splitting = [] self.splitting = []
## Half timestep in all directions ## Half timestep in all directions
# [self.splitting.append((i, 0.5)) for i in xrange(dim)] # [self.splitting.append((i, 0.5)) for i in xrange(dim)]
...@@ -30,12 +45,16 @@ class Splitting(DiscreteOperator): ...@@ -30,12 +45,16 @@ class Splitting(DiscreteOperator):
[self.splitting.append((i, 0.5)) for i in xrange(dim - 1)] [self.splitting.append((i, 0.5)) for i in xrange(dim - 1)]
self.splitting.append((dim - 1, 1.)) self.splitting.append((dim - 1, 1.))
[self.splitting.append((dim - 2 - i, 0.5)) for i in xrange(dim - 1)] [self.splitting.append((dim - 2 - i, 0.5)) for i in xrange(dim - 1)]
## Compute time detailed per directions and per operators
self.compute_time_details = [[0. for i in xrange(dim)] for op in self.operators] self.compute_time_details = [[0. for i in xrange(dim)] for op in self.operators]
self.name = "splitting" self.name = "splitting"
def apply(self, t, dt): def apply(self, t, dt):
""" """
Apply Remeshing operator. Apply Remeshing operator.
@param t : current time.
@param dt : time step.
""" """
c_time = 0. c_time = 0.
for split in self.splitting: for split in self.splitting:
...@@ -66,5 +85,5 @@ class Splitting(DiscreteOperator): ...@@ -66,5 +85,5 @@ class Splitting(DiscreteOperator):
if __name__ == "__main__": if __name__ == "__main__":
print __doc__ print __doc__
print "- Provided class : RemeshingDOp" print "- Provided class : Splitting"
print RemeshingDOp.__doc__ print RemeshingDOp.__doc__
# -*- coding: utf-8 -*-
""" """
@package parmepy.operator.Advection @package parmepy.operator.transport
Advection operator representation Transport operator representation.
""" """
from .continuous import ContinuousOperator from .continuous import ContinuousOperator
from .advection_d import Advection_P from .transport_d import Transport_d
class Transport(ContinuousOperator): class Transport(ContinuousOperator):
""" """
Advection operator representation Transport operator representation.
""" """
def __init__(self, velocity, scalar): def __init__(self, velocity, scalar):
""" """
Constructor. Create an Transport operator from given variables velocity and scalar.
Create an Advection operator from given variables velocity and scalar.
@param velocity ContinuousField : velocity variable. @param velocity : velocity variable.
@param scalar ContinuousField : scalar variable. @param scalar : scalar variable.
""" """
ContinuousOperator.__init__(self) ContinuousOperator.__init__(self)
## advection velocity variable ## Transport velocity
self.velocity = velocity self.velocity = velocity
## advected scalar variable ## Transported scalar
self.scalar = scalar self.scalar = scalar
self.needSplitting = True
self.addVariable([velocity, scalar]) self.addVariable([velocity, scalar])
def discretize(self, idVelocityD=0, idScalarD=0, result_position=None, result_scalar=None, method=None): def discretize(self, idVelocityD=0, idScalarD=0, result_position=None, result_scalar=None, method=None):
""" """
Advection operator discretization method. Transport operator discretization method.
Create an AdvectionDOp.AdvectionDOp from given specifications. Create an discrete Transport operator from given specifications.
@param spec : discretization specifications, not used. @param idVelocityD : Index of velocity discretisation to use.
@param idScalarD : Index of scalar discretisation to use.
@param result_position : result position.
@param result_scalar : result scalar.
@param method : the method to use.
""" """
if self.discreteOperator is None: if self.discreteOperator is None:
self.discreteOperator = Advection_P(self, idVelocityD, idScalarD, result_position, result_scalar) self.discreteOperator = Transport_d(self, idVelocityD, idScalarD, result_position, result_scalar, method)
def __str__(self): def __str__(self):
"""ToString method""" """ToString method"""
......
# -*- coding: utf-8 -*-
""" """
@package operator @package parmepy.operator.transport_d
Discrete advection representation
Discrete transport representation
""" """
from ..constants import * from ..constants import *
from .discrete import DiscreteOperator from .discrete import DiscreteOperator
import pyopencl as cl import pyopencl as cl
import numpy as np
import time import time
class Advection_P(DiscreteOperator): class Transport_d(DiscreteOperator):
""" """
Particle avection operator representation. Particle transport operator representation.
DiscreteOperator.DiscreteOperator specialization.
""" """
def __init__(self, advec, idVelocityD=0, idScalarD=0, result_position=None, result_scalar=None, method=None): def __init__(self, advec, idVelocityD=0, idScalarD=0, result_position=None, result_scalar=None, method=None):
""" """
Constructor. Create a Advection operator.
Create a Advection operator on a given continuous domain.
Work on a given scalar at a given velocity to produce scalar distribution at new positions. Work on a given scalar at a given velocity to produce scalar distribution at new positions.
@param advec : Advection operator. @param advec : Transport operator
@param idVelocityD : Index of velocity discretisation to use.
@param idScalarD : Index of scalar discretisation to use.
@param result_position : result position.
@param result_scalar : result scalar.
@param method : the method to use.
""" """
DiscreteOperator.__init__(self) DiscreteOperator.__init__(self)
## Velocity. ## Velocity.
...@@ -33,21 +36,37 @@ class Advection_P(DiscreteOperator): ...@@ -33,21 +36,37 @@ class Advection_P(DiscreteOperator):
self.res_position = result_position.discreteField[0] self.res_position = result_position.discreteField[0]
## Result scalar ## Result scalar
self.res_scalar = result_scalar.discreteField[0] self.res_scalar = result_scalar.discreteField[0]
## input fields
self.input = [self.velocity, self.scalar] self.input = [self.velocity, self.scalar]
## output fields
self.output = [self.res_position, self.res_scalar] self.output = [self.res_position, self.res_scalar]
self.numMethod = method self.method = method
self.needSplitting = True ## Previous splitting direction
self.old_splitting_direction = None self.old_splitting_direction = None
## Compute time detailed per directions
self.compute_time = [0., 0., 0.] self.compute_time = [0., 0., 0.]
## Compute time for copy detailed per directions
self.compute_time_copy = [0., 0., 0.] self.compute_time_copy = [0., 0., 0.]
## Compute time for transposition detailed per directions
self.compute_time_transpose = [0., 0., 0.] self.compute_time_transpose = [0., 0., 0.]
self.name = "advection" self.name = "advection"
def apply(self, t, dt, splittingDirection): def apply(self, t, dt, splittingDirection):
"""
Apply advection operator.
@param t : current time.
@param dt : time step.
@param splittingDirection : Direction of splitting.
Advection algorithm:
@li 1. Particle initialization : \n
- by copy scalar from grid to particles if previous splitting direction equals current splitting direction.\n
- by transposition of scalar from grid to particle.
@li 2. Particle advection :\n
- compute particle position in splitting direction as a scalar. Performs a RK2 resolution of dx_p/dt = a_p.
@li 3. Profile timings of OpenCL kernels.
"""
c_time, c_time_init = 0., 0. c_time, c_time_init = 0., 0.
if self.numMethod is not None: if self.numMethod is not None and self.init_transpose is not None and self.init_copy is not None:
# Particle init # Particle init
if (self.old_splitting_direction == splittingDirection) or self.old_splitting_direction is None: if (self.old_splitting_direction == splittingDirection) or self.old_splitting_direction is None:
evt_init = self.init_copy.launch(self.scalar.gpu_data, evt_init = self.init_copy.launch(self.scalar.gpu_data,
...@@ -73,6 +92,7 @@ class Advection_P(DiscreteOperator): ...@@ -73,6 +92,7 @@ class Advection_P(DiscreteOperator):
PARMES_REAL_GPU(self.scalar.domain.step[splittingDirection])) PARMES_REAL_GPU(self.scalar.domain.step[splittingDirection]))
for df in self.output: for df in self.output:
df.contains_data = False df.contains_data = False
# Get timpings from OpenCL events
self.numMethod.queue.finish() self.numMethod.queue.finish()
c_time_init = (evt_init.profile.end - evt_init.profile.start) * 1e-9 c_time_init = (evt_init.profile.end - evt_init.profile.start) * 1e-9
c_time = (evt.profile.end - evt.profile.start) * 1e-9 c_time = (evt.profile.end - evt.profile.start) * 1e-9
...@@ -97,5 +117,5 @@ class Advection_P(DiscreteOperator): ...@@ -97,5 +117,5 @@ class Advection_P(DiscreteOperator):
if __name__ == "__main__": if __name__ == "__main__":
print __doc__ print __doc__
print "- Provided class : AdvectionDOp" print "- Provided class : Transport_d"
print AdvectionDOp.__doc__ print AdvectionDOp.__doc__
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment