No commit message

No commit message

HySoP/hysop/physics/compute_forces.py

@@ -32,10 +32,14 @@ class Compute_forces(object):
         self.dim = topology.dim
         self.res = topology.mesh.resolution
         self.step = topology.mesh.size
-        self.coordMin = topology.mesh.origin
-        self.coordMax = (topology.mesh.end + 1) * topology.mesh.size + topology.domain.origin
+        self.coordMin = topology.mesh.origin + (topology.ghosts* self.step)
+        self.coordMax = (topology.mesh.end + 1) * self.step + topology.domain.origin - 2.* (topology.ghosts* self.step)
+        self.ghosts = topology.ghosts
+        self.periods = topology.periods
+        print 'periods, ', self.periods
         self.compute_control_box()
 
+        print 'res, dx, cMin, cMax, bmin, bMax', self.res, self.step, self.coordMin, self.coordMax, self.boxMin, self.boxMax
 
     def compute_control_box(self) :
         """
@@ -61,21 +65,23 @@ class Compute_forces(object):
 
         distMin = self.boxMin - self.coordMin
         distMax = self.boxMax - self.coordMax
+        print 'distMin, distMax', distMin, distMax
         for i in xrange(self.dim):
             if (distMin[i]>=0. and distMax[i]<=0.): 
             # the control volume is included inside the local domain
-                tmp_ind[i][0] = distMin[i] / self.step[i] # Down, West, South
-                tmp_ind[i][1] = self.res[i] + distMax[i] / self.step[i] - 1 # Up, East, North
+                tmp_ind[i][0] = distMin[i] // self.step[i] # Down, West, South
+                tmp_ind[i][1] = self.res[i] -2. * self.ghosts[i] + distMax[i] // self.step[i] -1 * np.invert(self.periods[i])# Up, East, North
+                print 'tmp_ind', tmp_ind[i][0], tmp_ind[i][1]
             if (distMin[i]>=0. and self.boxMin[i]<= self.coordMax[i] and distMax[i]>=0.): 
             # only min corner of control volume is included inside the local domain
-                tmp_ind[i][0] = distMin[i] / self.step[i] 
-                tmp_ind[i][1] = self.res[i] - 1
+                tmp_ind[i][0] = distMin[i] // self.step[i]
+                tmp_ind[i][1] = self.res[i] -2. * self.ghosts[i] - 1
             if (distMin[i]<=0. and self.boxMax[i]>= self.coordMin[i] and distMax[i]<=0.) : 
             # only max corner of control volume is included inside the local domain
-                tmp_ind[i][1] = self.res[i] + distMax[i]/self.step[i] - 1
+                tmp_ind[i][1] = self.res[i] -2.*self.ghosts[i] + distMax[i] // self.step[i] -1 * np.invert(self.periods[i])
             if (distMin[i]<=0. and distMax[i]>=0.): 
             # the local domain is included inside the control volume
-                tmp_ind[i][1] = self.res[i] - 1
+                tmp_ind[i][1] = self.res[i] - 2. * self.ghosts[i] - 1
 
             # otherwise, there is no volume control inside the local domain
 
@@ -86,11 +92,12 @@ class Compute_forces(object):
         ind[self.East] = tmp_ind[1][1]
         ind[self.South] = tmp_ind[2][0]
         ind[self.North] = tmp_ind[2][1]
+        print 'ind_Down, ind_South', ind[self.Down], ind[self.Up]
 
         # Remove last point to integrate properly ... ?????????
-        ind[self.Up] = ind[self.Up] - 1
-        ind[self.East] = ind[self.East] - 1
-        ind[self.North] = ind[self.North] - 1
+#        ind[self.Up] = ind[self.Up] - 1
+#        ind[self.East] = ind[self.East] - 1
+#        ind[self.North] = ind[self.North] - 1
     
         # Count the number of points on each face and inside the domain
         nbPoints = np.zeros(2 * self.dim, dtype=PARMES_INTEGER)
@@ -216,9 +223,9 @@ class Compute_forces(object):
         """
         localForce = np.zeros(self.dim, dtype=PARMES_REAL)
         force = np.zeros(self.dim, dtype=PARMES_REAL)
-        # computation of hessian matrices (d2_u1, d2_u2 and d2_u3) and jacobian matrix (d_u)
+        # computation of hessian matrices (d2_u1, d2_u2 and d2_u3) and jacobian matrix (d_u) (no ghost points)
         hessian1, hessian2, hessian3 = Fct2Op(vort, velo, choice = 'hessianU', topology=self.topo).apply(None,None)
-        self.jacobian, maxgersh = DifferentialOperator_d(vort, velo, choice='gradU', topology=self.topo).apply()
+        self.jacobian, maxgersh = DifferentialOperator_d(vort, velo, choice='gradV', topology=self.topo).apply()
         # computation of the components of Div(T) where T = 1/Re * (Nabla(U) + Nabla(U)T)
         self.a = 1. / Re * (2. * hessian1[0] + hessian1[4] + hessian2[1] + hessian1[8] + hessian3[2])
         self.b = 1. / Re * (hessian2[0] + hessian1[1] + 2. * hessian2[4] + hessian2[8] + hessian3[5])
@@ -254,13 +261,14 @@ class Compute_forces(object):
         int1 = np.zeros(self.dim, dtype=PARMES_REAL)
         # For all points in the box
         for ind in xrange (chi.shape[1]):
-            i = chi[0][ind]
-            j = chi[1][ind]
-            k = chi[2][ind]
+            # adjust indices for vector fields with ghost points
+            l = chi[0][ind] + self.ghosts[0]
+            m = chi[1][ind] + self.ghosts[1]
+            n = chi[2][ind] + self.ghosts[2]
             # coordinates of the current point
             coords = self.coordMin + chi[:,ind] * self.step
             # part1 of the force
-            int1 = int1 + np.vdot(coords,vort[i,j,k])
+            int1 = int1 + np.vdot(coords,vort[l,m,n])
 
         force = force + fact * (int1 - self.bufferForce) # 1st order time integration
         self.bufferForce = int1 # Save for next time step ...
@@ -284,20 +292,24 @@ class Compute_forces(object):
             i = chi[0][ind]
             j = chi[1][ind]
             k = chi[2][ind]
+            # adjust indices for vector fields with ghost points
+            l = chi[0][ind] + self.ghosts[0]
+            m = chi[1][ind] + self.ghosts[1]
+            n = chi[2][ind] + self.ghosts[2]
             # part1 = 1/2(velocity.velocity)n - (n.velocity)velocity - 1/(dim-1)((n.velocity)(coord X vorticity) + 1/(dim-1)(n.vorticity)(coord X velocity)
 
             # (velocity.velocity)
-            u_u = np.vdot(velo[i,j,k],velo[i,j,k])
+            u_u = np.vdot(velo[l,m,n],velo[l,m,n])
             # normal.velocity
-            n_u = np.vdot(velo[i,j,k],NormalVec[:])
+            n_u = np.vdot(velo[l,m,n],NormalVec[:])
             # normal.vorticity
-            n_w = np.vdot(vort[i,j,k],NormalVec[:])
+            n_w = np.vdot(vort[l,m,n],NormalVec[:])
             # coordinates of the current point
             coords = self.coordMin[:] + chi[:,ind] * self.step[:]
             # part1 of the force
-            int1 = int1 + 0.5 * u_u * NormalVec[:] - n_u * velo[i,j,k] \
-                - fact * n_u * np.cross(coords,vort[i,j,k]) \
-                + fact * n_w * np.cross(coords,velo[i,j,k])
+            int1 = int1 + 0.5 * u_u * NormalVec[:] - n_u * velo[l,m,n] \
+                - fact * n_u * np.cross(coords,vort[l,m,n]) \
+                + fact * n_w * np.cross(coords,velo[l,m,n])
 
             # part 2 of the force : part 2 = 1/(dim-1)(X^(n^Div(T)) + n.T
             if (chi.all() == self.chi_down.all()): # normalVec = (-1,0,0)