From 14015815d89e00d21aea4eceb9fb1f7b13cbaa75 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Franck=20P=C3=A9rignon?= <franck.perignon@imag.fr>
Date: Wed, 25 Sep 2013 08:09:51 +0000
Subject: [PATCH] Remove useless files
---
Examples/TaylorGreen3D.py | 18 +-
Examples/dataTG.py | 4 +-
Examples/postTaylor.py | 82 +-
HySoP/hysop/numerics/fct2op.py | 114 ---
HySoP/hysop/operator/discrete/density.py | 3 +-
HySoP/hysop/operator/discrete/multiphase.py | 2 -
.../hysop/operator/monitors/compute_forces.py | 6 +-
HySoP/unusedOrObsolet/differentialOperator.py | 698 +++++++++++++++++-
HySoP/unusedOrObsolet/synchronizeGhosts.py | 493 +++++++++++++
9 files changed, 1203 insertions(+), 217 deletions(-)
delete mode 100644 HySoP/hysop/numerics/fct2op.py
create mode 100644 HySoP/unusedOrObsolet/synchronizeGhosts.py
diff --git a/Examples/TaylorGreen3D.py b/Examples/TaylorGreen3D.py
index 09c38f9e3..f21b934f8 100755
--- a/Examples/TaylorGreen3D.py
+++ b/Examples/TaylorGreen3D.py
@@ -21,7 +21,7 @@ from parmepy.problem.navier_stokes import NSProblem
from parmepy.operator.monitors.printer import Printer
from parmepy.operator.monitors.energy_enstrophy import Energy_enstrophy
from dataTG import dim, nb, NBGHOSTS, ADVECTION_METHOD, VISCOSITY, \
- WITH_PROJ, OUTPUT_FREQ, FILENAME, simu
+ OUTPUT_FREQ, FILENAME, simu
from parmepy.constants import CONSERVATIVE
@@ -96,7 +96,6 @@ poisson = Poisson(velo, vorti,
vorti: nbElem})
# projection=WITH_PROJ)
-
# Bridges between the different topologies in order to
# redistribute data.
# 1 -Advection to stretching
@@ -110,7 +109,7 @@ distrStrPoiss = Redistribute([vorti, velo], stretch, poisson)
## Define the problem to solve
pb = NSProblem(operators=[advec, distrAdvStr, stretch,
distrStrPoiss,
- diffusion, poisson], # , distrPoissEnergy],
+ diffusion, poisson],
simulation=simu, dumpFreq=-1)
## Setting solver to Problem (only operators for computational tasks)
@@ -128,23 +127,16 @@ energy = Energy_enstrophy(velo, vorti,
## frequency=500,
## prefix='./res/TG_',
## ext='.vtk')
-
+vorti.setTopoInit(topofft)
+velo.setTopoInit(topofft)
pb.addMonitors([energy])
pb.setUp()
def run():
-## Solve problem
pb.solve()
-
-def printV(step):
- for vd in velo.discreteFields.values():
- print 'velo ', step, vd.norm()
- for vd in vorti.discreteFields.values():
- print 'vorti ', step, vd.norm()
-
-
+## Solve problem
from parmepy.mpi import MPI
print "Start computation ..."
time = MPI.Wtime()
diff --git a/Examples/dataTG.py b/Examples/dataTG.py
index e6efe1784..1e1a0031a 100644
--- a/Examples/dataTG.py
+++ b/Examples/dataTG.py
@@ -15,9 +15,9 @@ WITH_PROJ = True
# post treatment ...
OUTPUT_FREQ = 1
-FILENAME = './resp2_noproj/energy.dat'
+FILENAME = './TG/energy2.dat'
# simulation parameters
from parmepy.problem.simulation import Simulation
-simu = Simulation(tinit=0.0, tend=10., timeStep=0.01,
+simu = Simulation(tinit=0.0, tend=1., timeStep=0.002,
iterMax=1000000)
diff --git a/Examples/postTaylor.py b/Examples/postTaylor.py
index 0f4fe424b..2eedcdba0 100644
--- a/Examples/postTaylor.py
+++ b/Examples/postTaylor.py
@@ -1,56 +1,56 @@
import scitools.filetable as ft
import matplotlib.pyplot as plt
-file1 = open("resp1/energy.dat")
-ener = ft.read(file1)
-time = ener[:, 0]
-energy = ener[:, 1]
-enstrophy = ener[:, 2]
-#ratio = ener[:,3]
-#dedt = ener[:,4]
-#nuS = ener[:,5]
-#nuES = ener[:,6]
-
-size = energy.shape[0]
-import numpy as np
-em2 = np.zeros(size - 2)
-em1 = np.zeros(size - 2)
-em2[:] = energy[:-2]
-em1[:] = energy[1:-1]
-dt = np.zeros(size - 1)
-dt[:] = time[1:] - time[:-1]
-en = energy[2:]
-
-dedt_calc = - (3. * en[:] - 4. * em1[:] + em2[:]) / (2. * dt[1:])
+## file1 = open("resp1/energy.dat")
+## ener = ft.read(file1)
+## time = ener[:, 0]
+## energy = ener[:, 1]
+## enstrophy = ener[:, 2]
+## #ratio = ener[:,3]
+## #dedt = ener[:,4]
+## #nuS = ener[:,5]
+## #nuES = ener[:,6]
VISCOSITY = 1. / 1600.
-nuEff = dedt_calc / enstrophy[2:]
-ratio_calc = nuEff * 1600
-nuS_calc = VISCOSITY * enstrophy[2:]
-nuES_calc = nuEff * enstrophy[2:]
-
+## nuS_calc = VISCOSITY * enstrophy[2:]
# multiproc #
-file2 = open("resp2_noproj/energy.dat")
+file2 = open("resp2_RK2/energy.dat")
ener2 = ft.read(file2)
time2 = ener2[:, 0]
energy2 = ener2[:, 1]
enstrophy2 = ener2[:, 2]
+nuS_calc2 = enstrophy2[:] * VISCOSITY
-size = energy2.shape[0]
-em22 = np.zeros(size - 2)
-em12 = np.zeros(size - 2)
-em22[:] = energy2[:-2]
-em12[:] = energy2[1:-1]
-dt2 = np.zeros(size - 1)
-dt2[:] = time2[1:] - time2[:-1]
-en2 = energy2[2:]
-
-nuS_calc2 = VISCOSITY * enstrophy2[2:]
plt.ioff()
plt.figure(1)
-print time2[2:].shape
-print nuS_calc2.shape
-plt.plot(time[2:], nuS_calc)
-plt.plot(time2[2:], nuS_calc2, 'o--')
+#plt.plot(time, nuS_calc)
+plt.plot(time2, enstrophy2, 'o--')
plt.show()
+
+
+## size = energy.shape[0]
+## import numpy as np
+## em2 = np.zeros(size - 2)
+## em1 = np.zeros(size - 2)
+## em2[:] = energy[:-2]
+## em1[:] = energy[1:-1]
+## dt = np.zeros(size - 1)
+## dt[:] = time[1:] - time[:-1]
+## en = energy[2:]
+
+## dedt_calc = - (3. * en[:] - 4. * em1[:] + em2[:]) / (2. * dt[1:])
+## nuEff = dedt_calc / enstrophy[2:]
+## ratio_calc = nuEff * 1600
+## nuES_calc = nuEff * enstrophy[2:]
+
+
+## size = energy2.shape[0]
+## em22 = np.zeros(size - 2)
+## em12 = np.zeros(size - 2)
+## em22[:] = energy2[:-2]
+## em12[:] = energy2[1:-1]
+## dt2 = np.zeros(size - 1)
+## dt2[:] = time2[1:] - time2[:-1]
+## en2 = energy2[2:]
+
diff --git a/HySoP/hysop/numerics/fct2op.py b/HySoP/hysop/numerics/fct2op.py
deleted file mode 100644
index ee3ea47d0..000000000
--- a/HySoP/hysop/numerics/fct2op.py
+++ /dev/null
@@ -1,114 +0,0 @@
-# -*- coding: utf-8 -*-
-"""
-@file fct2op.py
-
-Convert differential operator for stretching into a function
-to Discrete time integration
-"""
-from parmepy.numerics.differentialOperator import DifferentialOperator
-import numpy as np
-
-
-class Fct2Op(object):
- """
- Fct2Op
- make the link between differentialOperator and
- the integration time scheme (Euler,RK,...)
- """
-
- def __init__(self, field1, field2, method='', choice='', topology=None):
- """
- Constructor.
- store field2 to work with it in the apply function
-
- @param field1 : vorticity field
- @param field2 : velocity field or grad(U)
- depending on 'choice' parameter
- @param choice : can be 'gradV' for GradU.w or 'divWU' for div(wu)
- """
- self.field1 = field1
- self.field2 = field2
- self.method = method
- self.choice = choice
- self.topology = topology
-
- def apply(self, t, y):
- """
- Apply operator.
- Fct2Op.apply is used like a def function
- """
-# self.synchro.apply()
- if self.choice == 'gradV': # field2 = gradU
- temp0 = y[0][...] * 0.
- temp1 = y[0][...] * 0.
- temp2 = y[0][...] * 0.
-
- ind0a = self.topology.mesh.local_start[0]
- ind0b = self.topology.mesh.local_end[0] + 1
- ind1a = self.topology.mesh.local_start[1]
- ind1b = self.topology.mesh.local_end[1] + 1
- ind2a = self.topology.mesh.local_start[2]
- ind2b = self.topology.mesh.local_end[2] + 1
- ind0 = np.arange(ind0a, ind0b)
- ind1 = np.arange(ind1a, ind1b)
- ind2 = np.arange(ind2a, ind2b)
-
- ## Perform Matrix/Scalar product on local mesh
- temp0[ind0a:ind0b, ind1a:ind1b, ind2a:ind2b] =\
- self.field2[0][ind0a:ind0b, ind1a:ind1b, ind2a:ind2b] *\
- y[0][ind0a:ind0b, ind1a:ind1b, ind2a:ind2b] +\
- self.field2[1][ind0a:ind0b, ind1a:ind1b, ind2a:ind2b] *\
- y[1][ind0a:ind0b, ind1a:ind1b, ind2a:ind2b] +\
- self.field2[2][ind0a:ind0b, ind1a:ind1b, ind2a:ind2b] *\
- y[2][ind0a:ind0b, ind1a:ind1b, ind2a:ind2b]
- temp1[ind0a:ind0b, ind1a:ind1b, ind2a:ind2b] =\
- self.field2[3][ind0a:ind0b, ind1a:ind1b, ind2a:ind2b] *\
- y[0][ind0a:ind0b, ind1a:ind1b, ind2a:ind2b] +\
- self.field2[4][ind0a:ind0b, ind1a:ind1b, ind2a:ind2b] *\
- y[1][ind0a:ind0b, ind1a:ind1b, ind2a:ind2b] +\
- self.field2[5][ind0a:ind0b, ind1a:ind1b, ind2a:ind2b] *\
- y[2][ind0a:ind0b, ind1a:ind1b, ind2a:ind2b]
- temp2[ind0a:ind0b, ind1a:ind1b, ind2a:ind2b] =\
- self.field2[6][ind0a:ind0b, ind1a:ind1b, ind2a:ind2b] *\
- y[0][ind0a:ind0b, ind1a:ind1b, ind2a:ind2b] +\
- self.field2[7][ind0a:ind0b, ind1a:ind1b, ind2a:ind2b] *\
- y[1][ind0a:ind0b, ind1a:ind1b, ind2a:ind2b] +\
- self.field2[8][ind0a:ind0b, ind1a:ind1b, ind2a:ind2b] *\
- y[2][ind0a:ind0b, ind1a:ind1b, ind2a:ind2b]
-# result = np.concatenate((
-# np.array([temp0]),np.array([temp1]),np.array([temp2])))
- return temp0, temp1, temp2
-
- if self.choice == 'divWU': # field2 = velocity field
- result = DifferentialOperator(
- y, self.field2, self.method, self.choice).__call__()
- return result[0], result[1], result[2]
-
- if self.choice == 'hessianU': # field2 = velocity field
- gradU, maxgersh = DifferentialOperator(
- self.field1, self.field2,
- self.method, choice='gradV').__call__()
- GradUline1 = np.concatenate(([gradU[0]], [gradU[1]], [gradU[2]]))
- GradUline2 = np.concatenate(([gradU[3]], [gradU[4]], [gradU[5]]))
- GradUline3 = np.concatenate(([gradU[6]], [gradU[7]], [gradU[8]]))
- result1, maxgersh = DifferentialOperator(
- self.field1, GradUline1,
- self.method, choice='gradV').__call__()
- result2, maxgersh = DifferentialOperator(
- self.field1, GradUline2,
- self.method, choice='gradV').__call__()
- result3, maxgersh = DifferentialOperator(
- self.field1, GradUline3,
- self.method, choice='gradV').__call__()
- return result1, result2, result3
- else:
- print 'choice', choice, 'in Fct2Op is not defined'
-
- def __str__(self):
- s = "fct2Op(DiscreteOperator). " + DiscreteOperator.__str__(self)
- return s
-
-if __name__ == "__main__":
- print __doc__
- print "- Provided class : Fct2Op"
- print Fct2Op.__doc__
diff --git a/HySoP/hysop/operator/discrete/density.py b/HySoP/hysop/operator/discrete/density.py
index f3c0bc56a..e42603e6b 100644
--- a/HySoP/hysop/operator/discrete/density.py
+++ b/HySoP/hysop/operator/discrete/density.py
@@ -4,7 +4,6 @@
Discrete MultiPhase Rot Grad P
"""
from discrete import DiscreteOperator
-from parmepy.numerics.differentialOperator import DifferentialOperator
from parmepy.constants import np, debug
from parmepy.tools.timers import timed_function
@@ -14,7 +13,7 @@ class DensityVisco_d(DiscreteOperator):
"""
@debug
- def __init__(self, density, viscosity,
+ def __init__(self, density, viscosity,
method=None, densityVal=None, viscoVal=None):
"""
Constructor.
diff --git a/HySoP/hysop/operator/discrete/multiphase.py b/HySoP/hysop/operator/discrete/multiphase.py
index 8e90160a6..2dc970dd7 100644
--- a/HySoP/hysop/operator/discrete/multiphase.py
+++ b/HySoP/hysop/operator/discrete/multiphase.py
@@ -4,8 +4,6 @@
Discrete MultiPhase Rot Grad P
"""
from discrete import DiscreteOperator
-from parmepy.numerics.differentialOperator import DifferentialOperator
-from parmepy.numerics.fct2op import Fct2Op
from parmepy.constants import np, debug
from parmepy.tools.timers import timed_function
diff --git a/HySoP/hysop/operator/monitors/compute_forces.py b/HySoP/hysop/operator/monitors/compute_forces.py
index 99323f889..723073315 100644
--- a/HySoP/hysop/operator/monitors/compute_forces.py
+++ b/HySoP/hysop/operator/monitors/compute_forces.py
@@ -3,10 +3,8 @@
@file compute_forces.py
Compute forces
"""
-from parmepy.constants import debug, np, PARMES_REAL, ORDER, \
-PARMES_INTEGER, PI, XDIR, YDIR, ZDIR
-from parmepy.mpi.topology import Cartesian
-from parmepy.numerics.fct2op import Fct2Op
+from parmepy.constants import debug, np, ORDER, \
+ PARMES_INTEGER, PI, XDIR, YDIR, ZDIR
from parmepy.operator.updateGhosts import UpdateGhosts
from parmepy.numerics.differential_operations import DivStressTensor3D
from parmepy.numerics.finite_differences import FD_C_4
diff --git a/HySoP/unusedOrObsolet/differentialOperator.py b/HySoP/unusedOrObsolet/differentialOperator.py
index fff6def42..7e377277e 100755
--- a/HySoP/unusedOrObsolet/differentialOperator.py
+++ b/HySoP/unusedOrObsolet/differentialOperator.py
@@ -1,63 +1,683 @@
-# -*- coding: utf-8 -*-
+# -*- codingind utf-8 -*-
"""
-@package parmepy.operator.Stretching
-
-Penalization operator representation
+@file differentialOperator.py
+Discrete stretching representation
"""
-from parmepy.operator.continuous import Operator
-from parmepy.operator.differentialOperator_d import DifferentialOperator_d
+import numpy as np
+from parmepy.constants import np, PARMES_REAL, ORDER
+from parmepy.numerics.method import NumMethod
+#from parmepy.tools.cpu_data_transfer import Synchronize
+from parmepy.tools.cpu_data_transfer_S import SynchronizeS
-class DifferentialOperator(Operator):
+class DifferentialOperator(NumMethod):
"""
- Differential operator
+ Operator representation.
+ DiscreteOperator.DiscreteOperator specialization.
"""
- def __init__(self, field1, field2, choice=None,
- resolutions=None, method='',
- **other_config):
+ def __init__(self, field1, field2, method='', choice=''):
"""
Constructor.
- Create a Differential operator from given velocity and vorticity variables.
+ Create a Stretching operator on a given continuous domain.
+ Work on a given field2 and field1 to compute the stretching term.
- @param field1 : velocity variable.
- @param field2 : vorticity variable.
- @param choice : choice of the differential operation in : curl , divWU , gradV, gradS
+ @param field1 : field n1 (vorticity).
+ @param field2 : field n2 (velocity or density).
+ @param method : name and order of the spacial discretization method.
+ @param choice : differential operator to discretize
"""
- Operator.__init__(self)
+ self.name = "Differential Operator Discretization"
self.field1 = field1
self.field2 = field2
self.choice = choice
- self.addVariable([field1, field2])
- self.resolutions = resolutions
self.method = method
- self.config = other_config
+ self.compute_time = 0.
+ self.topology = self.field2.topology
+ self.meshSize = self.topology.mesh.space_step
- def setUp(self):
- """
- Stretching operator discretization method.
- Create a discrete Stretching operator from given specifications.
+# def setUp(self):
+# pass
- @param field1 : Index of velocity discretisation to use.
- @param field2 : Index of vorticity discretisation to use.
- @param topology : topology.
+ def __call__(self, synchroOp=None):
+ """
+ Apply operator.
"""
- Operator.setUp(self)
- self.discreteOperator = DifferentialOperator_d(self, method=self.method,
- **self.config)
- self.discreteOperator.setUp()
+ ind0a = self.topology.mesh.local_start[0]
+ ind0b = self.topology.mesh.local_end[0] + 1
+ ind1a = self.topology.mesh.local_start[1]
+ ind1b = self.topology.mesh.local_end[1] + 1
+ ind2a = self.topology.mesh.local_start[2]
+ ind2b = self.topology.mesh.local_end[2] + 1
+ ind0 = np.arange(ind0a, ind0b)
+ ind1 = np.arange(ind1a, ind1b)
+ ind2 = np.arange(ind2a, ind2b)
+ if self.choice.find('divWU') >= 0:
+
+ # Ghosts synchronization
+# linenb = 0
+# if (self.topology.rank == 0):
+# time.sleep(2)
+# print 'Bligne 1', field1[:,linenb,linenb]
+# time.sleep(2)
+# print 'Bligne 2', field1[linenb,:,linenb]
+# time.sleep(2)
+# print 'Bligne 3', field1[linenb,linenb,:]
+# time.sleep(2)
+
+# if (self.topology.rank == 0):
+# time.sleep(2)
+# print 'Aligne 1', field1[:,linenb,linenb]
+# time.sleep(2)
+# print 'Aligne 2', field1[linenb,:,linenb]
+# time.sleep(2)
+# print 'Aligne 3', field1[linenb,linenb,:]
+# time.sleep(2)
+ OpSynchronize = Synchronize(self.topology)
+ OpSynchronize.apply(self.field2, self.field1)
+
+ synchroOp.apply()
+ if self.method.find('FD_order2') >= 0:
+ raise ValueError("2nd order scheme Not yet implemented")
+## X components of temp and result
+# 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][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][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][ind, ind]= temp1 + temp2 + temp3
+
+## Y components of temp and result
+# 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][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][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][ind, ind]= temp1 + temp2 + temp3
+
+## Z components of temp and result
+# 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][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][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])
+
+ else:
+# X components of temp and result
+ temp1 = self.field2[0][...] * 0.
+ temp1[ind0a:ind0b, ind1a:ind1b, ind2a:ind2b] = (
+ 1.0 * self.field1[0][ind0-2, ind1a:ind1b, ind2a:ind2b] *
+ self.field2[0][ind0-2, ind1a:ind1b, ind2a:ind2b] -
+ 8.0 * self.field1[0][ind0-1, ind1a:ind1b, ind2a:ind2b] *
+ self.field2[0][ind0-1, ind1a:ind1b, ind2a:ind2b] +
+ 8.0 * self.field1[0][ind0+1, ind1a:ind1b, ind2a:ind2b] *
+ self.field2[0][ind0+1, ind1a:ind1b, ind2a:ind2b] -
+ 1.0 * self.field1[0][ind0+2, ind1a:ind1b, ind2a:ind2b] *
+ 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.field1[1][ind0a:ind0b, ind1-2, ind2a:ind2b] *
+ self.field2[0][ind0a:ind0b, ind1-2, ind2a:ind2b] -
+ 8.0 * self.field1[1][ind0a:ind0b, ind1-1, ind2a:ind2b] *
+ self.field2[0][ind0a:ind0b, ind1-1, ind2a:ind2b] +
+ 8.0 * self.field1[1][ind0a:ind0b, ind1+1, ind2a:ind2b] *
+ self.field2[0][ind0a:ind0b, ind1+1, ind2a:ind2b] -
+ 1.0 * self.field1[1][ind0a:ind0b, ind1+2, ind2a:ind2b] *
+ self.field2[0][ind0a:ind0b, ind1+2, ind2a:ind2b]
+ ) / (12. * self.meshSize[1])
+
+ temp3 = self.field2[0][...] * 0.
+ temp3[ind0a:ind0b, ind1a:ind1b, ind2a:ind2b] = (
+ 1.0 * self.field1[2][ind0a:ind0b, ind1a:ind1b, ind2-2] *
+ self.field2[0][ind0a:ind0b, ind1a:ind1b, ind2-2] -
+ 8.0 * self.field1[2][ind0a:ind0b, ind1a:ind1b, ind2-1] *
+ self.field2[0][ind0a:ind0b, ind1a:ind1b, ind2-1] +
+ 8.0 * self.field1[2][ind0a:ind0b, ind1a:ind1b, ind2+1] *
+ self.field2[0][ind0a:ind0b, ind1a:ind1b, ind2+1] -
+ 1.0 * self.field1[2][ind0a:ind0b, ind1a:ind1b, ind2+2] *
+ self.field2[0][ind0a:ind0b, ind1a:ind1b, ind2+2]
+ ) / (12. * self.meshSize[2])
+
+# tmp1 = np.array([temp1 + temp2 + temp3])
+ tmp1 = temp1 + temp2 + temp3
+
+# Y components of temp and result
+ temp1 = self.field2[1][...] * 0.
+ temp1[ind0a:ind0b, ind1a:ind1b, ind2a:ind2b] = (
+ 1.0 * self.field1[0][ind0-2, ind1a:ind1b, ind2a:ind2b] *
+ self.field2[1][ind0-2, ind1a:ind1b, ind2a:ind2b] -
+ 8.0 * self.field1[0][ind0-1, ind1a:ind1b, ind2a:ind2b] *
+ self.field2[1][ind0-1, ind1a:ind1b, ind2a:ind2b] +
+ 8.0 * self.field1[0][ind0+1, ind1a:ind1b, ind2a:ind2b] *
+ self.field2[1][ind0+1, ind1a:ind1b, ind2a:ind2b] -
+ 1.0 * self.field1[0][ind0+2, ind1a:ind1b, ind2a:ind2b] *
+ self.field2[1][ind0+2, ind1a:ind1b, ind2a:ind2b]
+ ) / (12. * self.meshSize[0])
+
+ temp2 = self.field2[1][...] * 0.
+ temp2[ind0a:ind0b, ind1a:ind1b, ind2a:ind2b] = (
+ 1.0 * self.field1[1][ind0a:ind0b, ind1-2, ind2a:ind2b] *
+ self.field2[1][ind0a:ind0b, ind1-2, ind2a:ind2b] -
+ 8.0 * self.field1[1][ind0a:ind0b, ind1-1, ind2a:ind2b] *
+ self.field2[1][ind0a:ind0b, ind1-1, ind2a:ind2b] +
+ 8.0 * self.field1[1][ind0a:ind0b, ind1+1, ind2a:ind2b] *
+ self.field2[1][ind0a:ind0b, ind1+1, ind2a:ind2b] -
+ 1.0 * self.field1[1][ind0a:ind0b, ind1+2, ind2a:ind2b] *
+ self.field2[1][ind0a:ind0b, ind1+2, ind2a:ind2b]
+ ) / (12. * self.meshSize[1])
+
+ temp3 = self.field2[1][...] * 0.
+ temp3[ind0a:ind0b, ind1a:ind1b, ind2a:ind2b] = (
+ 1.0 * self.field1[2][ind0a:ind0b, ind1a:ind1b, ind2-2] *
+ self.field2[1][ind0a:ind0b, ind1a:ind1b, ind2-2] -
+ 8.0 * self.field1[2][ind0a:ind0b, ind1a:ind1b, ind2-1] *
+ self.field2[1][ind0a:ind0b, ind1a:ind1b, ind2-1] +
+ 8.0 * self.field1[2][ind0a:ind0b, ind1a:ind1b, ind2+1] *
+ self.field2[1][ind0a:ind0b, ind1a:ind1b, ind2+1] -
+ 1.0 * self.field1[2][ind0a:ind0b, ind1a:ind1b, ind2+2] *
+ self.field2[1][ind0a:ind0b, ind1a:ind1b, ind2+2]
+ ) / (12. * self.meshSize[2])
+
+# tmp2 = np.array([temp1 + temp2 + temp3])
+ tmp2 = temp1 + temp2 + temp3
+
+# Z components of temp and result
+ temp1 = self.field2[2][...] * 0.
+ temp1[ind0a:ind0b, ind1a:ind1b, ind2a:ind2b] = (
+ 1.0 * self.field1[0][ind0-2, ind1a:ind1b, ind2a:ind2b] *
+ self.field2[2][ind0-2, ind1a:ind1b, ind2a:ind2b] -
+ 8.0 * self.field1[0][ind0-1, ind1a:ind1b, ind2a:ind2b] *
+ self.field2[2][ind0-1, ind1a:ind1b, ind2a:ind2b] +
+ 8.0 * self.field1[0][ind0+1, ind1a:ind1b, ind2a:ind2b] *
+ self.field2[2][ind0+1, ind1a:ind1b, ind2a:ind2b] -
+ 1.0 * self.field1[0][ind0+2, ind1a:ind1b, ind2a:ind2b] *
+ self.field2[2][ind0+2, ind1a:ind1b, ind2a:ind2b]
+ ) / (12. * self.meshSize[0])
+
+ temp2 = self.field2[2][...] * 0.
+ temp2[ind0a:ind0b, ind1a:ind1b, ind2a:ind2b] = (
+ 1.0 * self.field1[1][ind0a:ind0b, ind1-2, ind2a:ind2b] *
+ self.field2[2][ind0a:ind0b, ind1-2, ind2a:ind2b] -
+ 8.0 * self.field1[1][ind0a:ind0b, ind1-1, ind2a:ind2b] *
+ self.field2[2][ind0a:ind0b, ind1-1, ind2a:ind2b] +
+ 8.0 * self.field1[1][ind0a:ind0b, ind1+1, ind2a:ind2b] *
+ self.field2[2][ind0a:ind0b, ind1+1, ind2a:ind2b] -
+ 1.0 * self.field1[1][ind0a:ind0b, ind1+2, ind2a:ind2b] *
+ self.field2[2][ind0a:ind0b, ind1+2, ind2a:ind2b]
+ ) / (12. * self.meshSize[1])
+
+ temp3 = self.field2[2][...] * 0.
+ temp3[ind0a:ind0b, ind1a:ind1b, ind2a:ind2b] = (
+ 1.0 * self.field1[2][ind0a:ind0b, ind1a:ind1b, ind2-2] *
+ self.field2[2][ind0a:ind0b, ind1a:ind1b, ind2-2] -
+ 8.0 * self.field1[2][ind0a:ind0b, ind1a:ind1b, ind2-1] *
+ self.field2[2][ind0a:ind0b, ind1a:ind1b, ind2-1] +
+ 8.0 * self.field1[2][ind0a:ind0b, ind1a:ind1b, ind2+1] *
+ self.field2[2][ind0a:ind0b, ind1a:ind1b, ind2+1] -
+ 1.0 * self.field1[2][ind0a:ind0b, ind1a:ind1b, ind2+2] *
+ self.field2[2][ind0a:ind0b, ind1a:ind1b, ind2+2]
+ ) / (12. * self.meshSize[2])
+
+# tmp3 = np.array([temp1 + temp2 + temp3])
+ tmp3 = temp1 + temp2 + temp3
+ result = np.concatenate((
+ np.array(tmp1), np.array(tmp2), np.array(tmp3)))
+# return result
+ return tmp1, tmp2, tmp3
+
+ elif self.choice.find('gradV') >= 0:
+
+ # Ghosts synchronization
+ synchroOp.apply()
+# OpSynchronize = Synchronize(self.topology)
+# OpSynchronize.apply(self.field2)
+ if self.method.find('FD_order2') >= 0:
+ raise ValueError("2nd order scheme Not yet implemented")
+
+ else:
+ maxArray = np.zeros(5, dtype=PARMES_REAL, order=ORDER)
+
+## Fourth order scheme
+# X components of temp and result
+# temp1 = np.zeros(
+# (self.resolution), dtype=PARMES_REAL, order=ORDER)
+ 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 = np.zeros(
+# (self.resolution), dtype=PARMES_REAL, order=ORDER)
+ 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 = np.zeros(
+# (self.resolution), dtype=PARMES_REAL, order=ORDER)
+ 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 = np.zeros(
+# (self.resolution), dtype=PARMES_REAL, order=ORDER)
+ 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 = np.zeros(
+# (self.resolution), dtype=PARMES_REAL, order=ORDER)
+ 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 = np.zeros(
+# (self.resolution), dtype=PARMES_REAL, order=ORDER)
+ 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(temp4)+abs(temp5)+abs(temp6))
+ maxadv2 = np.max(abs(temp5))
+
+# Z components of temp and result
+# temp7 = np.zeros(
+# (self.resolution), dtype=PARMES_REAL, order=ORDER)
+ 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 = np.zeros(
+# (self.resolution), dtype=PARMES_REAL, order=ORDER)
+ 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 = np.zeros(
+# (self.resolution), dtype=PARMES_REAL, order=ORDER)
+ 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(temp7)+abs(temp8)+abs(temp9))
+ maxadv3 = np.max(abs(temp9))
+
+ # grad(V) result
+ result = np.concatenate((
+ np.array([temp1, temp2, temp3]),
+ np.array([temp4, temp5, temp6]),
+ np.array([temp7, temp8, temp9])
+ ))
+ # div(V) result
+ divergence = np.array(temp1 + temp5 + temp9)
+
+ # maxima of partial derivatives : needed for stab conditions
+ # advection stab
+ maxArray[0] = max(maxstr1, maxstr2, maxstr3)
+ # stretching stab
+ maxArray[1] = max(maxadv1, maxadv2, maxadv3)
+
+ return result, maxArray, divergence
+
+ elif (self.choice.find('gradS') >= 0):
+
+ # Ghosts synchronization
+ #OpSynchronize = SynchronizeS(self.topology)
+ #OpSynchronize.apply(self.field2)
+ synchroOp.apply()
+
+ if self.method.find('FD_order2') >= 0:
+ raise ValueError("2nd order scheme Not yet implemented")
+
+ else:
+ # Fourth order scheme
+ temp1 = self.field2[...] * 0.
+ temp1[ind0a:ind0b, ind1a:ind1b, ind2a:ind2b] = (
+ 1.0 * self.field2[ind0-2, ind1a:ind1b, ind2a:ind2b] -
+ 8.0 * self.field2[ind0-1, ind1a:ind1b, ind2a:ind2b] +
+ 8.0 * self.field2[ind0+1, ind1a:ind1b, ind2a:ind2b] -
+ 1.0 * self.field2[ind0+2, ind1a:ind1b, ind2a:ind2b]
+ ) / (12. * self.meshSize[0])
+
+ temp2 = self.field2[...] * 0.
+ temp2[ind0a:ind0b, ind1a:ind1b, ind2a:ind2b] = (
+ 1.0 * self.field2[ind0a:ind0b, ind1-2, ind2a:ind2b] -
+ 8.0 * self.field2[ind0a:ind0b, ind1-1, ind2a:ind2b] +
+ 8.0 * self.field2[ind0a:ind0b, ind1+1, ind2a:ind2b] -
+ 1.0 * self.field2[ind0a:ind0b, ind1+2, ind2a:ind2b]
+ ) / (12. * self.meshSize[1])
+
+ temp3 = self.field2[...] * 0.
+ temp3[ind0a:ind0b, ind1a:ind1b, ind2a:ind2b] = (
+ 1.0 * self.field2[ind0a:ind0b, ind1a:ind1b, ind2-2] -
+ 8.0 * self.field2[ind0a:ind0b, ind1a:ind1b, ind2-1] +
+ 8.0 * self.field2[ind0a:ind0b, ind1a:ind1b, ind2+1] -
+ 1.0 * self.field2[ind0a:ind0b, ind1a:ind1b, ind2+2]
+ ) / (12. * self.meshSize[2])
+
+ result = np.array([temp1, temp2, temp3])
+
+ return result
+
+ elif self.choice.find('laplacianV') >= 0:
+
+ # Ghosts synchronization
+# OpSynchronize = Synchronize(self.topology)
+# OpSynchronize.apply(self.field2)
+
+ synchroOp.apply()
+ if self.method.find('FD_order2') >= 0:
+ raise ValueError("2nd order scheme Not yet implemented")
+
+ else:
+## Fourth order scheme for the laplacian operator
+# X components of temp and result
+ temp1 = self.field2[0][...] * 0.
+ temp1[ind0a:ind0b, ind1a:ind1b, ind2a:ind2b] = (
+ - 1. / 12. *
+ self.field2[0][ind0-2, ind1a:ind1b, ind2a:ind2b] +
+ 4. / 3. *
+ self.field2[0][ind0-1, ind1a:ind1b, ind2a:ind2b] +
+ - 5. / 2. *
+ self.field2[0][ind0+0, ind1a:ind1b, ind2a:ind2b] +
+ 4. / 3. *
+ self.field2[0][ind0+1, ind1a:ind1b, ind2a:ind2b] +
+ - 1. / 12. *
+ self.field2[0][ind0+2, ind1a:ind1b, ind2a:ind2b]
+ ) / (self.meshSize[0]**2)
+
+ temp2 = self.field2[0][...] * 0.
+ temp2[ind0a:ind0b, ind1a:ind1b, ind2a:ind2b] = (
+ -1./12. *
+ self.field2[0][ind0a:ind0b, ind1-2, ind2a:ind2b] +
+ 4. / 3. *
+ self.field2[0][ind0a:ind0b, ind1-1, ind2a:ind2b] +
+ -5. / 2. *
+ self.field2[0][ind0a:ind0b, ind1+0, ind2a:ind2b] +
+ 4. / 3. *
+ self.field2[0][ind0a:ind0b, ind1+1, ind2a:ind2b] +
+ -1. / 12. *
+ self.field2[0][ind0a:ind0b, ind1+2, ind2a:ind2b]
+ ) / (self.meshSize[1]**2)
+
+ temp3 = self.field2[0][...] * 0.
+ temp3[ind0a:ind0b, ind1a:ind1b, ind2a:ind2b] = (
+ -1. / 12. *
+ self.field2[0][ind0a:ind0b, ind1a:ind1b, ind2-2] +
+ 4. / 3. *
+ self.field2[0][ind0a:ind0b, ind1a:ind1b, ind2-1] +
+ -5. / 2. *
+ self.field2[0][ind0a:ind0b, ind1a:ind1b, ind2+0] +
+ 4. / 3. *
+ self.field2[0][ind0a:ind0b, ind1a:ind1b, ind2+1] +
+ -1. / 12. *
+ self.field2[0][ind0a:ind0b, ind1a:ind1b, ind2+2]
+ ) / (self.meshSize[2]**2)
+ temp1 = temp1 + temp2 + temp3
+
+# Y components of temp and result
+ temp4 = self.field2[1][...] * 0.
+ temp4[ind0a:ind0b, ind1a:ind1b, ind2a:ind2b] = (
+ -1. / 12. *
+ self.field2[1][ind0-2, ind1a:ind1b, ind2a:ind2b] +
+ 4. / 3. *
+ self.field2[1][ind0-1, ind1a:ind1b, ind2a:ind2b] +
+ - 5. / 2. *
+ self.field2[1][ind0+0, ind1a:ind1b, ind2a:ind2b] +
+ 4. / 3. *
+ self.field2[1][ind0+1, ind1a:ind1b, ind2a:ind2b] +
+ -1. / 12. *
+ self.field2[1][ind0+2, ind1a:ind1b, ind2a:ind2b]
+ ) / (self.meshSize[0]**2)
+
+ temp5 = self.field2[1][...] * 0.
+ temp5[ind0a:ind0b, ind1a:ind1b, ind2a:ind2b] = (
+ -1. / 12. *
+ self.field2[1][ind0a:ind0b, ind1-2, ind2a:ind2b] +
+ 4. / 3. *
+ self.field2[1][ind0a:ind0b, ind1-1, ind2a:ind2b] +
+ -5. / 2. *
+ self.field2[1][ind0a:ind0b, ind1+0, ind2a:ind2b] +
+ 4. / 3. *
+ self.field2[1][ind0a:ind0b, ind1+1, ind2a:ind2b] +
+ -1. / 12. *
+ self.field2[1][ind0a:ind0b, ind1+2, ind2a:ind2b]
+ ) / (self.meshSize[1]**2)
+
+ temp6 = self.field2[1][...] * 0.
+ temp6[ind0a:ind0b, ind1a:ind1b, ind2a:ind2b] = (
+ -1. / 12. *
+ self.field2[1][ind0a:ind0b, ind1a:ind1b, ind2-2] +
+ 4. / 3. *
+ self.field2[1][ind0a:ind0b, ind1a:ind1b, ind2-1] +
+ -5. / 2. *
+ self.field2[1][ind0a:ind0b, ind1a:ind1b, ind2+0] +
+ 4. / 3. *
+ self.field2[1][ind0a:ind0b, ind1a:ind1b, ind2+1] +
+ -1. / 12. *
+ self.field2[1][ind0a:ind0b, ind1a:ind1b, ind2+2]
+ ) / (self.meshSize[2]**2)
+ temp4 = temp4 + temp5 + temp6
+
+# Z components of temp and result
+ temp7 = self.field2[2][...] * 0.
+ temp7[ind0a:ind0b, ind1a:ind1b, ind2a:ind2b] = (
+ -1. / 12. *
+ self.field2[2][ind0-2, ind1a:ind1b, ind2a:ind2b] +
+ 4./3. *
+ self.field2[2][ind0-1, ind1a:ind1b, ind2a:ind2b] +
+ - 5. / 2. *
+ self.field2[2][ind0+0, ind1a:ind1b, ind2a:ind2b] +
+ 4. / 3. *
+ self.field2[2][ind0+1, ind1a:ind1b, ind2a:ind2b] +
+ -1. / 12. *
+ self.field2[2][ind0+2, ind1a:ind1b, ind2a:ind2b]
+ ) / (self.meshSize[0]**2)
+
+ temp8 = self.field2[2][...] * 0.
+ temp8[ind0a:ind0b, ind1a:ind1b, ind2a:ind2b] = (
+ -1. / 12. *
+ self.field2[2][ind0a:ind0b, ind1-2, ind2a:ind2b] +
+ 4. / 3. *
+ self.field2[2][ind0a:ind0b, ind1-1, ind2a:ind2b] +
+ -5. / 2. *
+ self.field2[2][ind0a:ind0b, ind1+0, ind2a:ind2b] +
+ 4. / 3. *
+ self.field2[2][ind0a:ind0b, ind1+1, ind2a:ind2b] +
+ -1. / 12. *
+ self.field2[2][ind0a:ind0b, ind1+2, ind2a:ind2b]
+ ) / (self.meshSize[1]**2)
+
+ temp9 = self.field2[2][...] * 0.
+ temp9[ind0a:ind0b, ind1a:ind1b, ind2a:ind2b] = (
+ -1. / 12. *
+ self.field2[2][ind0a:ind0b, ind1a:ind1b, ind2-2] +
+ 4. / 3. *
+ self.field2[2][ind0a:ind0b, ind1a:ind1b, ind2-1] +
+ -5. / 2. *
+ self.field2[2][ind0a:ind0b, ind1a:ind1b, ind2+0] +
+ 4. / 3. *
+ self.field2[2][ind0a:ind0b, ind1a:ind1b, ind2+1] +
+ -1. / 12. *
+ self.field2[2][ind0a:ind0b, ind1a:ind1b, ind2+2]
+ ) / (self.meshSize[2]**2)
+ temp7 = temp7 + temp8 + temp9
+
+ result = np.concatenate((
+ np.array([temp1]), np.array([temp4]), np.array([temp7])
+ ))
+
+ return result
+
+ elif self.choice.find('curl') >= 0:
+
+ # Ghosts synchronization
+# OpSynchronize = Synchronize(self.topology)
+# OpSynchronize.apply(self.field1)
+
+ synchroOp.apply()
+ if self.method.find('FD_order2') >= 0:
+ raise ValueError("2nd order scheme Not yet implemented")
+
+ else:
+ temp1 = self.field1[0][...] * 0.
+ temp2 = self.field1[0][...] * 0.
+
+ temp1[ind0a:ind0b, ind1a:ind1b, ind2a:ind2b] = (
+ 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[ind0a:ind0b, ind1a:ind1b, ind2a:ind2b] = (
+ 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])
+ res1 = temp1 - temp2
+
+ temp1[ind0a:ind0b, ind1a:ind1b, ind2a:ind2b] = (
+ 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[ind0a:ind0b, ind1a:ind1b, ind2a:ind2b] = (
+ 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])
+ res2 = temp1 - temp2
+
+ temp1[ind0a:ind0b, ind1a:ind1b, ind2a:ind2b] = (
+ 1.0 * self.field1[1][ind0-2, ind1a:ind1b, ind2a:ind2b] -
+ 8.0 * self.field1[1][ind0-1, ind1a:ind1b, ind2a:ind2b] +
+ 8.0 * self.field1[1][ind0+1, ind1a:ind1b, ind2a:ind2b] -
+ 1.0 * self.field1[1][ind0+2, ind1a:ind1b, ind2a:ind2b]
+ ) / (12. * self.meshSize[0])
+
+ temp2[ind0a:ind0b, ind1a:ind1b, ind2a:ind2b] = (
+ 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])
+ res3 = temp1 - temp2
+ result = np.concatenate((
+ np.array([res1]), np.array([res2]), np.array([res3])
+ ))
+ return result
- def __str__(self):
- """ToString method"""
- s = " Differencial operator (ContinuousOperator)"
- if self.discreteOperator is not None:
- s += "with the following discretization:\n"
- s += str(self.discreteOperator)
else:
- s += ". Not discretised"
- return s + "\n"
+ raise ValueError(
+ "Unknown differiential operator property " +
+ 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 ind ", self.total_time
+ print "\t Differential Operator ind", self.compute_time
+
+ def __str__(self):
+ s = "DifferentialOperator (DiscreteOperator). " +\
+ DiscreteOperator.__str__(self)
+ return s
if __name__ == "__main__":
print __doc__
- print "- Provided class : DifferentialOperator"
+ print "- Provided class ind DifferentialOperator"
print DifferentialOperator.__doc__
diff --git a/HySoP/unusedOrObsolet/synchronizeGhosts.py b/HySoP/unusedOrObsolet/synchronizeGhosts.py
new file mode 100644
index 000000000..0eff6764f
--- /dev/null
+++ b/HySoP/unusedOrObsolet/synchronizeGhosts.py
@@ -0,0 +1,493 @@
+"""
+@file operator/discrete/synchronizeGhosts.py
+
+Update ghost points for some fields defined on a specific topology.
+"""
+from parmepy.constants import debug, ORDER, PARMES_INTEGER, PARMES_REAL
+from parmepy.mpi.main_var import main_comm, main_rank
+from discrete import DiscreteOperator
+from parmepy.tools.timers import timed_function
+import numpy as np
+
+
+class SynchronizeGhosts_d(DiscreteOperator):
+ """
+ Ghost points synchronization.
+ """
+
+ @debug
+ def __init__(self, fieldslist, topology, transferMethod='greaterSend'):
+ """
+ Defines a way to send/recv values at ghosts points for a list
+ of fields, discretized on a given topology.
+
+ @param fieldslist : a list of fields
+ @param topology : the topology common to all fields.
+ @param transferMethod : which type of exchange is used.
+ """
+ DiscreteOperator.__init__(self, fieldslist, method=transferMethod,
+ name="Ghosts Synchronization")
+ self.input = fieldslist
+ self.output = fieldslist
+ self.compute_time = 0.
+ self.fieldslist = fieldslist
+ self.topology = topology
+ if (main_rank == 0):
+ print 'CUT SPACE', self.topology.dims
+ self.transferMethod = transferMethod
+
+ @debug
+ def setUp(self):
+
+ self.topoMPI = \
+ self.topology.comm.Create_cart(
+ self.topology.dims,
+ periods=self.topology.periods)
+ ghosts = self.topology.ghosts
+ resolution = self.topology.mesh.resolution
+ self.fieldslist = np.asarray([self.fieldslist])
+ if (self.transferMethod == 'SendRecv'):
+ # Allocation of space for reciev data
+ self.dataRecvX = np.zeros((
+ ghosts[0], resolution[1], resolution[2]),
+ dtype=PARMES_REAL, order=ORDER)
+ self.dataRecvY = np.zeros((
+ resolution[0], ghosts[1], resolution[2]),
+ dtype=PARMES_REAL, order=ORDER)
+ self.dataRecvZ = np.zeros((
+ resolution[0], resolution[1], ghosts[2]),
+ dtype=PARMES_REAL, order=ORDER)
+ self.SendRecvList = np.zeros(
+ self.topology.dim * 2,
+ order=ORDER, dtype=PARMES_INTEGER)
+ self.SendRecvList[0], self.SendRecvList[1] = \
+ self.topology.Shift(0, 1)
+ self.SendRecvList[2], self.SendRecvList[3] = \
+ self.topology.Shift(1, 1)
+ if (self.topology.dim > 2):
+ self.SendRecvList[4], self.SendRecvList[5] = \
+ self.topology.Shift(2, 1)
+
+ if (self.transferMethod == 'greaterSend'):
+ tab = self.topology.neighbours
+ self.tabSort = np.zeros(
+ [self.topology.dim * 2, 2],
+ order=ORDER, dtype=PARMES_INTEGER)
+ self.tabSort[:, 1] = range(self.topology.dim * 2)
+ self.tabSort[0, 0] = tab[0, 0]
+ self.tabSort[1, 0] = tab[1, 0]
+ self.tabSort[2, 0] = tab[0, 1]
+ self.tabSort[3, 0] = tab[1, 1]
+ if(self.topology.dim > 2):
+ self.tabSort[4, 0] = tab[0, 2]
+ self.tabSort[5, 0] = tab[1, 2]
+ tabSort1 = self.tabSort[self.tabSort[:, 0] <= main_rank]
+ tabSort2 = self.tabSort[self.tabSort[:, 0] > main_rank]
+ if (np.asarray(tabSort1[:, 0].shape) > 0):
+ tabSort1 = tabSort1[tabSort1[:, 0].argsort()]
+ if (np.asarray(tabSort2[:, 0].shape) > 0):
+ tabSort2 = tabSort2[tabSort2[:, 0].argsort()]
+ tabSort2 = self.TabInvSort(tabSort2)
+ self.tabSort = np.concatenate((tabSort1, tabSort2))
+
+ if ((self.transferMethod == 'greaterSend') or
+ (self.transferMethod == 'SendRecv')):
+ # Allocation of space for send and reciev data
+ self.dataSendX = np.zeros((
+ ghosts[0], resolution[1], resolution[2]),
+ dtype=PARMES_REAL, order=ORDER)
+ self.dataSendY = np.zeros((
+ resolution[0], ghosts[1], resolution[2]),
+ dtype=PARMES_REAL, order=ORDER)
+ self.dataSendZ = np.zeros((
+ resolution[0], resolution[1], ghosts[2]),
+ dtype=PARMES_REAL, order=ORDER)
+ # 1 - get discrete fields arrays
+ # 2 - get list of ghosts area
+ # 3 - get neighbours
+ # 4 - find what is to be send/recv, the size of the arrays and
+ # allocate buffers if required
+ # 5 - set send/recv order
+ #self.discreteOperator.setUp()
+
+ @debug
+ @timed_function
+ def apply(self, simulation=None):
+ self.resolution = self.topology.mesh.resolution
+ resolution = self.topology.mesh.resolution
+ ghosts = self.topology.ghosts
+ dim = self.topology.dim
+ # Do the send/recv as defined in setup.
+ #### if args is None: \todo : A verifier.
+# args = self.fieldslist
+# self.fieldslist = np.asarray([self.fieldslist])
+ if (self.transferMethod == 'SendRecv'):
+ for f in fieldslist:
+ # DOWN
+ if (main_rank != self.SendRecvList[0]):
+ self.dataSendZ = f[j][:, :, ghosts[2]:2. * ghosts[2]]
+ main_comm.Sendrecv(
+ self.dataSendZ,
+ dest=self.SendRecvList[0],
+ sendtag=22,
+ recvbuf=self.dataRecvZ,
+ source=main_rank,
+ recvtag=11, status=None)
+ f[j][:, :, resolution[2] - ghosts[2]: resolution[2]] = \
+ self.dataRecvZ
+ else:
+ f[j][:, :, resolution[2] - ghosts[2]: resolution[2]
+ ] = f[j][:, :, ghosts[2]:2 * ghosts[2]]
+ # UP
+ if (main_rank != self.SendRecvList[1]):
+ self.dataSendZ = np.array(np.copy(f[j][
+ :, :,
+ resolution[2] - ghosts[2]:resolution[2]],
+ order=ORDER, dtype=PARMES_REAL))
+ main_comm.Sendrecv(
+ self.dataSendZ,
+ dest=self.SendRecvList[1],
+ sendtag=11,
+ recvbuf=self.dataRecvZ,
+ source=main_rank,
+ recvtag=22, status=None)
+ else:
+ f[j][:, :, ghosts[2]:2 * ghosts[2]
+ ] = \
+ f[j][:, :, resolution[2] - ghosts[2]:resolution[2]]
+ # SOUTH
+ if (main_rank != self.SendRecvList[2]):
+ self.dataSendY = f[j][:, ghosts[1]:2. * ghosts[1], :]
+ main_comm.Sendrecv(
+ self.dataSendY,
+ dest=self.SendRecvList[2],
+ sendtag=44,
+ recvbuf=self.dataRecvY,
+ source=main_rank,
+ recvtag=33, status=None)
+ f[j][:, resolution[2] - ghosts[2]:resolution[2], :] \
+ = self.dataRecvY
+ else:
+ f[j][:, resolution[1] - ghosts[1]:resolution[1], :
+ ] = f[j][:, ghosts[1]:2 * ghosts[1], :]
+ # NORTH
+ if (main_rank != self.SendRecvList[3]):
+ self.dataSendY = np.array(np.copy(f[j][
+ :,
+ resolution[1] - ghosts[1]:resolution[1],
+ :],
+ order=ORDER, dtype=PARMES_REAL))
+ main_comm.Sendrecv(
+ self.dataSendY,
+ dest=self.SendRecvList[3],
+ sendtag=33,
+ recvbuf=self.dataRecvY,
+ source=main_rank,
+ recvtag=44, status=None)
+ else:
+ f[j][:, ghosts[1]:2 * ghosts[1], :] =\
+ f[j][:, resolution[1] - ghosts[1]:resolution[1], :]
+ # EAST
+ if (main_rank != self.SendRecvList[4]):
+ self.dataSendX = f[j][ghosts[0]:2. * ghosts[0], :, :]
+ main_comm.Sendrecv(
+ self.dataSendX,
+ dest=self.SendRecvList[4],
+ sendtag=66,
+ recvbuf=self.dataRecvX,
+ source=main_rank,
+ recvtag=55, status=None)
+ f[j][resolution[0] - ghosts[0]:resolution[0], :, :] \
+ = self.dataRecvX
+ else:
+ f[j][resolution[0] - ghosts[0]:resolution[0], :, :
+ ] = f[j][:, :, ghosts[2]:2 * ghosts[2]]
+ # WEST
+ if (main_rank != self.SendRecvList[5]):
+ self.dataSendX = np.array(np.copy(f[j][
+ resolution[0] - ghosts[0]:resolution[0],
+ :, :],
+ order=ORDER, dtype=PARMES_REAL))
+ main_comm.Sendrecv(
+ self.dataSendX,
+ dest=self.SendRecvList[5],
+ sendtag=55,
+ recvbuf=self.dataRecvX,
+ source=main_rank,
+ recvtag=66, status=None)
+ else:
+ f[j][ghosts[0]:2 * ghosts[0], :, :] =\
+ f[j][resolution[0] - ghosts[0]:resolution[0], :, :]
+
+ if (self.transferMethod == 'greaterSend'):
+ for f in self.fieldslist:
+# print 'tabsort', main_rank, self.tabSort
+ for i in xrange(self.tabSort[:, 0].shape[0]):
+# print 'Cours', main_rank, self.tabSort[i, :]
+ if (main_rank == self.tabSort[i, 0]):
+ ## Without communication
+ # WEST
+ if (self.tabSort[i, 1] == 0):
+ for j in xrange(dim):
+ #print '\nTEST', resolution, f[j][...].shape, '\t' , ghosts
+ f[j][resolution[0] - ghosts[0]: resolution[0],
+ :, :] =\
+ f[j][ghosts[0]:2 * ghosts[0], :, :]
+ # EAST
+ if (self.tabSort[i, 1] == 1):
+ for j in xrange(dim):
+ f[j][0:ghosts[0], :, :] =\
+ f[j][
+ resolution[0] - 2 *
+ ghosts[0]:resolution[0] - ghosts[0],
+ :, :]
+
+ # SOUTH
+ if (self.tabSort[i, 1] == 2):
+ for j in xrange(dim):
+ f[j][:, resolution[1] -
+ ghosts[1]:resolution[1], :] =\
+ f[j][:, ghosts[1]:2 * ghosts[1], :]
+ #NORTH
+ if(self.tabSort[i, 1] == 3):
+ for j in xrange(dim):
+ f[j][:, 0:ghosts[1], :] =\
+ f[j][
+ :,
+ resolution[1] - 2 *
+ ghosts[1]:resolution[1] - ghosts[1],
+ :]
+
+ # DOWN
+ if(self.tabSort[i, 1] == 4):
+ for j in xrange(dim):
+ f[j][:, :, resolution[2] -
+ ghosts[2]:resolution[2]] =\
+ f[j][:, :, ghosts[2]:2 * ghosts[2]]
+ # UP
+ if(self.tabSort[i, 1] == 5):
+ for j in xrange(dim):
+ f[j][:, :, 0:ghosts[2]] =\
+ f[j][
+ :, :,
+ resolution[2] - 2 *
+ ghosts[2]:resolution[2] - ghosts[2]]
+# print 'ok', main_rank, self.tabSort[i, :]
+ else:
+ ### PART COMMUNICATION
+ # WEST
+ if(self.tabSort[i, 1] == 0):
+ if (main_rank < self.tabSort[i, 0]):
+ self.WESTsend(self.tabSort[i, 0], f)
+ else:
+ self.WESTrecv(self.tabSort[i, 0], f)
+ # EST
+ if (self.tabSort[i, 1] == 1):
+ if (main_rank < self.tabSort[i, 0]):
+ self.EASTsend(self.tabSort[i, 0], f)
+ else:
+ self.EASTrecv(self.tabSort[i, 0], f)
+ # SOUTH
+ if(self.tabSort[i, 1] == 2):
+ if (main_rank < self.tabSort[i, 0]):
+ self.SOUTHsend(self.tabSort[i, 0], f)
+ else:
+ self.SOUTHrecv(self.tabSort[i, 0], f)
+ # NORTH
+ if (self.tabSort[i, 1] == 3):
+ if (main_rank < self.tabSort[i, 0]):
+ self.NORTHsend(self.tabSort[i, 0], f)
+ else:
+ self.NORTHrecv(self.tabSort[i, 0], f)
+ # DOWN
+ if (self.tabSort[i, 1] == 4):
+ if (main_rank < self.tabSort[i, 0]):
+ self.DOWNsend(self.tabSort[i, 0], f)
+ else:
+ self.DOWNrecv(self.tabSort[i, 0], f)
+ # UP
+ if (self.tabSort[i, 1] == 5):
+ if (main_rank < self.tabSort[i, 0]):
+ self.UPsend(self.tabSort[i, 0], f)
+ else:
+ self.UPrecv(self.tabSort[i, 0], f)
+# print 'ok', main_rank, self.tabSort[i, :]
+# return self.discreteOperator.apply(*args)
+
+ def TabInvSort(self, tab):
+ tmp = tab[0, 0]
+ deb = 0
+ tabcont = 0
+ tabfinal = np.copy(tab)
+ if tab[:, 0].shape[0] == 1:
+ return tab
+ for i in xrange(1, tab[:, 0].shape[0]):
+ if tmp == tab[i, 0]:
+ tabcont = tabcont + 1
+ else:
+ tmp = tab[i, 0]
+ if tabcont == 0:
+ tabfinal[deb, :] = tab[deb, :]
+ else:
+ tabNew = tab[deb:i, :]
+# for l in xrange(int(tabNew[:,1].shape[0] / 2)):
+# tmp = tabNew[2 * l, 1]
+# tabNew[2 * l, 1] = tabNew[2 * l + 1,1]
+# tabNew[2 * l + 1, 1] = tmp
+ tabNew = tabNew[np.invert(tabNew[:, 1].argsort())]
+ tabfinal[deb:i, :] = tabNew
+ deb = i
+ tabcont = 0
+ if deb == tab[:, 0].shape[0]:
+ tabfinal[deb, :] = tab[deb, :]
+ if tabcont > 0:
+ tabNew = tab[deb:tab[:, 0].shape[0], :]
+ tabNew = tabNew[np.invert(tabNew[:, 1].argsort())]
+ tabfinal[deb:tab[:, 0].shape[0], :] = tabNew
+ return tabfinal
+
+ def WESTsend(self, proc, listFields):
+ ghosts = self.topology.ghosts
+ resolution = self.resolution
+ for i in xrange(self.topology.dim):
+ self.dataSendX = np.array(np.copy(listFields[i][
+ ghosts[0]:ghosts[0] * 2, :, :]), order=ORDER)
+ main_comm.Ssend(self.dataSendX, dest=proc, tag=66)
+ main_comm.Recv(self.dataSendX, source=proc, tag=55)
+ listFields[i][0:ghosts[0], :, :] = self.dataSendX
+
+ def EASTsend(self, proc, listFields):
+ ghosts = self.topology.ghosts
+ resolution = self.resolution
+ for i in xrange(self.topology.dim):
+ self.dataSendX = np.array(np.copy(listFields[i][
+ resolution[0] - 2 * ghosts[0]:resolution[0] - ghosts[0],
+ :, :]),
+ order=ORDER)
+ main_comm.Ssend(self.dataSendX, dest=proc, tag=66)
+ main_comm.Recv(self.dataSendX, source=proc, tag=55)
+ listFields[i][resolution[0] - ghosts[0]:resolution[0], :, :] =\
+ self.dataSendX
+
+ def SOUTHsend(self, proc, listFields):
+ ghosts = self.topology.ghosts
+ resolution = self.resolution
+ for i in xrange(self.topology.dim):
+ self.dataSendY = np.array(np.copy(listFields[i][
+ :, ghosts[1]:2 * ghosts[1], :]), order=ORDER)
+ main_comm.Ssend(self.dataSendY, dest=proc, tag=44)
+ main_comm.Recv(self.dataSendY, source=proc, tag=33)
+ listFields[i][:, 0:ghosts[1], :] = self.dataSendY
+
+ def NORTHsend(self, proc, listFields):
+ ghosts = self.topology.ghosts
+ resolution = self.resolution
+ for i in xrange(self.topology.dim):
+ self.dataSendY = np.array(np.copy(listFields[i][
+ :,
+ resolution[1] - 2 * ghosts[1]:resolution[1] - ghosts[1], :]),
+ order=ORDER)
+ main_comm.Ssend(self.dataSendY, dest=proc, tag=33)
+ main_comm.Recv(self.dataSendY, source=proc, tag=44)
+ listFields[i][:, resolution[1] - ghosts[1]:resolution[1], :] =\
+ self.dataSendY
+
+ def DOWNsend(self, proc, listFields):
+ ghosts = self.topology.ghosts
+ resolution = self.resolution
+ for i in xrange(self.topology.dim):
+ self.dataSendZ = np.array(np.copy(listFields[i][
+ :, :, ghosts[2]:2 * ghosts[2]]), order=ORDER)
+ main_comm.Ssend(self.dataSendZ, dest=proc, tag=11)
+ main_comm.Recv(self.dataSendZ, source=proc, tag=22)
+ listFields[i][:, :, 0:ghosts[2]] = self.dataSendZ
+
+ def UPsend(self, proc, listFields):
+ ghosts = self.topology.ghosts
+ resolution = self.resolution
+ for i in xrange(self.topology.dim):
+ self.dataSendZ = np.array(np.copy(listFields[i][
+ :, :,
+ resolution[2] - 2 * ghosts[2]:resolution[2] - ghosts[2]
+ ]), order=ORDER)
+ main_comm.Ssend(self.dataSendZ, dest=proc, tag=22)
+ main_comm.Recv(self.dataSendZ, source=proc, tag=11)
+ listFields[i][
+ :, :,
+ resolution[2] - ghosts[2]:resolution[2]] =\
+ self.dataSendZ
+
+ def WESTrecv(self, proc, listFields):
+ ghosts = self.topology.ghosts
+ resolution = self.resolution
+ for i in xrange(self.topology.dim):
+ main_comm.Recv(self.dataSendX, source=proc, tag=66)
+ listFields[i][0:ghosts[0], :, :] = self.dataSendX
+ self.dataSendX = np.array(np.copy(listFields[i][
+ ghosts[0]:2 * ghosts[0], :, :]), order=ORDER)
+ main_comm.Ssend(self.dataSendX, dest=proc, tag=55)
+
+ def EASTrecv(self, proc, listFields):
+ ghosts = self.topology.ghosts
+ resolution = self.resolution
+ for i in xrange(self.topology.dim):
+ main_comm.Recv(self.dataSendX, source=proc, tag=66)
+ listFields[i][resolution[0] - ghosts[0]:resolution[0], :, :] =\
+ self.dataSendX
+ self.dataSendX = np.array(np.copy(listFields[i][
+ resolution[0] - 2 * ghosts[0]:resolution[0] - ghosts[0],
+ :, :]),
+ order=ORDER)
+ main_comm.Ssend(self.dataSendX, dest=proc, tag=55)
+
+ def SOUTHrecv(self, proc, listFields):
+ ghosts = self.topology.ghosts
+ resolution = self.resolution
+ for i in xrange(self.topology.dim):
+ main_comm.Recv(self.dataSendY, source=proc, tag=33)
+ listFields[i][:, 0:ghosts[1], :] = self.dataSendY
+ self.dataSendY = np.array(np.copy(listFields[i][
+ :, ghosts[1]:2 * ghosts[1], :]), order=ORDER)
+ main_comm.Ssend(self.dataSendY, dest=proc, tag=44)
+
+ def NORTHrecv(self, proc, listFields):
+ ghosts = self.topology.ghosts
+ resolution = self.resolution
+ for i in xrange(self.topology.dim):
+ main_comm.Recv(self.dataSendY, source=proc, tag=44)
+ listFields[i][:, resolution[1] - ghosts[1]:resolution[1], :] =\
+ self.dataSendY
+ self.dataSendY = np.array(np.copy(listFields[i][
+ :, resolution[1] - 2 * ghosts[1]:resolution[1] -
+ ghosts[1], :]),
+ order=ORDER)
+ main_comm.Ssend(self.dataSendY, dest=proc, tag=33)
+
+ def DOWNrecv(self, proc, listFields):
+ ghosts = self.topology.ghosts
+ resolution = self.resolution
+ for i in xrange(self.topology.dim):
+ main_comm.Recv(self.dataSendZ, source=proc, tag=22)
+ listFields[i][:, :, 0:ghosts[2]] = self.dataSendZ
+ data = np.array(np.copy(listFields[i][
+ :, :, ghosts[2]:ghosts[2] * 2]), order=ORDER)
+ main_comm.Ssend(self.dataSendZ, dest=proc, tag=11)
+
+ def UPrecv(self, proc, listFields):
+ ghosts = self.topology.ghosts
+ resolution = self.resolution
+ for i in xrange(self.topology.dim):
+ main_comm.Recv(self.dataSendZ, source=proc, tag=11)
+ listFields[i][:, :, resolution[2] - ghosts[2]:resolution[2]] =\
+ self.dataSendZ
+ self.dataSendZ = np.array(np.copy(listFields[i][
+ :, :,
+ resolution[2] - 2 * ghosts[2]:resolution[2] - ghosts[2]
+ ]), order=ORDER)
+ main_comm.Ssend(self.dataSendZ, dest=proc, tag=22)
+
+
+if (__name__ == "__main__"):
+ print __doc__
+ print "- Provided class : SynchronizeGhosts"
+ print SynchronizeGhosts_d.__doc__
--
GitLab