Update arrays operations in stretching time integration + momentally fix origin problem in printer

HySoP/hysop/numerics/differentialOperator.py

@@ -5,7 +5,7 @@ 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
+from parmepy.tools.cpu_data_transfer import Synchronize
 
 
 class DifferentialOperator(NumMethod):
@@ -20,8 +20,8 @@ class DifferentialOperator(NumMethod):
         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 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
         """
@@ -31,7 +31,7 @@ class DifferentialOperator(NumMethod):
         self.choice = choice
         self.method = method
         self.compute_time = 0.
-        self.topology = self.field2.topology
+        self.topology = self.field1.topology
         self.meshSize = self.topology.mesh.space_step
 
 #    def setUp(self):
@@ -271,14 +271,14 @@ class DifferentialOperator(NumMethod):
         elif self.choice.find('gradV') >= 0:
 
             # Ghosts synchronization
-#            OpSynchronize = Synchronize(self.topology)
-#            OpSynchronize.apply(self.field2)
+            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)
+                maxArray = np.zeros(5, dtype=PARMES_REAL, order=ORDER)
 
 ##              Fourth order scheme
 #               X components of temp and result
@@ -292,6 +292,8 @@ class DifferentialOperator(NumMethod):
                     1.0 * self.field2[0][ind0+2, ind1a:ind1b, ind2a:ind2b]
                 ) / (12. * self.meshSize[0])
 
+                max1 = np.max(abs(temp1))
+
 #               temp2 = np.zeros(
 #                   (self.resolution), dtype=PARMES_REAL, order=ORDER)
                 temp2 = self.field2[0][...] * 0.
@@ -302,6 +304,8 @@ class DifferentialOperator(NumMethod):
                     1.0 * self.field2[0][ind0a:ind0b, ind1+2, ind2a:ind2b]
                 ) / (12. * self.meshSize[1])
 
+                max2 = np.max(abs(temp2))
+
 #               temp3 = np.zeros(
 #                   (self.resolution), dtype=PARMES_REAL, order=ORDER)
                 temp3 = self.field2[0][...] * 0.
@@ -314,6 +318,7 @@ class DifferentialOperator(NumMethod):
 
                 maxstr1 = np.max(abs(temp1) + abs(temp2) + abs(temp3))
                 maxadv1 = np.max(abs(temp1))
+                max3 = np.max(abs(temp3))
 
 #               Y components of temp and result
 #               temp4 = np.zeros(
@@ -326,6 +331,8 @@ class DifferentialOperator(NumMethod):
                     1.0 * self.field2[1][ind0+2, ind1a:ind1b, ind2a:ind2b]
                 ) / (12. * self.meshSize[0])
 
+                max4 = np.max(abs(temp4))
+
 #               temp5 = np.zeros(
 #                   (self.resolution), dtype=PARMES_REAL, order=ORDER)
                 temp5 = self.field2[1][...] * 0.
@@ -336,6 +343,8 @@ class DifferentialOperator(NumMethod):
                     1.0 * self.field2[1][ind0a:ind0b, ind1+2, ind2a:ind2b]
                 ) / (12. * self.meshSize[1])
 
+                max5 = np.max(abs(temp5))
+
 #               temp6 = np.zeros(
 #                   (self.resolution), dtype=PARMES_REAL, order=ORDER)
                 temp6 = self.field2[1][...] * 0.
@@ -348,6 +357,7 @@ class DifferentialOperator(NumMethod):
 
                 maxstr2 = np.max(abs(temp4)+abs(temp5)+abs(temp6))
                 maxadv2 = np.max(abs(temp5))
+                max6 = np.max(abs(temp6))
 
 #               Z components of temp and result
 #               temp7 = np.zeros(
@@ -360,6 +370,8 @@ class DifferentialOperator(NumMethod):
                     1.0 * self.field2[2][ind0+2, ind1a:ind1b, ind2a:ind2b]
                 ) / (12. * self.meshSize[0])
 
+                max7 = np.max(abs(temp7))
+
 #               temp8 = np.zeros(
 #                   (self.resolution), dtype=PARMES_REAL, order=ORDER)
                 temp8 = self.field2[2][...] * 0.
@@ -370,6 +382,8 @@ class DifferentialOperator(NumMethod):
                     1.0 * self.field2[2][ind0a:ind0b, ind1+2, ind2a:ind2b]
                 ) / (12. * self.meshSize[1])
 
+                max8 = np.max(abs(temp8))
+
 #               temp9 = np.zeros(
 #                   (self.resolution), dtype=PARMES_REAL, order=ORDER)
                 temp9 = self.field2[2][...] * 0.
@@ -382,21 +396,40 @@ class DifferentialOperator(NumMethod):
 
                 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)
+                max9 = 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])
                 ))
-
-            return result, maxgersh
+                # div(V) result
+                divergence = np.array(temp1 + temp5 + temp9)
+
+                maxdiv = np.max(abs(divergence))
+                maxMax = max(max1, max2, max3,
+                             max4, max5, max6,
+                             max7, max8, max9)
+                maxproj = maxdiv / maxMax
+
+                # maxima of partial derivatives : needed for stab conditions
+                    # advection stab
+                maxArray[0] = max(maxstr1, maxstr2, maxstr3)
+                    # stretching stab
+                maxArray[1] = max(maxadv1, maxadv2, maxadv3)
+                    # solenoidal reprojection criterion
+                maxArray[2] = maxproj
+                maxArray[3] = maxdiv
+                maxArray[4] = maxMax
+
+            return result, maxArray, divergence
 
         elif (self.choice.find('gradS') >= 0):
 
             # Ghosts synchronization
-#            OpSynchronize = SynchronizeS(self.topology)
-#            OpSynchronize.apply(self.field2)
+            OpSynchronize = SynchronizeS(self.topology)
+            OpSynchronize.apply(self.field2)
 
             if self.method.find('FD_order2') >= 0:
                 raise ValueError("2nd order scheme Not yet implemented")
@@ -434,8 +467,8 @@ class DifferentialOperator(NumMethod):
         elif self.choice.find('laplacianV') >= 0:
 
             # Ghosts synchronization
-#            OpSynchronize = Synchronize(self.topology)
-#            OpSynchronize.apply(self.field2)
+            OpSynchronize = Synchronize(self.topology)
+            OpSynchronize.apply(self.field2)
 
             if self.method.find('FD_order2') >= 0:
                 raise ValueError("2nd order scheme Not yet implemented")
@@ -583,8 +616,8 @@ class DifferentialOperator(NumMethod):
         elif self.choice.find('curl') >= 0:
 
             # Ghosts synchronization
-#            OpSynchronize = Synchronize(self.topology)
-#            OpSynchronize.apply(self.field1)
+            OpSynchronize = Synchronize(self.topology)
+            OpSynchronize.apply(self.field1)
 
             if self.method.find('FD_order2') >= 0:
                 raise ValueError("2nd order scheme Not yet implemented")

HySoP/hysop/numerics/fct2op.py

@@ -71,8 +71,8 @@ class Fct2Op(object):
                                                             y[1][ind0a:ind0b,ind1a:ind1b,ind2a:ind2b] + \
                                                         self.field2[8][ind0a:ind0b,ind1a:ind1b,ind2a:ind2b]* \
                                                             y[2][ind0a:ind0b,ind1a:ind1b,ind2a:ind2b]
-            result = np.concatenate((np.array([temp0]),np.array([temp1]),np.array([temp2])))
-            return result
+#            result = np.concatenate((np.array([temp0]),np.array([temp1]),np.array([temp2])))
+            return temp0, temp1, temp2
 
         if self.choice == 'divWU' : # field2 = velocity field
             result = DifferentialOperator(y, self.field2, self.method, self.choice).__call__()

HySoP/hysop/numerics/integrators/euler.py

@@ -38,8 +38,14 @@ class Euler(ODESolver):
         @param t : current time.
         @param dt : time step.
 #        """
-        result = y + dt * np.asarray(f(t, y), dtype=PARMES_REAL, order=ORDER)
-        return result
+#        result = y + dt * np.asarray(f(t, y), dtype=PARMES_REAL, order=ORDER)
+#        return result
+
+        RHS = f(t, y)
+
+        return y[0] + dt * RHS[0], \
+               y[1] + dt * RHS[1], \
+               y[2] + dt * RHS[2]
 
     def __str__(self):
         """ToString method"""

HySoP/hysop/numerics/integrators/runge_kutta2.py

@@ -37,12 +37,22 @@ class RK2(ODESolver):
         @param y : function Y(t,space)
         """
 
-        # K1 = dt *f[t,y(t,x,y,z)]
-        K1 = dt * np.asarray(f(t, y), dtype=PARMES_REAL, order=ORDER)
-        # result = y(t,x,y,z)] + dt *f[t + dt/2 , y(t,x,y,z) + k1/2.]
-        result = y + dt * np.asarray(f(t + dt / 2., y + K1 / 2.),
-                                     dtype=PARMES_REAL, order=ORDER)
-        return result
+#        # K1 = dt *f[t,y(t,x,y,z)]
+#        K1 = dt * np.asarray(f(t, y), dtype=PARMES_REAL, order=ORDER)
+#        # result = y(t,x,y,z)] + dt *f[t + dt/2 , y(t,x,y,z) + k1/2.]
+
+#        result = y + dt * np.asarray(f(t + dt / 2., y + K1 / 2.),
+#                                     dtype=PARMES_REAL, order=ORDER)
+#        return result
+
+        RHS1 = f(t, y)
+        K1 = np.array([dt * RHS1[0], dt * RHS1[1], dt * RHS1[2]])
+
+        RHS2 = f(t + dt / 2., y + K1 / 2.)
+
+        return y[0] + dt * RHS2[0], \
+               y[1] + dt * RHS2[1], \
+               y[2] + dt * RHS2[2]
 
     def __str__(self):
         """ToString method"""

HySoP/hysop/numerics/integrators/runge_kutta3.py

@@ -36,16 +36,30 @@ class RK3(ODESolver):
         @param dt : time step.
         @param y : function Y(t,space)
         """
-        # K1 = dt *f[t,y(t,x,y,z)]
-        K1 = dt * np.asarray(f(t, y), dtype=PARMES_REAL, order=ORDER)
-        # result = y(t,x,y,z)] + dt *f[t + dt/2 , y(t,x,y,z) + k1/2.]
-        K2 = dt * np.asarray(f(t + dt / 3., y + K1[:, :, :, :] / 3.),
-                             dtype=PARMES_REAL, order=ORDER)
-        K3 = dt * np.asarray(f(t + dt * 2. / 3., y + K2[:, :, :, :] * 2. / 3.),
-                             dtype=PARMES_REAL, order=ORDER)
+#        # K1 = dt *f[t,y(t,x,y,z)]
+#        K1 = dt * np.asarray(f(t, y), dtype=PARMES_REAL, order=ORDER)
+#        # result = y(t,x,y,z)] + dt *f[t + dt/2 , y(t,x,y,z) + k1/2.]
+#        K2 = dt * np.asarray(f(t + dt / 3., y + K1[:, :, :, :] / 3.),
+#                             dtype=PARMES_REAL, order=ORDER)
+#        K3 = dt * np.asarray(f(t + dt * 2. / 3., y + K2[:, :, :, :] * 2. / 3.),
+#                             dtype=PARMES_REAL, order=ORDER)
+#        # result = y(t,x,y,z)] + 1/4 [ k1 + 3 k3]
+#        result = y + 1. / 4. * (K1 + 3. * K3)
+#        return result
+
+        RHS1 = f(t, y)
+        K1 = np.array([dt * RHS1[0], dt * RHS1[1], dt * RHS1[2]])
+
+        RHS2 = f(t + dt / 3., y + K1 / 3.)
+        K2 =  np.array([dt * RHS2[0], dt * RHS2[1], dt * RHS2[2]])
+
+        RHS3 = f(t + dt * 2. / 3., y + K2 * 2. / 3.)
+        K3 = np.array([dt * RHS3[0], dt * RHS3[1], dt * RHS3[2]])
+
         # result = y(t,x,y,z)] + 1/4 [ k1 + 3 k3]
-        result = y + 1. / 4. * (K1 + 3. * K3)
-        return result
+        return y[0] + 1. / 4. * (K1[0] + 3. * K3[0]), \
+               y[1] + 1. / 4. * (K1[1] + 3. * K3[1]), \
+               y[2] + 1. / 4. * (K1[2] + 3. * K3[2])
 
     def __str__(self):
         """ToString method"""

HySoP/hysop/numerics/integrators/runge_kutta4.py

@@ -38,18 +38,35 @@ class RK4(ODESolver):
         @param y : function Y(t,space)
         """
 
-        # K1 = dt *f[t, y(t, x, y, z)]
-        K1 = dt * np.asarray(f(t, y), dtype=PARMES_REAL, order=ORDER)
-        K2 = dt * np.asarray(f(t + dt / 2., y + K1[:, :, :, :] / 2.),
-                             dtype=PARMES_REAL, order=ORDER)
-        K3 = dt * np.asarray(f(t + dt / 2., y + K2[:, :, :, :] / 2.),
-                             dtype=PARMES_REAL, order=ORDER)
-        K4 = dt * np.asarray(f(t + dt, y + K3[:, :, :, :]),
-                             dtype=PARMES_REAL, order=ORDER)
+#        # K1 = dt *f[t, y(t, x, y, z)]
+#        K1 = dt * np.asarray(f(t, y), dtype=PARMES_REAL, order=ORDER)
+#        K2 = dt * np.asarray(f(t + dt / 2., y + K1[:, :, :, :] / 2.),
+#                             dtype=PARMES_REAL, order=ORDER)
+#        K3 = dt * np.asarray(f(t + dt / 2., y + K2[:, :, :, :] / 2.),
+#                             dtype=PARMES_REAL, order=ORDER)
+#        K4 = dt * np.asarray(f(t + dt, y + K3[:, :, :, :]),
+#                             dtype=PARMES_REAL, order=ORDER)
+#        # result = y(t, x, y, z)]  +  dt *f[t + dt/2 ,  y(t, x, y, z) + k1/2.]
+#        result = y + 1. / 6. * (K1[:, :, :, :] + 2. * K2[:, :, :, :] +
+#                                2. * K3[:, :, :, :] + K4[:, :, :, :])
+#        return result
+
+        RHS1 = f(t, y)
+        K1 = np.array([dt * RHS1[0], dt * RHS1[1], dt * RHS1[2]])
+
+        RHS2 = f(t + dt / 2., y + K1 / 1.)
+        K2 =  np.array([dt * RHS2[0], dt * RHS2[1], dt * RHS2[2]])
+
+        RHS3 = f(t + dt * 2., y + K2 * 2.)
+        K3 = np.array([dt * RHS3[0], dt * RHS3[1], dt * RHS3[2]])
+
+        RHS4 = f(t + dt, y + K3)
+        K4 = np.array([dt * RHS4[0], dt * RHS4[1], dt * RHS4[2]])
+
         # result = y(t, x, y, z)]  +  dt *f[t + dt/2 ,  y(t, x, y, z) + k1/2.]
-        result = y + 1. / 6. * (K1[:, :, :, :] + 2. * K2[:, :, :, :] +
-                                2. * K3[:, :, :, :] + K4[:, :, :, :])
-        return result
+        return y[0] + 1. / 6. * (K1[0] + 2. * K2[0] + 2. * K3[0] + K4[0]), \
+               y[1] + 1. / 6. * (K1[1] + 2. * K2[1] + 2. * K3[1] + K4[1]), \
+               y[2] + 1. / 6. * (K1[2] + 2. * K2[2] + 2. * K3[2] + K4[2])
 
     def __str__(self):
         """ToString method"""

HySoP/hysop/operator/discrete/stretching.py

@@ -60,9 +60,7 @@ class Stretching_d(DiscreteOperator):
         else:
             self.num_method = RK3
             self.cststretch = 2.5127
-        self.cststretch = self.cststretch / 2.
 
-    @debug
     @timed_function
     def apply(self, simulation):
         """
@@ -106,13 +104,13 @@ class Stretching_d(DiscreteOperator):
                     method='FD_order4',
                     choice='gradV')
 
-                gradResult, maxgersh = self.stretchOp()
+                gradResult, maxArray, divergence = self.stretchOp()
 
-                print "Maxgersh", maxgersh
-                print "dtAdvec", 0.5 / maxgersh[1]
+                print "maxadv, maxgersh", maxArray[1], maxArray[0]
+                print "dtAdvec", 0.5 / maxArray[1]
 
                 ## Determination of Stretching time step (stability condition)
-                self.ndt = self.stabilityTest(dt, maxgersh[0],
+                self.ndt = self.stabilityTest(dt, maxArray[0],
                                               self.num_method, self.cststretch)
                 self.dtstretch = dt / float(self.ndt)
 
@@ -127,9 +125,14 @@ class Stretching_d(DiscreteOperator):
             t = simulation.time
             for i in range(self.ndt):
                 methodInt = self.num_method(self.function)
-                self.vorticity.data[...] = methodInt(self.function.apply,
-                                                     t, self.dtstretch,
-                                                     self.vorticity.data)
+#                self.vorticity.data[...] = methodInt(self.function.apply,
+#                                                     t, self.dtstretch,
+#                                                     self.vorticity.data)
+                self.vorticity.data[0][...], \
+                self.vorticity.data[1][...], \
+                self.vorticity.data[2][...] = methodInt(self.function.apply,
+                                                        t, self.dtstretch,
+                                                        self.vorticity.data)
 #            print 'norm omega', LA.norm(self.vorticity.data)
 
     def stabilityTest(self, dt, maxi, num_method, cststretch):

HySoP/hysop/operator/monitors/printer.py

@@ -117,7 +117,8 @@ class Printer(Monitoring):
             print "Print to file " + filename
             ## VTK output \todo: Need fix in 2D, getting an IOError.
             try:
-                evtk.imageToVTK(filename, pointData=self._build_vtk_dict())
+                orig = self.variables[0].discreteFields.values()[0].topology.mesh.global_start
+                evtk.imageToVTK(filename, origin=(orig[0], orig[1], orig[2]), pointData=self._build_vtk_dict())
             except IOError:
                 if self.variables[0].domain.dimension == 2:
                     ## Standard output