From d19cc0cdc5babb1e1b1e4b0ba854505c3a7901b3 Mon Sep 17 00:00:00 2001
From: Chloe Mimeau <Chloe.Mimeau@imag.fr>
Date: Thu, 31 Jan 2013 15:08:37 +0000
Subject: [PATCH] about forces
---
HySoP/hysop/physics/compute_forces.py | 137 ++++++++----------
HySoP/hysop/test/test_operator/test_Forces.py | 14 +-
2 files changed, 68 insertions(+), 83 deletions(-)
diff --git a/HySoP/hysop/physics/compute_forces.py b/HySoP/hysop/physics/compute_forces.py
index ec9ebfa39..d49d92e2d 100644
--- a/HySoP/hysop/physics/compute_forces.py
+++ b/HySoP/hysop/physics/compute_forces.py
@@ -33,14 +33,12 @@ class Compute_forces(object):
self.res = topology.mesh.resolution
self.step = topology.mesh.size
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.coordMax = topology.mesh.end * self.step + topology.domain.origin - 2.* (topology.ghosts* self.step)
self.ghosts = topology.ghosts
self.periods = topology.periods
- print 'periods, ', self.periods
+# print 'res, dx, cMin, cMax, bmin, bMax', self.res, self.step, self.coordMin, self.coordMax, self.boxMin, self.boxMax
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) :
"""
Compute indicator functions for the control box (including the obstacle)
@@ -65,20 +63,18 @@ 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] -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]
+ tmp_ind[i][1] = self.res[i] -2. * self.ghosts[i] + distMax[i] // self.step[i] -1# Up, East, North
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] -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] -2.*self.ghosts[i] + distMax[i] // self.step[i] -1 * np.invert(self.periods[i])
+ tmp_ind[i][1] = self.res[i] -2.*self.ghosts[i] + distMax[i] // self.step[i] -1
if (distMin[i]<=0. and distMax[i]>=0.):
# the local domain is included inside the control volume
tmp_ind[i][1] = self.res[i] - 2. * self.ghosts[i] - 1
@@ -92,7 +88,7 @@ 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]
+# print 'ind_Down, ind_South', ind[self.Down], ind[self.Up]
# Remove last point to integrate properly ... ?????????
# ind[self.Up] = ind[self.Up] - 1
@@ -108,6 +104,8 @@ class Compute_forces(object):
nbPoints[self.South] = (ind[self.East] - ind[self.West] + 1) * (ind[self.Up] - ind[self.Down] + 1)
nbPoints[self.North] = (ind[self.East] - ind[self.West] + 1) * (ind[self.Up] - ind[self.Down] + 1)
+# print 'nb points plan Down', nbPoints[self.Down]
+
# if there is no volume control inside the local domain => nbPoints=0 => forces are not computed
if (distMin[0]<0. and distMax[0]<0. and self.boxMax[0]< self.coordMin[0]):
nbPoints[self.Down] = 0
@@ -128,94 +126,84 @@ class Compute_forces(object):
boxDim[2] = ind[self.North] - ind[self.South] + 1
nbPointsBox = boxDim[0] * boxDim[1] * boxDim[2]
- tmp_box = np.zeros([self.dim,nbPointsBox], dtype=PARMES_INTEGER)
- coords = np.zeros(self.dim, dtype=PARMES_INTEGER)
+# print 'nb points box', nbPointsBox
+
+ tmp_box = np.zeros([nbPointsBox,self.dim], dtype=PARMES_INTEGER)
+ coords = np.zeros(self.dim, dtype=PARMES_REAL)
count_box = 0
if(boxDim.all() >0):
- for k in xrange (ind[self.South],ind[self.North]): ## enlever les boucles !!
+ for k in xrange (ind[self.South],ind[self.North]+1): ## enlever les boucles !!
coords[2] = self.coordMin[2] + k * self.step[2]
- for j in xrange (ind[self.West],ind[self.East]):
+ for j in xrange (ind[self.West],ind[self.East]+1):
coords[1] = self.coordMin[1] + j * self.step[1]
- for i in xrange (ind[self.Down],ind[self.Up]):
+ for i in xrange (ind[self.Down],ind[self.Up]+1):
coords[0] = self.coordMin[0] + i * self.step[0]
dist = np.vdot(coords - self.obstacle.center,coords - self.obstacle.center) - self.obstacle.radius ** 2
if(dist >=0.0): # We are on or outside the sphere ...
+ tmp_box[count_box] = [i, j, k]
count_box = count_box + 1
- tmp_box[0][count_box] = i
- tmp_box[1][count_box] = j
- tmp_box[2][count_box] = k
+
+# print 'tmp_box', tmp_box, tmp_box.shape, count_box
# local indices of points in the control box which are outside the obstacle
- self.chi_box = tmp_box[:][0:count_box]
+ self.chi_box = np.zeros([count_box,self.dim], dtype=PARMES_INTEGER)
+ self.chi_box = tmp_box[0:count_box,:] ### PBM !!!
- self.chi_up = np.zeros([self.dim,nbPoints[self.Up]], dtype=PARMES_INTEGER)
- self.chi_down = np.zeros([self.dim,nbPoints[self.Down]], dtype=PARMES_INTEGER)
- self.chi_east = np.zeros([self.dim,nbPoints[self.East]], dtype=PARMES_INTEGER)
- self.chi_west = np.zeros([self.dim,nbPoints[self.West]], dtype=PARMES_INTEGER)
- self.chi_north = np.zeros([self.dim,nbPoints[self.North]], dtype=PARMES_INTEGER)
- self.chi_south = np.zeros([self.dim,nbPoints[self.South]], dtype=PARMES_INTEGER)
+ self.chi_up = np.zeros([nbPoints[self.Up],self.dim], dtype=PARMES_INTEGER)
+ self.chi_down = np.zeros([nbPoints[self.Down],self.dim], dtype=PARMES_INTEGER)
+ self.chi_east = np.zeros([nbPoints[self.East],self.dim], dtype=PARMES_INTEGER)
+ self.chi_west = np.zeros([nbPoints[self.West],self.dim], dtype=PARMES_INTEGER)
+ self.chi_north = np.zeros([nbPoints[self.North],self.dim], dtype=PARMES_INTEGER)
+ self.chi_south = np.zeros([nbPoints[self.South],self.dim], dtype=PARMES_INTEGER)
# local indices of points located on the surfaces of the control box
if(nbPoints[self.South]!=0.): # South boundary is in the domain
count = 0
- self.chi_south[2][:] = ind[self.South]
- for j in xrange (ind[self.West],ind[self.East]):
- for i in xrange (ind[self.Down],ind[self.Up]):
- self.chi_south[0][count] = i
- self.chi_south[1][count] = j
+ for j in xrange (ind[self.West],ind[self.East] + 1):
+ for i in xrange (ind[self.Down],ind[self.Up] + 1):
+ self.chi_south[count] = [i, j, ind[self.South]]
count = count + 1
if(nbPoints[self.East]!=0.): # East boundary is in the domain
count = 0
- self.chi_east[1][:] = ind[self.East]
- for k in xrange (ind[self.South],ind[self.North]):
- for i in xrange (ind[self.Down],ind[self.Up]):
- self.chi_east[0][count] = i
- self.chi_east[2][count] = k
+ for k in xrange (ind[self.South],ind[self.North] + 1):
+ for i in xrange (ind[self.Down],ind[self.Up] + 1):
+ self.chi_east[count] = [i, ind[self.East], k]
count = count + 1
if(nbPoints[self.Down]!=0.): # Down boundary is in the domain
count = 0
- self.chi_down[0][:] = ind[self.Down]
- for k in xrange (ind[self.South],ind[self.North]):
- for j in xrange (ind[self.West],ind[self.East]):
- self.chi_down[1][count] = j
- self.chi_down[2][count] = k
+ for k in xrange (ind[self.South],ind[self.North] + 1):
+ for j in xrange (ind[self.West],ind[self.East] + 1):
+ self.chi_down[count] = [ind[self.Down], j, k]
count = count + 1
if(nbPoints[self.North]!=0.): # North boundary is in the domain
count = 0
- self.chi_north[2][:] = ind[self.North]
- for j in xrange (ind[self.West],ind[self.East]):
- for i in xrange (ind[self.Down],ind[self.Up]):
- self.chi_north[0][count] = i
- self.chi_north[1][count] = j
+ for j in xrange (ind[self.West],ind[self.East] + 1):
+ for i in xrange (ind[self.Down],ind[self.Up] + 1):
+ self.chi_north[count] = [i, j, ind[self.North]]
count = count + 1
if(nbPoints[self.West]!=0.): # West boundary is in the domain
count = 0
- self.chi_west[1][:] = ind[self.West]
- for k in xrange (ind[self.South],ind[self.North]):
- for i in xrange (ind[self.Down],ind[self.Up]):
- self.chi_west[0][count] = i
- self.chi_west[2][count] = k
+ for k in xrange (ind[self.South],ind[self.North] + 1):
+ for i in xrange (ind[self.Down],ind[self.Up] + 1):
+ self.chi_west[count] = [i, ind[self.West], k]
count = count + 1
if(nbPoints[self.Up]!=0.): # Up boundary is in the domain
count = 0
- self.chi_up[0][:] = ind[self.Up]
- for k in xrange (ind[self.South],ind[self.North]):
- for j in xrange (ind[self.West],ind[self.East]):
- self.chi_up[1][count] = j
- self.chi_up[2][count] = k
+ for k in xrange (ind[self.South],ind[self.North] + 1):
+ for j in xrange (ind[self.West],ind[self.East] + 1):
+ self.chi_up[count] = [ind[self.Up], j, k]
count = count + 1
self.bufferForce = 0.
- # Compute coef used to compute the drag coefficient
+ # Compute coef used to determine drag coefficient : cD = coef.Fx = 2.Fx/rho.uinf**2.D
uinf = 1.
self.coef = 2. / (uinf ** 2 * PI * self.obstacle.radius ** 2)
-
def apply(self, t, dt, velo, vort, Re):
"""
Computation of the drag according to the "impulsion" formula presented by
@@ -223,7 +211,7 @@ 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) (no ghost points)
+ # computation of hessian matrices (d2_u1, d2_u2 and d2_u3) and jacobian matrix (d_u)
hessian1, hessian2, hessian3 = Fct2Op(vort, velo, choice = 'hessianU', topology=self.topo).apply(None,None)
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)
@@ -262,13 +250,13 @@ class Compute_forces(object):
# For all points in the box
for ind in xrange (chi.shape[1]):
# 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]
+ i = chi[ind][0] + self.ghosts[0]
+ j = chi[ind][1] + self.ghosts[1]
+ k = chi[ind][2] + self.ghosts[2]
# coordinates of the current point
- coords = self.coordMin + chi[:,ind] * self.step
+ coords = self.coordMin + chi[ind] * self.step
# part1 of the force
- int1 = int1 + np.vdot(coords,vort[l,m,n])
+ int1 = int1 + np.vdot(coords,vort[i,j,k])
force = force + fact * (int1 - self.bufferForce) # 1st order time integration
self.bufferForce = int1 # Save for next time step ...
@@ -288,28 +276,25 @@ class Compute_forces(object):
X_n_DivT = np.zeros(self.dim, dtype=PARMES_REAL)
n_T = np.zeros(self.dim, dtype=PARMES_REAL)
# For each point of the current plane
- for ind in xrange (chi.shape[1]):
- i = chi[0][ind]
- j = chi[1][ind]
- k = chi[2][ind]
+ for ind in xrange (chi.shape[0]):
# 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]
+ i = chi[ind][0] + self.ghosts[0]
+ j = chi[ind][1] + self.ghosts[1]
+ k = chi[ind][2] + 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[l,m,n],velo[l,m,n])
+ u_u = np.vdot(velo[i,j,k],velo[i,j,k])
# normal.velocity
- n_u = np.vdot(velo[l,m,n],NormalVec[:])
+ n_u = np.vdot(velo[i,j,k],NormalVec[:])
# normal.vorticity
- n_w = np.vdot(vort[l,m,n],NormalVec[:])
+ n_w = np.vdot(vort[i,j,k],NormalVec[:])
# coordinates of the current point
- coords = self.coordMin[:] + chi[:,ind] * self.step[:]
+ coords = self.coordMin[:] + chi[ind] * self.step[:]
# part1 of the force
- 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])
+ 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])
# 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)
diff --git a/HySoP/hysop/test/test_operator/test_Forces.py b/HySoP/hysop/test/test_operator/test_Forces.py
index 169bd80ff..d78333fef 100755
--- a/HySoP/hysop/test/test_operator/test_Forces.py
+++ b/HySoP/hysop/test/test_operator/test_Forces.py
@@ -54,9 +54,9 @@ class test_Forces(unittest.TestCase):
def testComputeForces(self):
# Parameters
- nb = 32
+ nb = 11
timeStep = 0.09
- finalTime = 0.18
+ finalTime = 0.36
self.t = 0.
t0 = time.time()
@@ -73,7 +73,7 @@ class test_Forces(unittest.TestCase):
vorti = pp.AnalyticalField(domain=box, formula=self.vorticite, name='Vorticity', vector=True)
## Solver creation (discretisation of objects is done in solver initialisation)
- topo3D = pp.CartesianTopology(domain=box, resolution=[nb, nb, nb], dim=3, ghosts=[2.,2.,2.])
+ topo3D = pp.CartesianTopology(domain=box, resolution=[nb, nb, nb], dim=3, ghosts=[2,2,2])
## Fields discretization
vorti.discretize(topo3D)
@@ -83,16 +83,16 @@ class test_Forces(unittest.TestCase):
# Forces computation
Re = 200.
- noca = Compute_forces(topo3D, sphere, boxMin= [0.25, 0.25, 0.25], boxMax=[0.75, 0.75, 0.75])
+ noca = Compute_forces(topo3D, sphere, boxMin= [0.2, 0.2, 0.2], boxMax=[0.8, 0.8, 0.8])
if (topo3D.rank == 0):
f = open('./parmepy/test/test_operator/NocaForces.dat', 'w')
while (self.t <= finalTime):
nocares = noca.apply(self.t, timeStep, velo.discreteField[velo._fieldId], vorti.discreteField[vorti._fieldId], Re)
if (topo3D.rank == 0):
-# f.write(t, force, force / self.coef, '\n')
- f.write(str(nocares))
- f.write('\n')
+ # print time and forces values in the following order : time, cX, cY, cZ
+ f.write("%s %s %s %s\n" % (self.t, nocares[0], nocares[1], nocares[2]))
+
self.t = self.t + timeStep
--
GitLab