diff --git a/HySoP/hysop/physics/compute_forces.py b/HySoP/hysop/physics/compute_forces.py
index 1751883e50dc27e75e6a30ac7bb365ddef723f96..ec9ebfa397611801b63e1e1271abcd4191a18310 100644
--- a/HySoP/hysop/physics/compute_forces.py
+++ b/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)