From bb745d01f721d1df8b138fb884191b8d0cbc846b Mon Sep 17 00:00:00 2001
From: Jean-Matthieu Etancelin <jean-matthieu.etancelin@imag.fr>
Date: Tue, 10 Jul 2012 16:06:21 +0000
Subject: [PATCH] parmepy get back pure python solver workflow. (Need pure
python code for advection and remeshing. Need GPU solver.)
---
Examples/mainJM.py | 70 +++------
HySoP/hysop/ParticularSolvers/basic.py | 145 ++++++++---------
.../interpolation/interpolation.py | 2 -
.../hysop/ParticularSolvers/remesh/method.py | 1 -
HySoP/hysop/__init__.py.in | 30 +++-
HySoP/hysop/constants.py | 4 +-
HySoP/hysop/domain/box.py | 8 +-
HySoP/hysop/domain/grid.py | 1 +
HySoP/hysop/fields/analytical.py | 28 ++--
HySoP/hysop/fields/continuous.py | 17 +-
HySoP/hysop/fields/discrete.py | 50 +++---
HySoP/hysop/fields/topology.py | 31 ++--
HySoP/hysop/operator/RemeshingDOp.py | 113 --------------
HySoP/hysop/operator/Transport.py | 2 +-
HySoP/hysop/operator/advection_d.py | 29 ++--
HySoP/hysop/operator/continuous.py | 3 +
HySoP/hysop/operator/discrete.py | 35 ++++-
HySoP/hysop/operator/remeshing.py | 105 +++++++++++++
HySoP/hysop/operator/splitting.py | 55 +++++++
HySoP/hysop/problem/problem.py | 36 +++--
HySoP/hysop/tools/__init__.py | 0
HySoP/hysop/tools/printer.py | 147 +++++++++---------
HySoP/hysop/tools/timer.py | 43 +++++
HySoP/setup.py.in | 36 +++--
24 files changed, 565 insertions(+), 426 deletions(-)
delete mode 100644 HySoP/hysop/operator/RemeshingDOp.py
create mode 100644 HySoP/hysop/operator/remeshing.py
create mode 100644 HySoP/hysop/operator/splitting.py
create mode 100644 HySoP/hysop/tools/__init__.py
create mode 100644 HySoP/hysop/tools/timer.py
diff --git a/Examples/mainJM.py b/Examples/mainJM.py
index 31370b43b..4723fbfec 100644
--- a/Examples/mainJM.py
+++ b/Examples/mainJM.py
@@ -1,25 +1,29 @@
# -*- coding: utf-8 -*-
import time
import parmepy as pp
+from math import *
-def vitesse(x):
- return [1., 1., 1.]
+def vitesse(x, y, z):
+ vx = 1. + x
+ vy = - x * y
+ vz = x * y * z + 10.
+ return vx, vy, vz
-def scalaire(x):
- if x[0] < 0.5 and x[1] < 0.5 and x[2] < 0.5:
+def scalaire(x, y, z):
+ if x < 0.5 and y < 0.5 and z < 0.5:
return 1.
else:
return 0.
def run():
- # Parameter
- # timeStep = 0.02
- # FinalTime = 1.
- # outputFilePrefix = './res/RK2_'
- # outputModulo = 0
+ # Parameters
+ timeStep = 0.02
+ finalTime = 1.
+ outputFilePrefix = './res/RK2_'
+ outputModulo = 0
t0 = time.time()
@@ -27,10 +31,10 @@ def run():
box = pp.Box(3, length=[1., 1., 1.], origin=[0., 0., 0.])
## Fields
- scal = pp.ContinuousField(domain=box, name='Scalar')
- velo = pp.ContinuousField(domain=box, name='Velocity', vector=True)
+ scal = pp.AnalyticalField(domain=box, formula=scalaire, name='Scalar')
+ velo = pp.AnalyticalField(domain=box, formula=vitesse, name='Velocity', vector=True)
- ## Operator
+ ## Operators
advec = pp.Transport(velo, scal)
## Solver creation (discretisation of objects is done in solver initialisation)
@@ -40,49 +44,21 @@ def run():
pb = pp.Problem(topo3D, [advec])
## Setting solver to Problem
- pb.setSolver(pp.ParticularSolver)
+ pb.setSolver(pp.ParticularSolver, finalTime, timeStep)
# ## Setting io to Problem
- # pb.setSolver(pp.Printer(...))
+ pb.setIO(pp.Printer(fields=[scal, velo], frequency=outputModulo, outputPrefix=outputFilePrefix))
t1 = time.time()
- # ## Solve problem
- # pb.solve(T=FinalTime, dt=timeStep)
+ ## Solve problem
+ pb.solve()
- # tf = time.time()
+ tf = time.time()
- print pb
- print "\n TODO : -terminer l'initialisation du solver (opérateurs) \n\t -Ajouter les opérateurs annexes (timer et io) au probleme \n\t -Ecrire méthode de résolution"
-
- # print "Total time : ", tf - t0, "sec (CPU)"
+ print "Total time : ", tf - t0, "sec (CPU)"
print "Init time : ", t1 - t0, "sec (CPU)"
- # print "Solving time : ", tf - t1, "sec (CPU)"
-
-
- # InterpolationMethod = pp.Linear(grid=box.discreteDomain, gvalues=velo.discreteVariable)
- # RemeshingMethod = pp.M6Prime(grid=box.discreteDomain, gvalues=scal.discreteVariable)
- # ODESolver = pp.RK2(InterpolationMethod, box.discreteDomain.applyConditions)
- # ## cpu
- # #solver = pp.ParticularSolver(cTspPb, ODESolver, InterpolationMethod, RemeshingMethod)
- # ## gpu
- # solver = pp.GPUParticularSolver(cTspPb,
- # ODESolver, InterpolationMethod, RemeshingMethod,
- # device_type='gpu', src='./examples/mainJM_kernels.cl'
- # )
- # #solver = pp.GPUParticularSolver_GLRender(cTspPb, ODESolver, InterpolationMethod, RemeshingMethod, device_type='gpu', dt=timeStep)
- # cTspPb.setSolver(solver)
- # cTspPb.setPrinter(pp.Printer(frequency=outputModulo, outputPrefix=outputFilePrefix, problem=cTspPb))
-
- # t1 = time.time()
-
- # cTspPb.solve(T=FinalTime, dt=timeStep)
-
- # tf = time.time()
-
- # print "Total time : ", tf - t0, "sec (CPU)"
- # print "Init time : ", t1 - t0, "sec (CPU)"
- # print "Solving time : ", tf - t1, "sec (CPU)"
+ print "Solving time : ", tf - t1, "sec (CPU)"
if __name__ == "__main__":
diff --git a/HySoP/hysop/ParticularSolvers/basic.py b/HySoP/hysop/ParticularSolvers/basic.py
index d6f7189b1..16ffa8a14 100644
--- a/HySoP/hysop/ParticularSolvers/basic.py
+++ b/HySoP/hysop/ParticularSolvers/basic.py
@@ -5,6 +5,11 @@ Particular solver description.
"""
from solver import Solver
from ..fields.continuous import ContinuousField
+from ..operator.Transport import Transport
+from ..operator.remeshing import Remeshing
+from ..operator.splitting import Splitting
+from ..tools.timer import Timer
+from ..tools.printer import Printer
class ParticularSolver(Solver):
@@ -12,76 +17,41 @@ class ParticularSolver(Solver):
Particular solver description.
Link with differents numericals methods used.
"""
- def __init__(self, problem,
- ODESolver = None,
- InterpolationMethod = None,
- RemeshingMethod = None,
- splittingMethod='strang'):
+ def __init__(self, problem, t_end, dt,
+ ODESolver=None,
+ InterpolationMethod=None,
+ RemeshingMethod=None,
+ splittingMethod='strang',
+ timer=None,
+ io=None):
"""
Constructor.
Create a solver description.
"""
Solver.__init__(self, problem)
self.problem = problem
- self.p_position = ContinuousField(domain=problem.domains[0], name='Particle Position')
- self.p_scalar = ContinuousField(domain=problem.domains[0], name='Particle Scalar')
- self.problem.addVariable([self.p_position, self.p_scalar])
- # self.discreteProblem = problem.discreteProblem
- # ## ODE Solver method.
- # self.ODESolver = ODESolver
- # ## Interpolation method.
- # self.InterpolationMethod = InterpolationMethod
- # ## Remeshing method.
- # self.RemeshingMethod = RemeshingMethod
- # ## Splitting Method
- # self.splittingMethod = splittingMethod
- # self.dim = self.discreteProblem.domains[0].dimension
- # self.grid = self.discreteProblem.domains[0]
- # self.gvelo = self.discreteProblem.advecOp.velocity
- # self.gscal = self.discreteProblem.advecOp.scalar
- # self.parts = ParticleField(self.grid)
- # self.ppos = DiscreteVariable(domain=self.parts, dimension=self.gvelo.dimension, name="Particle position")
- # self.parts.setPositions(self.ppos)
- # self.pvelo = DiscreteVariable(domain=self.parts, dimension=self.gvelo.dimension, name="Particle velocity")
- # self.pscal = DiscreteVariable(domain=self.parts, dimension=self.gscal.dimension, scalar=True, name="Particle scalar")
- # #self.pbloc = DiscreteVariable(domain=self.parts, dimension=1, integer=True, initialValue=lambda x: 0, scalar=True)
- # #self.ptag = DiscreteVariable(domain=self.parts, dimension=1, integer=True, initialValue=lambda x: 0, scalar=True)
- # # Interpolation : particles positions, grid velocity -> particles velocity
- # interpol = InterpolationDOp(position=self.ppos,
- # velocityField=self.discreteProblem.advecOp.velocity,
- # resvelo=self.pvelo,
- # method=self.InterpolationMethod)
- # if isinstance(self.RemeshingMethod, Lambda2):
- # # Tag particles : particles positions, particles velocities -> particles bloc, particles Tag
- # tag = TagDOp(partPositions=self.ppos,
- # partVelocities=self.pvelo,
- # grid=self.discreteProblem.domains[0],
- # resbloc=self.pbloc,
- # restag=self.ptag)
- # # Remeshing : particles positions, particles scalar, particles tag, particles bloc -> grid scalar
- # self.remesh = RemeshingDOp(partPositions=self.ppos,
- # partScalar=self.pscal,
- # resscal=self.discreteProblem.advecOp.scalar,
- # method=self.RemeshingMethod)
- # # Advection : grid position, grid velocity, grid scalar -> particles position, particles scalar
- # self.discreteProblem.advecOp.setResultVariable(respos=self.ppos, resscal=self.pscal)
- # self.discreteProblem.advecOp.setMethod(self.ODESolver)
- # if isinstance(self.RemeshingMethod, Lambda2):
- # self.discreteProblem.addOperator([interpol, tag, self.remesh])
- # else:
- # self.discreteProblem.addOperator([self.remesh])
-
- # # Splitting Step TODO:Transformer le splitting en opérateur d'opérateur
- # splitting = []
- # if self.splittingMethod == 'strang':
- # [splitting.append((i, 0.5)) for i in xrange(self.dim - 1)]
- # splitting.append((self.dim - 1, 1.))
- # [splitting.append((self.dim - 2 - i, 0.5)) for i in xrange(self.dim - 1)]
- # else:
- # raise ValueError("Unknown splitting method : ", self.splittingMethod)
- # self.discreteProblem.operator_list = []
- # [self.discreteProblem.operator_list.append((op, 1.0, 0)) for op in self.discreteProblem.operators if not op.is_splittable]
- # [[self.discreteProblem.operator_list.append((op, i[1], i[0])) for op in self.discreteProblem.operators if op.is_splittable] for i in splitting]
+ self.advection = None
+ ## ODE Solver method.
+ self.ODESolver = ODESolver
+ ## Interpolation method.
+ self.InterpolationMethod = InterpolationMethod
+ ## Remeshing method.
+ self.RemeshingMethod = RemeshingMethod
+ ## Splitting Method
+ self.splittingMethod = splittingMethod
+ if timer is None:
+ self.problem.setTimer(Timer(t_end, dt))
+ else:
+ self.problem.setTimer(timer)
+ if io is None:
+ self.problem.setIO(Printer())
+ else:
+ self.problem.setIO(io)
+ for op in self.problem.operators:
+ if isinstance(op, Transport):
+ self.advection = op
+ if self.advection is None:
+ raise ValueError("Particular Solver : Cannot create from a problem with no Transport operator")
def initialize(self):
"""
@@ -89,25 +59,34 @@ class ParticularSolver(Solver):
@param discreteProblem : Problem to initialize.
"""
- # ### Create grid.points TODO: A enlever
- # self.grid.points = np.empty_like(self.gvelo.values, dtype=dtype_real)
- # if self.dim == 1:
- # for x in xrange(self.grid.elementNumber[0]):
- # self.grid.points[x] = self.grid.axes[0][x]
- # if self.dim == 2:
- # for x in xrange(self.grid.elementNumber[0]):
- # for y in xrange(self.grid.elementNumber[1]):
- # self.grid.points[x, y][0] = self.grid.axes[0][x]
- # self.grid.points[x, y][1] = self.grid.axes[1][y]
- # if self.dim == 3:
- # for x in xrange(self.grid.elementNumber[0]):
- # for y in xrange(self.grid.elementNumber[1]):
- # for z in xrange(self.grid.elementNumber[2]):
- # self.grid.points[x, y, z][0] = self.grid.axes[0][x]
- # self.grid.points[x, y, z][1] = self.grid.axes[1][y]
- # self.grid.points[x, y, z][2] = self.grid.axes[2][z]
- # self.gvelo.init()
- # self.gscal.init()
+ self.p_position = ContinuousField(domain=self.problem.domains[0], name='Particle Position')
+ self.p_scalar = ContinuousField(domain=self.problem.domains[0], name='Particle Scalar')
+ self.problem.addVariable([self.p_position, self.p_scalar])
+ ## Discretise domains
+ for d in self.problem.domains:
+ d.discretize(self.problem.topology.resolution)
+ ## Discretise fields and initialize
+ for v in self.problem.variables:
+ v.discretize(self.problem.topology)
+ v.initialize()
+ ## Discretise operators
+ for op in self.problem.operators:
+ if op is self.advection:
+ op.discretize(result_position=self.p_position, result_scalar=self.p_scalar, method=self.ODESolver)
+ else:
+ op.discretize()
+ ## Create remeshing operator
+ remesh = Remeshing(partPositions=self.p_position,
+ partScalar=self.p_scalar,
+ resscal=self.advection.scalar,
+ method=self.RemeshingMethod)
+ for op in self.problem.operators:
+ if op.discreteOperator.needSplitting:
+ if op is self.advection:
+ index = self.problem.operators.index(op)
+ self.problem.operators[index] = Splitting([op, remesh])
+ else:
+ self.problem.operators[index] = Splitting([op])
def __str__(self):
"""ToString method"""
diff --git a/HySoP/hysop/ParticularSolvers/interpolation/interpolation.py b/HySoP/hysop/ParticularSolvers/interpolation/interpolation.py
index d33dd4c05..308ca5ae3 100644
--- a/HySoP/hysop/ParticularSolvers/interpolation/interpolation.py
+++ b/HySoP/hysop/ParticularSolvers/interpolation/interpolation.py
@@ -4,8 +4,6 @@
InterpolationMethod interface.
"""
-providedClass = "InterpolationMethod"
-
class InterpolationMethod:
"""
diff --git a/HySoP/hysop/ParticularSolvers/remesh/method.py b/HySoP/hysop/ParticularSolvers/remesh/method.py
index 5b023b992..ade97b75e 100644
--- a/HySoP/hysop/ParticularSolvers/remesh/method.py
+++ b/HySoP/hysop/ParticularSolvers/remesh/method.py
@@ -3,7 +3,6 @@
@package Utils
RemeshingMethod interface.
"""
-from ..Param import *
class RemeshingMethod:
diff --git a/HySoP/hysop/__init__.py.in b/HySoP/hysop/__init__.py.in
index 86afe6a33..9bf57ac7e 100755
--- a/HySoP/hysop/__init__.py.in
+++ b/HySoP/hysop/__init__.py.in
@@ -4,26 +4,26 @@
Python package dedicated to flow simulation using particular methods on hybrid architectures (MPI-GPU)
"""
-__version__=1.00
-__all__=['Box','CartesianTopology','ScalarField']
+__version__ = 1.00
+__all__ = ['Box', 'CartesianTopology', 'ScalarField']
import os
import site
import mpi4py.MPI as MPI
-if(os.environ.has_key('LD_LIBRARY_PATH')):
- os.environ['LD_LIBRARY_PATH'] = os.environ['LD_LIBRARY_PATH']+':'+"@CMAKE_INSTALL_PREFIX@"
+if 'LD_LIBRARY_PATH' in os.environ:
+ os.environ['LD_LIBRARY_PATH'] = os.environ['LD_LIBRARY_PATH'] + ':' + "@CMAKE_INSTALL_PREFIX@"
else:
os.environ['LD_LIBRARY_PATH'] = "@CMAKE_INSTALL_PREFIX@"
rank_world = MPI.COMM_WORLD.Get_rank()
if(rank_world == 0):
- print "Starting @PACKAGE_NAME@ version "+str(__version__)+".\n"
+ print "Starting @PACKAGE_NAME@ version " + str(__version__) + ".\n"
import constants
import domain.box
-## Box-type physical domain
+## Box-type physical domain
Box = domain.box.Box
## Cartesian grid
@@ -35,5 +35,23 @@ CartesianTopology = fields.topology.CartesianTopology
## Fields
import fields.discrete
import fields.continuous
+import fields.analytical
ScalarField = fields.discrete.ScalarField
ContinuousField = fields.continuous.ContinuousField
+AnalyticalField = fields.analytical.AnalyticalField
+
+## Operators
+import operator.Transport
+Transport = operator.Transport.Transport
+
+## Problem
+import problem.problem
+Problem = problem.problem.Problem
+
+## Solver
+import ParticularSolvers.basic
+ParticularSolver = ParticularSolvers.basic.ParticularSolver
+
+## Tools
+import tools.printer
+Printer = tools.printer.Printer
diff --git a/HySoP/hysop/constants.py b/HySoP/hysop/constants.py
index 2a0ab9d49..1af66a9bc 100644
--- a/HySoP/hysop/constants.py
+++ b/HySoP/hysop/constants.py
@@ -10,7 +10,7 @@ import math
#
PI = math.pi
-# Set default type for real and integer numbers
+# Set default type for real and integer numbers
PARMES_REAL = np.float64
PARMES_INTEGER = np.uint32
PARMES_REAL_GPU = np.float32
@@ -25,3 +25,5 @@ YDIR = 1
ZDIR = 2
## Tag for periodic boundary conditions
PERIODIC = 0
+## Directions string
+S_DIR = ["_X", "_Y", "_Z"]
diff --git a/HySoP/hysop/domain/box.py b/HySoP/hysop/domain/box.py
index ed0528797..c2ce204a7 100644
--- a/HySoP/hysop/domain/box.py
+++ b/HySoP/hysop/domain/box.py
@@ -1,5 +1,5 @@
# -*- coding: utf-8 -*-
-from ..constants import PERIODIC
+from ..constants import *
from .domain import Domain
from .grid import CartesianGrid
import numpy as np
@@ -25,13 +25,13 @@ class Box(Domain):
raise ValueError("Box parameters inconsistents dimensions")
Domain.__init__(self, dimension)
## Box length. length = max - min.
- self.length = np.array(length)
+ self.length = np.asarray(length, dtype=PARMES_REAL)
## Minimum Box position.
- self.origin = np.array(origin)
+ self.origin = np.asarray(origin, dtype=PARMES_REAL)
## Maximum Box position.
self.max = self.origin + self.length
# Boundary conditions type :
- self.boundaries = np.zeros((self.dimension))
+ self.boundaries = np.zeros((self.dimension), dtype=PARMES_INTEGER)
self.boundaries[:] = PERIODIC
self.discreteDomain = None
diff --git a/HySoP/hysop/domain/grid.py b/HySoP/hysop/domain/grid.py
index fe13e3149..6a32f40e4 100644
--- a/HySoP/hysop/domain/grid.py
+++ b/HySoP/hysop/domain/grid.py
@@ -16,6 +16,7 @@ class CartesianGrid(DiscreteDomain):
self.origin = domain.origin
self.length = domain.length
self.resolution = resolution
+ self.size = self.length / self.resolution
## self.elementNumber = spec
## self.elementSize = box.length / self.elementNumber
diff --git a/HySoP/hysop/fields/analytical.py b/HySoP/hysop/fields/analytical.py
index 7e6451064..d4ec1738f 100644
--- a/HySoP/hysop/fields/analytical.py
+++ b/HySoP/hysop/fields/analytical.py
@@ -4,7 +4,8 @@
"""
from .continuous import ContinuousField
-from .discrete import DiscreteField
+from .discrete import ScalarField
+from .discrete import VectorField
class AnalyticalField(ContinuousField):
@@ -13,7 +14,7 @@ class AnalyticalField(ContinuousField):
"""
- def __init__(self, domain, name="", formula=None):
+ def __init__(self, domain, formula, name="", vector=False):
"""
Constructor.
Create an AnalyticalField from a formula.
@@ -22,29 +23,24 @@ class AnalyticalField(ContinuousField):
@param dimension : variable dimension.
@param formula : function defining the variable.
"""
- ContinuousField.__init__(self, domain, name)
+ ContinuousField.__init__(self, domain, name, vector)
## Analytic formula.
self.formula = formula
- def discretize(self, topology=None):
- """
- AnalyticalField discretization method.
- Create an DiscreteField from given topology.
- """
- self.__fieldId += 1
- self.discreteField.append(ScalarField(self, topology, name=self.name + "_D_" + str(self.__fieldId),
- idFromParent=self.__fieldId))
- self.discreteField[self.__fieldId].initialize(self.formula)
- return self.discreteField[self.__fieldId], self.__fieldId
-
- def value(self, pos):
+ 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)
+ return self.formula(*pos)
+
+ def initialize(self):
+ if self._fieldId == -1:
+ raise ValueError("Cannot initialise analytical field non discretized.")
+ for dF in self.discreteField:
+ dF.initialize(self.formula)
def __str__(self):
"""ToString method"""
diff --git a/HySoP/hysop/fields/continuous.py b/HySoP/hysop/fields/continuous.py
index 284d4d404..7bfb54eb6 100644
--- a/HySoP/hysop/fields/continuous.py
+++ b/HySoP/hysop/fields/continuous.py
@@ -27,7 +27,7 @@ class ContinuousField(object):
## Name of this field
self.name = name
# number of different discretizations
- self.__fieldId = -1
+ self._fieldId = -1
## list of the various discretizations of this field
self.discreteField = []
## is field is a vector field
@@ -38,19 +38,22 @@ class ContinuousField(object):
discretization of the field on a topology
"""
- self.__fieldId += 1
+ self._fieldId += 1
print "Discretisation de ", self.name,
if self.vector:
print "Vector Field"
self.discreteField.append(VectorField(self, topology,
- name=self.name + "_D_" + str(self.__fieldId),
- idFromParent=self.__fieldId))
+ name=self.name + "_D_" + str(self._fieldId),
+ idFromParent=self._fieldId))
else:
print "Scalar Field"
self.discreteField.append(ScalarField(self, topology,
- name=self.name + "_D_" + str(self.__fieldId),
- idFromParent=self.__fieldId))
- return self.discreteField[self.__fieldId], self.__fieldId
+ name=self.name + "_D_" + str(self._fieldId),
+ idFromParent=self._fieldId))
+ return self.discreteField[self._fieldId], self._fieldId
+
+ def initialize(self):
+ pass
def __str__(self):
""" Field info display """
diff --git a/HySoP/hysop/fields/discrete.py b/HySoP/hysop/fields/discrete.py
index 8f05381c2..bccf1be47 100644
--- a/HySoP/hysop/fields/discrete.py
+++ b/HySoP/hysop/fields/discrete.py
@@ -27,6 +27,7 @@ class ScalarField(object):
self.data = np.zeros((resolution), dtype=PARMES_REAL, order=ORDER)
self.__parentField = parent
self.__idFromParent = idFromParent
+ self.domain = self.__parentField.domain.discreteDomain
def __str__(self):
""" Display class information """
@@ -52,32 +53,20 @@ class ScalarField(object):
def initialize(self, formula=None):
"""
Initialize values with given formula
- \todo Ecrire l'initialisation efficace en python.
"""
- raise NotImplementedError("Not yet implemented")
-
- ## def init(self):
- ## if not self.initialValue is None:
- ## if self.domain.dimension == 1:
- ## for x in xrange(self.shape[0]):
- ## self.values[x] = self.initialValue([self.domain.axes[0][x]])
- ## self.contains_data = True
- ## elif self.domain.dimension == 2:
- ## for x in xrange(self.shape[0]):
- ## for y in xrange(self.shape[1]):
- ## self.values[x, y] = self.initialValue([self.domain.axes[0][x], self.domain.axes[1][y]])
- ## self.contains_data = True
- ## elif self.domain.dimension == 3:
- ## for x in xrange(self.shape[0]):
- ## for y in xrange(self.shape[1]):
- ## for z in xrange(self.shape[2]):
- ## self.values[x, y, z] = self.initialValue([self.domain.axes[0][x], self.domain.axes[1][y], self.domain.axes[2][z]])
- ## self.contains_data = True
- ## else:
- ## raise ValueError("Creating DiscreteVariable with incorrect domain dimension.")
+ v_formula = np.vectorize(formula)
+ if self.dimension == 3:
+ self.data = v_formula(self.topology.mesh.coords[0],
+ self.topology.mesh.coords[1],
+ self.topology.mesh.coords[2])
+ elif self.dimension == 2:
+ self.data = v_formula(self.topology.mesh.coords[0],
+ self.topology.mesh.coords[1])
+ else:
+ self.data = v_formula(self.topology.mesh.coords[0])
def tofile(self, filename):
- evtk.imageToVTK(filename, pointData={"scalar": self.data})
+ evtk.imageToVTK(filename, pointData={self.name: self.data})
class VectorField(object):
@@ -99,6 +88,7 @@ class VectorField(object):
self.data = [np.zeros((resolution), dtype=PARMES_REAL, order=ORDER) for d in xrange(self.dimension)]
self.__parentField = parent
self.__idFromParent = idFromParent
+ self.domain = self.__parentField.domain.discreteDomain
def __str__(self):
""" Display class information """
@@ -141,12 +131,20 @@ class VectorField(object):
else:
self.data[i[-1]].__getitem__(tuple(i[0:-1]))
- def initialize(self, formula):
+ def initialize(self, formula=None):
"""
Initialize values with given formula
- \todo Ecrire l'initialisation efficace en python.
"""
- raise NotImplementedError("Not yet implemented")
+ #print self.topology.mesh.coords
+ if formula is not None:
+ v_formula = np.vectorize(formula)
+ if self.dimension == 3:
+ self.data[0], self.data[1], self.data[2] = v_formula(self.topology.mesh.coords[0],
+ self.topology.mesh.coords[1],
+ self.topology.mesh.coords[2])
+ elif self.dimension == 2:
+ self.data[0], self.data[1] = v_formula(self.topology.mesh.coords[0],
+ self.topology.mesh.coords[1])
if __name__ == "__main__":
print __doc__
diff --git a/HySoP/hysop/fields/topology.py b/HySoP/hysop/fields/topology.py
index 0b74994e2..41448182e 100644
--- a/HySoP/hysop/fields/topology.py
+++ b/HySoP/hysop/fields/topology.py
@@ -26,11 +26,11 @@ class CartesianTopology(object):
nbprocs = comm.Get_size()
if(dims is None):
assert(dim is not None)
- self.dims = np.asarray(MPI.Compute_dims(nbprocs, dim))
+ self.dims = np.asarray(MPI.Compute_dims(nbprocs, dim), dtype=PARMES_INTEGER)
# Topology dimension
self.dim = dim
else:
- self.dims = np.asarray(dims)
+ self.dims = np.asarray(dims, dtype=PARMES_INTEGER)
self.dim = self.dims.size
# Define the topology
# default period status is periodic
@@ -53,9 +53,9 @@ class CartesianTopology(object):
self.south, self.north = self.topo.Shift(ZDIR, 1)
# ghost layer (default = 0)
if(ghosts is None):
- self.ghosts = np.zeros((self.domain.dimension))
+ self.ghosts = np.zeros((self.domain.dimension), dtype=PARMES_INTEGER)
# Global resolution
- self.resolution = np.asarray(resolution)
+ self.resolution = np.asarray(resolution, dtype=PARMES_INTEGER)
nbpoints = self.resolution[:self.dim] // self.dims
remainingPoints = self.resolution[:self.dim] % self.dims
localResolution = self.resolution.copy()
@@ -63,9 +63,9 @@ class CartesianTopology(object):
for i in range(self.dim):
if(self.coords[i] == self.dims[i] - 1):
localResolution[i] = nbpoints[i] + remainingPoints[i]
- start = np.zeros((domain.dimension))
+ start = np.zeros((domain.dimension), dtype=PARMES_INTEGER)
start[:self.dim] = self.coords * nbpoints
- self.mesh = LocalMesh(self.rank, resolution=localResolution, start=start, ghosts=self.ghosts)
+ self.mesh = LocalMesh(self.rank, resolution=localResolution, start=start, dom_size=self.domain.length/self.resolution, ghosts=self.ghosts)
def __str__(self):
""" Topology info display """
@@ -87,8 +87,8 @@ class FFTTopology(object):
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)
+ self.dims = np.ones((dim + 1), dtype=PARMES_INTEGER)
+ self.dims[0:dim] = np.asarray(dims, dtype=PARMES_INTEGER)
# Define the topology
self.topo = comm.Create_cart(self.dims, periods=periods[0:dim + 1])
# Size of the topology
@@ -117,7 +117,7 @@ class FFTTopology(object):
class LocalMesh(object):
"""
"""
- def __init__(self, rank, resolution, start, ghosts=np.zeros((3))):
+ def __init__(self, rank, resolution, start, dom_size, ghosts=np.zeros((3))):
# Current process rank in the topology that owns this mesh
self.rank = rank
# Local resolution
@@ -127,6 +127,19 @@ class LocalMesh(object):
# index of the upper point (in each dir), global mesh
self.end = self.start + self.resolution - 1
# Mesh step size
+ self.size = dom_size
+ # Mesh coordinates
+ coord_start = self.start * self.size
+ coord_end = (self.end + 1) * self.size
+ if self.start.size == 3:
+ self.coords = np.ogrid[coord_start[0]:coord_end[0]:self.size[0],
+ coord_start[1]:coord_end[1]:self.size[1],
+ coord_start[2]:coord_end[2]:self.size[2]]
+ elif self.start.size == 2:
+ self.coords = np.ogrid[coord_start[0]:coord_end[0]:self.size[0],
+ coord_start[1]:coord_end[1]:self.size[1]]
+ else:
+ self.coords = np.ogrid[coord_start[0]:coord_end[0]:self.size[0]]
def __str__(self):
""" Local mesh display """
diff --git a/HySoP/hysop/operator/RemeshingDOp.py b/HySoP/hysop/operator/RemeshingDOp.py
deleted file mode 100644
index 6942158c4..000000000
--- a/HySoP/hysop/operator/RemeshingDOp.py
+++ /dev/null
@@ -1,113 +0,0 @@
-# -*- coding: utf-8 -*-
-"""
-@package Operator
-Operator representation
-"""
-from ..Param import *
-from DiscreteOperator import DiscreteOperator
-import pyopencl as cl
-import time
-
-
-class RemeshingDOp(DiscreteOperator):
- """
- Remeshing operator representation.
- DiscreteOperator.DiscreteOperator specialization.
- """
-
- def __init__(self, partPositions, partScalar, resscal, method, partBloc=None, partTag=None):
- """
- Constructor.
- 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 partScalar : particles scalar.
- @param partBloc : particles bloc number.
- @param partTag : particle tag number.
- @param resscal DiscreteVariable.DiscreteVariable : new grid scalar values.
- """
- DiscreteOperator.__init__(self)
- ## Particles positions.
- self.ppos = partPositions
- ## Particles scalar values.
- self.pscal = partScalar
- ## Particles bloc
- self.pbloc = partBloc
- ## Particles tag
- self.ptag = partTag
- ## Result scalar
- self.resultScalar = resscal
- self.addVariable([self.ppos, self.pscal, self.resultScalar])
- if not self.pbloc is None:
- self.addVariable([self.pbloc])
- if not self.ptag is None:
- self.addVariable([self.ptag])
- self.numMethod = method
- self.gpu_kernel_name = "remeshing"
- self.compute_time = [0., 0., 0.]
- self.queued_time = 0.
- self.submit_time = 0.
- self.total_flop = 0
- self.total_bytes_accessed = 0
-
- def setResultVariable(self, resscal):
- """
- Set results variables
-
- @param resscal DiscreteVariable.DiscreteVariable : new grid scalar values.
- """
- ## Result scalar
- self.resultScalar = resscal
- self.addVariable([self.resultScalar])
-
- def applyOperator(self, t, dt, splittingDirection):
- """
- Apply Remeshing operator.
- """
- c_time = time.time()
- if self.gpu_kernel is None:
- c_time = time.time()
- ## ------------------- NumPy Optimisation
- self.resultScalar.values[...] = 0.
- self.numMethod.remesh(self.ppos.values, self.pscal.values, splittingDirection)
- ## -------------------
- self.compute_time[splittingDirection] += (time.time() - c_time)
- self.total_flop += self.numMethod.nb_flop * np.prod(self.resultScalar.values.shape)
- else:
- l = list(self.gpu_shape)
- nb_ligne = l.pop(splittingDirection)
- local_wg = list(l)
- if len(l) == 2:
- local_wg[0] = 1
- local_wg = tuple(local_wg)
- else:
- local_wg = None
- evt = self.gpu_kernel.remeshing(self.gpu_queue, tuple(l), local_wg,
- self.ppos.gpu_mem_object, self.pscal.gpu_mem_object, self.resultScalar.gpu_mem_object,
- dtype_real_gpu(t), dtype_real_gpu(dt), dtype_integer(splittingDirection), dtype_integer(nb_ligne),
- self.resultScalar.domain.min.astype(dtype_real_gpu),
- self.resultScalar.domain.max.astype(dtype_real_gpu),
- self.resultScalar.domain.elementNumber.astype(dtype_integer),
- self.resultScalar.domain.elementSize.astype(dtype_real_gpu)
- )
- # self.gpu_queue.finish()
- # cl.enqueue_copy(self.gpu_queue, self.resultScalar.values, self.resultScalar.gpu_mem_object)
- # print self.resultScalar.values
- self.gpu_queue.finish()
- self.queued_time += (evt.profile.submit - evt.profile.queued) * 1e-9
- self.submit_time += (evt.profile.start - evt.profile.submit) * 1e-9
- self.compute_time[splittingDirection] += (evt.profile.end - evt.profile.start) * 1e-9
- self.total_flop += (75 * nb_ligne) * np.prod(l)
- self.total_bytes_accessed += ((2 + 1) * nb_ligne) * np.prod(l) * 4 # 2 float lus + 1 float écrits
- #print "Remeshing M'6:", self.compute_time * 1e-9
-
- def __str__(self):
- """ToString method"""
- return "RemeshingDOp (DiscreteOperator)"
-
-
-if __name__ == "__main__":
- print __doc__
- print "- Provided class : RemeshingDOp"
- print RemeshingDOp.__doc__
diff --git a/HySoP/hysop/operator/Transport.py b/HySoP/hysop/operator/Transport.py
index 00541de94..8bf58be17 100644
--- a/HySoP/hysop/operator/Transport.py
+++ b/HySoP/hysop/operator/Transport.py
@@ -28,7 +28,7 @@ class Transport(ContinuousOperator):
self.scalar = scalar
self.addVariable([velocity, scalar])
- def discretize(self, idVelocityD=0, idScalarD=0, result_position=None, result_scalar=None):
+ def discretize(self, idVelocityD=0, idScalarD=0, result_position=None, result_scalar=None, method=None):
"""
Advection operator discretization method.
Create an AdvectionDOp.AdvectionDOp from given specifications.
diff --git a/HySoP/hysop/operator/advection_d.py b/HySoP/hysop/operator/advection_d.py
index 94084acd9..9cb78ff52 100644
--- a/HySoP/hysop/operator/advection_d.py
+++ b/HySoP/hysop/operator/advection_d.py
@@ -14,7 +14,7 @@ class Advection_P(DiscreteOperator):
DiscreteOperator.DiscreteOperator specialization.
"""
- def __init__(self, advec, idVelocityD=0, idScalarD=0, result_position=None, result_scalar=None):
+ def __init__(self, advec, idVelocityD=0, idScalarD=0, result_position=None, result_scalar=None, method=None):
"""
Constructor.
Create a Advection operator on a given continuous domain.
@@ -23,18 +23,20 @@ class Advection_P(DiscreteOperator):
@param advec : Advection operator.
"""
DiscreteOperator.__init__(self)
- ## Velocity. \todo utiliser un vectorField
+ ## Velocity.
self.velocity = advec.velocity.discreteField[idVelocityD]
## Transported scalar.
self.scalar = advec.scalar.discreteField[idScalarD]
## Result position
- self.res_position = result_position
+ self.res_position = result_position.discreteField[0]
## Result scalar
- self.res_scalar = result_scalar
+ self.res_scalar = result_scalar.discreteField[0]
## input fields
self.input = [self.velocity, self.scalar]
## output fields
self.output = [self.res_position, self.res_scalar]
+ self.numMethod = method
+ self.needSplitting = True
# self.gpu_kernel_name = "advection"
# self.compute_time = [0., 0., 0.]
# self.submit_time = 0.
@@ -42,9 +44,14 @@ class Advection_P(DiscreteOperator):
# self.total_flop = 0
# self.total_bytes_accessed = 0
- def applyOperator(self, t, dt, splittingDirection):
- raise NotImplementedMethod("Not yet implemented")
- # """
+ def apply(self, t, dt, splittingDirection):
+ c_time = time.time()
+ print "Advection operator step, direction :", splittingDirection,
+ c_time = (time.time() - c_time)
+ self.compute_time += c_time
+ if self.numMethod is None:
+ print "nothing to do"
+ # """
# Apply advection operator.
# @param t : current time.
@@ -118,13 +125,7 @@ class Advection_P(DiscreteOperator):
# #print "Advection:", self.compute_time * 1e-9
def __str__(self):
- """ToString method"""
- s = "Advection_P (DiscreteOperator). Applying on following fields : \n Input fields : \n"
- for f in self.input:
- s += str(f) + "\n"
- s += "Output fields : \n"
- for f in self.output:
- s += str(f) + "\n"
+ s = "Advection_P (DiscreteOperator). " + DiscreteOperator.__str__(self)
return s
if __name__ == "__main__":
diff --git a/HySoP/hysop/operator/continuous.py b/HySoP/hysop/operator/continuous.py
index ed2b750a0..a1083989a 100644
--- a/HySoP/hysop/operator/continuous.py
+++ b/HySoP/hysop/operator/continuous.py
@@ -36,6 +36,9 @@ class ContinuousOperator:
if self.variables.count(v) == 0:
self.variables.append(v)
+ def apply(self, *args):
+ self.discreteOperator.apply(*args)
+
@abstractmethod
def discretize(self, spec=None):
"""
diff --git a/HySoP/hysop/operator/discrete.py b/HySoP/hysop/operator/discrete.py
index 2578a9755..301674518 100644
--- a/HySoP/hysop/operator/discrete.py
+++ b/HySoP/hysop/operator/discrete.py
@@ -19,13 +19,16 @@ class DiscreteOperator:
Constructor.
Create an empty discrete operator.
"""
- ## Application domains of the variables.
- self.domains = []
+ self.input, self.output = None, None
+ ## variables
+ self.variables = []
## Operator numerical method.
self.numMethod = []
- self.is_splittable = True
+ self.needSplitting = False
self.gpu_kernel = None
self.gpu_kernel_name = ""
+ self.discreteOperator = self
+ self.compute_time = 0.
def setMethod(self, method):
"""
@@ -35,8 +38,19 @@ class DiscreteOperator:
"""
self.numMethod = method
+ def addVariable(self, cVariable):
+ """
+ Add an continuous variables to the operator.
+ Also add variables' domains to the operator.
+
+ @param cVariable ContinuousVariable.ContinuousVariable : variables to add.
+ """
+ for v in cVariable:
+ if self.variables.count(v) == 0:
+ self.variables.append(v)
+
@abstractmethod
- def applyOperator(self):
+ def apply(self):
"""
Abstract method, apply operaton on a variable.
Must be implemented by sub-class.
@@ -45,6 +59,19 @@ class DiscreteOperator:
"""
raise NotImplementedError("Need to override method in a subclass of " + providedClass)
+ def __str__(self):
+ """ToString method"""
+ s = "Applying on following fields : \n"
+ if self.input is not None:
+ s += "Input fields : \n"
+ for f in self.input:
+ s += str(f) + "\n"
+ if self.output is not None:
+ s += "Output fields : \n"
+ for f in self.output:
+ s += str(f) + "\n"
+ return s
+
if __name__ == "__main__":
print __doc__
print "- Provided class : DiscreteOperator"
diff --git a/HySoP/hysop/operator/remeshing.py b/HySoP/hysop/operator/remeshing.py
new file mode 100644
index 000000000..c49c03753
--- /dev/null
+++ b/HySoP/hysop/operator/remeshing.py
@@ -0,0 +1,105 @@
+# -*- coding: utf-8 -*-
+"""
+@package Operator
+Operator representation
+"""
+from .discrete import DiscreteOperator
+import pyopencl as cl
+import time
+
+
+class Remeshing(DiscreteOperator):
+ """
+ Remeshing operator representation.
+ DiscreteOperator.DiscreteOperator specialization.
+ """
+
+ def __init__(self, partPositions, partScalar, resscal, method):
+ """
+ Constructor.
+ 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 partScalar : particles scalar.
+ @param partBloc : particles bloc number.
+ @param partTag : particle tag number.
+ @param resscal DiscreteVariable.DiscreteVariable : new grid scalar values.
+ """
+ DiscreteOperator.__init__(self)
+ ## Particles positions.
+ self.ppos = partPositions
+ ## Particles scalar values.
+ self.pscal = partScalar
+ ## Result scalar
+ self.resultScalar = resscal
+ self.addVariable([self.ppos, self.pscal, self.resultScalar])
+ self.numMethod = method
+ ## input fields
+ self.input = [self.ppos, self.pscal]
+ ## output fields
+ self.output = [self.resultScalar]
+ self.needSplitting = True
+ # self.gpu_kernel_name = "remeshing"
+ # self.compute_time = [0., 0., 0.]
+ # self.queued_time = 0.
+ # self.submit_time = 0.
+ # self.total_flop = 0
+ # self.total_bytes_accessed = 0
+
+ def apply(self, t, dt, splittingDirection):
+ """
+ Apply Remeshing operator.
+ """
+ c_time = time.time()
+ print "Remeshing operator step, direction :", splittingDirection,
+ c_time = (time.time() - c_time)
+ self.compute_time += c_time
+ if self.numMethod is None:
+ print "nothing to do"
+
+ # c_time = time.time()
+ # if self.gpu_kernel is None:
+ # c_time = time.time()
+ # ## ------------------- NumPy Optimisation
+ # self.resultScalar.values[...] = 0.
+ # self.numMethod.remesh(self.ppos.values, self.pscal.values, splittingDirection)
+ # ## -------------------
+ # self.compute_time[splittingDirection] += (time.time() - c_time)
+ # self.total_flop += self.numMethod.nb_flop * np.prod(self.resultScalar.values.shape)
+ # else:
+ # l = list(self.gpu_shape)
+ # nb_ligne = l.pop(splittingDirection)
+ # local_wg = list(l)
+ # if len(l) == 2:
+ # local_wg[0] = 1
+ # local_wg = tuple(local_wg)
+ # else:
+ # local_wg = None
+ # evt = self.gpu_kernel.remeshing(self.gpu_queue, tuple(l), local_wg,
+ # self.ppos.gpu_mem_object, self.pscal.gpu_mem_object, self.resultScalar.gpu_mem_object,
+ # dtype_real_gpu(t), dtype_real_gpu(dt), dtype_integer(splittingDirection), dtype_integer(nb_ligne),
+ # self.resultScalar.domain.min.astype(dtype_real_gpu),
+ # self.resultScalar.domain.max.astype(dtype_real_gpu),
+ # self.resultScalar.domain.elementNumber.astype(dtype_integer),
+ # self.resultScalar.domain.elementSize.astype(dtype_real_gpu)
+ # )
+ # # self.gpu_queue.finish()
+ # # cl.enqueue_copy(self.gpu_queue, self.resultScalar.values, self.resultScalar.gpu_mem_object)
+ # # print self.resultScalar.values
+ # self.gpu_queue.finish()
+ # self.queued_time += (evt.profile.submit - evt.profile.queued) * 1e-9
+ # self.submit_time += (evt.profile.start - evt.profile.submit) * 1e-9
+ # self.compute_time[splittingDirection] += (evt.profile.end - evt.profile.start) * 1e-9
+ # self.total_flop += (75 * nb_ligne) * np.prod(l)
+ # self.total_bytes_accessed += ((2 + 1) * nb_ligne) * np.prod(l) * 4 # 2 float lus + 1 float écrits
+ # #print "Remeshing M'6:", self.compute_time * 1e-9
+
+ def __str__(self):
+ s = "Advection_P (DiscreteOperator). " + DiscreteOperator.__str__(self)
+ return s
+
+if __name__ == "__main__":
+ print __doc__
+ print "- Provided class : RemeshingDOp"
+ print RemeshingDOp.__doc__
diff --git a/HySoP/hysop/operator/splitting.py b/HySoP/hysop/operator/splitting.py
new file mode 100644
index 000000000..82b6657c1
--- /dev/null
+++ b/HySoP/hysop/operator/splitting.py
@@ -0,0 +1,55 @@
+# -*- coding: utf-8 -*-
+"""
+@package Operator
+Operator representation
+"""
+from .discrete import DiscreteOperator
+import time
+
+
+class Splitting(DiscreteOperator):
+ """
+ Splitting operator representation.
+ DiscreteOperator.DiscreteOperator specialization.
+ """
+
+ def __init__(self, operators=[], dim=3):
+ """
+ Constructor.
+ 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.
+
+ """
+ DiscreteOperator.__init__(self)
+ self.operators = operators
+ self.splitting = []
+ [self.splitting.append((i, 0.5)) for i in xrange(dim - 1)]
+ self.splitting.append((dim - 1, 1.))
+ [self.splitting.append((dim - 2 - i, 0.5)) for i in xrange(dim - 1)]
+
+ def apply(self, t, dt):
+ """
+ Apply Remeshing operator.
+ """
+ print "Splitting operator step :"
+ c_time = time.time()
+ for split in self.splitting:
+ print "direction : " + str(split[0]) + " dt*" + str(split[1])
+ for op in self.operators:
+ op.apply(t, split[1] * dt, split[0])
+ c_time = (time.time() - c_time)
+ self.compute_time += c_time
+
+ def __str__(self):
+ s = "Splitting (DiscreteOperator). Splitting steps : \n"
+ for split in self.splitting:
+ s += "direction : " + str(split[0]) + " dt=" + str(split[1]) + "\n"
+ s += "Concerned operators : \n"
+ for op in self.operators:
+ s += str(op)
+ return s
+
+if __name__ == "__main__":
+ print __doc__
+ print "- Provided class : RemeshingDOp"
+ print RemeshingDOp.__doc__
diff --git a/HySoP/hysop/problem/problem.py b/HySoP/hysop/problem/problem.py
index 2aa412159..74e40d065 100644
--- a/HySoP/hysop/problem/problem.py
+++ b/HySoP/hysop/problem/problem.py
@@ -25,18 +25,27 @@ class Problem():
self.addDomain(v.domain)
self.topology = topology
self.solver = None
+ self.timer = None
+ self.io = None
- def setSolver(self, solver):
- self.solver = solver(self)
- ## Discretise domains
- for d in self.domains:
- d.discretize(self.topology.resolution)
- ## Discretise fields
- for v in self.variables:
- v.discretize(self.topology)
- ## Discretise operators
- for op in self.operators:
- op.discretize(result_position=self.solver.p_position, result_scalar=self.solver.p_scalar)
+ def setTimer(self, t):
+ self.timer = t
+
+ def setIO(self, io):
+ self.io = io
+
+ def setSolver(self, solver, t_end, dt):
+ self.solver = solver(self, t_end, dt)
+ self.solver.initialize()
+
+ def solve(self):
+ print"\n\n Start solving ..."
+ while not self.timer.end:
+ print "\n==== Iteration : " + str(self.timer.ite + 1) + "\t t=" + str(self.timer.t + self.timer.dt) + " ===="
+ for op in self.operators:
+ op.apply(self.timer.t, self.timer.dt)
+ self.timer.step()
+ self.io.step()
def addVariable(self, cVariable):
"""
@@ -67,6 +76,11 @@ class Problem():
if self.domains.count(cDomain) == 0:
self.domains.append(cDomain)
+ def addOperator(self, op, index=None):
+ if index is None:
+ if not op in self.operators:
+ self.operators.append(op)
+
def __str__(self):
"""ToString method"""
s = "Problem based on\n"
diff --git a/HySoP/hysop/tools/__init__.py b/HySoP/hysop/tools/__init__.py
new file mode 100644
index 000000000..e69de29bb
diff --git a/HySoP/hysop/tools/printer.py b/HySoP/hysop/tools/printer.py
index 8fbb07ea1..0fc402187 100644
--- a/HySoP/hysop/tools/printer.py
+++ b/HySoP/hysop/tools/printer.py
@@ -3,8 +3,8 @@
@package Utils
Results Printer
"""
-from ..Param import *
-from ..Utils.GPUParticularSolver import GPUParticularSolver
+from ..constants import *
+import evtk.hl as evtk
class Printer():
@@ -12,7 +12,7 @@ class Printer():
Results printer.
"""
- def __init__(self, frequency=0, outputPrefix='./out', problem=None):
+ def __init__(self, fields=[], frequency=0, outputPrefix='./out_'):
"""
Constructor.
Create a results printer from a ContinuousProblem.ContinuousProblem with a given output prefix (relative path) and an output frequency.
@@ -23,81 +23,88 @@ class Printer():
"""
self.freq = frequency
self.outputPrefix = outputPrefix
- self.problem = problem
+ self.fields = fields
self.ite = 0
if self.freq != 0:
- self.printStep = self._printStep
+ self.step = self._printStep
else:
- self.printStep = self._passStep
- self.skip = 1
- self.grid = self.problem.discreteProblem.advecOp.velocity.domain
- self.gvelo = self.problem.discreteProblem.advecOp.velocity.values
- self.ppos = self.problem.discreteProblem.solver.ppos.values
- self.pscal = self.problem.discreteProblem.solver.pscal.values
- self.gscal = self.problem.discreteProblem.advecOp.scalar.values
- self.printStep()
+ self.step = self._passStep
+ self.__dict_data = {}
+ for f in self.fields:
+ for df in f.discreteField:
+ if f.vector:
+ for d in xrange(len(df.data)):
+ self.__dict_data[df.name + S_DIR[d]] = df.data[d]
+ else:
+ self.__dict_data[df.name] = df.data
+ self.step()
def _passStep(self):
pass
def _printStep(self):
- self.skip -= 1
- if self.skip == 0:
- if isinstance(self.problem.discreteProblem.solver, GPUParticularSolver) and self.ite > 0:
- self.problem.discreteProblem.solver.collect_data()
- print "Writing file : " + self.outputPrefix + "results_{0:05d}.dat".format(self.ite)
- self.skip = self.freq
- f = open(self.outputPrefix + "results_{0:05d}.dat".format(self.ite), 'w')
- dim = self.grid.dimension
- f.write("# iteration : {0}, time = {1}".format(self.ite, self.problem.discreteProblem.t))
- f.write("# gid({0}) gpos({1}) gvelo({2}) gscal(1) ppos({3}) pvelo({4}) pscal(1) pbloc(1) ptag(1)\n".format(dim, dim, dim, dim, dim))
- if self.grid.dimension == 1:
- for x in xrange(self.grid.elementNumber[0]):
- f.write(" {0:6d}".format(x))
- f.write(" {:8.12}".format(self.grid.axes[0][x]))
- f.write(" {:8.12}".format(self.gvelo[x][0]))
- f.write(" {:8.12}".format(self.gscal[x]))
- f.write(" {:8.12}".format(self.ppos[x][0]))
- f.write(" {:8.12}".format(self.pscal[x]))
- f.write("\n")
- f.write("\n")
- if self.grid.dimension == 2:
- for x in xrange(self.grid.elementNumber[0]):
- for y in xrange(self.grid.elementNumber[1]):
- f.write(" {0:6d}".format(x))
- f.write(" {0:6d}".format(y))
- f.write(" {:8.12}".format(self.grid.axes[0][x]))
- f.write(" {:8.12}".format(self.grid.axes[1][y]))
- f.write(" {:8.12}".format(self.gvelo[x, y][0]))
- f.write(" {:8.12}".format(self.gvelo[x, y][1]))
- f.write(" {:8.12}".format(self.gscal[x, y]))
- #f.write(" {:8.12}".format(self.ppos[x, y][0]))
- #f.write(" {:8.12}".format(self.ppos[x, y][1]))
- f.write(" {:8.12}".format(self.pscal[x, y]))
- f.write("\n")
- f.write("\n")
- if self.grid.dimension == 3:
- for x in xrange(self.grid.elementNumber[0]):
- for y in xrange(self.grid.elementNumber[1]):
- for z in xrange(self.grid.elementNumber[2]):
- f.write(" {0:6d}".format(x))
- f.write(" {0:6d}".format(y))
- f.write(" {0:6d}".format(z))
- f.write(" {:8.12}".format(self.grid.axes[0][x]))
- f.write(" {:8.12}".format(self.grid.axes[1][y]))
- f.write(" {:8.12}".format(self.grid.axes[2][z]))
- f.write(" {:8.12}".format(self.gvelo[x, y, z][0]))
- f.write(" {:8.12}".format(self.gvelo[x, y, z][1]))
- f.write(" {:8.12}".format(self.gvelo[x, y, z][2]))
- f.write(" {:8.12}".format(self.gscal[x, y, z]))
- f.write(" {:8.12}".format(self.ppos[x, y, z][0]))
- f.write(" {:8.12}".format(self.ppos[x, y, z][1]))
- f.write(" {:8.12}".format(self.ppos[x, y, z][2]))
- f.write(" {:8.12}".format(self.pscal[x, y, z]))
- f.write("\n")
- f.write("\n")
- f.close()
- self.ite += 1
+ if (self.ite % self.freq) == 0:
+ filename = self.outputPrefix + "results_{0:05d}.dat".format(self.ite)
+ print "Print file : " + filename
+ evtk.imageToVTK(filename, pointData=self.__dict_data)
+ self.ite += 1
+ # self.skip -= 1
+ # if self.skip == 0:
+ # if isinstance(self.problem.discreteProblem.solver, GPUParticularSolver) and self.ite > 0:
+ # self.problem.discreteProblem.solver.collect_data()
+ # print "Writing file : " + self.outputPrefix + "results_{0:05d}.dat".format(self.ite)
+ # self.skip = self.freq
+ # f = open(self.outputPrefix + "results_{0:05d}.dat".format(self.ite), 'w')
+ # dim = self.grid.dimension
+ # f.write("# iteration : {0}, time = {1}".format(self.ite, self.problem.discreteProblem.t))
+ # f.write("# gid({0}) gpos({1}) gvelo({2}) gscal(1) ppos({3}) pvelo({4}) pscal(1) pbloc(1) ptag(1)\n".format(dim, dim, dim, dim, dim))
+ # if self.grid.dimension == 1:
+ # for x in xrange(self.grid.elementNumber[0]):
+ # f.write(" {0:6d}".format(x))
+ # f.write(" {:8.12}".format(self.grid.axes[0][x]))
+ # f.write(" {:8.12}".format(self.gvelo[x][0]))
+ # f.write(" {:8.12}".format(self.gscal[x]))
+ # f.write(" {:8.12}".format(self.ppos[x][0]))
+ # f.write(" {:8.12}".format(self.pscal[x]))
+ # f.write("\n")
+ # f.write("\n")
+ # if self.grid.dimension == 2:
+ # for x in xrange(self.grid.elementNumber[0]):
+ # for y in xrange(self.grid.elementNumber[1]):
+ # f.write(" {0:6d}".format(x))
+ # f.write(" {0:6d}".format(y))
+ # f.write(" {:8.12}".format(self.grid.axes[0][x]))
+ # f.write(" {:8.12}".format(self.grid.axes[1][y]))
+ # f.write(" {:8.12}".format(self.gvelo[x, y][0]))
+ # f.write(" {:8.12}".format(self.gvelo[x, y][1]))
+ # f.write(" {:8.12}".format(self.gscal[x, y]))
+ # #f.write(" {:8.12}".format(self.ppos[x, y][0]))
+ # #f.write(" {:8.12}".format(self.ppos[x, y][1]))
+ # f.write(" {:8.12}".format(self.pscal[x, y]))
+ # f.write("\n")
+ # f.write("\n")
+ # if self.grid.dimension == 3:
+ # for x in xrange(self.grid.elementNumber[0]):
+ # for y in xrange(self.grid.elementNumber[1]):
+ # for z in xrange(self.grid.elementNumber[2]):
+ # f.write(" {0:6d}".format(x))
+ # f.write(" {0:6d}".format(y))
+ # f.write(" {0:6d}".format(z))
+ # f.write(" {:8.12}".format(self.grid.axes[0][x]))
+ # f.write(" {:8.12}".format(self.grid.axes[1][y]))
+ # f.write(" {:8.12}".format(self.grid.axes[2][z]))
+ # f.write(" {:8.12}".format(self.gvelo[x, y, z][0]))
+ # f.write(" {:8.12}".format(self.gvelo[x, y, z][1]))
+ # f.write(" {:8.12}".format(self.gvelo[x, y, z][2]))
+ # f.write(" {:8.12}".format(self.gscal[x, y, z]))
+ # f.write(" {:8.12}".format(self.ppos[x, y, z][0]))
+ # f.write(" {:8.12}".format(self.ppos[x, y, z][1]))
+ # f.write(" {:8.12}".format(self.ppos[x, y, z][2]))
+ # f.write(" {:8.12}".format(self.pscal[x, y, z]))
+ # f.write("\n")
+ # f.write("\n")
+ # f.close()
+ # self.ite += 1
if __name__ == "__main__":
diff --git a/HySoP/hysop/tools/timer.py b/HySoP/hysop/tools/timer.py
new file mode 100644
index 000000000..8d8c0ba9f
--- /dev/null
+++ b/HySoP/hysop/tools/timer.py
@@ -0,0 +1,43 @@
+# -*- coding: utf-8 -*-
+"""
+@package Operator
+Operator representation
+"""
+import time
+
+
+class Timer:
+ """
+ Timer. Manage time steps
+ DiscreteOperator.DiscreteOperator specialization.
+ """
+
+ def __init__(self, t_end, dt, t_init=0.):
+ """
+ Constructor.
+
+ """
+ self.t_end = t_end
+ self.dt = dt
+ self.t = t_init
+ self.end = False
+ self.ite = 0
+
+ def step(self):
+ """
+ Apply Remeshing operator.
+ """
+ if self.t < self.t_end:
+ self.ite += 1
+ self.t += self.dt
+ if self.t >= self.t_end:
+ self.end = True
+
+ def __str__(self):
+ s = "Timer (DiscreteOperator). "
+ return s
+
+if __name__ == "__main__":
+ print __doc__
+ print "- Provided class : RemeshingDOp"
+ print RemeshingDOp.__doc__
diff --git a/HySoP/setup.py.in b/HySoP/setup.py.in
index 2d2fc4387..2e931b6cd 100644
--- a/HySoP/setup.py.in
+++ b/HySoP/setup.py.in
@@ -5,7 +5,8 @@ setup.py file for @PYPACKAGE_NAME@
"""
from numpy.distutils.core import setup, Extension
-import os, glob
+import os
+import glob
import sys
# List of all directories for sources
#os.listdir(os.path.join('ParMePy','New_ParmePy'))
@@ -16,7 +17,10 @@ name = '@PYPACKAGE_NAME@'
packages = ['parmepy',
'parmepy.domain',
'parmepy.fields',
- 'parmepy.operator']
+ 'parmepy.operator',
+ 'parmepy.problem',
+ 'parmepy.ParticularSolvers',
+ 'parmepy.tools']
# 'examples']
# 'new_ParMePy',
# 'new_ParMePy/Domain',
@@ -28,15 +32,20 @@ packages = ['parmepy',
# Enable this to get debug info
# DISTUTILS_DEBUG=1
-parmes_dir=['@CMAKE_BINARY_DIR@/Modules']
-parmes_libdir=['@CMAKE_BINARY_DIR@/src']
-parmeslib=['@PARMES_LIBRARY_NAME@']
-f2py_options=['--no-lower']
+parmes_dir = ['@CMAKE_BINARY_DIR@/Modules']
+parmes_libdir = ['@CMAKE_BINARY_DIR@/src']
+parmeslib = ['@PARMES_LIBRARY_NAME@']
+f2py_options = ['--no-lower']
scales_src = glob.glob('@CMAKE_SOURCE_DIR@/parmepy/scales2py/*.f90')
-scalesModule=Extension(name='parmepy.scales2py',f2py_options=f2py_options,sources=scales_src,include_dirs=parmes_dir,
- library_dirs=parmes_libdir,libraries=parmeslib,define_macros = [('F2PY_REPORT_ON_ARRAY_COPY','1')])
+scalesModule = Extension(name='parmepy.scales2py',
+ f2py_options=f2py_options,
+ sources=scales_src,
+ include_dirs=parmes_dir,
+ library_dirs=parmes_libdir,
+ libraries=parmeslib,
+ define_macros=[('F2PY_REPORT_ON_ARRAY_COPY', '1')])
ext_modules = [scalesModule]
print ext_modules
@@ -47,8 +56,13 @@ withfftw = "@WITH_FFTW@"
if(withfftw is "ON"):
fortran_src = glob.glob('@CMAKE_SOURCE_DIR@/parmepy/parmesfftw2py/*.f90')
# sys.path.extend('config_fc -DF2PY_REPORT_ON_ARRAY_COPY'.split())
- parpyModule=Extension(name='parmepy.fftw2py',f2py_options=f2py_options,sources=fortran_src,include_dirs=parmes_dir,
- library_dirs=parmes_libdir,libraries=parmeslib,define_macros = [('F2PY_REPORT_ON_ARRAY_COPY','1')])
+ parpyModule = Extension(name='parmepy.fftw2py',
+ f2py_options=f2py_options,
+ sources=fortran_src,
+ include_dirs=parmes_dir,
+ library_dirs=parmes_libdir,
+ libraries=parmeslib,
+ define_macros=[('F2PY_REPORT_ON_ARRAY_COPY', '1')])
ext_modules.append(parpyModule)
print ext_modules
@@ -59,7 +73,7 @@ setup(name=name,
author_email='parmes-devel@lists.forge.imag.fr',
url='https://forge.imag.fr/projects/parmes/',
license='GNU public license',
- package_dir = {'': '@CMAKE_SOURCE_DIR@'},
+ package_dir={'': '@CMAKE_SOURCE_DIR@'},
ext_modules=ext_modules,
packages=packages
#data_files=[('new_ParMePy/Utils', ['./new_ParMePy/Utils/gpu_src.cl'])]
--
GitLab