From df65bb850e98c1f9fd30e963a58f454b103931d4 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Franck=20P=C3=A9rignon?= <franck.perignon@imag.fr>
Date: Wed, 27 Feb 2013 16:29:00 +0000
Subject: [PATCH] updates for new topology (not tested...) + pep8
---
HySoP/hysop/__init__.py.in | 1 +
HySoP/hysop/mpi/topology.py | 27 +++---
HySoP/hysop/operator/stretching.py | 6 +-
HySoP/hysop/operator/stretching_d.py | 106 ++++++++++++------------
HySoP/hysop/particular_solvers/basic.py | 48 +++++++----
5 files changed, 102 insertions(+), 86 deletions(-)
diff --git a/HySoP/hysop/__init__.py.in b/HySoP/hysop/__init__.py.in
index b243b260e..5f72aec27 100755
--- a/HySoP/hysop/__init__.py.in
+++ b/HySoP/hysop/__init__.py.in
@@ -10,6 +10,7 @@ __all__ = ['Box', 'CartesianTopology']#, 'ScalarField']
# MPI
if("@USE_MPI@" is "ON"):
+ #import parmepy.mpi as mpi
from mpi import topology
CartesianTopology = topology.CartesianTopology
if(mpi.main_rank == 0):
diff --git a/HySoP/hysop/mpi/topology.py b/HySoP/hysop/mpi/topology.py
index 27df9147b..571594a7a 100644
--- a/HySoP/hysop/mpi/topology.py
+++ b/HySoP/hysop/mpi/topology.py
@@ -387,40 +387,39 @@ class SubMesh(object):
## Topology that creates (and owns) this mesh
self._topology = topo
## Local resolution of this mesh, including ghost points
- self._resolution = resolution
+ self.resolution = resolution
## index of the lowest point of this mesh
## (in each dir) in the global mesh
- self._global_start = g_start
+ self.global_start = g_start
## index of the upper point (in each dir), global mesh
- self._global_end = self._global_start + self._resolution - 1
+ self.global_end = self.global_start + self.resolution - 1
## Mesh step size
- self._space_step_size = topo.domain.length / \
+ self.space_step_size = topo.domain.length / \
(topo.globalMeshResolution - 1)
## Mesh local indices, only for "computed" points
## (i.e. excluding ghosts)
- self._local_start = topo.ghosts.copy()
- self._local_end = self._resolution.copy() - topo.ghosts[:] - 1
+ self.local_start = topo.ghosts.copy()
+ self.local_end = self.resolution.copy() - topo.ghosts[:] - 1
## Coordinates of the "lowest" point of this mesh (including ghost)
## Warning FP : this strongly depends on where is the origin
## of the domain, of the type of boundary conditions
## and if origin is on ghosts or "real" points.
- self._origin = topo.domain.origin.copy()
- self._origin[:] += self._space_step_size * \
- (self._global_start[:] - topo.ghosts[:])
+ self.origin = topo.domain.origin.copy()
+ self.origin[:] += self.space_step_size * \
+ (self.global_start[:] - topo.ghosts[:])
def __str__(self):
""" Sub mesh display """
s = 'Coords:' + str(self._topology.coords[:])
- s += ' Sub mesh resolution:' + str(self._resolution) + '\n'
- s += 'Starting point global indices :' + str(self._global_start) + '\n'
+ s += ' Sub mesh resolution:' + str(self.resolution) + '\n'
+ s += 'Starting point global indices :' + str(self.global_start) + '\n'
s += 'Local indices for "computed" points (excluding ghosts) :\n'
- s += ' start = ' + str(self._local_start)
- s += ', end = ' + str(self._local_end) + '\n'
+ s += ' start = ' + str(self.local_start)
+ s += ', end = ' + str(self.local_end) + '\n'
return s
if __name__ == "__main__":
print "This module defines the following classes:"
print "- CartesianTopology: ", CartesianTopology.__doc__
- print "- FFTTopology: ", FFTTopology.__doc__
print "- LocalMesh: ", LocalMesh.__doc__
diff --git a/HySoP/hysop/operator/stretching.py b/HySoP/hysop/operator/stretching.py
index c7e7d64c5..835252cc8 100755
--- a/HySoP/hysop/operator/stretching.py
+++ b/HySoP/hysop/operator/stretching.py
@@ -4,9 +4,9 @@
Penalization operator representation
"""
-from .continuous import ContinuousOperator
-from .stretching_d import Stretching_d
-from ..particular_solvers.integrator.runge_kutta3 import RK3
+from continuous import ContinuousOperator
+from stretching_d import Stretching_d
+from parmepy.integrator import RK3
class Stretching(ContinuousOperator):
diff --git a/HySoP/hysop/operator/stretching_d.py b/HySoP/hysop/operator/stretching_d.py
index 0b5d957e2..bb9d28862 100755
--- a/HySoP/hysop/operator/stretching_d.py
+++ b/HySoP/hysop/operator/stretching_d.py
@@ -1,23 +1,15 @@
# -*- coding: utf-8 -*-
"""
-@package operator
+@package parmepy.operator.stretching_d
+
Discrete stretching representation
"""
-from ..constants import *
-from .discrete import DiscreteOperator
-from .differentialOperator import DifferentialOperator
-from .differentialOperator_d import DifferentialOperator_d
-from .fct2op import Fct2Op
-from ..particular_solvers.integrator.euler import Euler
-from ..particular_solvers.integrator.runge_kutta2 import RK2
-from ..particular_solvers.integrator.runge_kutta3 import RK3
-from ..particular_solvers.integrator.runge_kutta4 import RK4
-from ..particular_solvers.integrator.integrator import ODESolver
-from ..tools.cpu_data_transfer import Synchronize
-import numpy as np
-import numpy.testing as npt
+
+from discrete import DiscreteOperator
+from differentialOperator_d import DifferentialOperator_d
+from fct2op import Fct2Op
+import parmepy.integrator as integrator
import time
-import sys
class Stretching_d(DiscreteOperator):
@@ -26,7 +18,8 @@ class Stretching_d(DiscreteOperator):
DiscreteOperator.DiscreteOperator specialization.
"""
- def __init__(self, stretch, topology=None, idVelocityD=0, idVorticityD=0, propertyOp='divConservation', method=None):
+ def __init__(self, stretch, topology=None, idVelocityD=0, idVorticityD=0,
+ propertyOp='divConservation', method=None):
"""
Constructor.
Create a Stretching operator on a given continuous domain.
@@ -46,7 +39,6 @@ class Stretching_d(DiscreteOperator):
self.method = method
## self.name = "stretching"
self.topology = topology
- self.OpSynchronize = Synchronize(self.topology)
def apply(self, t, dt):
"""
@@ -56,8 +48,10 @@ class Stretching_d(DiscreteOperator):
@param dt : time step.
Stretching algorithm:
- @li 1. discretization of the divergence operator (div wu) or w grad(u) : \n
- @li 2. Time integration. Performs a method choose to resolove the equation dw/dt = div(wu) or dw/dt = w grad(u).
+ @li 1. discretization of the divergence operator (div wu)
+ or w grad(u) : \n
+ @li 2. Time integration. Performs a method choose to resolove
+ the equation dw/dt = div(wu) or dw/dt = w grad(u).
-define the function of ODEsolver as the result of the first step
-integrate with the chosen method (RK4. RK2. euler)
"""
@@ -65,43 +59,48 @@ class Stretching_d(DiscreteOperator):
if self.method is not None:
- if self.propertyOp == 'divConservation' :
+ if self.propertyOp == 'divConservation':
## Space discretisation of grad(u)
- self.stretchOp = DifferentialOperator_d(self.vorticity, self.velocity, choice='gradV', topology=self.topology)
+ self.stretchOp = DifferentialOperator_d(
+ self.vorticity, self.velocity, choice='gradV',
+ topology=self.topology)
+
gradResult, maxgersh = self.stretchOp.apply()
## Determination of Stretching time step (stability condition)
self.ndt = self.stabilityTest(dt, maxgersh[0], self.method)
- self.dtstretch = dt/float(self.ndt)
+ self.dtstretch = dt / float(self.ndt)
print "Maxgersh", maxgersh
print "dtStretch, ndt", self.dtstretch, self.ndt
- print "dtAdvec", 1./maxgersh[1]
+ print "dtAdvec", 0.5 / maxgersh[1]
## fct2Op
- self.fonction = Fct2Op(self.vorticity.data, gradResult,choice = 'gradV', topology=self.topology)
+ self.fonction = Fct2Op(self.vorticity.data,
+ gradResult, choice='gradV',
+ topology=self.topology)
- elif self.propertyOp == 'massConservation' :
+ elif self.propertyOp == 'massConservation':
## Determination of Stretching time step (stability condition)
self.ndt = 1
self.dtstretch = dt
## fct2Op
- self.fonction = Fct2Op(self.vorticity, self.velocity,choice = 'divWU', topology=self.topology)
+ self.fonction = Fct2Op(self.vorticity,
+ self.velocity, choice='divWU',
+ topology=self.topology)
else:
- raise ValueError("Unknown Stretching Operator property: " + str(self.propertyOp))
-
- # Time integration
- for i in range (self.ndt) :
- methodInt=self.method(self.fonction)
- self.vorticity.data = methodInt.integrate(self.fonction.apply, t, self.dtstretch, self.vorticity.data)
-
-
- # Ghosts synchronization
-# self.OpSynchronize.apply([self.vorticity.data, self.velocity.data])
+ raise ValueError("Unknown Stretching Operator property: "
+ + str(self.propertyOp))
+ # Time integration
+ for i in range(self.ndt):
+ methodInt = self.method(self.fonction)
+ self.vorticity.data = methodInt.integrate(self.fonction.apply,
+ t, self.dtstretch,
+ self.vorticity.data)
self.compute_time = time.time() - self.compute_time
self.total_time += self.compute_time
@@ -110,25 +109,27 @@ class Stretching_d(DiscreteOperator):
def stabilityTest(self, dt, maxi, method):
dtstretch = dt
- cststretch =0.
- if method == Euler :
+ cststretch = 0.
+ if method == integrator.Euler:
cststretch = 2.0
- dtstretch = min(dtstretch, cststretch/maxi)
- if method == RK2 :
+ dtstretch = min(dtstretch, cststretch / maxi)
+ if method == integrator.RK2:
cststretch = 2.0
- dtstretch = min(dtstretch, cststretch/maxi)
- if method == RK3 :
+ dtstretch = min(dtstretch, cststretch / maxi)
+ if method == integrator.RK3:
cststretch = 2.5127
- dtstretch = min(dtstretch, cststretch/maxi)
- if method == RK4 :
+ dtstretch = min(dtstretch, cststretch / maxi)
+ if method == integrator.RK4:
cststretch = 2.7853
- dtstretch = min(dtstretch, cststretch/maxi)
- if dt == dtstretch :
- print "dt, ndt , dtstretch, upperbound", dt, int(dt/dtstretch), dt/float(int(dt/dtstretch)), cststretch/maxi
- return int(dt/dtstretch)
- else :
- print "dt, ndt , dtstretch, upperbound", dt, int(dt/dtstretch)+1, dt/float(int(dt/dtstretch)+1), cststretch/maxi
- return int(dt/dtstretch) +1
+ dtstretch = min(dtstretch, cststretch / maxi)
+ if dt == dtstretch:
+ print "dt, ndt , dtstretch, upperbound", dt, int(dt / dtstretch),\
+ dt / float(int(dt / dtstretch)), cststretch / maxi
+ return int(dt / dtstretch)
+ else:
+ print "dt, ndt , dtstretch, upperbound", dt, int(dt / dtstretch)\
+ + 1, dt / float(int(dt / dtstretch) + 1), cststretch / maxi
+ return int(dt / dtstretch) + 1
def printComputeTime(self):
self.timings_info[0] = "\"Stretching total\""
@@ -137,7 +138,8 @@ class Stretching_d(DiscreteOperator):
print "Time of the last Stretching iteration :", self.compute_time
def __str__(self):
- s = "Stretching_d (DiscreteOperator). " + DiscreteOperator.__str__(self)
+ s = "Stretching_d (DiscreteOperator). " +\
+ DiscreteOperator.__str__(self)
return s
if __name__ == "__main__":
diff --git a/HySoP/hysop/particular_solvers/basic.py b/HySoP/hysop/particular_solvers/basic.py
index c18d9bc62..b793f481f 100644
--- a/HySoP/hysop/particular_solvers/basic.py
+++ b/HySoP/hysop/particular_solvers/basic.py
@@ -17,6 +17,7 @@ from parmepy.tools.timer import Timer
from parmepy.tools.printer import Printer
from parmepy.integrator import RK3
+
class ParticleSolver(Solver):
"""
Particular solver description.
@@ -46,7 +47,8 @@ class ParticleSolver(Solver):
Solver.__init__(self, problem)
## Problem to solve
self.problem = problem
- ## Advection operator of the problem. Advection need special treatment in particle methods.
+ ## Advection operator of the problem. Advection need special
+ ## treatment in particle methods.
self.advection = None
self.stretch = None
self.penal = None
@@ -87,19 +89,23 @@ class ParticleSolver(Solver):
#if isinstance(op, Velocity)
#self.velocity = op
#if self.advection is None:
- # raise ValueError("Particular Solver : Cannot create from a problem with no Transport operator")
+ # raise ValueError("Particular Solver :
+ # Cannot create from a problem with no Transport operator")
def initialize(self):
"""
Solver initialisation.
Initialize a particle method solver regarding the problem.
"""
- self.p_position = ContinuousField(domain=self.problem.domains[0], name='ParticlePosition')
- self.p_scalar = ContinuousField(domain=self.problem.domains[0], name='ParticleScalar')
+ self.p_position = ContinuousField(domain=self.problem.domains[0],
+ name='ParticlePosition')
+ self.p_scalar = ContinuousField(domain=self.problem.domains[0],
+ name='ParticleScalar')
self.problem.addVariable([self.p_position, self.p_scalar])
## Discretise domains
+ ## Note FP : this is probably useless???
for d in self.problem.domains:
- d.discretize(self.problem.topology.resolution)
+ d.discretize(self.problem.topology.globalMeshResolution)
## Discretise fields and initialize
for v in self.problem.variables:
v.discretize(self.problem.topology)
@@ -110,9 +116,12 @@ class ParticleSolver(Solver):
## 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)
+ op.discretize(result_position=self.p_position,
+ result_scalar=self.p_scalar,
+ method=self.ODESolver)
if op is self.stretch:
- op.discretize(topology=self.problem.topology, method=self.ODESolver)
+ op.discretize(topology=self.problem.topology,
+ method=self.ODESolver)
if op is self.penal:
op.obstacle.discretize(self.problem.topology)
op.discretize()
@@ -124,22 +133,27 @@ class ParticleSolver(Solver):
op.discretize(topology=self.problem.topology)
## Velocity is only program on gpu for instance
#elif op is self.velocity:
- # op.discretize(result_position=self.p_position, result_scalar=self.p_scalar, method=self.ODESolver)
+ # op.discretize(result_position=self.p_position,
+ # result_scalar=self.p_scalar, method=self.ODESolver)
else:
op.discretize()
if self.advection is not None and not self.advection.include_remesh:
## Create remeshing operator
- self.remeshing = Remeshing(partPositions=self.advection.discreteOperator.res_position,
- partScalar=self.advection.discreteOperator.res_scalar,
- resscal=self.advection.discreteOperator.scalar,
- method=self.RemeshingMethod)
+ self.remeshing = Remeshing(
+ partPositions=self.advection.discreteOperator.res_position,
+ partScalar=self.advection.discreteOperator.res_scalar,
+ resscal=self.advection.discreteOperator.scalar,
+ method=self.RemeshingMethod)
+
for op in self.problem.operators:
if op.needSplitting:
index = self.problem.operators.index(op)
if op is self.advection and not self.advection.include_remesh:
- self.problem.operators[index] = Splitting([op, self.remeshing], dim=self.problem.topology.dim)
+ self.problem.operators[index] = Splitting(
+ [op, self.remeshing], dim=self.problem.topology.dim)
else:
- self.problem.operators[index] = Splitting([op], dim=self.problem.topology.dim)
+ self.problem.operators[index] = Splitting(
+ [op], dim=self.problem.topology.dim)
self.isInitialized = True
def end(self):
@@ -147,10 +161,10 @@ class ParticleSolver(Solver):
def __str__(self):
"""ToString method"""
- s = " Particular solver "
+ s = " Particle solver "
return s
if __name__ == "__main__":
print __doc__
- print "- Provided class : ParticularSolver"
- print ParticularSolver.__doc__
+ print "- Provided class : ParticleSolver"
+ print ParticleSolver.__doc__
--
GitLab