diff --git a/Examples/NavierStokes3d.py b/Examples/NavierStokes3d.py
index df497b387f102c3aadbd73ac4a369ed47b88dded..5344afa60474d1d4e7cdeae0e198462b4519a4fc 100755
--- a/Examples/NavierStokes3d.py
+++ b/Examples/NavierStokes3d.py
@@ -136,6 +136,7 @@ penal = pp.Penalization(velocity, vorticity, sphere, lambd)
topofft = pp.CartesianTopology(domain=myDomain3d,
resolution=resolution3d, dim=3)
+##,ghost=2
navierStokes = pp.Problem(topofft, [penal, stretch])
@@ -167,6 +168,8 @@ omega_z = scales.solve_advection(timeStep, vx, vy, vz, omega_z)
# solve stretching
navierStokes.solve()
+print 'vorticity', vorticity.discreteField[0][0] , vorticity.discreteField[0][1] , vorticity.discreteField[0][2]
+
# solve diffusion
nudt = 0.0001
diff --git a/HySoP/hysop/domain/topology.py b/HySoP/hysop/domain/topology.py
index f74b923fe5f2b0f465ab4b925f4fb56e2e5be493..edd7163186ac547280bead41f70d0f981c086797 100644
--- a/HySoP/hysop/domain/topology.py
+++ b/HySoP/hysop/domain/topology.py
@@ -26,6 +26,7 @@ class CartesianTopology(object):
if(dims is None):
assert(dim is not None)
self.dims = np.asarray(MPI.Compute_dims(nbprocs, dim), dtype=PARMES_INTEGER)
+ print 'dims = ' , dims
# Topology dimension
self.dim = dim
else:
@@ -50,9 +51,14 @@ class CartesianTopology(object):
self.west, self.east = self.topo.Shift(YDIR, 1)
if(self.dim > 2):
self.south, self.north = self.topo.Shift(ZDIR, 1)
+ self.tabSort = np.array([[self.up,0] , [self.down,1] , [self.east,2], [self.west,3],
+ [self.north,4],[self.south,5]])
+ self.tabSort=self.tabSort[self.tabSort[:,0].argsort()]
# ghost layer (default = 0)
if(ghosts is None):
self.ghosts = np.zeros((self.domain.dimension), dtype=PARMES_INTEGER)
+ else:
+ self.ghosts=np.asarray(ghosts, dtype=PARMES_INTEGER)
# Global resolution
self.resolution = np.asarray(resolution, dtype=PARMES_INTEGER)
nbpoints = self.resolution[:self.dim] // self.dims
diff --git a/HySoP/hysop/operator/differentialOperator.py b/HySoP/hysop/operator/differentialOperator.py
index 46888ddefef1b33645e70869ab25fc4885ebe0de..076bcdb3b9c9d741dcb8667fe3d8f141d167940b 100755
--- a/HySoP/hysop/operator/differentialOperator.py
+++ b/HySoP/hysop/operator/differentialOperator.py
@@ -28,18 +28,17 @@ class DifferentialOperator(ContinuousOperator):
self.choice = choice
self.addVariable([field1, field2])
- def discretize(self, field1=None, field2=None, resolution=None, meshSize=None):
+ def discretize(self, field1=None, field2=None, topology=None):
"""
Stretching operator discretization method.
Create a discrete Stretching operator from given specifications.
@param field1 : Index of velocity discretisation to use.
@param field2 : Index of vorticity discretisation to use.
- @param resolution : local resolution.
- @param meshSize : mesh size.
+ @param topology : topology.
"""
if self.discreteOperator is None:
- self.discreteOperator = DifferentialOperator_d(field1, field2, self.choice, resolution, meshSize)
+ self.discreteOperator = DifferentialOperator_d(field1, field2, self.choice, topology)
def __str__(self):
"""ToString method"""
diff --git a/HySoP/hysop/operator/differentialOperator_d.py b/HySoP/hysop/operator/differentialOperator_d.py
index 8955cbf12088a5e3d44ded8a084f26908245cf78..0ad4bd0b9f38899713a7b88a52b8da2c22c561d0 100755
--- a/HySoP/hysop/operator/differentialOperator_d.py
+++ b/HySoP/hysop/operator/differentialOperator_d.py
@@ -1,10 +1,11 @@
-# -*- coding: utf-8 -*-
+# -*- codingind utf-8 -*-
"""
@package operator
Discrete stretching representation
"""
from ..constants import *
from .discrete import DiscreteOperator
+from ..tools.cpu_data_transfer import Synchronize
import numpy as np
import sys
import time
@@ -15,122 +16,135 @@ class DifferentialOperator_d(DiscreteOperator):
DiscreteOperator.DiscreteOperator specialization.
"""
- def __init__(self, field1, field2, choice=None, resolution=None, meshSize=None):
+ def __init__(self, field1, field2, choice=None, topology=None):
"""
Constructor.
Create a Stretching operator on a given continuous domain.
Work on a given field2 and field1 to compute the stretching term.
- @param stretch : Stretching operator.
+ @param stretch ind Stretching operator.
"""
DiscreteOperator.__init__(self)
+ ## Topology of the obstacle
+ if(topology is not None):
+ self.topology = topology
+ else:
+ raise NotImplementedError()
self.field1 = field1
self.field2 = field2
self.choice = choice
- self.resolution = resolution
- self.meshSize = meshSize
+ self.resolution = self.topology.mesh.resolution
+ self.meshSize = self.topology.mesh.size
self.compute_time = 0.
def apply(self):
"""
Apply operator.
"""
-# if False in [(self.field1[i][...].shape == self.field2[i][...].shape) for i in xrange(self.field2.topology.dim)] :
+# if False in [(self.field1[i][ind,ind].shape == self.field2[i][ind,ind].shape) for i in xrange(self.field2.topology.dim)] ind
# print "Error, not the same topology for field1 and field2"
# sys.exit(0)
-# else:
+# elseind
# mesh_size = self.field2.topology.mesh.size
# resolution = self.field2.topology.mesh.resolution
- ind=np.arange(self.resolution[0])
+ ind0a = self.topology.ghosts[0]
+ ind0b=self.resolution[0]-self.topology.ghosts[0]
+ ind1a = self.topology.ghosts[1]
+ ind1b=self.resolution[1]-self.topology.ghosts[1]
+ ind2a = self.topology.ghosts[2]
+ ind2b=self.resolution[2]-self.topology.ghosts[2]
+ ind0=np.arange(ind0a,ind0b)
+ ind1=np.arange(ind1a,ind1b)
+ ind2=np.arange(ind2a,ind2b)
if self.choice == 'divWU':
## second order scheme
## X components of temp and result
-# temp1 = ( self.field1[0][np.roll(ind,-1),...] * self.field2[0][np.roll(ind,-1),...] -\
-# self.field1[0][np.roll(ind,+1),...] * self.field2[0][np.roll(ind,+1),...]) / (2. * self.meshSize[0])
+# temp1 = ( self.field1[0][ind+1,ind,ind] * self.field2[0][ind+1,ind,ind] -\
+# self.field1[0][ind-1,ind,ind] * self.field2[0][ind-1,ind,ind]) / (2. * self.meshSize[0])
-# temp2 = (self.field1[1][:,np.roll(ind,-1),:] * self.field2[0][:,np.roll(ind,-1),:] -\
-# self.field1[1][:,np.roll(ind,+1),:] * self.field2[0][:,np.roll(ind,+1),:] ) / (2. * self.meshSize[1])
+# temp2 = (self.field1[1][ind,ind+1,ind] * self.field2[0][ind,ind+1,ind] -\
+# self.field1[1][ind,ind-1,ind] * self.field2[0][ind,ind-1,ind] ) / (2. * self.meshSize[1])
#
-# temp3 = (self.field1[2][...,np.roll(ind,-1)] * self.field2[0][...,np.roll(ind,-1)] -\
-# self.field1[2][...,np.roll(ind,+1)] * self.field2[0][...,np.roll(ind,+1)] ) / (2. * self.meshSize[2])
+# temp3 = (self.field1[2][ind,ind,ind+1] * self.field2[0][ind,ind,ind+1] -\
+# self.field1[2][ind,ind,ind-1] * self.field2[0][ind,ind,ind-1] ) / (2. * self.meshSize[2])
-# self.result[0][...]= temp1 + temp2 + temp3
+# self.result[0][ind,ind]= temp1 + temp2 + temp3
## Y components of temp and result
-# temp1 = (self.field1[0][np.roll(ind,-1),...] * self.field2[1][np.roll(ind,-1),...] -\
-# self.field1[0][np.roll(ind,+1),...] * self.field2[1][np.roll(ind,+1),...] ) / (2. * self.meshSize[0])
+# temp1 = (self.field1[0][ind+1,ind,ind] * self.field2[1][ind+1,ind,ind] -\
+# self.field1[0][ind-1,ind,ind] * self.field2[1][ind-1,ind,ind] ) / (2. * self.meshSize[0])
-# temp2 = (self.field1[1][:,np.roll(ind,-1),:] * self.field2[1][:,np.roll(ind,-1),:] -\
-# self.field1[1][:,np.roll(ind,+1),:] * self.field2[1][:,np.roll(ind,+1),:] ) / (2. * self.meshSize[1])
+# temp2 = (self.field1[1][ind,ind+1,ind] * self.field2[1][ind,ind+1,ind] -\
+# self.field1[1][ind,ind-1,ind] * self.field2[1][ind,ind-1,ind] ) / (2. * self.meshSize[1])
#
-# temp3 = (self.field1[2][...,np.roll(ind,-1)] * self.field2[1][...,np.roll(ind,-1)] -\
-# self.field1[2][...,np.roll(ind,+1)] * self.field2[1][...,np.roll(ind,+1)] ) / (2. * self.meshSize[2])
+# temp3 = (self.field1[2][ind,ind,ind+1] * self.field2[1][ind,ind,ind+1] -\
+# self.field1[2][ind,ind,ind-1] * self.field2[1][ind,ind,ind-1] ) / (2. * self.meshSize[2])
-# self.result[1][...]= temp1 + temp2 +temp3
+# self.result[1][ind,ind]= temp1 + temp2 +temp3
## Z components of temp and result
-# temp1 = (self.field1[0][np.roll(ind,-1),...] * self.field2[2][np.roll(ind,-1),...] -\
-# self.field1[0][np.roll(ind,+1),...] * self.field2[2][np.roll(ind,+1),...] ) / (2. * self.meshSize[0])
+# temp1 = (self.field1[0][ind+1,ind,ind] * self.field2[2][ind+1,ind,ind] -\
+# self.field1[0][ind-1,ind,ind] * self.field2[2][ind-1,ind,ind] ) / (2. * self.meshSize[0])
-# temp2 = (self.field1[1][:,np.roll(ind,-1),:] * self.field2[2][:,np.roll(ind,-1),:] -\
-# self.field1[1][:,np.roll(ind,+1),:] * self.field2[2][:,np.roll(ind,+1),:] ) / (2. * self.meshSize[1])
+# temp2 = (self.field1[1][ind,ind+1,ind] * self.field2[2][ind,ind+1,ind] -\
+# self.field1[1][ind,ind-1,ind] * self.field2[2][ind,ind-1,ind] ) / (2. * self.meshSize[1])
#
-# temp3 = (self.field1[2][...,np.roll(ind,-1)] * self.field2[2][...,np.roll(ind,-1)] -\
-# self.field1[2][...,np.roll(ind,+1)] * self.field2[2][...,np.roll(ind,+1)]) / (2. * self.meshSize[2])
+# temp3 = (self.field1[2][ind,ind,ind+1] * self.field2[2][ind,ind,ind+1] -\
+# self.field1[2][ind,ind,ind-1] * self.field2[2][ind,ind,ind-1]) / (2. * self.meshSize[2])
# X components of temp and result
- temp1 = (1.0 * self.field1[0][np.roll(ind,+2),...] * self.field2[0][np.roll(ind,+2),...] -
- 8.0 * self.field1[0][np.roll(ind,+1),...] * self.field2[0][np.roll(ind,+1),...] +\
- 8.0 * self.field1[0][np.roll(ind,-1),...] * self.field2[0][np.roll(ind,-1),...] -\
- 1.0 * self.field1[0][np.roll(ind,-2),...] * self.field2[0][np.roll(ind,-2),...]) / (12. * self.meshSize[0])
-
- temp2 = (1.0 * self.field1[1][:,np.roll(ind,+2),:] * self.field2[0][:,np.roll(ind,+2),:] -\
- 8.0 * self.field1[1][:,np.roll(ind,+1),:] * self.field2[0][:,np.roll(ind,+1),:] +\
- 8.0 * self.field1[1][:,np.roll(ind,-1),:] * self.field2[0][:,np.roll(ind,-1),:] -\
- 1.0 * self.field1[1][:,np.roll(ind,-2),:] * self.field2[0][:,np.roll(ind,-2),:]) / (12. * self.meshSize[1])
+ temp1 = (1.0 * self.field1[0][ind-2,ind,ind] * self.field2[0][ind-2,ind,ind] -
+ 8.0 * self.field1[0][ind-1,ind,ind] * self.field2[0][ind-1,ind,ind] +\
+ 8.0 * self.field1[0][ind+1,ind,ind] * self.field2[0][ind+1,ind,ind] -\
+ 1.0 * self.field1[0][ind+2,ind,ind] * self.field2[0][ind+2,ind,ind]) / (12. * self.meshSize[0])
+
+ temp2 = (1.0 * self.field1[1][ind,ind-2,ind] * self.field2[0][ind,ind-2,ind] -\
+ 8.0 * self.field1[1][ind,ind-1,ind] * self.field2[0][ind,ind-1,ind] +\
+ 8.0 * self.field1[1][ind,ind+1,ind] * self.field2[0][ind,ind+1,ind] -\
+ 1.0 * self.field1[1][ind,ind+2,ind] * self.field2[0][ind,ind+2,ind]) / (12. * self.meshSize[1])
- temp3 = (1.0 * self.field1[2][...,np.roll(ind,+2)] * self.field2[0][...,np.roll(ind,+2)] -\
- 8.0 * self.field1[2][...,np.roll(ind,+1)] * self.field2[0][...,np.roll(ind,+1)] +\
- 8.0 * self.field1[2][...,np.roll(ind,-1)] * self.field2[0][...,np.roll(ind,-1)] -\
- 1.0 * self.field1[2][...,np.roll(ind,-2)] * self.field2[0][...,np.roll(ind,-2)]) / (12. * self.meshSize[2])
+ temp3 = (1.0 * self.field1[2][ind,ind,ind-2] * self.field2[0][ind,ind,ind-2] -\
+ 8.0 * self.field1[2][ind,ind,ind-1] * self.field2[0][ind,ind,ind-1] +\
+ 8.0 * self.field1[2][ind,ind,ind+1] * self.field2[0][ind,ind,ind+1] -\
+ 1.0 * self.field1[2][ind,ind,ind+2] * self.field2[0][ind,ind,ind+2]) / (12. * self.meshSize[2])
tmp1 = np.array([temp1 + temp2 + temp3])
# Y components of temp and result
- temp1 = (1.0 * self.field1[0][np.roll(ind,+2),...] * self.field2[1][np.roll(ind,+2),...] -
- 8.0 * self.field1[0][np.roll(ind,+1),...] * self.field2[1][np.roll(ind,+1),...] +\
- 8.0 * self.field1[0][np.roll(ind,-1),...] * self.field2[1][np.roll(ind,-1),...] -\
- 1.0 * self.field1[0][np.roll(ind,-2),...] * self.field2[1][np.roll(ind,-2),...]) / (12. * self.meshSize[0])
-
- temp2 = (1.0 * self.field1[1][:,np.roll(ind,+2),:] * self.field2[1][:,np.roll(ind,+2),:] -\
- 8.0 * self.field1[1][:,np.roll(ind,+1),:] * self.field2[1][:,np.roll(ind,+1),:] +\
- 8.0 * self.field1[1][:,np.roll(ind,-1),:] * self.field2[1][:,np.roll(ind,-1),:] -\
- 1.0 * self.field1[1][:,np.roll(ind,-2),:] * self.field2[1][:,np.roll(ind,-2),:]) / (12. * self.meshSize[1])
+ temp1 = (1.0 * self.field1[0][ind-2,ind,ind] * self.field2[1][ind-2,ind,ind] -
+ 8.0 * self.field1[0][ind-1,ind,ind] * self.field2[1][ind-1,ind,ind] +\
+ 8.0 * self.field1[0][ind+1,ind,ind] * self.field2[1][ind+1,ind,ind] -\
+ 1.0 * self.field1[0][ind+2,ind,ind] * self.field2[1][ind+2,ind,ind]) / (12. * self.meshSize[0])
+
+ temp2 = (1.0 * self.field1[1][ind,ind-2,ind] * self.field2[1][ind,ind-2,ind] -\
+ 8.0 * self.field1[1][ind,ind-1,ind] * self.field2[1][ind,ind-1,ind] +\
+ 8.0 * self.field1[1][ind,ind+1,ind] * self.field2[1][ind,ind+1,ind] -\
+ 1.0 * self.field1[1][ind,ind+2,ind] * self.field2[1][ind,ind+2,ind]) / (12. * self.meshSize[1])
- temp3 = (1.0 * self.field1[2][...,np.roll(ind,+2)] * self.field2[1][...,np.roll(ind,+2)] -\
- 8.0 * self.field1[2][...,np.roll(ind,+1)] * self.field2[1][...,np.roll(ind,+1)] +\
- 8.0 * self.field1[2][...,np.roll(ind,-1)] * self.field2[1][...,np.roll(ind,-1)] -\
- 1.0 * self.field1[2][...,np.roll(ind,-2)] * self.field2[1][...,np.roll(ind,-2)]) / (12. * self.meshSize[2])
+ temp3 = (1.0 * self.field1[2][ind,ind,ind-2] * self.field2[1][ind,ind,ind-2] -\
+ 8.0 * self.field1[2][ind,ind,ind-1] * self.field2[1][ind,ind,ind-1] +\
+ 8.0 * self.field1[2][ind,ind,ind+1] * self.field2[1][ind,ind,ind+1] -\
+ 1.0 * self.field1[2][ind,ind,ind+2] * self.field2[1][ind,ind,ind+2]) / (12. * self.meshSize[2])
tmp2 = np.array([temp1 + temp2 + temp3])
# Z components of temp and result
- temp1 = (1.0 * self.field1[0][np.roll(ind,+2),...] * self.field2[2][np.roll(ind,+2),...] -
- 8.0 * self.field1[0][np.roll(ind,+1),...] * self.field2[2][np.roll(ind,+1),...] +\
- 8.0 * self.field1[0][np.roll(ind,-1),...] * self.field2[2][np.roll(ind,-1),...] -\
- 1.0 * self.field1[0][np.roll(ind,-2),...] * self.field2[2][np.roll(ind,-2),...]) / (12. * self.meshSize[0])
-
- temp2 = (1.0 * self.field1[1][:,np.roll(ind,+2),:] * self.field2[2][:,np.roll(ind,+2),:] -\
- 8.0 * self.field1[1][:,np.roll(ind,+1),:] * self.field2[2][:,np.roll(ind,+1),:] +\
- 8.0 * self.field1[1][:,np.roll(ind,-1),:] * self.field2[2][:,np.roll(ind,-1),:] -\
- 1.0 * self.field1[1][:,np.roll(ind,-2),:] * self.field2[2][:,np.roll(ind,-2),:]) / (12. * self.meshSize[1])
+ temp1 = (1.0 * self.field1[0][ind-2,ind,ind] * self.field2[2][ind-2,ind,ind] -
+ 8.0 * self.field1[0][ind-1,ind,ind] * self.field2[2][ind-1,ind,ind] +\
+ 8.0 * self.field1[0][ind+1,ind,ind] * self.field2[2][ind+1,ind,ind] -\
+ 1.0 * self.field1[0][ind+2,ind,ind] * self.field2[2][ind+2,ind,ind]) / (12. * self.meshSize[0])
+
+ temp2 = (1.0 * self.field1[1][ind,ind-2,ind] * self.field2[2][ind,ind-2,ind] -\
+ 8.0 * self.field1[1][ind,ind-1,ind] * self.field2[2][ind,ind-1,ind] +\
+ 8.0 * self.field1[1][ind,ind+1,ind] * self.field2[2][ind,ind+1,ind] -\
+ 1.0 * self.field1[1][ind,ind+2,ind] * self.field2[2][ind,ind+2,ind]) / (12. * self.meshSize[1])
- temp3 = (1.0 * self.field1[2][...,np.roll(ind,+2)] * self.field2[2][...,np.roll(ind,+2)] -\
- 8.0 * self.field1[2][...,np.roll(ind,+1)] * self.field2[2][...,np.roll(ind,+1)] +\
- 8.0 * self.field1[2][...,np.roll(ind,-1)] * self.field2[2][...,np.roll(ind,-1)] -\
- 1.0 * self.field1[2][...,np.roll(ind,-2)] * self.field2[2][...,np.roll(ind,-2)]) / (12. * self.meshSize[2])
+ temp3 = (1.0 * self.field1[2][ind,ind,ind-2] * self.field2[2][ind,ind,ind-2] -\
+ 8.0 * self.field1[2][ind,ind,ind-1] * self.field2[2][ind,ind,ind-1] +\
+ 8.0 * self.field1[2][ind,ind,ind+1] * self.field2[2][ind,ind,ind+1] -\
+ 1.0 * self.field1[2][ind,ind,ind+2] * self.field2[2][ind,ind,ind+2]) / (12. * self.meshSize[2])
tmp3 = np.array([temp1 + temp2 + temp3])
result = np.concatenate((np.array([tmp1]), np.array([tmp2]), np.array([tmp3])))
@@ -138,118 +152,134 @@ class DifferentialOperator_d(DiscreteOperator):
elif self.choice == 'gradU':
+ # Ghosts synchronization
+ OpSynchronize = Synchronize(self.topology)
+ print 'Fields avant', self.field2
+ OpSynchronize.apply(self.field2)
+
+
maxgersh = np.zeros(2, dtype=PARMES_REAL, order=ORDER)
+
## Fourth order scheme
# X components of temp and result
- temp1 = (1.0 * self.field2[0][np.roll(ind,+2),...] -
- 8.0 * self.field2[0][np.roll(ind,+1),...] +\
- 8.0 * self.field2[0][np.roll(ind,-1),...] -\
- 1.0 * self.field2[0][np.roll(ind,-2),...]) / (12. * self.meshSize[0])
-
- temp2 = (1.0 * self.field2[0][:,np.roll(ind,+2),:] -\
- 8.0 * self.field2[0][:,np.roll(ind,+1),:] +\
- 8.0 * self.field2[0][:,np.roll(ind,-1),:] -\
- 1.0 * self.field2[0][:,np.roll(ind,-2),:]) / (12. * self.meshSize[1])
+ temp1 = self.field2[0][...]*0.
+ temp1[ind0a:ind0b,ind1a:ind1b,ind2a:ind2b] = (1.0 * self.field2[0][ind0-2,ind1a:ind1b,ind2a:ind2b] -
+ 8.0 * self.field2[0][ind0-1,ind1a:ind1b,ind2a:ind2b] +\
+ 8.0 * self.field2[0][ind0+1,ind1a:ind1b,ind2a:ind2b] -\
+ 1.0 * self.field2[0][ind0+2,ind1a:ind1b,ind2a:ind2b]) / (12. * self.meshSize[0])
+
+ temp2 = self.field2[0][...]*0.
+ temp2[ind0a:ind0b,ind1a:ind1b,ind2a:ind2b] = (1.0 * self.field2[0][ind0a:ind0b,ind1-2,ind2a:ind2b] -\
+ 8.0 * self.field2[0][ind0a:ind0b,ind1-1,ind2a:ind2b] +\
+ 8.0 * self.field2[0][ind0a:ind0b,ind1+1,ind2a:ind2b] -\
+ 1.0 * self.field2[0][ind0a:ind0b,ind1+2,ind2a:ind2b]) / (12. * self.meshSize[1])
- temp3 = (1.0 * self.field2[0][...,np.roll(ind,+2)] -\
- 8.0 * self.field2[0][...,np.roll(ind,+1)] +\
- 8.0 * self.field2[0][...,np.roll(ind,-1)] -\
- 1.0 * self.field2[0][...,np.roll(ind,-2)]) / (12. * self.meshSize[2])
+ temp3 = self.field2[0][...]*0.
+ temp3[ind0a:ind0b,ind1a:ind1b,ind2a:ind2b] = (1.0 * self.field2[0][ind0a:ind0b,ind1a:ind1b,ind2-2] -\
+ 8.0 * self.field2[0][ind0a:ind0b,ind1a:ind1b,ind2-1] +\
+ 8.0 * self.field2[0][ind0a:ind0b,ind1a:ind1b,ind2+1] -\
+ 1.0 * self.field2[0][ind0a:ind0b,ind1a:ind1b,ind2+2]) / (12. * self.meshSize[2])
maxstr1= np.max(abs(temp1)+abs(temp2)+abs(temp3))
maxadv1= np.max(abs(temp1))
# Y components of temp and result
- temp4 = (1.0 * self.field2[1][np.roll(ind,+2),...] -
- 8.0 * self.field2[1][np.roll(ind,+1),...] +\
- 8.0 * self.field2[1][np.roll(ind,-1),...] -\
- 1.0 * self.field2[1][np.roll(ind,-2),...]) / (12. * self.meshSize[0])
-
- temp5 = (1.0 * self.field2[1][:,np.roll(ind,+2),:] -\
- 8.0 * self.field2[1][:,np.roll(ind,+1),:] +\
- 8.0 * self.field2[1][:,np.roll(ind,-1),:] -\
- 1.0 * self.field2[1][:,np.roll(ind,-2),:]) / (12. * self.meshSize[1])
+ temp4 = self.field2[1][...]*0.
+ temp4[ind0a:ind0b,ind1a:ind1b,ind2a:ind2b] = (1.0 * self.field2[1][ind0-2,ind1a:ind1b,ind2a:ind2b] -
+ 8.0 * self.field2[1][ind0-1,ind1a:ind1b,ind2a:ind2b] +\
+ 8.0 * self.field2[1][ind0+1,ind1a:ind1b,ind2a:ind2b] -\
+ 1.0 * self.field2[1][ind0+2,ind1a:ind1b,ind2a:ind2b]) / (12. * self.meshSize[0])
+
+ temp5 = self.field2[1][...]*0.
+ temp5[ind0a:ind0b,ind1a:ind1b,ind2a:ind2b] = (1.0 * self.field2[1][ind0a:ind0b,ind1-2,ind2a:ind2b] -\
+ 8.0 * self.field2[1][ind0a:ind0b,ind1-1,ind2a:ind2b] +\
+ 8.0 * self.field2[1][ind0a:ind0b,ind1+1,ind2a:ind2b] -\
+ 1.0 * self.field2[1][ind0a:ind0b,ind1+2,ind2a:ind2b]) / (12. * self.meshSize[1])
- temp6 = (1.0 * self.field2[1][...,np.roll(ind,+2)] -\
- 8.0 * self.field2[1][...,np.roll(ind,+1)] +\
- 8.0 * self.field2[1][...,np.roll(ind,-1)] -\
- 1.0 * self.field2[1][...,np.roll(ind,-2)]) / (12. * self.meshSize[2])
+ temp6 = self.field2[1][...]*0.
+ temp6[ind0a:ind0b,ind1a:ind1b,ind2a:ind2b] = (1.0 * self.field2[1][ind0a:ind0b,ind1a:ind1b,ind2-2] -\
+ 8.0 * self.field2[1][ind0a:ind0b,ind1a:ind1b,ind2-1] +\
+ 8.0 * self.field2[1][ind0a:ind0b,ind1a:ind1b,ind2+1] -\
+ 1.0 * self.field2[1][ind0a:ind0b,ind1a:ind1b,ind2+2]) / (12. * self.meshSize[2])
maxstr2= np.max(abs(temp1)+abs(temp2)+abs(temp3))
maxadv2= np.max(abs(temp2))
# Z components of temp and result
- temp7 = (1.0 * self.field2[2][np.roll(ind,+2),...] -
- 8.0 * self.field2[2][np.roll(ind,+1),...] +\
- 8.0 * self.field2[2][np.roll(ind,-1),...] -\
- 1.0 * self.field2[2][np.roll(ind,-2),...]) / (12. * self.meshSize[0])
-
- temp8 = (1.0 * self.field2[2][:,np.roll(ind,+2),:] -\
- 8.0 * self.field2[2][:,np.roll(ind,+1),:] +\
- 8.0 * self.field2[2][:,np.roll(ind,-1),:] -\
- 1.0 * self.field2[2][:,np.roll(ind,-2),:]) / (12. * self.meshSize[1])
+ temp7 = self.field2[2][...]*0.
+ temp7[ind0a:ind0b,ind1a:ind1b,ind2a:ind2b] = (1.0 * self.field2[2][ind0-2,ind1a:ind1b,ind2a:ind2b] -
+ 8.0 * self.field2[2][ind0-1,ind1a:ind1b,ind2a:ind2b] +\
+ 8.0 * self.field2[2][ind0+1,ind1a:ind1b,ind2a:ind2b] -\
+ 1.0 * self.field2[2][ind0+2,ind1a:ind1b,ind2a:ind2b]) / (12. * self.meshSize[0])
+
+ temp8 = self.field2[2][...]*0.
+ temp8[ind0a:ind0b,ind1a:ind1b,ind2a:ind2b] = (1.0 * self.field2[2][ind0a:ind0b,ind1-2,ind2a:ind2b] -\
+ 8.0 * self.field2[2][ind0a:ind0b,ind1-1,ind2a:ind2b] +\
+ 8.0 * self.field2[2][ind0a:ind0b,ind1+1,ind2a:ind2b] -\
+ 1.0 * self.field2[2][ind0a:ind0b,ind1+2,ind2a:ind2b]) / (12. * self.meshSize[1])
- temp9 = (1.0 * self.field2[2][...,np.roll(ind,+2)] -\
- 8.0 * self.field2[2][...,np.roll(ind,+1)] +\
- 8.0 * self.field2[2][...,np.roll(ind,-1)] -\
- 1.0 * self.field2[2][...,np.roll(ind,-2)]) / (12. * self.meshSize[2])
+ temp9 = self.field2[2][...]*0.
+ temp9[ind0a:ind0b,ind1a:ind1b,ind2a:ind2b] = (1.0 * self.field2[2][ind0a:ind0b,ind1a:ind1b,ind2-2] -\
+ 8.0 * self.field2[2][ind0a:ind0b,ind1a:ind1b,ind2-1] +\
+ 8.0 * self.field2[2][ind0a:ind0b,ind1a:ind1b,ind2+1] -\
+ 1.0 * self.field2[2][ind0a:ind0b,ind1a:ind1b,ind2+2]) / (12. * self.meshSize[2])
maxstr3= np.max(abs(temp1)+abs(temp2)+abs(temp3))
maxadv3= np.max(abs(temp3))
maxgersh[0] = max(maxstr1,maxstr2,maxstr3)
maxgersh[1] = max(maxadv1,maxadv2,maxadv3)
result = np.concatenate((np.array([temp1, temp2, temp3]),np.array([temp4, temp5, temp6]),np.array([temp7, temp8, temp9])))
+
return result, maxgersh
elif self.choice == 'curl':
- temp1 = (1.0 * self.field1[2][:,np.roll(ind,+2),:] -
- 8.0 * self.field1[2][:,np.roll(ind,+1),:] +\
- 8.0 * self.field1[2][:,np.roll(ind,-1),:] -\
- 1.0 * self.field1[2][:,np.roll(ind,-2),:] ) / (12. * self.meshSize[1])
+ temp1 = (1.0 * self.field1[2][ind0a:ind0b,ind1-2,ind2a:ind2b] -
+ 8.0 * self.field1[2][ind0a:ind0b,ind1-1,ind2a:ind2b] +\
+ 8.0 * self.field1[2][ind0a:ind0b,ind1+1,ind2a:ind2b] -\
+ 1.0 * self.field1[2][ind0a:ind0b,ind1+2,ind2a:ind2b] ) / (12. * self.meshSize[1])
- temp2 = (1.0 * self.field1[1][...,np.roll(ind,+2)] -\
- 8.0 * self.field1[1][...,np.roll(ind,+1)] +\
- 8.0 * self.field1[1][...,np.roll(ind,-1)] -\
- 1.0 * self.field1[1][...,np.roll(ind,-2)] ) / (12. * self.meshSize[2])
+ temp2 = (1.0 * self.field1[1][ind0a:ind0b,ind1a:ind1b,ind2-2] -\
+ 8.0 * self.field1[1][ind0a:ind0b,ind1a:ind1b,ind2-1] +\
+ 8.0 * self.field1[1][ind0a:ind0b,ind1a:ind1b,ind2+1] -\
+ 1.0 * self.field1[1][ind0a:ind0b,ind1a:ind1b,ind2+2] ) / (12. * self.meshSize[2])
tmp1 = temp1 - temp2
- temp1 = (1.0 * self.field1[0][...,np.roll(ind,+2)] -\
- 8.0 * self.field1[0][...,np.roll(ind,+1)] +\
- 8.0 * self.field1[0][...,np.roll(ind,-1)] -\
- 1.0 * self.field1[0][...,np.roll(ind,-2)] ) / (12. * self.meshSize[2])
+ temp1 = (1.0 * self.field1[0][ind0a:ind0b,ind1a:ind1b,ind2-2] -\
+ 8.0 * self.field1[0][ind0a:ind0b,ind1a:ind1b,ind2-1] +\
+ 8.0 * self.field1[0][ind0a:ind0b,ind1a:ind1b,ind2+1] -\
+ 1.0 * self.field1[0][ind0a:ind0b,ind1a:ind1b,ind2+2] ) / (12. * self.meshSize[2])
- temp2 = (1.0 * self.field1[2][np.roll(ind,+2),...] -
- 8.0 * self.field1[2][np.roll(ind,+1),...] +\
- 8.0 * self.field1[2][np.roll(ind,-1),...] -\
- 1.0 * self.field1[2][np.roll(ind,-2),...] ) / (12. * self.meshSize[0])
+ temp2 = (1.0 * self.field1[2][ind0-2,ind1a:ind1b,ind2a:ind2b] -
+ 8.0 * self.field1[2][ind0-1,ind1a:ind1b,ind2a:ind2b] +\
+ 8.0 * self.field1[2][ind0+1,ind1a:ind1b,ind2a:ind2b] -\
+ 1.0 * self.field1[2][ind0+2,ind1a:ind1b,ind2a:ind2b] ) / (12. * self.meshSize[0])
tmp2 = temp1 - temp2
- temp1 = (1.0 * self.field1[1][np.roll(ind,+2),...] -
- 8.0 * self.field1[1][np.roll(ind,+1),...] +\
- 8.0 * self.field1[1][np.roll(ind,-1),...] -\
- 1.0 * self.field1[1][np.roll(ind,-2),...] ) / (12. * self.meshSize[0])
+ temp1 = (1.0 * self.field1[1][ind-2,ind1a:ind1b,ind2a:ind2b] -
+ 8.0 * self.field1[1][ind-1,ind1a:ind1b,ind2a:ind2b] +\
+ 8.0 * self.field1[1][ind+1,ind1a:ind1b,ind2a:ind2b] -\
+ 1.0 * self.field1[1][ind+2,ind1a:ind1b,ind2a:ind2b] ) / (12. * self.meshSize[0])
- temp2 = (1.0 * self.field1[0][:,np.roll(ind,+2),:] -\
- 8.0 * self.field1[0][:,np.roll(ind,+1),:] +\
- 8.0 * self.field1[0][:,np.roll(ind,-1),:] -\
- 1.0 * self.field1[0][:,np.roll(ind,-2),:] ) / (12. * self.meshSize[1])
+ temp2 = (1.0 * self.field1[0][ind0a:ind0b,ind1-2,ind2a:ind2b] -\
+ 8.0 * self.field1[0][ind0a:ind0b,ind1-1,ind2a:ind2b] +\
+ 8.0 * self.field1[0][ind0a:ind0b,ind1+1,ind2a:ind2b] -\
+ 1.0 * self.field1[0][ind0a:ind0b,ind1+2,ind2a:ind2b] ) / (12. * self.meshSize[1])
tmp3 = temp1 - temp2
result = np.concatenate((np.array([tmp1]), np.array([tmp2]), np.array([tmp3])))
- print 'curly' , result
return result
else:
- raise ValueError("Unknown differiential operator: " + str(self.choice))
+ raise ValueError("Unknown differiential operatorind " + str(self.choice))
def printComputeTime(self):
self.timings_info[0] = "\"Differential Operator total\""
self.timings_info[1] = str(self.total_time)
self.timings_info[1] += str(self.compute_time[0]) + " " + str(self.compute_time[1]) + " " + str(self.compute_time[2]) + " "
- print "Differential Operator total time : ", self.total_time
- print "\t Differential Operator :", self.compute_time
+ print "Differential Operator total time ind ", self.total_time
+ print "\t Differential Operator ind", self.compute_time
def __str__(self):
@@ -258,5 +288,5 @@ class DifferentialOperator_d(DiscreteOperator):
if __name__ == "__main__":
print __doc__
- print "- Provided class : DifferentialOperator_d"
+ print "- Provided class ind DifferentialOperator_d"
print DifferentialOperator_d.__doc__
diff --git a/HySoP/hysop/operator/fct2op.py b/HySoP/hysop/operator/fct2op.py
index 485c57d9789250b0b9fe6cf0629dba33a1624976..fe085a5e5dea5e4950c93b64bf6db7565f58d863 100644
--- a/HySoP/hysop/operator/fct2op.py
+++ b/HySoP/hysop/operator/fct2op.py
@@ -15,13 +15,13 @@ class Fct2Op(object):
make the link between differentialOperator_d and the integration time scheme (Euler,RK,...)
"""
- def __init__(self, field1, field2, choice=None, resolution = None , meshSize = None):
+ def __init__(self, field1, field2, choice=None, topology=None):
"""
Constructor.
store field2 to work with it in the apply function
@param field1 : vorticity field
- @param field2 : velocity field
+ @param field2 : velocity field or grad(U) depending on 'choice' parameter
@param choice : can be 'gradU' for w.Du or divWU for div(wu)
@param resolution : resolution of the integration domain (array of number of points)
@param meshSize : array of the space step
@@ -29,8 +29,9 @@ class Fct2Op(object):
self.field1 = field1
self.field2 = field2
self.choice = choice
- self.resolution = resolution
- self.meshSize = meshSize
+ self.topology = topology
+ self.resolution = self.topology.mesh.resolution
+ self.meshSize = self.topology.mesh.size
def apply(self, t , y):
"""
@@ -46,7 +47,8 @@ class Fct2Op(object):
result = np.concatenate((np.array([temp0]),np.array([temp1]),np.array([temp2])))
return result
if self.choice == 'divWU' :
- result = DifferentialOperator_d(y, self.field2, self.choice, self.resolution , self.meshSize).apply()
+ print 'taille field22:', self.field2.shape
+ result = DifferentialOperator_d(y, self.field2, self.choice, self.topology).apply()
return result
else :
print 'choice', choice , 'is not define'
diff --git a/HySoP/hysop/operator/stretching.py b/HySoP/hysop/operator/stretching.py
index 55d8fc1105b23663b8e157c61040224039937c74..c7e7d64c595cb45a6766e2eabfa1ffe41a829be1 100755
--- a/HySoP/hysop/operator/stretching.py
+++ b/HySoP/hysop/operator/stretching.py
@@ -29,7 +29,7 @@ class Stretching(ContinuousOperator):
self.vorticity = vorticity
self.addVariable([velocity, vorticity])
- def discretize(self, idVelocityD=0, idVorticityD=0, propertyOp='divConservation', method=RK3):
+ def discretize(self, topology=None, idVelocityD=0, idVorticityD=0, propertyOp='divConservation', method=RK3):
"""
Stretching operator discretization method.
Create a discrete Stretching operator from given specifications.
@@ -40,7 +40,7 @@ class Stretching(ContinuousOperator):
@param method : the method to use. (Euler, RK, ..) default : RK3
"""
if self.discreteOperator is None:
- self.discreteOperator = Stretching_d(self, idVelocityD, idVorticityD, propertyOp, method)
+ self.discreteOperator = Stretching_d(self, topology, idVelocityD, idVorticityD, propertyOp, method)
def __str__(self):
"""ToString method"""
diff --git a/HySoP/hysop/operator/stretching_d.py b/HySoP/hysop/operator/stretching_d.py
index ada0a269cd833715331aa7d9c993011a40948afb..3e107c5d34798c2267a9c1b497692fd32f9b7449 100755
--- a/HySoP/hysop/operator/stretching_d.py
+++ b/HySoP/hysop/operator/stretching_d.py
@@ -13,6 +13,7 @@ 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
import time
@@ -25,7 +26,7 @@ class Stretching_d(DiscreteOperator):
DiscreteOperator.DiscreteOperator specialization.
"""
- def __init__(self, stretch, 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.
@@ -44,6 +45,8 @@ class Stretching_d(DiscreteOperator):
self.propertyOp = propertyOp
self.method = method
## self.name = "stretching"
+ self.topology = topology
+ self.OpSynchronize = Synchronize(self.topology)
def apply(self, t, dt):
"""
@@ -66,7 +69,7 @@ class Stretching_d(DiscreteOperator):
if self.propertyOp == 'divConservation' :
## Space discretisation of grad(u)
- self.stretchOp = DifferentialOperator_d(self.vorticity, self.velocity, choice='gradU', resolution=self.vorticity.topology.mesh.resolution, meshSize=self.vorticity.topology.mesh.size)
+ self.stretchOp = DifferentialOperator_d(self.vorticity, self.velocity, choice='gradU', topology=self.topology)
gradResult, maxgersh = self.stretchOp.apply()
## Determination of Stretching time step (stability condition)
@@ -78,7 +81,7 @@ class Stretching_d(DiscreteOperator):
print "dtAdvec", 1./maxgersh[1]
## fct2Op
- self.fonction = Fct2Op(self.vorticity.data, gradResult,choice = 'gradU', resolution=self.vorticity.topology.mesh.resolution, meshSize=self.vorticity.topology.mesh.size)
+ self.fonction = Fct2Op(self.vorticity.data, gradResult,choice = 'gradU', topology=self.topology)
elif self.propertyOp == 'massConservation' :
## Determination of Stretching time step (stability condition)
@@ -86,7 +89,7 @@ class Stretching_d(DiscreteOperator):
self.dtstretch = dt
## fct2Op
- self.fonction = Fct2Op(self.vorticity, self.velocity,choice = 'divWU', resolution=self.vorticity.topology.mesh.resolution, meshSize=self.vorticity.topology.mesh.size)
+ self.fonction = Fct2Op(self.vorticity, self.velocity,choice = 'divWU', topology=self.topology)
else:
raise ValueError("Unknown Stretching Operator property: " + str(self.propertyOp))
@@ -95,10 +98,15 @@ class Stretching_d(DiscreteOperator):
for i in range (self.ndt) :
methodInt=self.method(self.fonction)
self.res_vorticity = methodInt.integrate(self.fonction.apply, t, self.dtstretch, self.vorticity.data)
-
self.vorticity.data= np.copy(self.res_vorticity)
+
+ # Ghosts synchronization
+# self.OpSynchronize.apply([self.vorticity.data, self.velocity.data])
+
+
self.compute_time = time.time() - self.compute_time
self.total_time += self.compute_time
+
return self.compute_time
def stabilityTest(self, dt, maxi, method):
diff --git a/HySoP/hysop/particular_solvers/basic.py b/HySoP/hysop/particular_solvers/basic.py
index fb62e489083d60e271ae2eca2442da3deb64ecd0..de31f57adbfff24e6e0f119c3aa6702494b977ff 100644
--- a/HySoP/hysop/particular_solvers/basic.py
+++ b/HySoP/hysop/particular_solvers/basic.py
@@ -105,7 +105,7 @@ class ParticleSolver(Solver):
if op is self.advection:
op.discretize(result_position=self.p_position, result_scalar=self.p_scalar, method=self.ODESolver)
if op is self.stretch:
- op.discretize(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()
diff --git a/HySoP/hysop/particular_solvers/integrator/euler.py b/HySoP/hysop/particular_solvers/integrator/euler.py
index 0a656e695ee4e5e257d1e3d9912493d3f748736d..819c94712b6720acafad7c229a7bd90191f74d54 100755
--- a/HySoP/hysop/particular_solvers/integrator/euler.py
+++ b/HySoP/hysop/particular_solvers/integrator/euler.py
@@ -38,7 +38,7 @@ class Euler(ODESolver):
@param t : current time.
@param dt : time step.
# """
- result = y[:,:,:,:] + dt * np.asarray(f(t,y[:,:,:,:]))
+ result = y + dt * np.asarray(f(t,y))
return result
diff --git a/HySoP/hysop/physics/__init__.py b/HySoP/hysop/physics/__init__.py
new file mode 100644
index 0000000000000000000000000000000000000000..36aa39913183557e0a4b87007a07d1e1bfba8cc9
--- /dev/null
+++ b/HySoP/hysop/physics/__init__.py
@@ -0,0 +1,6 @@
+"""
+@package physics
+
+All about physical quantities
+
+"""
diff --git a/HySoP/hysop/physics/compute_forces.py b/HySoP/hysop/physics/compute_forces.py
new file mode 100644
index 0000000000000000000000000000000000000000..bc6dfc4ca67c4dadc8974228380637ffd1a43a40
--- /dev/null
+++ b/HySoP/hysop/physics/compute_forces.py
@@ -0,0 +1,114 @@
+# -*- coding: utf-8 -*-
+"""
+@package physics
+Compute forces
+"""
+from ..constants import *
+import numpy as np
+import mpi4py.MPI as MPI
+import time
+
+
+class Compute_forces(object):
+ """
+ Compute the forces according the Noca s formula
+ """
+
+ def __init__(self, topology, obstacle, boxMin= None, boxMax=None):
+ """
+ Constructor.
+ @param topology : Local topology
+ @param obstacle : obstacle
+ @param boxMin : Global minimum coordinates of the control volume
+ @param boxMax : Global maximum coordinates of the control volume
+ """
+ self.topo = topology
+ self.obstacle = obstacle
+ self.boxMin = boxMin
+ self.boxMax = boxMax
+ self.dim = topology.dim
+ self.res = topology.mesh.resolution
+ self.step = topology.mesh.size
+ self.coordMin = topology.mesh.coord_start
+ self.coordMax = topology.mesh.coord_end
+ compute_control_box()
+
+
+ def compute_control_box(self) :
+ normal=np.zeros([self.dim*2, self.dim])
+ self.Up = 0
+ self.Down = 1
+ self.East = 2
+ self.West = 3
+ self.North = 4
+ self.South = 5
+ normal[self.Up][0] = 1.0
+ normal[self.Down][0] = -1.0
+ normal[self.East][1] = 1.0
+ normal[self.West][1] = -1.0
+ normal[self.North][2] = 1.0
+ normal[self.South][2] = -1.0
+
+ tmp_ind=np.zeros([self.dim, 2])
+
+ distMin = self.boxMin - self.coordMin
+ distMax = self.boxMax - self.coordMax
+ for i in xrange(self.dim):
+ if (distMin[i]>=0. .and. distMax[i]<=0.): # the control volume is included inside the local domain
+ tmp_ind[i][0] = self.coordMin[i] + distMin[i]/self.step[i] # Up, East, North
+ tmp_ind[i][1] = self.coordMax[i] - distMax[i]/self.step[i] # Down, West, South
+ if (distMin[i]>=0. .and. distMax[i]>=0.): # only min corner of control volume is included inside the local domain
+ tmp_ind[i][0] = self.coordMin[i] + distMin[i]/self.step[i]
+ tmp_ind[i][1] = self.res[i]
+ if (distMin[i]<=0. .and. distMax[i]<=0.): # only max corner of control volume is included inside the local domain
+ tmp_ind[i][1] = self.coordMax[i] - distMax[i]/self.step[i]
+ if (distMin[i]<=0. .and. distMax[i]>=0.): # no control volume in the local domain
+ tmp_ind[i][1] = self.res[i]
+
+ ind = np.zeros(self.dim*2)
+ ind[Up] = tmp_ind[0][0]
+ ind[Down] = tmp_ind[0][1]
+ ind[East] = tmp_ind[1][0]
+ ind[West] = tmp_ind[1][1]
+ ind[North] = tmp_ind[2][0]
+ ind[South] = tmp_ind[2][1]
+
+ # Remove last point to integrate properly ...
+ ind(2:2*dime:2)=ind(2:2*dime:2)-1
+
+ # Count the number of points on each face and inside the domain
+ nbPoints=np.zeros(2*self.dim)
+ if(isMinIn(1)):
+ nbPoints[Down]=(ind[East]-ind[West]+1)*(ind[North]-ind[South]+1)
+ if(isMaxIn(1)):
+ nbPoints[Up]=(ind[East]-ind[West]+1)*(ind[North]-ind[South]+1)
+ if(isMinIn(3)):
+ nbPoints[South]=(ind[East]-ind[West]+1)*(ind[Up]-ind[Down]+1)
+ if(isMaxIn(3)):
+ nbPoints[North]=(ind[East]-ind[West]+1)*(ind[Up]-ind[Down]+1)
+ if(isMinIn(2)):
+ nbPoints[East]=(ind[Up]-ind[Down]+1)*(ind[North]-ind[South]+1)
+ if(isMaxIn(2)):
+ nbPoints[West]=(ind[Up]-ind[Down]+1)*(ind[North]-ind[South]+1)
+
+ boxDim(c_X)=ind(Up)-ind(Down)+1
+ boxDim(c_Y)=ind(East)-ind(West)+1
+ boxDim(c_Z)=ind(North)-ind(South)+1
+ nbPointsBox = boxDim(c_X)*boxDim(c_Y)*boxDim(c_Z)
+
+
+
+
+ def apply(self, *listFields):
+
+
+
+
+ def __str__(self):
+ s = "Compute_forces. "
+ return s
+
+if __name__ == "__main__":
+ print __doc__
+ print "- Provided class : Compute_forces"
+ print Compute_forces.__doc__
diff --git a/HySoP/hysop/test/test_operator/test_Stretching.py b/HySoP/hysop/test/test_operator/test_Stretching.py
index b1588d92b54c9f8d5a903ff1641c595608f9d27f..1ee8a2e576b769bf2975ec8f063cceb84b5a0a6b 100755
--- a/HySoP/hysop/test/test_operator/test_Stretching.py
+++ b/HySoP/hysop/test/test_operator/test_Stretching.py
@@ -57,7 +57,7 @@ class test_Stretching(unittest.TestCase):
def testOperatorStretching(self):
# Parameters
- nb = 32
+ nb = 16
timeStep = 0.09
finalTime = 0.09
self.t = 0.
@@ -74,14 +74,14 @@ class test_Stretching(unittest.TestCase):
stretch = pp.Stretching(velo,vorti)
## Solver creation (discretisation of objects is done in solver initialisation)
- topo3D = pp.CartesianTopology(domain=box, resolution=[nb, nb, nb], dim=3)
+ topo3D = pp.CartesianTopology(domain=box, resolution=[nb, nb, nb], dim=3, ghosts=[2.,2.,2.])
##Problem
pb = pp.Problem(topo3D, [stretch])
## Setting solver to Problem
pb.setSolver(finalTime, timeStep, solver_type='basic')
- pb.solver.ODESolver= RK4# RK3# RK2# Euler#
+ pb.solver.ODESolver= RK3#RK2#Euler#RK4#
pb.initSolver()
t1 = time.time()
diff --git a/HySoP/hysop/tools/cpu_data_transfer.py b/HySoP/hysop/tools/cpu_data_transfer.py
new file mode 100644
index 0000000000000000000000000000000000000000..825c7020cbd5ac28f7e07399e5d74baf4193e38d
--- /dev/null
+++ b/HySoP/hysop/tools/cpu_data_transfer.py
@@ -0,0 +1,354 @@
+# -*- coding: utf-8 -*-
+"""
+@package tools
+Communicator MPI to synchronize ghost
+"""
+from ..constants import *
+import numpy as np
+import mpi4py.MPI as MPI
+import time
+
+
+class Synchronize(object):
+ """
+ Synchronization of ghosts values on domain
+ """
+
+ def __init__(self, topology):
+ """
+ Constructor.
+ @param topology : Local topology
+ """
+ self.topo = topology
+ self.total_time = 0.
+ self.size= topology.mesh.resolution
+
+
+ def apply(self, *listFields):
+
+# self.topo.tabSort = np.array([[self.up,0] , [self.down,1] , [self.east,2], [self.west,3],
+# [self.north,4],[self.south,5]])
+ self.compute_time = time.time()
+
+ for f in listFields:
+# f[:] = np.asarray(ff[:])
+ for i in xrange(self.topo.tabSort[:,0].shape[0]) :
+ print 'rang, destinataire, orientation', self.topo.rank, self.topo.tabSort[i,0], self.topo.tabSort[i,1]
+ if ( self.topo.rank == self.topo.tabSort[i,0]):
+ # UP DOWN
+ if(self.topo.tabSort[i,1] == 0):
+ f[0][self.topo.mesh.resolution[0]-self.topo.ghosts[0]:self.topo.mesh.resolution[0],:,:] =\
+ f[0][self.topo.ghosts[0]:2*self.topo.ghosts[0],:,:]
+
+ if(self.topo.tabSort[i,1] == 1):
+ f[0][0:self.topo.ghosts[0],:,:] =\
+ f[0][self.topo.mesh.resolution[0]-2*self.topo.ghosts[0]:self.topo.mesh.resolution[0]-self.topo.ghosts[0],:,:]
+
+ # WEST EAST
+ if(self.topo.tabSort[i,1] == 2):
+ f[0][:,self.topo.mesh.resolution[1]-self.topo.ghosts[1]:self.topo.mesh.resolution[1],:] =\
+ f[0][:,self.topo.ghosts[1]:2*self.topo.ghosts[1],:]
+
+ if(self.topo.tabSort[i,1] == 3):
+ f[0][:,0:self.topo.ghosts[1],:] =\
+ f[0][:,self.topo.mesh.resolution[1]-2*self.topo.ghosts[1]:self.topo.mesh.resolution[1]-self.topo.ghosts[1],:]
+
+ # NORTH SOUTH
+ if(self.topo.tabSort[i,1] == 4):
+# print 'test', f[:][:,:,:].shape , 'copie',\
+# f[:,:,:,self.topo.ghosts[2]:(2*self.topo.ghosts[2])].shape, '222' , self.topo.mesh.resolution[2]-self.topo.ghosts[2],self.topo.mesh.resolution[2]
+# print 'resol' , self.topo.mesh.resolution , ' ghost++', self.topo.ghosts
+
+ f[0][:,:,self.topo.mesh.resolution[2]-self.topo.ghosts[2]:self.topo.mesh.resolution[2]] =\
+ f[0][:,:,self.topo.ghosts[2]:2*self.topo.ghosts[2]]
+
+ if(self.topo.tabSort[i,1] == 5):
+ f[0][:,:,0:self.topo.ghosts[2]] =\
+ f[0][:,:,self.topo.mesh.resolution[2]-2*self.topo.ghosts[2]:self.topo.mesh.resolution[2]-self.topo.ghosts[2]]
+ else :
+# if(self.topo.tabSort[i,1] == 'up'):
+ if(self.topo.tabSort[i,1] == 0):
+ if( self.topo.rank < self.topo.tabSort[i,0]) :
+ self.UPsend(self.topo.tabSort[i,0], f)
+ else :
+ self.UPrecv(self.topo.tabSort[i,0], f)
+# if(self.topo.tabSort[i,1] == 'down'):
+ if(self.topo.tabSort[i,1] == 1):
+ if( self.topo.rank < self.topo.tabSort[i,0]) :
+ self.DOWNsend(self.topo.tabSort[i,0], f)
+ else :
+ self.DOWNrecv(self.topo.tabSort[i,0], f)
+# if(self.topo.tabSort[i,1] == 'east'):
+ if(self.topo.tabSort[i,1] == 2):
+ if( self.topo.rank < self.topo.tabSort[i,0]) :
+ self.EASTsend(self.topo.tabSort[i,0], f)
+ else :
+ self.EASTrecv(self.topo.tabSort[i,0], f)
+# if(self.topo.tabSort[i,1] == 'west'):
+ if(self.topo.tabSort[i,1] == 3):
+ if( self.topo.rank < self.topo.tabSort[i,0]) :
+ self.WESTsend(self.topo.tabSort[i,0], f)
+ else :
+ self.WESTrecv(self.topo.tabSort[i,0], f)
+# if(self.topo.tabSort[i,1] == 'north'):
+ if(self.topo.tabSort[i,1] == 4):
+ if( self.topo.rank < self.topo.tabSort[i,0]) :
+ self.NORTHsend(self.topo.tabSort[i,0], f)
+ else :
+ self.NORTHrecv(self.topo.tabSort[i,0], f)
+# if(self.topo.tabSort[i,1] == 'south'):
+ if(self.topo.tabSort[i,1] == 5):
+ if( self.topo.rank < self.topo.tabSort[i,0]) :
+ self.SOUTHsend(self.topo.tabSort[i,0], f)
+ else :
+ self.SOUTHrecv(self.topo.tabSort[i,0], f)
+
+
+ self.compute_time = time.time() - self.compute_time
+ self.total_time += self.compute_time
+
+ return self.compute_time
+
+
+ def UPsend(self, proc, listFields):
+ comm = MPI.COMM_WORLD
+ comm.Send(listFields[0][ \
+ self.topo.mesh.resolution[0]-2*self.topo.ghosts[0]:self.topo.mesh.resolution[0]-self.topo.ghosts[0],:,:]
+ , dest=proc, tag=77)
+ data = listFields[0][ \
+ self.topo.mesh.resolution[0]-self.topo.ghosts[0]:self.topo.mesh.resolution[0],:,:]
+ comm.Recv(data, source=proc, tag=77)
+ listFields[0][self.topo.mesh.resolution[0]-self.topo.ghosts[0]:self.topo.mesh.resolution[0],:,:] = data
+
+ def DOWNsend(self, proc, listFields):
+ comm = MPI.COMM_WORLD
+ comm.Send(listFields[0][ \
+ self.topo.ghosts[0]:2*self.topo.ghosts[0],:,:]
+ , dest=proc, tag=77)
+ data = listFields[0][ \
+ 0:self.topo.ghosts[0],:,:]
+ comm.Recv(data, source=proc, tag=77)
+ listFields[0][0:self.topo.ghosts[0],:,:] = data
+
+ def EASTsend(self, proc, listFields):
+ comm = MPI.COMM_WORLD
+ subsizes = np.asarray(listFields[0][ \
+ :,self.topo.mesh.resolution[1]-2*self.topo.ghosts[1]:self.topo.mesh.resolution[1]-self.topo.ghosts[1],:].shape)
+ startsSend = np.asarray([0,self.topo.mesh.resolution[1]-2*self.topo.ghosts[1],0])
+ print '===============SEND===================='
+ print '==================================='
+ print 'subsizes',subsizes, 'sizes',self.size
+ print 'startsSend de EAST',startsSend
+ print 'list flags' , listFields[0][ \
+ :,self.topo.mesh.resolution[1]-2*self.topo.ghosts[1]:self.topo.mesh.resolution[1]-self.topo.ghosts[1],:].flags
+ print 'test dstack', np.dstack(listFields[0][ \
+ :,self.topo.mesh.resolution[1]-2*self.topo.ghosts[1]:self.topo.mesh.resolution[1]-self.topo.ghosts[1],:]).flags
+ print '==================================='
+ EASTtypeSend = MPI.FLOAT.Create_subarray(self.size, subsizes, startsSend, order=MPI.ORDER_F) ##attention, checker le MPI type
+ EASTtypeSend.Commit()
+# for i in xrange(self.topo.dim) :
+ i=0
+ print '******************SEND***********************'
+ print '*****************************************'
+ print 'visutype', EASTtypeSend.Get_envelope()
+ print '*****************************************'
+ print '*****************************************'
+ print 'taille LFS' , listFields[i][ \
+ :,self.topo.mesh.resolution[1]-2*self.topo.ghosts[1]:self.topo.mesh.resolution[1]-self.topo.ghosts[1],:].shape , 'recv', listFields[i][ \
+ :,self.topo.mesh.resolution[1]-self.topo.ghosts[1]:self.topo.mesh.resolution[1],:].shape
+ comm.Send([ np.dstack(listFields[i][ \
+ :,self.topo.mesh.resolution[1]-2*self.topo.ghosts[1]:self.topo.mesh.resolution[1]-self.topo.ghosts[1],:]), EASTtypeSend ]
+ , dest=proc, tag=33)
+
+ startsRecv = [0,self.topo.mesh.resolution[1]-self.topo.ghosts[1],0]
+ print '===============SEND===================='
+ print '==================================='
+ print 'subsizes',subsizes, 'sizes',self.size
+ print 'startsRecv de EAST',startsRecv
+ print '==================================='
+ EASTtypeRecv = MPI.FLOAT.Create_subarray(self.size, subsizes, startsRecv, order=MPI.ORDER_F)
+ EASTtypeRecv.Commit()
+ for i in xrange(self.topo.dim) :
+ data = listFields[i][ \
+ :,self.topo.mesh.resolution[1]-self.topo.ghosts[1]:self.topo.mesh.resolution[1],:]
+ comm.Recv([np.dstack(data), EASTtypeRecv] , source=proc, tag=33)
+ listFields[i][:,self.topo.mesh.resolution[1]-self.topo.ghosts[1]:self.topo.mesh.resolution[1],:] = data
+
+ def WESTsend(self, proc, listFields):
+ comm = MPI.COMM_WORLD
+ subsizes = np.asarray(listFields[ \
+ :,self.topo.ghosts[1]:2*self.topo.ghosts[1],:].shape)
+ startsSend = np.asarray([0,self.topo.ghosts[1],0])
+ print '===============SEND===================='
+ print '==================================='
+ print 'subsizes',subsizes, 'sizes',self.size
+ print 'startsSend de WEST',startsSend
+ print '==================================='
+ WESTtypeSend = MPI.FLOAT.Create_subarray(self.size, subsizes, startsSend, order=MPI.ORDER_F)
+ WESTtypeSend.Commit()
+ for i in xrange(self.topo.dim) :
+ print 'taille LF' , listFields[i][ \
+ :,self.topo.ghosts[1]:2*self.topo.ghosts[1],:].shape
+ comm.Send([ np.dstack(listFields[i][ \
+ :,self.topo.ghosts[1]:2*self.topo.ghosts[1],:]),WESTtypeSend ]
+ , dest=proc, tag=33)
+
+ startsRecv = [0,0,0]
+ print '==============SEND====================='
+ print '==================================='
+ print 'subsizes',subsizes, 'sizes',self.size
+ print 'startsRecv de WEST',startsRecv
+ print '==================================='
+ WESTtypeRecv = MPI.FLOAT.Create_subarray(self.size, subsizes, startsRecv, order=MPI.ORDER_F)
+ WESTtypeRecv.Commit()
+ for i in xrange(self.topo.dim) :
+ data = listFields[i][ \
+ :,0:self.topo.ghosts[1],:]
+ comm.Recv([np.dstack(data), WESTtypeRecv], source=proc, tag=33)
+ listFields[i][:,0:self.topo.ghosts[1],:] = data
+
+ def NORTHsend(self, proc, listFields):
+ comm = MPI.COMM_WORLD
+ comm.Send(listFields[0][ \
+ :,:,self.topo.mesh.resolution[2]-2*self.topo.ghosts[2]:self.topo.mesh.resolution[2]-self.topo.ghosts[2]]
+ , dest=proc, tag=55)
+ data = listFields[0][ \
+ :,:,self.topo.mesh.resolution[2]-self.topo.ghosts[2]:self.topo.mesh.resolution[2]]
+ comm.Recv(data, source=proc, tag=55)
+ listFields[0][:,:,self.topo.mesh.resolution[2]-self.topo.ghosts[2]:self.topo.mesh.resolution[2]] = data
+
+ def SOUTHsend(self, proc, listFields):
+ comm = MPI.COMM_WORLD
+ comm.Send(listFields[0][ \
+ :,:,self.topo.ghosts[2]:2*self.topo.ghosts[2]]
+ , dest=proc, tag=55)
+ data = listFields[0][ \
+ :,:,0:self.topo.ghosts[2]]
+ comm.Recv(data, source=proc, tag=55)
+ listFields[0][:,:,0:self.topo.ghosts[2]] = data
+
+ def UPrecv(self, proc, listFields):
+ comm = MPI.COMM_WORLD
+ data = listFields[0][ \
+ self.topo.mesh.resolution[0]-self.topo.ghosts[0]:self.topo.mesh.resolution[0],:,:]
+ comm.Recv(data, source=proc, tag=77)
+ listFields[0][ \
+ self.topo.mesh.resolution[0]-self.topo.ghosts[0]:self.topo.mesh.resolution[0],:,:] =data
+ comm.Send(listFields[0][
+ self.topo.mesh.resolution[0]-2*self.topo.ghosts[0]:self.topo.mesh.resolution[0]-self.topo.ghosts[0],:,:]
+ , dest=proc, tag=77)
+
+ def DOWNrecv(self, proc, listFields):
+ comm = MPI.COMM_WORLD
+ data = listFields[0][ \
+ 0:self.topo.ghosts[0],:,:]
+ comm.Recv(data, source=proc, tag=77)
+ listFields[0][ \
+ 0:self.topo.ghosts[0],:,:] =data
+ comm.Send(listFields[0][ \
+ self.topo.ghosts[0]:2*self.topo.ghosts[0],:,:]
+ , dest=proc, tag=77)
+
+ def EASTrecv(self, proc, listFields):
+ comm = MPI.COMM_WORLD
+ subsizes = np.asarray(listFields[0][ \
+ :,self.topo.mesh.resolution[1]-self.topo.ghosts[1]:self.topo.mesh.resolution[1],:].shape)
+ startsRecv = np.asarray([0,self.topo.mesh.resolution[1]-self.topo.ghosts[1],0])
+ print '==============RECV coucou====================='
+ print '==================================='
+## print 'ploupli', dir(listFields),listFields.shape
+ print 'subsizes',subsizes, 'sizes',self.size
+ print 'startsRecv de EAST',startsRecv
+ print '==================================='
+ EASTtypeRecv = MPI.FLOAT.Create_subarray(self.size, subsizes, startsRecv, order=MPI.ORDER_F) ###attention, plus de ORDER donc par defaut : C !!
+ EASTtypeRecv.Commit()
+ print '**************RECV***************************'
+ print '*****************************************'
+ print 'visutype', EASTtypeRecv.Get_contents()
+ print '*****************************************'
+ print '*****************************************'
+# for i in xrange(self.topo.dim) :
+## print 'i=',i,listFields.shape
+ i=0
+ data = listFields[i][:,self.size[1]-self.topo.ghosts[1]:self.size[1],:]
+ comm.Recv([np.dstack(data), EASTtypeRecv], source=proc, tag=33)
+ listFields[i][ \
+ :,self.topo.mesh.resolution[1]-self.topo.ghosts[1]:self.topo.mesh.resolution[1],:] =data
+
+ startsSend = [0,self.topo.mesh.resolution[1]-2*self.topo.ghosts[1],0]
+ print '================RECV==================='
+ print '==================================='
+ print 'subsizes',subsizes, 'sizes',self.size
+ print 'startsSend de EAST',startsSend
+ print '==================================='
+ EASTtypeSend = MPI.FLOAT.Create_subarray(self.size, subsizes, startsSend, order=MPI.ORDER_F)
+ EASTtypeSend.Commit()
+ for i in xrange(self.topo.dim) :
+ comm.Send([np.dstack(listFields[i][ \
+ :,self.topo.mesh.resolution[1]-2*self.topo.ghosts[1]:self.topo.mesh.resolution[1]-self.topo.ghosts[1],:]), EASTtypeSend]
+ , dest=proc, tag=33)
+
+ def WESTrecv(self, proc, listFields):
+ comm = MPI.COMM_WORLD
+ subsizes = np.asarray(listFields[0][ \
+ :,0:self.topo.ghosts[1],:].shape)
+ startsRecv = np.asarray([0,0,0])
+ print '===============RECV===================='
+ print '==================================='
+ print 'subsizes',subsizes, 'sizes',self.size
+ print 'startsRecv de WEST',startsRecv
+ print '==================================='
+ WESTtypeRecv= MPI.FLOAT.Create_subarray(self.size, subsizes, startsRecv, order=MPI.ORDER_F)
+ WESTtypeRecv.Commit()
+ for i in xrange(self.topo.dim) :
+ data=listFields[i][:,0:self.topo.ghost[1],:]
+ comm.Recv([np.dstack(data), WESTtypeRecv], source = proc, tag =33)
+ listFields[i][:,0:self.topo.ghost[1],:] = data
+
+ startsSend = [0,self.topo.ghosts[1],0]
+ print '================RECV==================='
+ print '==================================='
+ print 'subsizes',subsizes, 'sizes',self.size
+ print 'startsSend de WEST',startsSend
+ print '==================================='
+ WESTtypeSend= MPI.FLOAT.Create_subarray(self.size, subsizes, startsSend, order=MPI.ORDER_F)
+ WESTtypeSend.Commit()
+ for i in xrange(self.topo.dim) :
+ comm.Send([np.dstack(listFields[i][:,self.topo.ghost[1]:2*self.topo.ghost[1],:]), WESTtypeSend], dest = proc, tag =33)
+
+ def NORTHrecv(self, proc, listFields):
+ comm = MPI.COMM_WORLD
+ data = listFields[0][ \
+ :,:,self.topo.mesh.resolution[2]-self.topo.ghosts[2]:self.topo.mesh.resolution[2]]
+ comm.Recv(data, source=proc, tag=55)
+ listFields[0][ \
+ :,:,self.topo.mesh.resolution[2]-self.topo.ghosts[2]:self.topo.mesh.resolution[2]] =data
+ comm.Send(listFields[0][ \
+ :,:,self.topo.mesh.resolution[2]-2*self.topo.ghosts[2]:self.topo.mesh.resolution[2]-self.topo.ghosts[2]]
+ , dest=proc, tag=55)
+
+ def SOUTHrecv(self, proc, listFields):
+ comm = MPI.COMM_WORLD
+ data = listFields[0][ \
+ :,:,0:self.topo.ghosts[2]]
+ comm.Recv(data, source=proc, tag=55)
+ listFields[0][ \
+ :,:,0:self.topo.ghosts[2]] =data
+ comm.Send(listFields[0][ \
+ :,:,self.topo.ghosts[2]:2*self.topo.ghosts[2]]
+ , dest=proc, tag=55)
+
+ def printComputeTime(self):
+ self.timings_info[0] = "\"Synchronization total\""
+ self.timings_info[1] = str(self.total_time)
+ print "Synchronization total time : ", self.total_time
+ print "Time of the last Synchronization iteration :", self.compute_time
+
+ def __str__(self):
+ s = "Synchronize. "
+ return s
+
+if __name__ == "__main__":
+ print __doc__
+ print "- Provided class : Synchronize"
+ print Synchronize.__doc__