Parallel velocity correction fixed

Examples/NS_bluff_bodies.py

@@ -43,13 +43,13 @@ sin = np.sin
 
 ## Domain
 #box = pp.Box(dim, length=[1., 1., 1.], origin=[0., 0., 0.])
-box = pp.Box(dim, length=[2.0 * pi, pi, pi], origin=[0., 0., 0.])
+box = pp.Box(dim, length=[4.0 * pi, 2.0 * pi, 2.0 * pi], origin=[-pi, -pi, -pi])
 #box = pp.Box(dim, length=[12., 10., 10.], origin=[-4., -5., -5.])
 #box = pp.Box(dim, length=[2.0 * pi, 2.0 * pi, 2.0 * pi])
 
 ## Global resolution
 #nbElem = [nb] * dim
-nbElem = [65, 33, 33]
+nbElem = [129, 65, 65]
 
 # Upstream flow velocity
 uinf = 1.0
@@ -68,6 +68,9 @@ def computeVel(res, x, y, z, t):
     res[0][...] = - ((uinf * x) / module) * (1 - (RADIUS ** 3 / module ** 3))
     res[1][...] = ((uinf * y) / module) * (1 + (RADIUS ** 3 / (2. * module ** 2)))
     res[2][...] = 0.
+#    res[0][...] = uinf
+#    res[1][...] = 0.
+#    res[2][...] = 0.
     return res
 
 # Function to compute initial vorticity
@@ -152,15 +155,6 @@ poisson = Poisson(velo, vorti,
                                vorti: nbElem},
                   projection=PROJ)
 
-correction = VelocityCorrection(velo, vorti,
-                                resolutions={velo: nbElem,
-                                             vorti: nbElem},
-                                uinf=uinf)
-
-penal = Penalization(velo, [sphere], 
-                     coeff=[1e8],
-                     resolutions={velo: nbElem})
-
 curl = Curl(velo, vorti,
             resolutions={velo: nbElem,
                          vorti: nbElem},
@@ -182,7 +176,7 @@ distrCurlAdv = Redistribute([vorti, velo], curl, advec)
 
 ## Simulation with fixed time step
 simu = Simulation(tinit=0.0,
-                  tend=5., timeStep=0.005,
+                  tend=10., timeStep=0.005,
                   iterMax=1000000)
 
 ## Define the problem to solve
@@ -196,20 +190,17 @@ simu = Simulation(tinit=0.0,
 #                          distrStrDiff,
 #                          diffusion, 
 #                          poisson],
-##                          correction,
-##                          penal],
 #               simulation=simu, dumpFreq=-1)
 
 pb = NSProblem(operators=[curl,
-                          distrAdvStr_velo, # Redistribute while advec
+                          distrCurlAdv,
+                          distrCurlStr_velo, # Redistribute while advec
                           advec,             # <--------------------:
                           distrAdvStr_vorti, 
                           stretch, 
                           distrStrDiff,
                           diffusion, 
-                          poisson,
-                          correction,
-                          penal],
+                          poisson],
                simulation=simu, dumpFreq=-1)
 
 ## Setting solver to Problem (only operators for computational tasks)
@@ -218,6 +209,18 @@ pb.pre_setUp()
 ## Topology without ghost points
 topofft = poisson.discreteFields[poisson.vorticity].topology
 
+## Add operators related to the problem and depending on FFT topology
+correction = VelocityCorrection(velo, vorti,
+                                resolutions={velo: nbElem,
+                                             vorti: nbElem},
+                                uinf=uinf,
+                                topo=topofft)
+
+penal = Penalization(velo, [sphere],
+                     coeff=[1e8],
+                     topo=topofft,
+                     resolutions={velo: nbElem})
+
 ## Diagnostics/monitors related to the problem
 
 forces = Forces(velo, vorti, topo, 
@@ -226,7 +229,7 @@ forces = Forces(velo, vorti, topo,
                 frequency=1, prefix='./res/Noca_sphere.dat')
 
 printer = Printer(variables=[velo, vorti],
-                  topo=topo,
+                  topo=topofft,
                   frequency=100,
                   prefix='./res/NS_sphere',
                   formattype=VTK)
@@ -234,14 +237,18 @@ printer = Printer(variables=[velo, vorti],
 energy = Energy_enstrophy(velo, vorti,
                           topo=topofft,
                           viscosity=VISCOSITY,
-                          isNormalized=False,
+                          isNormalized=True,
                           frequency=OUTPUT_FREQ,
                           prefix=FILENAME)
 
 ## Add monitors and setting solver to Problem
 vorti.setTopoInit(topofft)
 velo.setTopoInit(topofft)
-pb.addMonitors([energy, printer])
+pb.addMonitors([correction, penal, energy, printer])
+penal.discretize()
+correction.discretize()
+correction.setUp()
+penal.setUp()
 pb.setUp()
 
 penal.apply(simu)

Examples/dataNS_bb.py

@@ -50,11 +50,11 @@ PROJ = [False, 10]
 ## pi constant
 pi = np.pi
 # Obstacle description
-OBST_Ox = np.pi / 2.0
-OBST_Oy = np.pi / 2.0
-OBST_Oz = np.pi / 2.0
+OBST_Ox = 0.
+OBST_Oy = 0.
+OBST_Oz = 0.
 RADIUS = 0.5
 
 # post treatment ...
 OUTPUT_FREQ = 1
-FILENAME = './res/energy_NS_potFlow_res2_2.dat'
+FILENAME = './res/energy_NS_sphere_FFTcurl_8cpu.dat'

HySoP/hysop/operator/discrete/velocity_correction.py

@@ -50,16 +50,15 @@ class VelocityCorrection_D(DiscreteOperator):
 
     def setUp(self):
         spaceStep = self.velocity.topology.mesh.space_step
-        self.subLength = self.velocity.topology.mesh.end - \
-                         self.velocity.topology.mesh.origin
-        self.coeff_mean = np.prod(spaceStep) / np.prod(self.subLength)
+        lengths = self.velocity.topology.domain.length
+        self.coeff_mean = np.prod(spaceStep) / np.prod(lengths)
         self.x0_coord = self.velocity.topology.domain.origin[XDIR]
         self.x_coord = self.velocity.topology.mesh.coords[XDIR]
         if self.case is '2D':
-            self.surf = self.subLength[YDIR]
+            self.surf = lengths[YDIR]
             self.dvol = spaceStep[YDIR]
         elif self.case is '3D':
-            self.surf = self.subLength[YDIR] * self.subLength[ZDIR]
+            self.surf = lengths[YDIR] * lengths[ZDIR]
             self.dvol = spaceStep[YDIR] * spaceStep[ZDIR]
         else:
             raise ValueError("invalid problem dimension")
@@ -69,11 +68,12 @@ class VelocityCorrection_D(DiscreteOperator):
 
     @timed_function
     def apply(self, simulation):
-        # Calling for requirements completion
+        ## Calling for requirements completion
         DiscreteOperator.apply(self, simulation)
         ctime = MPI.Wtime()
 
-        # Computation of the flowrates evaluated from current (ie non corrected) velocity
+        ## Computation of the flowrates evaluated from 
+        ## current (ie non corrected) velocity
         local_comput_flowrates = npw.zeros((3))
         for i in xrange(self.velocity.nbComponents):
             local_comput_flowrates[i] = np.sum(self.velocity[i][...])
@@ -85,26 +85,45 @@ class VelocityCorrection_D(DiscreteOperator):
         for i in xrange(self.velocity.nbComponents):
             self.comput_flowrates[i] = recvbuff[i] * self.dvol
 
-        # Correction of the X-velocity component
+        ## Correction of the X-velocity component
         self.velocity[XDIR][...] += self.uinf - \
                                     self.comput_flowrates[XDIR] / self.surf
 
-        # Correction of the other velocity components
+        ## Correction of the other velocity components
         if self.case is '2D':
-            vort_mean = np.sum(self.vorticity[0][...]) * self.coeff_mean
+            # Computation of vorticity mean (in space)
+            local_vort_mean = np.sum(self.vorticity[0][...])
+            recvbuff = 0.0
+            self.velocity.topology.topo.Allreduce(local_vort_mean,
+                                                  recvbuff)
+            vort_mean = recvbuff * self.coeff_mean
+
+            # Correction of the Y-velocity component
             self.velocity[YDIR][...] += vort_mean * \
                                         (self.x_coord - self.x0_coord) - \
-                                        self.comput_flowrates[YDIR] / self.surf
+                                        self.comput_flowrates[YDIR] \
+                                        / self.surf
 
         elif self.case is '3D':
-            vort_Y_mean = np.sum(self.vorticity[YDIR][...]) * self.coeff_mean
-            vort_Z_mean = np.sum(self.vorticity[ZDIR][...]) * self.coeff_mean
-            self.velocity[YDIR][...] += vort_Z_mean * \
+            # Computation of Y and Z-vorticity means (in space)
+            local_vort_Y_mean = np.sum(self.vorticity[YDIR][...])
+            local_vort_Z_mean = np.sum(self.vorticity[ZDIR][...])
+
+            sendbuff = npw.zeros((2))
+            recvbuff = npw.zeros((2))
+            sendbuff[:] = [local_vort_Y_mean, local_vort_Z_mean]
+            self.velocity.topology.topo.Allreduce(sendbuff, recvbuff)
+            vort_mean = recvbuff * self.coeff_mean
+
+            # Correction of the Y and Z-velocity components
+            self.velocity[YDIR][...] += vort_mean[1] * \
                                         (self.x_coord - self.x0_coord) - \
-                                        self.comput_flowrates[YDIR] / self.surf
-            self.velocity[ZDIR][...] += vort_Y_mean * \
+                                        self.comput_flowrates[YDIR] \
+                                        / self.surf
+            self.velocity[ZDIR][...] += vort_mean[0] * \
                                         (self.x_coord - self.x0_coord) - \
-                                        self.comput_flowrates[ZDIR] / self.surf
+                                        self.comput_flowrates[ZDIR] \
+                                        / self.surf
 
         else:
             raise ValueError("invalid problem dimension")