diff --git a/HySoP/hysop/numerics/differentialOperator.py b/HySoP/hysop/numerics/differentialOperator.py
new file mode 100755
index 0000000000000000000000000000000000000000..5693647139d5142c33e14b41e40de0cb07545696
--- /dev/null
+++ b/HySoP/hysop/numerics/differentialOperator.py
@@ -0,0 +1,460 @@
+# -*- codingind utf-8 -*-
+"""
+@package operator
+Discrete stretching representation
+"""
+from parmepy.constants import np, PARMES_REAL, ORDER
+from parmepy.numerics.method import NumMethod
+from parmepy.tools.cpu_data_transfer import Synchronize
+import sys
+import time
+
+
+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 n°1 (vorticity).
+ @param field2 : field n°2 (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):
+ """
+ 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)
+
+ 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 = (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 = (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 = (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])
+
+# Y components of temp and result
+ temp1 = (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 = (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 = (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])
+
+# Z components of temp and result
+ temp1 = (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 = (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 = (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])
+ result = np.concatenate((np.array(tmp1), np.array(tmp2), np.array(tmp3)))
+ return result
+
+ elif self.choice.find('gradV') >= 0:
+
+ # Ghosts synchronization
+ OpSynchronize = Synchronize(self.topology)
+ OpSynchronize.apply(self.field2)
+
+ if self.method.find('FD_order2') >= 0:
+ raise ValueError("2nd order scheme Not yet implemented")
+
+ else:
+ maxgersh = np.zeros(2, 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))
+ maxgersh[0] = max(maxstr1,maxstr2,maxstr3)
+ maxgersh[1] = max(maxadv1,maxadv2,maxadv3)
+ result = np.concatenate((np.array([temp1, temp2, temp3]),np.array([temp4, temp5, temp6]),np.array([temp7, temp8, temp9])))
+
+ return result, maxgersh
+
+
+ elif self.choice.find('gradS') >= 0:
+
+ # Ghosts synchronization
+ OpSynchronize = SynchronizeS(self.topology)
+ OpSynchronize.apply(self.field2)
+
+ 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)
+
+ 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)
+
+ 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_d (DiscreteOperator). " + DiscreteOperator.__str__(self)
+ return s
+
+if __name__ == "__main__":
+ print __doc__
+ print "- Provided class ind DifferentialOperator_d"
+ print DifferentialOperator_d.__doc__
diff --git a/HySoP/hysop/operator/fct2op.py b/HySoP/hysop/numerics/fct2op.py
similarity index 67%
rename from HySoP/hysop/operator/fct2op.py
rename to HySoP/hysop/numerics/fct2op.py
index 5fac6363eec8b56d7df102da2e2e54e7037e13b1..f6390def06013d1abe67824f8b5e11bef5dca838 100644
--- a/HySoP/hysop/operator/fct2op.py
+++ b/HySoP/hysop/numerics/fct2op.py
@@ -4,34 +4,30 @@
Convert differential operator for stretching into a function to Discrete time integration
"""
-from discrete import DiscreteOperator
-from differentialOperator_d import DifferentialOperator_d
+from parmepy.numerics.differentialOperator import DifferentialOperator
import numpy as np
class Fct2Op(object):
"""
Fct2Op
- make the link between differentialOperator_d and the integration time scheme (Euler,RK,...)
+ make the link between differentialOperator and the integration time scheme (Euler,RK,...)
"""
- def __init__(self, field1, field2, choice=None, topology=None):
+ 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 w.Du or divWU for div(wu)
- @param resolution : resolution of the integration domain (array of number of points)
- @param meshSize : array of the space step
+ @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
- self.resolution = self.topology.mesh.resolution
- self.meshSize = self.topology.mesh.size
def apply(self, t , y):
"""
@@ -43,16 +39,18 @@ class Fct2Op(object):
temp0 = y[0][...] * 0.
temp1 = y[0][...] * 0.
temp2 = y[0][...] * 0.
- ind0a = self.topology.ghosts[0]
- ind0b=self.resolution[0]-self.topology.ghosts[0]
- ind1a = self.topology.ghosts[1]
- ind1b=self.resolution[1]-self.topology.ghosts[1]
- ind2a = self.topology.ghosts[2]
- ind2b=self.resolution[2]-self.topology.ghosts[2]
- ind0=np.arange(ind0a,ind0b)
- ind1=np.arange(ind1a,ind1b)
- ind2=np.arange(ind2a,ind2b)
- #Make Scalar product
+
+ 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]* \
@@ -73,21 +71,22 @@ class Fct2Op(object):
y[2][ind0a:ind0b,ind1a:ind1b,ind2a:ind2b]
result = np.concatenate((np.array([temp0]),np.array([temp1]),np.array([temp2])))
return result
+
if self.choice == 'divWU' : # field2 = velocity field
- print 'taille field22:', self.field2.shape
- result = DifferentialOperator_d(y, self.field2, self.choice, self.topology).apply()
+ result = DifferentialOperator(y, self.field2, self.method, self.choice).__call__()
return result
+
if self.choice == 'hessianU' : # field2 = velocity field
- gradU, maxgersh = DifferentialOperator_d(self.field1, self.field2, choice='gradV', topology=self.topology).apply()
+ gradU, maxgersh = DifferentialOperator(self.field1, self.field2, choice='gradV', topology=self.topology).__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_d(self.field1, GradUline1, choice='gradV', topology=self.topology).apply()
- result2, maxgersh = DifferentialOperator_d(self.field1, GradUline2, choice='gradV', topology=self.topology).apply()
- result3, maxgersh = DifferentialOperator_d(self.field1, GradUline3, choice='gradV', topology=self.topology).apply()
+ result1, maxgersh = DifferentialOperator(self.field1, GradUline1, choice='gradV', topology=self.topology).__call__()
+ result2, maxgersh = DifferentialOperator(self.field1, GradUline2, choice='gradV', topology=self.topology).__call__()
+ result3, maxgersh = DifferentialOperator(self.field1, GradUline3, choice='gradV', topology=self.topology).__call__()
return result1, result2, result3
else :
- print 'choice', choice , 'is not defined'
+ print 'choice', choice , 'in Fct2Op is not defined'
def __str__(self):
diff --git a/HySoP/hysop/operator/continuous.py b/HySoP/hysop/operator/continuous.py
index 1ce5ca8b869df44205a351921fddf22f6020d61b..bbb873794a63d2c58f8ab78f21da79e4e8c794b1 100644
--- a/HySoP/hysop/operator/continuous.py
+++ b/HySoP/hysop/operator/continuous.py
@@ -23,7 +23,7 @@ class Operator(object):
## Domain
self.domain = self.variables[0].domain
## Fields discretes id
- self.dicreteFieldId = {}
+ self.discreteFieldId = {}
## Operator discretization.
self.discreteOperator = None
## Timings informations
diff --git a/HySoP/hysop/operator/penalization.py b/HySoP/hysop/operator/penalization.py
index 4f45229b391b5df7c6cf8671921a7f47b3042a13..bb8c171b179a4437f368cb046ae7da28b2d4f75b 100644
--- a/HySoP/hysop/operator/penalization.py
+++ b/HySoP/hysop/operator/penalization.py
@@ -5,6 +5,7 @@
Penalization operator representation
"""
from parmepy.operator.continuous import Operator
+from parmepy.mpi.topology import Cartesian
from parmepy.operator.penalization_d import Penalization_d
@@ -42,6 +43,14 @@ class Penalization(Operator):
"""
Operator.setUp(self)
+ print self.method
+ #variables discretization
+ for v in self.variables:
+ topoId, topo = self.domain.addTopology(
+ Cartesian(domain=self.domain, dim=self.domain.dimension,
+ globalMeshResolution=self.resolutions[v]))
+ df, dfId = v.discretize(topo)
+ self.discreteFieldId[v] = dfId
self.discreteOperator = Penalization_d(self, method=self.method,
**self.config)
self.discreteOperator.setUp()
diff --git a/HySoP/hysop/operator/penalization_d.py b/HySoP/hysop/operator/penalization_d.py
index 2aae7651b18ad98d768608b4414c2db57bc35d1b..e0b636de643aab12f39f228699c7ee6270814a27 100644
--- a/HySoP/hysop/operator/penalization_d.py
+++ b/HySoP/hysop/operator/penalization_d.py
@@ -4,8 +4,7 @@
Discrete penalization representation
"""
from parmepy.operator.discrete import DiscreteOperator
-from parmepy.operator.differentialOperator_d import DifferentialOperator_d
-#from parmepy.physics.compute_forces import Compute_forces
+from parmepy.numerics.differentialOperator import DifferentialOperator
from parmepy.constants import np
import time
@@ -24,31 +23,23 @@ class Penalization_d(DiscreteOperator):
@param penal : Penalization operator.
"""
- DiscreteOperator.__init__(self)
+ DiscreteOperator.__init__(self, penal)
## Continous penalization operator
self.penal = penal
## Velocity.
- self.velocity = self.penal.velocity.getDiscreteField(
- self.penal.resolutions[self.penal.velocity])
+ self.velocity = self.penal.velocity.discreteField[
+ self.penal.discreteFieldId[self.penal.velocity]]
## Vorticity.
- self.vorticity = self.penal.vorticity.getDiscreteField(
- self.penal.resolutions[self.penal.vorticity])
- ## Obstacle represented by characteristic function Chi
+ self.vorticity = self.penal.vorticity.discreteField[
+ self.penal.discreteFieldId[self.penal.vorticity]]
+ ## Obstacle
self.obstacle = penal.obstacle
-# ## Topology
- self.topology = self.obstacle.topology
+ self.obstacle.setUp()
+ ## Method
self.method = method
- ## input field
-# self.input = [self.velocity]
-# ## output field
-# self.output = [self.res_velocity]
- self.obstacle.chiFunctions()
-# self.noca = Compute_forces(self.topology, self.obstacle, boxMin= [0.2, 0.2, 0.2], boxMax=[0.8, 0.8, 0.8])
-# if (self.topology.rank == 0):
-# self.f = open('./res/NocaForces.dat', 'w')
def setUp(self):
- pass
+ self.obstacle.definedObst.setUp(self.velocity.topology)
def apply(self, t, dt, ite):
@@ -56,32 +47,38 @@ class Penalization_d(DiscreteOperator):
coef = 1.0 / (1.0 + dt * self.penal.lambd[:])
- for k in self.obstacle.chiBoundary[:]:
- #-------- attempt loop suppression :
-# self.velocity[0][np.asarray(self.obstacle.chiBoundary[:])] = coef[2] * self.velocity[0][np.asarray(self.obstacle.chiBoundary[:])]
- self.velocity[0][k[0],k[1],k[2]] = coef[2] * self.velocity[0][k[0],k[1],k[2]]
- self.velocity[1][k[0],k[1],k[2]] = coef[2] * self.velocity[1][k[0],k[1],k[2]]
- self.velocity[2][k[0],k[1],k[2]] = coef[2] * self.velocity[2][k[0],k[1],k[2]]
+ b = self.obstacle.definedObst.chiBoundary
+ s = self.obstacle.definedObst.chiSolid
+ p = self.obstacle.definedObst.chiPorous
- for k in self.obstacle.chiPorous[:]:
- self.velocity[0][k[0],k[1],k[2]] = coef[1] * self.velocity[0][k[0],k[1],k[2]]
- self.velocity[1][k[0],k[1],k[2]] = coef[1] * self.velocity[1][k[0],k[1],k[2]]
- self.velocity[2][k[0],k[1],k[2]] = coef[1] * self.velocity[2][k[0],k[1],k[2]]
+ if self.velocity.dimension == 2 :
+ self.velocity[0][b[0], b[1]] = coef[0] * self.velocity[0][b[0], b[1]]
+ self.velocity[1][b[0], b[1]] = coef[0] * self.velocity[1][b[0], b[1]]
- for k in self.obstacle.chiSolid[:]:
- self.velocity[0][k[0],k[1],k[2]] = coef[2] * self.velocity[0][k[0],k[1],k[2]]
- self.velocity[1][k[0],k[1],k[2]] = coef[2] * self.velocity[1][k[0],k[1],k[2]]
- self.velocity[2][k[0],k[1],k[2]] = coef[2] * self.velocity[2][k[0],k[1],k[2]]
+ self.velocity[0][s[0], s[1]] = coef[0] * self.velocity[0][s[0], s[1]]
+ self.velocity[1][s[0], s[1]] = coef[0] * self.velocity[1][s[0], s[1]]
+
+ self.velocity[0][p[0], p[1]] = coef[1] * self.velocity[0][p[0], p[1]]
+ self.velocity[1][p[0], p[1]] = coef[1] * self.velocity[1][p[0], p[1]]
+
+ if self.velocity.dimension == 3 :
+ self.velocity[0][b[0], b[1], b[2]] = coef[0] * self.velocity[0][b[0], b[1], b[2]]
+ self.velocity[1][b[0], b[1], b[2]] = coef[0] * self.velocity[1][b[0], b[1], b[2]]
+ self.velocity[2][b[0], b[1], b[2]] = coef[0] * self.velocity[2][b[0], b[1], b[2]]
+
+
+ self.velocity[0][s[0], s[1], s[2]] = coef[0] * self.velocity[0][s[0], s[1], s[2]]
+ self.velocity[1][s[0], s[1], s[2]] = coef[0] * self.velocity[1][s[0], s[1], s[2]]
+ self.velocity[2][s[0], s[1], s[2]] = coef[0] * self.velocity[2][s[0], s[1], s[2]]
+
+
+ self.velocity[0][p[0], p[1], p[2]] = coef[1] * self.velocity[0][p[0], p[1], p[2]]
+ self.velocity[1][p[0], p[1], p[2]] = coef[1] * self.velocity[1][p[0], p[1], p[2]]
+ self.velocity[2][p[0], p[1], p[2]] = coef[1] * self.velocity[2][p[0], p[1], p[2]]
#-------- vorticity penalization (curl on grid) :
- self.vorticity.data = DifferentialOperator_d(self.velocity, self.vorticity, choice='curl', topology=self.topology).apply()
+# self.vorticity.data = DifferentialOperator(self.velocity, self.vorticity, choice='curl', topology=self.velocity.topology).apply()
- #---------Computation of forces exerced on the obstacle (according to Noca's formula)
-# Re = 200.
-# nocares = self.noca.apply(t, dt, self.velocity, self.vorticity, Re)
-# if (self.topology.rank == 0):
-# # print time and forces values in the following order : time, cX, cY, cZ
-# self.f.write("%s %s %s %s\n" % (t, nocares[0], nocares[1], nocares[2]))
self.compute_time = time.time() - self.compute_time
self.total_time += self.compute_time
@@ -89,8 +86,8 @@ class Penalization_d(DiscreteOperator):
return self.compute_time
def printComputeTime(self):
- self.timings_info[0] = "\"Penalization total\""
- self.timings_info[1] = str(self.total_time)
+# self.timings_info[0] = "\"Penalization total\""
+# self.timings_info[1] = str(self.total_time)
print "Penalization total time : ", self.total_time
print "Time of the last Penalization iteration :", self.compute_time
diff --git a/HySoP/hysop/operator/stretching.py b/HySoP/hysop/operator/stretching.py
index 0ce2fda7f911be93c2477dfd2144b61405b10340..6337a152b2861f8498e98845c5d7e70e5ad12fda 100755
--- a/HySoP/hysop/operator/stretching.py
+++ b/HySoP/hysop/operator/stretching.py
@@ -7,6 +7,7 @@ Stretching operator representation
"""
from parmepy.operator.continuous import Operator
+from parmepy.mpi.topology import Cartesian
from parmepy.operator.stretching_d import Stretching_d
@@ -16,8 +17,7 @@ class Stretching(Operator):
"""
def __init__(self, velocity, vorticity,
- resolutions=None, method='RK3',
- propertyOp='divConservation',
+ resolutions=None, method='FD_order4 RK3',
**other_config):
"""
Constructor.
@@ -43,11 +43,29 @@ class Stretching(Operator):
@param idVelocityD : Index of velocity discretisation to use.
@param idVorticityD : Index of vorticity discretisation to use.
- @param propertyOp : choice of the property garantied by
- the stretching Op ( divConservation or massConservation)
- @param method : the method to use. (Euler, RK, ..) default : RK3
+ @param propertyOp : property garantied by the stretching Op
+ (divConservation, massConservation) default : divConservation
+ @param method : the method to use :
+ space : (FD_order2, FD_order4, ..) default : FD_order4
+ and time : (Euler, RK, ..) default : RK3
"""
Operator.setUp(self)
+ print self.method
+ nbGhosts = 0
+ if self.method.find('FD_order2') >= 0:
+ nbGhosts = 1
+ elif self.method.find('FD_order4') >= 0:
+ nbGhosts = 2
+ else:
+ nbGhosts = 2
+ #variables discretization
+ for v in self.variables:
+ topoId, topo = self.domain.addTopology(
+ Cartesian(domain=self.domain, dim=self.domain.dimension,
+ globalMeshResolution=self.resolutions[v],
+ ghosts=[nbGhosts, nbGhosts, nbGhosts]))
+ df, dfId = v.discretize(topo)
+ self.discreteFieldId[v] = dfId
self.discreteOperator = Stretching_d(self, method=self.method,
**self.config)
self.discreteOperator.setUp()
diff --git a/HySoP/hysop/operator/stretching_d.py b/HySoP/hysop/operator/stretching_d.py
index ee99e4bff966b18a455f581bc743d560044c3950..ccdd63286a1509729c8b8edb3b780cbb6954b449 100755
--- a/HySoP/hysop/operator/stretching_d.py
+++ b/HySoP/hysop/operator/stretching_d.py
@@ -6,8 +6,8 @@ Discrete stretching representation
"""
from parmepy.operator.discrete import DiscreteOperator
-from parmepy.operator.differentialOperator_d import DifferentialOperator_d
-from fct2op import Fct2Op
+from parmepy.numerics.differentialOperator import DifferentialOperator
+from parmepy.numerics.fct2op import Fct2Op
from parmepy.numerics.integrators.euler import Euler
from parmepy.numerics.integrators.runge_kutta2 import RK2
from parmepy.numerics.integrators.runge_kutta3 import RK3
@@ -21,34 +21,46 @@ class Stretching_d(DiscreteOperator):
DiscreteOperator.DiscreteOperator specialization.
"""
- def __init__(self, stretch, method='RK3', propertyOp='divConservation'):
+ def __init__(self, stretch, method='FD_order4 RK3', propertyOp='divConservation'):
"""
Constructor.
Create a Stretching operator on a given continuous domain.
Work on a given velocity and vorticity to compute the stretching term.
@param stretch : Stretching operator.
+ @param method : Method to use
+ @param propertyOp : Operator property to satisfy (divConservation or massConservation)
"""
- DiscreteOperator.__init__(self)
+ DiscreteOperator.__init__(self, stretch)
self.stretch = stretch
## Velocity.
- self.velocity = self.stretch.velocity.getDiscreteField(
- self.stretch.resolutions[self.stretch.velocity])
+ self.velocity = self.stretch.velocity.discreteField[
+ self.stretch.discreteFieldId[self.stretch.velocity]]
## Vorticity.
- self.vorticity = self.stretch.vorticity.getDiscreteField(
- self.stretch.resolutions[self.stretch.vorticity])
+ self.vorticity = self.stretch.vorticity.discreteField[
+ self.stretch.discreteFieldId[self.stretch.vorticity]]
# ## input fields
# self.input = [self.velocity, self.vorticity]
self.propertyOp = propertyOp
self.method = method
def setUp(self):
- if self.method == 'RK3':
+ self.cststretch = 0.
+ if self.method.find('Euler') >= 0:
+ self.num_method = Euler
+ self.cststretch = 2.0
+ 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.name = "stretching"
- #on recupere la topologie des variables discretes
- self.topology = self.velocity.topology
- self.topology = self.vorticity.topology
+ self.cststretch = 2.5127
+ elif self.method.find('RK4') >= 0:
+ self.num_method = RK4
+ self.cststretch = 2.7853
+ else:
+ self.num_method = RK3
+ self.cststretch = 2.5127
def apply(self, t, dt, ite):
"""
@@ -56,6 +68,7 @@ class Stretching_d(DiscreteOperator):
@param t : current time.
@param dt : time step.
+ @param ite : current iteration.
Stretching algorithm:
@li 1. discretization of the divergence operator (div wu)
@@ -69,69 +82,61 @@ class Stretching_d(DiscreteOperator):
if self.num_method is not None:
- if self.propertyOp == 'divConservation':
+ if (self.propertyOp.find('divConservation') < 0
+ and self.propertyOp.find('massConservation') < 0):
+ print 'No or wrong stretching operator property specified, use default divConservation'
+ self.propertyOp = 'divConservation'
- ## Space discretisation of grad(u)
- self.stretchOp = DifferentialOperator_d(
- self.vorticity, self.velocity, choice='gradV',
- topology=self.topology)
+ if self.propertyOp == 'massConservation':
+ ## Determination of Stretching time step (stability condition)
+ self.ndt = 1
+ self.dtstretch = dt
- gradResult, maxgersh = self.stretchOp.apply()
+ ## fct2Op
+ self.function = Fct2Op(self.vorticity,
+ self.velocity,
+ choice='divWU',
+ topology=self.velocity.topology)
- ## Determination of Stretching time step (stability condition)
- self.ndt = self.stabilityTest(dt, maxgersh[0], self.num_method)
- self.dtstretch = dt / float(self.ndt)
+ else: # self.propertyOp == 'divConservation'
- print "Maxgersh", maxgersh
- print "dtStretch, ndt", self.dtstretch, self.ndt
- print "dtAdvec", 0.5/maxgersh[1]
+ ## Space discretization of grad(u)
+ self.stretchOp = DifferentialOperator(
+ self.vorticity, self.velocity,
+ method='FD_order4',
+ choice='gradV')
- ## fct2Op
- self.fonction = Fct2Op(self.vorticity.data,
- gradResult, choice='gradV',
- topology=self.topology)
+ gradResult, maxgersh = self.stretchOp()
+
+ print "Maxgersh", maxgersh
+ print "dtAdvec", 0.5 / maxgersh[1]
- elif self.propertyOp == 'massConservation':
## Determination of Stretching time step (stability condition)
- self.ndt = 1
- self.dtstretch = dt
+ self.ndt = self.stabilityTest(dt, maxgersh[0], self.num_method, self.cststretch)
+ self.dtstretch = dt / float(self.ndt)
## fct2Op
- self.fonction = Fct2Op(self.vorticity,
- self.velocity, choice='divWU',
- topology=self.topology)
-
- else:
- raise ValueError("Unknown Stretching Operator property: "
- + str(self.propertyOp))
+ self.function = Fct2Op(self.vorticity.data,
+ gradResult,
+ method='FD_order4',
+ choice='gradV',
+ topology=self.vorticity.topology)
# Time integration
for i in range(self.ndt):
- methodInt = self.num_method(self.fonction)
- self.vorticity.data = methodInt.integrate(self.fonction.apply,
- t, self.dtstretch,
- self.vorticity.data)
+ methodInt = self.num_method(self.function)
+ self.vorticity.data = methodInt(self.function.apply,
+ t, self.dtstretch,
+ self.vorticity.data)
self.compute_time = time.time() - self.compute_time
self.total_time += self.compute_time
return self.compute_time
- def stabilityTest(self, dt, maxi, num_method):
+ def stabilityTest(self, dt, maxi, num_method, cststretch):
dtstretch = dt
- cststretch = 0.
- if isinstance(num_method, Euler):
- cststretch = 2.0
- dtstretch = min(dtstretch, cststretch / maxi)
- if isinstance(num_method, RK2):
- cststretch = 2.0
- dtstretch = min(dtstretch, cststretch / maxi)
- if isinstance(num_method, RK3):
- cststretch = 2.5127
- dtstretch = min(dtstretch, cststretch / maxi)
- if isinstance(num_method, RK4):
- cststretch = 2.7853
- dtstretch = min(dtstretch, cststretch / maxi)
+ dtstretch = min(dtstretch, cststretch / maxi)
if dt == dtstretch:
print "dt, ndt , dtstretch, upperbound", dt, int(dt / dtstretch),\
dt / float(int(dt / dtstretch)), cststretch / maxi
@@ -142,8 +147,8 @@ class Stretching_d(DiscreteOperator):
return int(dt / dtstretch) + 1
def printComputeTime(self):
- self.timings_info[0] = "\"Stretching total\""
- self.timings_info[1] = str(self.total_time)
+# self.timings_info[0] = "\"Stretching total\""
+# self.timings_info[1] = str(self.total_time)
print "Stretching total time : ", self.total_time
print "Time of the last Stretching iteration :", self.compute_time
diff --git a/HySoP/hysop/operator/differentialOperator.py b/HySoP/hysop/unusedOrObsolet/differentialOperator.py
similarity index 100%
rename from HySoP/hysop/operator/differentialOperator.py
rename to HySoP/hysop/unusedOrObsolet/differentialOperator.py
diff --git a/HySoP/hysop/operator/differentialOperator_d.py b/HySoP/hysop/unusedOrObsolet/differentialOperator_d.py
similarity index 96%
rename from HySoP/hysop/operator/differentialOperator_d.py
rename to HySoP/hysop/unusedOrObsolet/differentialOperator_d.py
index 72ef8ddb2fb99c48a662eed1e25f92964bc7320a..1cf369491e454d6ab73425702c671a5930ce2208 100755
--- a/HySoP/hysop/operator/differentialOperator_d.py
+++ b/HySoP/hysop/unusedOrObsolet/differentialOperator_d.py
@@ -35,29 +35,21 @@ class DifferentialOperator_d(DiscreteOperator):
#La topologie est recuperee des variables discretes
self.topology = self.field1.topology
self.topology = self.field2.topology
- self.resolution = self.topology.mesh.resolution
- self.meshSize = self.topology.mesh.size
-
+ self.meshSize = self.topology.mesh.space_step_size
def apply(self):
"""
Apply operator.
"""
-# if False in [(self.field1[i][ind,ind].shape == self.field2[i][ind,ind].shape) for i in xrange(self.field2.topology.dim)] ind
-# print "Error, not the same topology for field1 and field2"
-# sys.exit(0)
-# elseind
-# mesh_size = self.field2.topology.mesh.size
-# resolution = self.field2.topology.mesh.resolution
- ind0a = self.topology.ghosts[0]
- ind0b=self.resolution[0]-self.topology.ghosts[0]
- ind1a = self.topology.ghosts[1]
- ind1b=self.resolution[1]-self.topology.ghosts[1]
- ind2a = self.topology.ghosts[2]
- ind2b=self.resolution[2]-self.topology.ghosts[2]
- ind0=np.arange(ind0a,ind0b)
- ind1=np.arange(ind1a,ind1b)
- ind2=np.arange(ind2a,ind2b)
+ 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 == 'divWU':
## second order scheme