Skip to content
Snippets Groups Projects
Commit d19cc0cd authored by Chloe Mimeau's avatar Chloe Mimeau
Browse files

about forces

parent 689ad021
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
HySoP/hysop/physics/compute_forces.py
+ 61
− 76
View file @ d19cc0cd
......@@ -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)
......
HySoP/hysop/test/test_operator/test_Forces.py
+ 7
− 7
View file @ d19cc0cd
......@@ -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
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment