diff --git a/HySoP/hysop/fields/continuous.py b/HySoP/hysop/fields/continuous.py
index 70c3c66b01d2011d4d01c9f1c8d561306b422941..33c8d421102aaa1085c741444381b9d38d6ab448 100644
--- a/HySoP/hysop/fields/continuous.py
+++ b/HySoP/hysop/fields/continuous.py
@@ -128,10 +128,10 @@ class Field(object):
"""
if self.topoInit is None:
self.discreteFields.values()[0].initialize(
- self._formula,*(args + self.extraParameters))
+ self._formula, *(args + self.extraParameters))
else:
self.discreteFields[self.topoInit].initialize(
- self._formula,*(args + self.extraParameters))
+ self._formula, *(args + self.extraParameters))
def value(self, *pos):
"""
diff --git a/HySoP/hysop/numerics/differentialOperator.py b/HySoP/hysop/numerics/differentialOperator.py
deleted file mode 100755
index 7e377277eff801ff7e757265342eeb24d21e31d6..0000000000000000000000000000000000000000
--- a/HySoP/hysop/numerics/differentialOperator.py
+++ /dev/null
@@ -1,683 +0,0 @@
-# -*- codingind utf-8 -*-
-"""
-@file differentialOperator.py
-Discrete stretching representation
-"""
-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(NumMethod):
- """
- Operator representation.
- DiscreteOperator.DiscreteOperator specialization.
- """
-
- def __init__(self, field1, field2, method='', choice=''):
- """
- Constructor.
- Create a Stretching operator on a given continuous domain.
- Work on a given field2 and field1 to compute the stretching term.
-
- @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
- """
- self.name = "Differential Operator Discretization"
- self.field1 = field1
- self.field2 = field2
- self.choice = choice
- self.method = method
- self.compute_time = 0.
- self.topology = self.field2.topology
- self.meshSize = self.topology.mesh.space_step
-
-# def setUp(self):
-# pass
-
- def __call__(self, synchroOp=None):
- """
- Apply operator.
- """
- 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
-
- else:
- 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 ind DifferentialOperator"
- print DifferentialOperator.__doc__
diff --git a/HySoP/hysop/numerics/differential_operations.py b/HySoP/hysop/numerics/differential_operations.py
index ac936b6a11399f3772865e8830bd09617da3c843..3ceb9ce43ddb1a282a0c46ed4070b5d6986682ff 100755
--- a/HySoP/hysop/numerics/differential_operations.py
+++ b/HySoP/hysop/numerics/differential_operations.py
@@ -28,63 +28,103 @@ class Curl(DifferentialOperation):
"""
Computes \f$ \nabla \ times (f) \f$, f being a vector field.
"""
- def __init__(self, topo, method=None):
+ def __init__(self, topo, method='FD_C_4', memshape=None):
+
+ if method.find('FD_order2') >= 0:
+ raise ValueError("2nd order scheme Not yet implemented")
+
+ elif method.find('FD_C_4_work') >= 0:
+ # - 4th ordered FD,
+ # - work space provided by used during function call,
+ # i.e. memshape not required,
+ # - fd.compute method
+
+ # declare and create fd scheme
+ self.fd_scheme = FD_C_4(topo.mesh.space_step)
+ # connect to fd function call
+ self.fcall = self.FDCentral4_with_work
- DifferentialOperation.__init__(self, topo)
- if method is not None and method.find('FD_order2') >= 0:
- raise ValueError("2nd order scheme Not yet implemented")
else:
- assert (topo.ghosts >= 2).all(),\
- 'you need a ghost layer for FD4 scheme.'
- self.fcall = self.FDCentral4
+ # - 4th ordered FD,
+ # - local work space => memshape required,
+ # - fd.compute method
+ if memshape is None:
+ raise ValueError("memshape required for this method")
+ self._worklength = 1
+ # declare and create fd scheme
self.fd_scheme = FD_C_4(topo.mesh.space_step)
- self.indices = topo.mesh.iCompute
- self.fd_scheme.computeIndices(self.indices)
- self.nbComponents = topo.domain.dimension
+ # connect to finite differences function ...
+ self.fcall = self.FDCentral4
+ # allocate local workspace
+ self._work = [npw.zeros(memshape)
+ for i in xrange(self._worklength)]
- def __call__(self, variable, data=None):
- return self.fcall(variable, data)
+ self._indices = topo.mesh.iCompute
+ # initialize fd scheme
+ self.fd_scheme.computeIndices(self._indices)
+ # check ghosts ... must be updated when
+ # other fd schemes will be implemented.
+ assert (topo.ghosts >= 2).all(),\
+ 'you need a ghost layer for FD4 scheme.'
- def FDCentral4(self, variable, data=None):
+ def __call__(self, variable, result, work=None):
+ return self.fcall(variable, result, work)
+
+ def FDCentral4(self, variable, result, work=None):
"""
- @param variable : a list of numpy arrays
+ @param variable : input vector field
+ @param result : list of numpy arrays to save result
"""
- # Get slices for subarrays ...
- # init (if required) return data
- if data is None:
- data = []
- for i in xrange(self.nbComponents):
- data.append(variable[i].copy())
- #data[i][...] = 0.0
-
- #data[0][...] = 0.0
- #data[1][...] = 0.0
- #data[2][...] = 0.0
-
- #-- d/dy vz -- in data[XDIR]
- self.fd_scheme.compute(variable[ZDIR], YDIR, data[XDIR])
- # -- -d/dz vy -- in data[YDIR]
- work = data[YDIR]
- self.fd_scheme.compute(variable[YDIR], ZDIR, work)
- # data_x = d/dy vz - d/dz vy
- data[XDIR][self.indices] -= work[self.indices]
-
- #-- d/dz vx -- in data[YDIR]
- self.fd_scheme.compute(variable[XDIR], ZDIR, data[YDIR])
- # -- -d/dx vz -- in data[ZDIR]
- work = data[ZDIR]
- self.fd_scheme.compute(variable[ZDIR], XDIR, work)
-
- # data_y = d/dz vx - d/dx vz
- data[YDIR][self.indices] -= work[self.indices]
-
- # - d/dy vx in data[ZDIR]
- self.fd_scheme.compute(variable[XDIR], YDIR, data[ZDIR])
- data[ZDIR][self.indices] *= -1.
- # add d/dx vy, add to data[ZDIR]
- self.fd_scheme.compute_and_add(variable[YDIR], XDIR, data[ZDIR])
+ #-- d/dy vz -- in result[XDIR]
+ self.fd_scheme.compute(variable[ZDIR], YDIR, result[XDIR])
+ # -- -d/dz vy -- in work
+ self.fd_scheme.compute(variable[YDIR], ZDIR, self._work[0])
+ # result_x = d/dy vz - d/dz vy
+ result[XDIR][self._indices] -= self._work[0][self._indices]
+
+ #-- d/dz vx -- in result[YDIR]
+ self.fd_scheme.compute(variable[XDIR], ZDIR, result[YDIR])
+ # -- -d/dx vz -- in work
+ self.fd_scheme.compute(variable[ZDIR], XDIR, self._work[0])
+ # result_y = d/dz vx - d/dx vz
+ result[YDIR][self._indices] -= self._work[0][self._indices]
+
+ #-- d/dx vy in result[ZDIR]
+ self.fd_scheme.compute(variable[YDIR], XDIR, result[ZDIR])
+ # result_z = d/dx vy - d/dy vx
+ self.fd_scheme.compute(variable[XDIR], YDIR, self._work[0])
+ result[ZDIR][self._indices] -= self._work[0][self._indices]
- return data
+ return result
+
+ def FDCentral4_with_work(self, variable, result, work):
+ """
+ @param variable : input vector field
+ @param result : list of numpy arrays to save result
+ """
+ assert len(work) >= 1
+ assert work[0].shape == variable[0].shape
+ #-- d/dy vz -- in result[XDIR]
+ self.fd_scheme.compute(variable[ZDIR], YDIR, result[XDIR])
+ # -- -d/dz vy -- in work
+ self.fd_scheme.compute(variable[YDIR], ZDIR, work[0])
+ # result_x = d/dy vz - d/dz vy
+ result[XDIR][self._indices] -= work[0][self._indices]
+
+ #-- d/dz vx -- in result[YDIR]
+ self.fd_scheme.compute(variable[XDIR], ZDIR, result[YDIR])
+ # -- -d/dx vz -- in work
+ self.fd_scheme.compute(variable[ZDIR], XDIR, work[0])
+ # result_y = d/dz vx - d/dx vz
+ result[YDIR][self._indices] -= work[0][self._indices]
+
+ #-- d/dx vy in result[ZDIR]
+ self.fd_scheme.compute(variable[YDIR], XDIR, result[ZDIR])
+ # result_z = d/dx vy - d/dy vx
+ self.fd_scheme.compute(variable[XDIR], YDIR, work[0])
+ result[ZDIR][self._indices] -= work[0][self._indices]
+
+ return result
class DivV(DifferentialOperation):
@@ -133,7 +173,7 @@ class DivV(DifferentialOperation):
# - fd.compute method
# declare and create fd scheme
- self.fd_scheme = FD_C_4((topo.mesh.space_step))
+ self.fd_scheme = FD_C_4(topo.mesh.space_step)
# connect to fd function call
self.fcall = self.FDCentral4_with_work
@@ -144,7 +184,7 @@ class DivV(DifferentialOperation):
# - fd.compute_and_add method
# declare and create fd scheme
- self.fd_scheme = FD_C_4((topo.mesh.space_step))
+ self.fd_scheme = FD_C_4(topo.mesh.space_step)
# connect to fd function call
self.fcall = self.FDCentral4_CAA_with_work
@@ -157,14 +197,14 @@ class DivV(DifferentialOperation):
# Max nb of components of input fields.
self.nbComponents = len(memshape)
- self._WORKLENGTH = 1
+ self._worklength = 1
# declare and create fd scheme
- self.fd_scheme = FD_C_4((topo.mesh.space_step))
+ self.fd_scheme = FD_C_4(topo.mesh.space_step)
# connect to finite differences function ...
self.fcall = self.FDCentral4_CAA
# allocate local workspace
self._work = [npw.zeros(memshape)
- for i in xrange(self._WORKLENGTH)]
+ for i in xrange(self._worklength)]
else:
# - 4th ordered FD,
@@ -176,16 +216,16 @@ class DivV(DifferentialOperation):
# Max nb of components of input fields.
self.nbComponents = len(memshape)
if self.nbComponents == 1: # 1D case
- self._WORKLENGTH = 1
+ self._worklength = 1
else: # 2 and 3D cases
- self._WORKLENGTH = 2
+ self._worklength = 2
# declare and create fd scheme
- self.fd_scheme = FD_C_4((topo.mesh.space_step))
+ self.fd_scheme = FD_C_4(topo.mesh.space_step)
# connect to finite differences function ...
self.fcall = self.FDCentral4
# allocate local workspace
self._work = [npw.zeros(memshape)
- for i in xrange(self._WORKLENGTH)]
+ for i in xrange(self._worklength)]
self._indices = topo.mesh.iCompute
# initialize fd scheme
@@ -354,8 +394,8 @@ class DivT(DifferentialOperation):
each component with the same size as var1[i].
"""
- def __init__(self, topo, method=None, memshape=None):
- if method is not None and method.find('FD_order2') >= 0:
+ def __init__(self, topo, method='FD_C_4', memshape=None):
+ if method.find('FD_order2') >= 0:
raise ValueError("2nd order scheme Not yet implemented")
elif method.find('FD_C_4_work') >= 0:
# - 4th ordered FD,
@@ -390,10 +430,10 @@ class DivT(DifferentialOperation):
# Max nb of components of input fields.
self.nbComponents = len(memshape)
- self._WORKLENGTH = 1
+ self._worklength = 1
# allocate local workspace
self._work = [npw.zeros(memshape)
- for i in xrange(self._WORKLENGTH)]
+ for i in xrange(self._worklength)]
# connect call function
self.fcall = self.FDCentral4
@@ -412,12 +452,12 @@ class DivT(DifferentialOperation):
# Max nb of components of input fields.
self.nbComponents = len(memshape)
if self.nbComponents == 1: # 1D case
- self._WORKLENGTH = 1
+ self._worklength = 1
else:
- self._WORKLENGTH = 2
+ self._worklength = 2
# allocate local workspace
self._work = [npw.zeros(memshape)
- for i in xrange(self._WORKLENGTH)]
+ for i in xrange(self._worklength)]
# connect call function
self.fcall = self.FDCentral4
@@ -467,7 +507,7 @@ class GradS(DifferentialOperation):
raise ValueError("2nd order scheme Not yet implemented")
elif method.find('FD_C_4') >= 0:
self.fcall = self.FDCentral4
- self.fd_scheme = FD_C_4((topo.mesh.space_step))
+ self.fd_scheme = FD_C_4(topo.mesh.space_step)
else:
raise ValueError("FD scheme Not yet implemented")
@@ -530,41 +570,91 @@ class GradVxW(DifferentialOperation):
Computes \f$ [\nabla(V)][W] \f$, with
\f$V\f$ and \f$W\f$ some vector fields.
"""
- def __init__(self, topo, method='FD_C_4'):
+ def __init__(self, topo, method='FD_C_4', memshape=None):
+
+ # dimension of the fields
+ self._dim = topo.domain.dimension
+
+ if self._dim == 1:
+ self._worklength = 1
+ elif self._dim == 2:
+ self._worklength = 2
+ else:
+ self._worklength = 2
+
if method.find('FD_order2') >= 0:
raise ValueError("2nd order scheme Not yet implemented")
+ elif method.find('FD_C_4_diag_work') >= 0:
+ # - 4th ordered FD
+ # - work space provided by used during function call,
+ # - fd.compute method
+ # - performs diagnostcis computation
+
+ # declare and create fd scheme
+ self.fd_scheme = FD_C_4(topo.mesh.space_step)
+ # connect to fd function call
+ self.fcall = self.FDCentral4_work_diag
+
elif method.find('FD_C_4_diag') >= 0:
# - 4th ordered FD, with diagnostics in output
- # - work space provided by used during function call,
+ # - use local work space
# - fd.compute method
# declare and create fd scheme
- self.fd_scheme = FD_C_4((topo.mesh.space_step))
+ self.fd_scheme = FD_C_4(topo.mesh.space_step)
# connect to fd function call
self.fcall = self.FDCentral4_diag
+ if memshape is None:
+ raise ValueError("memshape required for this method")
+ # allocate local workspace
+ self._work = [npw.zeros(memshape)
+ for i in xrange(self._worklength)]
else:
- self.fd_scheme = FD_C_4((topo.mesh.space_step))
+ # - 4th ordered FD
+ # - use local work space
+ # - fd.compute method
+ self.fd_scheme = FD_C_4(topo.mesh.space_step)
self.fcall = self.FDCentral4
+ if memshape is None:
+ raise ValueError("memshape required for this method")
+ # allocate local workspace
+ self._work = [npw.zeros(memshape)
+ for i in xrange(self._worklength)]
self._indices = topo.mesh.iCompute
self.fd_scheme.computeIndices(self._indices)
assert (topo.ghosts >= 2).all(),\
'you need a ghost layer for FD4 scheme.'
- # dimension of the fields
- self._dim = topo.domain.dimension
+ def __call__(self, var1, var2, result, diagnostics=None, work=None):
+ return self.fcall(var1, var2, result, diagnostics, work)
- if self._dim == 1:
- self._worklength = 1
- elif self._dim == 2:
- self._worklength = 2
- else:
- self._worklength = 2
-
- def __call__(self, var1, var2, result, work, diagnostics=None):
- return self.fcall(var1, var2, result, work, diagnostics)
+ def FDCentral4_diag(self, var1, var2, result, diagnostics, work=None):
+ """
+ """
+ assert len(result) == self._dim
+ nbc = len(var1)
+ diagnostics[:] = 0.0
+ for comp in xrange(nbc):
+ result[comp][...] = 0.0
+ self._work[1][...] = 0.0
+ for cdir in xrange(self._dim):
+ # self._work = d/dcdir (var1_comp)
+ self.fd_scheme.compute(var1[comp], cdir, self._work[0])
+ # some diagnostics ...
+ if cdir == comp:
+ tmp = np.max(abs(self._work[0]))
+ diagnostics[0] = max(diagnostics[0], tmp)
+ self._work[1][...] += abs(self._work[0])
+
+ # compute self._work = self._work.var2[cdir]
+ self._work[0][...] = self._work[0] * var2[cdir]
+ # sum to obtain nabla(var_comp) . var2, saved into result[comp]
+ np.add(result[comp], self._work[0], result[comp])
+ diagnostics[1] = max(diagnostics[1], np.max(self._work[1]))
+ return result, diagnostics
- def FDCentral4_diag(self, var1, var2, result, work, diagnostics):
+ def FDCentral4_work_diag(self, var1, var2, result, diagnostics, work):
"""
"""
assert len(result) == self._dim
@@ -590,21 +680,20 @@ class GradVxW(DifferentialOperation):
diagnostics[1] = max(diagnostics[1], np.max(work[1]))
return result, diagnostics
- def FDCentral4(self, var1, var2, result, work, diagnostics=None):
+ def FDCentral4(self, var1, var2, result, diagnostics=None, work=None):
"""
"""
assert len(result) == self._dim
- assert len(work) == self._worklength
nbc = len(var1)
for comp in xrange(nbc):
result[comp][...] = 0.0
for cdir in xrange(self._dim):
- # work = d/dcdir (var1_comp)
- self.fd_scheme.compute(var1[comp], cdir, work[0])
- # compute work = work.var2[cdir]
- work[0][...] = work[0] * var2[cdir]
+ # self._work = d/dcdir (var1_comp)
+ self.fd_scheme.compute(var1[comp], cdir, self._work[0])
+ # compute self._work = self._work.var2[cdir]
+ self._work[0][...] = self._work[0] * var2[cdir]
# sum to obtain nabla(var_comp) . var2, saved into result[comp]
- np.add(result[comp], work[0], result[comp])
+ np.add(result[comp], self._work[0], result[comp])
return result
diff --git a/HySoP/hysop/numerics/tests/test_diffOp.py b/HySoP/hysop/numerics/tests/test_diffOp.py
index 12f29c6638896383c50fd251fb718a628377a15d..b0f160fc60fe552dbb13f2d32d51095b8be3e7d5 100755
--- a/HySoP/hysop/numerics/tests/test_diffOp.py
+++ b/HySoP/hysop/numerics/tests/test_diffOp.py
@@ -1,250 +1,123 @@
# -*- coding: utf-8 -*-
import parmepy as pp
import numpy as np
-from parmepy.constants import PARMES_REAL, ORDER
-from parmepy.domain.box import Box
from parmepy.fields.continuous import Field
from parmepy.mpi.topology import Cartesian
-from parmepy.numerics.differential_operations import DivT, GradVxW, GradS, Curl, DivStressTensor3D
+from parmepy.numerics.differential_operations import DivT, GradVxW, Curl
import parmepy.tools.numpywrappers as npw
-import numpy.testing as npt
import math as m
+pi = m.pi
+cos = m.cos
+sin = m.sin
-def testDivWV():
- """Basic test for conservative formulation of the stretching term using FD scheme
- (comparison with periodic analytical solution based on Taylor Green benchmark)
- """
- # Domain and topologies definitions
- nb = 65
- Lx = Ly = Lz = 2. * np.pi
- box = Box(dimension=3, length=[Lx, Ly, Lz],
- origin=[0., 0., 0.])
- resol = [nb, nb, nb]
- nbghosts = 2
- ghosts = np.ones((box.dimension)) * nbghosts
- topoG = Cartesian(box, box.dimension, resol,
- ghosts=ghosts)
- topoNoG = Cartesian(box, box.dimension, resol)
-
- # Velocity, Vorticity and result fields
- velo = Field(domain=box, formula=vitesse,
- name='Velocity', isVector=True)
- velo_d = velo.discretize(topoG)
- velo.initialize()
-
- vorti = Field(domain=box, formula=vorticite,
- name='Vorticity', isVector=True)
- vorti_d = vorti.discretize(topoG)
- vorti.initialize()
-
- result = np.zeros_like(vorti_d.data)
-
- # Analytical field
- analyt = Field(domain=box, formula=analyticalDivWV,
- name='Analytical', isVector=True)
- analyt_d = analyt.discretize(topoNoG)
- analyt.initialize()
-
- # Div operator
- FD_method = 'FD_order4'
- StretchOp = DivT(topoG, FD_method, memshape=velo_d.data[0].shape)
- result = StretchOp(velo_d.data, vorti_d.data, result)
-
- # Numerical VS analytical
- ind = topoG.mesh.iCompute
- errX = errY = errZ = (Lx / (nb - 1)) ** 4
- print 'err = O(h**4) =', errX
- print '=============='
- print 'div(WU) test'
- print '=============='
- print 'Is the numerical solution equal to the analytical one at the order 4 ? :'
- print np.allclose(analyt_d[0], result[0][ind], rtol=errX)
- print np.allclose(analyt_d[1], result[1][ind], rtol=errY)
- print np.allclose(analyt_d[2], result[2][ind], rtol=errZ)
-
-
-def testGradVxW():
- """Basic test for non conservative formulation of the stretching term using FD scheme
- (comparison with periodic analytical solution based on Taylor Green benchmark)
- """
- # Domain and topologies definitions
- nb = 64
- Lx = Ly = Lz = 2. * np.pi
- box = Box(dimension=3, length=[Lx, Ly, Lz],
- origin=[0., 0., 0.])
- resol = [nb, nb, nb]
- nbghosts = 2
- ghosts = np.ones((box.dimension)) * nbghosts
- topoG = Cartesian(box, box.dimension, resol,
- ghosts=ghosts)
- topoNoG = Cartesian(box, box.dimension, resol)
-
- # Velocity, Vorticity and result fields
- velo = Field(domain=box, formula=vitesse,
- name='Velocity', isVector=True)
- velo_d = velo.discretize(topoG)
- velo.initialize()
-
- vorti = Field(domain=box, formula=vorticite,
- name='Vorticity', isVector=True)
- vorti_d = vorti.discretize(topoG)
- vorti.initialize()
-
- result = np.zeros_like(vorti_d.data)
-
- # Analytical field
- analyt = Field(domain=box, formula=analyticalGradVxW,
- name='Analytical', isVector=True)
- analyt_d = analyt.discretize(topoNoG)
- analyt.initialize()
-
- # Grad operator
- FD_method = 'FD_order4'
- StretchOp = GradVxW(topoG, FD_method)
- result, diag = StretchOp(velo_d.data, vorti_d.data)
-
- # Numerical VS analytical
- ind = topoG.mesh.iCompute
- errX = errY = errZ = (Lx / (nb - 1)) ** 4
- print 'err = O(h**4) =', errX
- print '=============='
- print 'grad(U).W test'
- print '=============='
- print 'Is the numerical solution equal to the analytical one at the order 4 ? :'
- print np.allclose(analyt_d[0], result[0][ind], rtol=errX)
- print np.allclose(analyt_d[1], result[1][ind], rtol=errY)
- print np.allclose(analyt_d[2], result[2][ind], rtol=errZ)
-
-def testCurl():
- """Basic test for Curl operator using FD scheme
- (comparison with periodic analytical solution)
- """
- # Domain and topologies definitions
- nb = 64
- Lx = Ly = Lz = 2. * np.pi
- box = Box(dimension=3, length=[Lx, Ly, Lz],
- origin=[0., 0., 0.])
- resol = [nb, nb, nb]
- nbghosts = 2
- ghosts = np.ones((box.dimension)) * nbghosts
- topoG = Cartesian(box, box.dimension, resol,
- ghosts=ghosts)
- topoNoG = Cartesian(box, box.dimension, resol)
-
- # Velocity and result fields
- velo = Field(domain=box, formula=vitesse,
- name='Velocity', isVector=True)
- velo_d = velo.discretize(topoG)
- velo.initialize()
- resCurl = np.zeros_like(velo_d.data)
-
- # Analytical field
- analyt = Field(domain=box, formula=vorticite,
- name='analytical', isVector=True)
- analyt_d = analyt.discretize(topoNoG)
- analyt.initialize()
-
- # Curl operator
- FD_method = 'FD_order4'
- CurlOp = Curl(topoG, FD_method)
- resCurl = CurlOp(velo_d.data)
-
- # Numerical VS analytical
- ind = topoG.mesh.iCompute
- errX = errY = errZ = (Lx / (nb - 1)) ** 4
- print 'err = O(h**4) =', errX
- print '=============='
- print 'curl(U) test'
- print '=============='
- print 'Is the numerical solution equal to the analytical one at the order 4 ? :'
- print np.allclose(analyt_d[0], resCurl[0][ind], rtol=errX)
- print np.allclose(analyt_d[1], resCurl[1][ind], rtol=errY)
- print np.allclose(analyt_d[2], resCurl[2][ind], rtol=errZ)
-
-def testDivStressTensor():
- """Basic test for divergence of stress tensor T
- (comparison with periodic analytical solution based on Taylor Green benchmark)
- """
- # Domain and topologies definitions
- nb = 64
- Lx = Ly = Lz = 2. * np.pi
- box = Box(dimension=3, length=[Lx, Ly, Lz],
- origin=[0., 0., 0.])
- resol = [nb, nb, nb]
- nbghosts = 2
- ghosts = np.ones((box.dimension)) * nbghosts
- topoG = Cartesian(box, box.dimension, resol,
- ghosts=ghosts)
- topoNoG = Cartesian(box, box.dimension, resol)
-
- # Velocity, work and result fields
- velo = Field(domain=box, formula=vitesse,
- name='Velocity', isVector=True)
- velo_d = velo.discretize(topoG)
- velo.initialize()
-
- result = npw.zeros_like(velo_d.data)
- work = npw.zeros_like(velo_d.data)
-
- # Analytical field
- analyt = Field(domain=box, formula=analyticalDivStressTensor,
- name='Analytical', isVector=True)
- analyt_d = analyt.discretize(topoNoG)
- analyt.initialize()
-
- # Grad operator
- FD_method = 'FD_order4'
- ind = topoG.mesh.iCompute
- divT = DivStressTensor3D(topoG, FD_method)
- result = divT(velo_d.data, result, work, ind)
-
- # Numerical VS analytical
- errX = errY = errZ = (Lx / (nb - 1)) ** 4
- print 'err = O(h**4) =', errX
- print '=============='
- print 'div(T) test'
- print '=============='
- print 'Is the numerical solution equal to the analytical one at the order 4 ? :'
- print np.allclose(analyt_d[0], result[0][ind], rtol=errX)
- print np.allclose(analyt_d[1], result[1][ind], rtol=errY)
- print np.allclose(analyt_d[2], result[2][ind], rtol=errZ)
-
-
-def vitesse(x, y, z):
- vx = m.sin(x) * m.cos(y) * m.cos(z)
- vy = - m.cos(x) * m.sin(y) * m.cos(z)
+## Function to compute TG velocity
+def computeVel(x, y, z):
+ vx = sin(x) * cos(y) * cos(z)
+ vy = - cos(x) * sin(y) * cos(z)
vz = 0.
return vx, vy, vz
-def vorticite(x, y, z):
- wx = - m.cos(x) * m.sin(y) * m.sin(z)
- wy = - m.sin(x) * m.cos(y) * m.sin(z)
- wz = 2. * m.sin(x) * m.sin(y) * m.cos(z)
+
+## Function to compute reference vorticity
+def computeVort(x, y, z):
+ wx = - cos(x) * sin(y) * sin(z)
+ wy = - sin(x) * cos(y) * sin(z)
+ wz = 2. * sin(x) * sin(y) * cos(z)
return wx, wy, wz
+
def analyticalDivWV(x, y, z):
- ax = - m.sin(y) * m.cos(y) * m.sin(z) * m.cos(z) * (- m.sin(x) * m.sin(x) + m.cos(x) * m.cos(x)) + \
- 2. * m.cos(x) * m.cos(x) * m.sin(z) * m.cos(z) * m.sin(y) * m.cos(y)
- ay = - 2. * m.cos(y) * m.cos(y) * m.sin(z) * m.cos(z) * m.sin(x) * m.cos(x) + \
- m.sin(x) * m.cos(x) * m.sin(z) * m.cos(z) * (m.cos(y) * m.cos(y) - m.sin(y) * m.sin(y))
+ ax = - sin(y) * cos(y) * sin(z) * cos(z) * \
+ (- sin(x) * sin(x) + cos(x) * cos(x)) + \
+ 2. * cos(x) * cos(x) * sin(z) * cos(z) * sin(y) * cos(y)
+
+ ay = - 2. * cos(y) * cos(y) * sin(z) * cos(z) * sin(x) * cos(x) + \
+ sin(x) * cos(x) * sin(z) * cos(z) * (cos(y) * cos(y)
+ - sin(y) * sin(y))
+
az = 0.
return ax, ay, az
+
def analyticalGradVxW(x, y, z):
- ax = - m.sin(y) * m.cos(y) * m.sin(z) * m.cos(z)
- ay = m.sin(x) * m.cos(x) * m.sin(z) * m.cos(z)
+ ax = - sin(y) * cos(y) * sin(z) * cos(z)
+ ay = sin(x) * cos(x) * sin(z) * cos(z)
az = 0.
return ax, ay, az
+
def analyticalDivStressTensor(x, y, z):
- ax = - 3. * m.sin(x) * m.cos(y) * m.cos(z)
- ay = 2. * m.sin(x) * m.cos(y) * m.cos(z) - m.sin(x) * m.sin(y) * m.sin(z)
+ ax = - 3. * sin(x) * cos(y) * cos(z)
+ ay = 2. * sin(x) * cos(y) * cos(z) - sin(x) * sin(y) * sin(z)
az = 0.
return ax, ay, az
-if __name__ == "__main__":
-# testDivWV()
-# testGradVxW()
-# testCurl()
- testDivStressTensor()
+
+nb = 65
+box = pp.Box(3, length=[2.0 * pi, 2.0 * pi, 2.0 * pi])
+nbElem = [nb] * 3
+topo = Cartesian(box, box.dimension, nbElem,
+ ghosts=[2, 2, 2])
+velo = Field(domain=box, formula=computeVel,
+ name='Velocity', isVector=True)
+vorti = Field(domain=box, formula=computeVort,
+ name='Vorticity', isVector=True)
+vorti.setTopoInit(topo)
+velo.discretize(topo)
+vorti.discretize(topo)
+velo.initialize()
+vorti.initialize()
+wd = vorti.discreteFields[topo]
+vd = velo.discreteFields[topo]
+ind = topo.mesh.iCompute
+
+
+def testCurl():
+ memshape = vd.data[0].shape
+ result = [npw.zeros(memshape) for i in xrange(3)]
+ curlOp = Curl(topo, memshape=memshape)
+ result = curlOp(vd.data, result)
+
+ for i in xrange(3):
+ assert np.allclose(wd.data[i][ind], result[i][ind])
+
+
+def testDivWV():
+ # Reference field
+ ref = Field(domain=box, formula=analyticalDivWV,
+ name='Analytical', isVector=True)
+ ref.discretize(topo)
+ ref.initialize()
+ refd = ref.discreteFields[topo]
+ memshape = vd.data[0].shape
+ # Div operator
+ divOp = DivT(topo, memshape=memshape)
+ result = [npw.zeros(memshape) for i in xrange(3)]
+ result = divOp(vd.data, wd.data, result)
+
+ # Numerical VS analytical
+ Lx = box.length[0]
+ errX = (Lx / (nb - 1)) ** 4
+ for i in xrange(3):
+ assert np.allclose(refd[i][ind], result[i][ind], rtol=errX)
+
+
+def testGradVxW():
+ # Reference field
+ ref = Field(domain=box, formula=analyticalGradVxW,
+ name='Analytical', isVector=True)
+ ref.discretize(topo)
+ ref.initialize()
+ refd = ref.discreteFields[topo]
+ memshape = vd.data[0].shape
+
+ gradOp = GradVxW(topo, memshape=memshape)
+ result = [npw.zeros(memshape) for i in xrange(3)]
+ result = gradOp(vd.data, wd.data, result)
+
+ # Numerical VS analytical
+ Lx = box.length[0]
+ errX = (Lx / (nb - 1)) ** 4
+ for i in xrange(3):
+ assert np.allclose(refd[i][ind], result[i][ind], rtol=errX)
diff --git a/HySoP/hysop/operator/curl.py b/HySoP/hysop/operator/curl.py
deleted file mode 100644
index 11a35350331f820d1d5bae210b9b453255fa0d68..0000000000000000000000000000000000000000
--- a/HySoP/hysop/operator/curl.py
+++ /dev/null
@@ -1,94 +0,0 @@
-"""
-@file operator/curl.py
-
-Curl of a field.
-"""
-from parmepy.constants import debug, np
-from parmepy.operator.continuous import Operator
-from parmepy.operator.discrete.curl import CurlFD
-from parmepy.operator.updateGhosts import UpdateGhosts
-
-
-class Curl(Operator):
- """
- Computes \f$ outVar = \nabla \times inVar \f$
- """
-
- @debug
- def __init__(self, invar, outvar, resolution, method=None, ghosts=None,
- topo=None):
- """
- Create an operator to compute the curl of a field.
- @param invar : the input vector field
- @param outvar : curl of input vector field
- @param resolution : resolution (global) of the fields
- @param method : method used (default = finite difference, 4th order)
- @param topo : a predefined topology to discretize invar/outvar
- """
-
- Operator.__init__(self, [invar, outvar], method, ghosts=ghosts,
- topo=topo)
-
- ## Dict of resolutions (one per variable)
- self.resolution = resolution
- self.output = [outvar]
- self.input = [invar]
- self.invar = invar
- self.outvar = outvar
- if self.method is None:
- self.method = 'FD_C4'
- ## An operator to update ghost points of input var.
- self.ghostsOperator = None
-
- @debug
- def setUp(self):
- if self.method.find('FD_order2') >= 0:
- nbGhosts = 1
- elif self.method.find('FD_C4') >= 0:
- nbGhosts = 2
- else:
- nbGhosts = 0
-
- # get (or create) the topology
- if self._predefinedTopo is not None:
- assert (self._predefinedTopo.ghosts >= nbGhosts).all()
- topo = self._predefinedTopo
-
- else:
- # same topo for all variables
- if self.ghosts is not None:
- self.ghosts[self.ghosts < nbGhosts] = nbGhosts
- else:
- self.ghosts = np.ones(self.domain.dimension) * nbGhosts
-
- # default topology constructor, with global resolution
- # and domain
- topo = self.domain.getOrCreateTopology(self.domain.dimension,
- self.resolution,
- ghosts=self.ghosts)
- # Now, discretisation of the variables
- for v in self.variables:
- self.discreteFields[v] = v.discretize(topo)
-
- # Create a discrete operator, according to the chosen method
- # (only finite differences at the time).
- self.discreteOperator = CurlFD(self.discreteFields[self.invar],
- self.discreteFields[self.outvar],
- method=self.method)
- self.discreteOperator.setUp()
- # ghost op
- self.ghostsOperator = UpdateGhosts(self.input, topo=topo)
-
- self._isUpToDate = True
-
- def apply(self, simulation=None):
- # update ghost points, if required
- if self.ghostsOperator is not None:
- self.ghostsOperator.apply()
- # computation ...
- self.discreteOperator.apply(simulation)
-
-if __name__ == "__main__":
- print __doc__
- print "- Provided class : Curl"
- print Curl.__doc__
diff --git a/HySoP/hysop/operator/differential.py b/HySoP/hysop/operator/differential.py
new file mode 100644
index 0000000000000000000000000000000000000000..aa8044b145daa4af4e7ca82b9822097c8f3929c5
--- /dev/null
+++ b/HySoP/hysop/operator/differential.py
@@ -0,0 +1,162 @@
+"""
+@file operator/differential.py
+
+Differential operators
+"""
+from parmepy.constants import debug, np
+from parmepy.operator.continuous import Operator
+from parmepy.operator.discrete.differential import CurlFD, GradFD
+from parmepy.operator.updateGhosts import UpdateGhosts
+from abc import ABCMeta, abstractmethod
+
+
+class Differential(Operator):
+ __metaclass__ = ABCMeta
+
+ ## @debug
+ ## def __new__(cls, *args, **kw):
+ ## return object.__new__(cls, *args, **kw)
+
+ @debug
+ @abstractmethod
+ def __init__(self, invar, outvar, resolutions, method=None,
+ topo=None, ghosts=None, ghostUpdate=None, name=None):
+
+ Operator.__init__(self, [invar, outvar], method, topo=topo,
+ ghosts=ghosts, name=name)
+ ## input variable
+ self.invar = invar
+ ## Curl of input
+ self.outvar = outvar
+ ## Grid resolution for each variable (dictionnary)
+ self.resolutions = resolutions
+ if self.method is None:
+ self.method = 'FD_C_4'
+ ## If ghostupdate == True, the operator will manage
+ ## ghost points update. Else, this must be done somewhere else ...
+ ## Default = yes
+ if ghostUpdate is None:
+ self.ghostUpdate = True
+ else:
+ self.ghostUpdate = ghostUpdate
+ self.output = [outvar]
+ self.input = [invar]
+
+ def pre_setUp(self):
+ if self.method.find('FD_order2') >= 0:
+ nbGhosts = 1
+ elif self.method.find('FD_C_4') >= 0:
+ nbGhosts = 2
+ else:
+ nbGhosts = 0
+
+ # get (or create) the topology
+ if self._predefinedTopo is not None:
+ assert (self._predefinedTopo.ghosts >= nbGhosts).all()
+ topo = self._predefinedTopo
+ for v in self.variables:
+ self.discreteFields[v] = v.discretize(topo)
+
+ else:
+ # same topo for all variables
+ if self.ghosts is not None:
+ self.ghosts[self.ghosts < nbGhosts] = nbGhosts
+ else:
+ self.ghosts = np.ones(self.domain.dimension) * nbGhosts
+
+ # default topology constructor, with global resolution
+ # and domain
+ for v in self.variables:
+ topo = self.domain.getOrCreateTopology(self.domain.dimension,
+ self.resolutions[v],
+ ghosts=self.ghosts)
+ self.discreteFields[v] = v.discretize(topo)
+
+ assert self.discreteFields[self.invar].topology == \
+ self.discreteFields[self.outvar].topology, \
+ 'Operator not yet implemented for multiple resolutions.'
+
+ if self.ghostUpdate:
+ topov = self.discreteFields[self.invar].topology
+ self.synchronize = UpdateGhosts(self.invar, topo=topov)
+ self.synchronize.setUp()
+
+ def apply(self, simulation=None):
+ if self.ghostUpdate:
+ self.synchronize.apply()
+ # computation ...
+ self.discreteOperator.apply(simulation)
+
+
+class Curl(Differential):
+ """
+ Computes \f$ outVar = \nabla inVar \f$
+ """
+
+ @debug
+ def __init__(self, invar, outvar, resolutions, method=None,
+ topo=None, ghosts=None, ghostUpdate=None):
+ """
+ Create an operator to compute the curl of a field.
+ @param invar : the input vector field
+ @param outvar : curl of input vector field
+ @param resolution : resolution (global) of the fields
+ @param method : method used (default = finite difference, 4th order)
+ @param topo : a predefined topology to discretize invar/outvar
+ @param ghostUpdate : true if this operator has to update its ghost
+ """
+
+ Differential.__init__(self, invar, outvar, resolutions, method=method,
+ topo=topo, ghosts=ghosts,
+ ghostUpdate=ghostUpdate, name='Curl')
+
+ @debug
+ def setUp(self):
+ self.pre_setUp()
+ # Create a discrete operator, according to the chosen method
+ # (only finite differences at the time).
+ self.discreteOperator = CurlFD(self.discreteFields[self.invar],
+ self.discreteFields[self.outvar],
+ method=self.method)
+ self.discreteOperator.setUp()
+ self._isUpToDate = True
+
+
+class Grad(Differential):
+ """
+ Computes \f$ outVar = \nabla inVar \f$
+ """
+
+ @debug
+ def __init__(self, invar, outvar, resolutions, method=None,
+ topo=None, ghosts=None, ghostUpdate=None):
+ """
+ Create an operator to compute the curl of a field.
+ @param invar : the input vector field
+ @param outvar : curl of input vector field
+ @param resolution : resolution (global) of the fields
+ @param method : method used (default = finite difference, 4th order)
+ @param topo : a predefined topology to discretize invar/outvar
+ @param ghostUpdate : true if this operator has to update its ghost
+ """
+ Differential.__init__(self, invar, outvar, resolutions, method=method,
+ topo=topo, ghosts=ghosts,
+ ghostUpdate=ghostUpdate, name='Grad')
+
+ @debug
+ def setUp(self):
+ self.pre_setUp()
+ # Create a discrete operator, according to the chosen method
+ # (only finite differences at the time).
+ self.discreteOperator = GradFD(self.discreteFields[self.invar],
+ self.discreteFields[self.outvar],
+ method=self.method)
+
+ self.discreteOperator.setUp()
+
+ self._isUpToDate = True
+
+if __name__ == "__main__":
+ print "- Provided class : Curl, Grad"
+ print Curl.__doc__
+ print Grad.__doc__
diff --git a/HySoP/hysop/operator/discrete/curl.py b/HySoP/hysop/operator/discrete/curl.py
deleted file mode 100644
index 78e2e991b709c130dba6289faf1ef342544446c3..0000000000000000000000000000000000000000
--- a/HySoP/hysop/operator/discrete/curl.py
+++ /dev/null
@@ -1,43 +0,0 @@
-# -*- coding: utf-8 -*-
-"""
-@file discrete/curl.py
-
-Discretisation of the curl operator
-"""
-from parmepy.constants import debug
-from discrete import DiscreteOperator
-from parmepy.tools.timers import timed_function
-from parmepy.numerics.differential_operations import Curl
-
-
-class CurlFD(DiscreteOperator):
- """
- Compute the curl of a discrete field.
- """
-
- @debug
- def __init__(self, invar, outvar, method=None):
- """
- @param[in] invar : input field
- @param[in,out] outvar : curl of the input field
- @param[in] method : how to compute the curl? Default = finite
- differences, 4th order.
- """
- self.invar = invar
- self.outvar = outvar
- DiscreteOperator.__init__(self, [self.invar, self.outvar], method)
- self.input = [self.invar]
- self.output = [self.outvar]
- self.curlFunction = Curl(self.invar.topology, self.method)
-
- @debug
- @timed_function
- def apply(self, simulation=None):
-
- self.curlFunction(self.invar, self.outvar)
-
-
-if __name__ == "__main__":
- print __doc__
- print "- Provided class : CurlFD"
- print CurlFD.__doc__
diff --git a/HySoP/hysop/operator/discrete/differential.py b/HySoP/hysop/operator/discrete/differential.py
new file mode 100644
index 0000000000000000000000000000000000000000..82874f9408dbed18402886fa2f61879f39c7839d
--- /dev/null
+++ b/HySoP/hysop/operator/discrete/differential.py
@@ -0,0 +1,101 @@
+# -*- coding: utf-8 -*-
+"""
+@file discrete/differential.py
+
+Discretisation of the differential operators (curl, grad ...)
+"""
+from parmepy.constants import debug
+from discrete import DiscreteOperator
+from parmepy.tools.timers import timed_function
+from parmepy.numerics.differential_operations import Curl, GradV
+import parmepy.tools.numpywrappers as npw
+from abc import ABCMeta, abstractmethod
+
+
+class Differential(DiscreteOperator):
+ __metaclass__ = ABCMeta
+
+ ## @debug
+ ## def __new__(cls, *args, **kw):
+ ## return object.__new__(cls, *args, **kw)
+
+ @debug
+ @abstractmethod
+ def __init__(self, invar, outvar, method='FD_C_4', name=None):
+ """
+ @param[in] invar : input field
+ @param[in,out] outvar : Grad of the input field.
+ Warning : must be a field with dim * invar.nbcomponents,
+ dim being the domain dimension.
+ @param[in] method : how to compute the grad? Default = finite
+ differences, 4th order.
+ """
+ self.invar = invar
+ self.outvar = outvar
+ DiscreteOperator.__init__(self, [self.invar, self.outvar],
+ method=method, name=name)
+ self.input = [self.invar]
+ self.output = [self.outvar]
+ self._dim = self.invar.topology.domain.dimension
+
+
+class CurlFD(Differential):
+ """
+ Compute the curl of a discrete field, using finite differences.
+ """
+
+ @debug
+ def __init__(self, invar, outvar, method=None):
+ """
+ @param[in] invar : input field
+ @param[in,out] outvar : Curl of the input field.
+ Warning : must be a field with dim * invar.nbcomponents,
+ dim being the domain dimension.
+ @param[in] method : how to compute the curl? Default = finite
+ differences, 4th order.
+ """
+ Differential.__init__(self, invar, outvar, method='FD_C_4_work',
+ name='Curl')
+
+ assert self.outvar.nbComponents == self.invar.nbComponents
+ self._function = Curl(self.invar.topology, self.method)
+ memshape = self.invar.data[0].shape
+ self._work = [npw.zeros(memshape)]
+
+ @debug
+ @timed_function
+ def apply(self, simulation=None):
+
+ self.outvar.data = self._function(self.invar.data,
+ self.outvar.data, self._work)
+
+
+class GradFD(Differential):
+ """
+ Compute the grad of a discrete field, using finite differences.
+ """
+
+ @debug
+ def __init__(self, invar, outvar, method='FD_C_4'):
+ """
+ @param[in] invar : input field
+ @param[in,out] outvar : Grad of the input field.
+ Warning : must be a field with dim * invar.nbcomponents,
+ dim being the domain dimension.
+ @param[in] method : how to compute the grad? Default = finite
+ differences, 4th order.
+ """
+ Differential.__init__(self, invar, outvar, method=method, name='Grad')
+ assert self.outvar.nbComponents == self._dim * self.invar.nbComponents
+ self._function = GradV(self.invar.topology, self.method)
+
+ @debug
+ @timed_function
+ def apply(self, simulation=None):
+
+ self.outvar.data = self._function(self.invar.data, self.outvar.data)
+
+if __name__ == "__main__":
+ print __doc__
+ print "- Provided class : CurlFD"
+ print CurlFD.__doc__
diff --git a/HySoP/hysop/operator/discrete/stretching.py b/HySoP/hysop/operator/discrete/stretching.py
index 95c8a1612acf9064ed798460e923cbe27bdfd150..96c2ca7622ce2086b37dea28a6e7917b4a40189b 100755
--- a/HySoP/hysop/operator/discrete/stretching.py
+++ b/HySoP/hysop/operator/discrete/stretching.py
@@ -46,7 +46,6 @@ class ConservativeStretching(DiscreteOperator):
assert self.velocity.topology == self.vorticity.topology,\
'Multiresolution case not yet implemented.'
- # We must have a topo for this operator ...
# A function to compute the gradient of a vector field
self.divFunc = DivT(self.velocity.topology, method=method,
memshape=self.velocity.data[0].shape)
@@ -56,36 +55,35 @@ class ConservativeStretching(DiscreteOperator):
result = self.divFunc(y, self.velocity.data, result)
return result
- shape = self.velocity.data[0].shape
+ memshape = self.velocity.data[0].shape
self.nbc = 3 # Stretching only in 3D and for vector fields.
- # Todo : check fd call with numpy and/or parmes fields.
+
# stability constant
self.cststretch = 0.
# set time integration method
if self.method.find('Euler') >= 0:
self.num_method = Euler
self.cststretch = 2.0
- self._workspace = [npw.zeros(shape)
+ self._workspace = [npw.zeros(memshape)
for i in xrange(2 * self.nbc)]
elif self.method.find('RK2') >= 0:
self.num_method = RK2
self.cststretch = 2.0
- self._workspace = [npw.zeros(shape)
+ self._workspace = [npw.zeros(memshape)
for i in xrange(2 * self.nbc)]
elif self.method.find('RK4') >= 0:
self.num_method = RK4
self.cststretch = 2.7853
- self._workspace = [npw.zeros(shape)
+ self._workspace = [npw.zeros(memshape)
for i in xrange(3 * self.nbc)]
else: # self.method.find('RK3') >= 0:
self.num_method = RK3
self.cststretch = 2.5127
- self._workspace = [npw.zeros(shape)
+ self._workspace = [npw.zeros(memshape)
for i in xrange(3 * self.nbc)]
# Create the time integrator
self.timeIntegrator = self.num_method(self.nbc, rhs, optim=WITH_GUESS)
# Prepare buffers used for send/recv operations in time integrator
- memshape = self.velocity.data[0].shape
self.timeIntegrator.setUp(self.velocity.topology, memshape)
@timed_function
@@ -107,7 +105,7 @@ class ConservativeStretching(DiscreteOperator):
self.vorticity.data = self.timeIntegrator(t, self.vorticity.data, dt,
work=self._workspace,
result=self.vorticity.data)
-
+
class GradStretching(DiscreteOperator):
"""
@@ -116,7 +114,7 @@ class GradStretching(DiscreteOperator):
"""
@debug
- def __init__(self, velocity, vorticity, method=None):
+ def __init__(self, velocity, vorticity, method='FD_C_4'):
"""
@param velocity : discrete field
@param vorticity : discrete field
@@ -138,52 +136,54 @@ class GradStretching(DiscreteOperator):
assert self.velocity.topology == self.vorticity.topology,\
'Multiresolution case not yet implemented.'
- # We must have a topo for this operator ...
- # A function to compute the gradient of a vector field
- self.graduv = GradVxW(self.velocity.topology, method=method)
- print 'grad UV'
+ memshape = self.velocity.data[0].shape
+ # A function to compute the gradient of a vector field:
+ # grad(V) x vorticity with some diagnostics as output.
+ self.graduv = GradVxW(self.velocity.topology,
+ method='FD_C_4_diag', memshape=memshape)
# Time integration scheme.
def rhs(t, y, result):
- result, diagnostics =\
- self.graduv(self.velocity.data, y)
- #result, diag = self.graduv(self.velocity.data, y)
+ result, self.diagnostics =\
+ self.graduv(self.velocity.data, y, result, self.diagnostics)
return result
- # Todo : check fd call with numpy and/or parmes fields.
+ self.nbc = 3 # Stretching only in 3D and for vector fields.
+
# stability constant
self.cststretch = 0.
# set time integration method
if self.method.find('Euler') >= 0:
self.num_method = Euler
self.cststretch = 2.0
+ self._workspace = [npw.zeros(memshape)
+ for i in xrange(2 * self.nbc)]
elif self.method.find('RK2') >= 0:
self.num_method = RK2
self.cststretch = 2.0
- elif self.method.find('RK3') >= 0:
- self.num_method = RK3
- self.cststretch = 2.5127
+ self._workspace = [npw.zeros(memshape)
+ for i in xrange(2 * self.nbc)]
elif self.method.find('RK4') >= 0:
self.num_method = RK4
self.cststretch = 2.7853
- else:
+ self._workspace = [npw.zeros(memshape)
+ for i in xrange(3 * self.nbc)]
+ else: # self.method.find('RK3') >= 0:
self.num_method = RK3
self.cststretch = 2.5127
- shape = self.velocity.data[0].shape
- self.nbc = 3 # Stretching only in 3D and for vector fields.
- self._workspace = [npw.zeros(shape)
- for i in xrange(3 * self.nbc)]
+ self._workspace = [npw.zeros(memshape)
+ for i in xrange(3 * self.nbc)]
- ## time integrator
- self.timeIntegrator = self.num_method(self.vorticity.nbComponents, rhs,
- optim=WITH_GUESS)
- ## a vector of float with diagnostics computed
- ## from GradVxW (max div ...)
- self.diagnostics = npw.ones((5))
+ # Create the time integrator
+ self.timeIntegrator = self.num_method(self.nbc, rhs, optim=WITH_GUESS)
+ # Prepare buffers used for send/recv operations in time integrator
+ self.timeIntegrator.setUp(self.velocity.topology, memshape)
+
+ ## a vector to save diagnostics computed from GradVxW (max div ...)
+ self.diagnostics = npw.ones((2))
@timed_function
def apply(self, simulation):
- ind = self.velocity.topology.mesh.iCompute
if simulation is None:
raise ValueError("Missing simulation value for computation.")
@@ -191,30 +191,22 @@ class GradStretching(DiscreteOperator):
dt = simulation.timeStep
# current time
t = simulation.time
- #self.diagnostics = self.graduv(self.velocity, self.vorticity, self.diagnostics)[1]
+
# Compute the number of required subcycles
ndt, subdt = self.checkStability(dt)
- # Call time integrator
- print "time steps ...", dt, subdt
- print "pre apply Grad...", self.vorticity.data[0].max()
-
assert sum(subdt) == dt
+
+ # Init workspace with a first evaluation of the
+ # rhs of the integrator
for i in xrange(ndt):
- print "Update workspace/Grad, pre", self._workspace[0][ind].max()
- self._workspace[:self.vorticity.nbComponents] = \
+ self._workspace[:self.nbc] = \
self.timeIntegrator.f(t, self.vorticity.data,
- self._workspace
- [:self.vorticity.nbComponents])
- print "Update workspace/Grad, post", self._workspace[0][ind].max()
+ self._workspace[:self.nbc])
# perform integration and save result in-place
- self.vorticity.data = self.timeIntegrator(
- t, self.vorticity.data,
- subdt[i],
- work=self._workspace,
- result=self.vorticity.data)
- #self.vorticity.data[0], self.vorticity.data[1],\
- # self.vorticity.data[2] =\
- # self.timeIntegrator(t, self.vorticity.data, subdt[i])
+ self.vorticity.data = \
+ self.timeIntegrator(t, self.vorticity.data, subdt[i],
+ work=self._workspace,
+ result=self.vorticity.data)
def checkStability(self, dt):
"""
diff --git a/HySoP/hysop/operator/discrete/synchronizeGhosts.py b/HySoP/hysop/operator/discrete/synchronizeGhosts.py
deleted file mode 100644
index 0eff6764f7709d4d2a2d1e546845794a850854a2..0000000000000000000000000000000000000000
--- a/HySoP/hysop/operator/discrete/synchronizeGhosts.py
+++ /dev/null
@@ -1,493 +0,0 @@
-"""
-@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__
diff --git a/HySoP/hysop/operator/energy_enstrophy.py b/HySoP/hysop/operator/energy_enstrophy.py
index ad2a8623a05a221f528ba7298907df9eb87d432e..b12115b85a7b12ade43836477dbe9381a6a3c564 100644
--- a/HySoP/hysop/operator/energy_enstrophy.py
+++ b/HySoP/hysop/operator/energy_enstrophy.py
@@ -4,7 +4,7 @@
Compute Energy and Enstrophy
"""
import numpy as np
-from parmepy.constants import debug # , PARMES_MPI_REAL
+from parmepy.constants import debug, XDIR # , PARMES_MPI_REAL
from parmepy.operator.monitors.monitoring import Monitoring
from parmepy.mpi import MPI
from parmepy.tools.timers import timed_function
@@ -113,8 +113,9 @@ class Energy_enstrophy(Monitoring):
# get the list of computation points (no ghosts)
nbc = vd.nbComponents
# Integrate (locally) velocity ** 2
+ self._work[...] = vd[XDIR][self.ind] ** 2
[np.add(self._work[...], vd[i][self.ind] ** 2, self._work[...])
- for i in xrange(nbc)]
+ for i in xrange(1, nbc)]
local_energy = np.sum(self._work)
# --- Enstrophy computation ---
diff --git a/HySoP/hysop/operator/stretching.py b/HySoP/hysop/operator/stretching.py
index 956bc5f39ee787ffd1eaefb38c243bf6eb7dd4f5..51151730b2200af8b760c051e1f0bf8b1ad4beb0 100755
--- a/HySoP/hysop/operator/stretching.py
+++ b/HySoP/hysop/operator/stretching.py
@@ -36,6 +36,8 @@ class Stretching(Operator):
@param topo : a predefined topology to discretize velocity/vorticity
@param ghosts : number of ghosts points. Default depends on the method.
Autom. computed if not set.
+ @param ghostUpdate : true if this operator has to update its ghost
+ points.
"""
Operator.__init__(self, [velocity, vorticity], method, topo=topo,
ghosts=ghosts, name="Stretching")
@@ -56,8 +58,9 @@ class Stretching(Operator):
self.formulation = formulation
else:
self.formulation = CONSERVATIVE
- ## Do this operator deals with ghost points update?
- ## default = yes
+ ## If ghostupdate == True, the operator will manage
+ ## ghost points update. Else, this must be done somewhere else ...
+ ## Default = yes
if ghostUpdate is None:
self.ghostUpdate = True
else:
@@ -122,9 +125,8 @@ class Stretching(Operator):
def apply(self, simulation=None):
# update ghost points, if required, for velocity and vorticity
- #if self.synchronize is not None:
- # print 'synchr ...'
- self.synchronize.apply()
+ if self.synchronize is not None:
+ self.synchronize.apply()
# computation ...
self.discreteOperator.apply(simulation)
#if self.synchronize is not None:
diff --git a/HySoP/hysop/operator/synchronizeGhosts.py b/HySoP/hysop/operator/synchronizeGhosts.py
deleted file mode 100644
index 652a8e99ab68e4a9e9594e4ed5cc6a386210f69d..0000000000000000000000000000000000000000
--- a/HySoP/hysop/operator/synchronizeGhosts.py
+++ /dev/null
@@ -1,78 +0,0 @@
-"""
-@file operator/synchronizeGhosts.py
-
-Update ghost points for some fields defined on a specific topology.
-"""
-from parmepy.operator.continuous import Operator
-from parmepy.operator.discrete.synchronizeGhosts import SynchronizeGhosts_d
-from parmepy.constants import debug
-from parmepy.fields.vector import VectorField
-from parmepy.fields.scalar import ScalarField
-
-class SynchronizeGhosts(Operator):
- """
- Ghost points synchronization.
- """
-
- @debug
- def __init__(self, fieldslist, resolution=None,
- method='', topology=None, **other_config):
- """
- 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
- """
- if isinstance(fieldslist, list):
- Operator.__init__(self, fieldslist, method,
- name="Ghost Synchronization")
- else:
- Operator.__init__(self, [fieldslist], method,
- name="Ghost Synchronization")
- self.fieldslist = fieldslist
- self.resolution = resolution
- self.method = method
- self.topology = topology
- self.config = other_config
-
- @debug
- def setUp(self, ghosts=None):
- """
- SynchroGhost operator discretization method.
- Create an discrete SynchroGhost operator from given specifications.
-
- """
- Operator.setUp(self)
- for v in self.variables:
- # the topology for v ...
- if self.topology is not None:
- topo = self.topology
- else:
- if ghosts is None:
- topo = self.domain.getOrCreateTopology(
- self.domain.dimension,
- self.resolution[v])
- else:
- topo = self.domain.getOrCreateTopology(
- self.domain.dimension,
- self.resolution[v],
- ghosts=ghosts)
- topo.ghosts = ghosts
- # ... and the corresponding discrete field
- if v.__class__ is VectorField or v.__class__ is ScalarField:
- self.discreteOperator = \
- SynchronizeGhosts_d(v, topo, self.method)
- else:
- self.discreteFields[v] = v.discretize(topo)
- self.discreteOperator = \
- SynchronizeGhosts_d(self.discreteFields[v],
- topo, self.method)
-
- self.discreteOperator.setUp()
- self._isUpToDate = True
-
-
-if (__name__ == "__main__"):
- print __doc__
- print "- Provided class : SynchronizeGhosts"
- print SynchronizeGhosts.__doc__