Skip to content
Snippets Groups Projects
Commit 77dc152e authored by Eric Dalissier's avatar Eric Dalissier
Browse files

No commit message 

No commit message
parent 52048ea0
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
HySoP/hysop/operator/fct2op.py
+ 11
− 2
View file @ 77dc152e
......@@ -39,17 +39,26 @@ class Fct2Op(object):
We used Fct2Op.apply like a def function
"""
if self.choice == 'gradU' :
if self.choice == 'gradU' : # field2 = gradU
#Make Scalar product
temp0 = self.field2[0][:,:,:]* y[0][:,:,:] +self.field2[1][:,:,:]* y[1][:,:,:]+self.field2[2][:,:,:]* y[2][:,:,:]
temp1 = self.field2[3][:,:,:]* y[0][:,:,:] +self.field2[4][:,:,:]* y[1][:,:,:]+self.field2[5][:,:,:]* y[2][:,:,:]
temp2 = self.field2[6][:,:,:]* y[0][:,:,:] +self.field2[7][:,:,:]* y[1][:,:,:]+self.field2[8][:,:,:]* y[2][:,:,:]
result = np.concatenate((np.array([temp0]),np.array([temp1]),np.array([temp2])))
return result
if self.choice == 'divWU' :
if self.choice == 'divWU' : # field2 = velocity field
print 'taille field22:', self.field2.shape
result = DifferentialOperator_d(y, self.field2, self.choice, self.topology).apply()
return result
if self.choice == 'hessianU' : # field2 = velocity field
gradU, maxgersh = DifferentialOperator_d(self.field1, self.field2, choice='gradU', topology=self.topology).apply()
GradUline1 = np.concatenate((gradU[0],gradU[1],gradU[2]))
GradUline2 = np.concatenate((gradU[3],gradU[4],gradU[5]))
GradUline3 = np.concatenate((gradU[6],gradU[7],gradU[8]))
result1, maxgersh = DifferentialOperator_d(self.field1, GradUline1, choice='gradU', topology=self.topology).apply()
result2, maxgersh = DifferentialOperator_d(self.field1, GradUline2, choice='gradU', topology=self.topology).apply()
result3, maxgersh = DifferentialOperator_d(self.field1, GradUline3, choice='gradU', topology=self.topology).apply()
return result1, result2, result3
else :
print 'choice', choice , 'is not define'
......
HySoP/hysop/physics/compute_forces.py
+ 56
− 27
View file @ 77dc152e
......@@ -5,6 +5,7 @@ Compute forces
"""
from ..constants import *
import numpy as np
from ..operator.fct2op import Fct2Op
import mpi4py.MPI as MPI
import time
......@@ -38,7 +39,7 @@ class Compute_forces(object):
"""
Compute indicator functions for the control box (including the obstacle)
"""
self.normal=np.zeros([self.dim*2, self.dim])
self.normal=np.zeros([self.dim*2, self.dim], dtype=PARMES_REAL)
self.Up = 0
self.Down = 1
self.East = 2
......@@ -52,7 +53,7 @@ class Compute_forces(object):
self.normal[self.North][2] = 1.0
self.normal[self.South][2] = -1.0
tmp_ind = np.zeros([self.dim, 2])
tmp_ind = np.zeros([self.dim, 2], dtype=PARMES_INTEGER)
distMin = self.boxMin - self.coordMin
distMax = self.boxMax - self.coordMax
......@@ -74,7 +75,7 @@ class Compute_forces(object):
# otherwise, there is no volume control inside the local domain
ind = np.zeros(self.dim*2)
ind = np.zeros(self.dim*2, dtype=PARMES_INTEGER)
ind[self.Down] = tmp_ind[0][0]
ind[self.Up] = tmp_ind[0][1]
ind[self.West] = tmp_ind[1][0]
......@@ -88,7 +89,7 @@ class Compute_forces(object):
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)
nbPoints = np.zeros(2*self.dim, dtype=PARMES_INTEGER)
nbPoints[self.Down] = (ind[self.East]-ind[self.West]+1)*(ind[self.North]-ind[self.South]+1)
nbPoints[self.Up] = (ind[self.East]-ind[self.West]+1)*(ind[self.North]-ind[self.South]+1)
nbPoints[self.West] = (ind[self.Up]-ind[self.Down]+1)*(ind[self.North]-ind[self.South]+1)
......@@ -110,13 +111,13 @@ class Compute_forces(object):
if (distMin[2]>0. and distMax[2]>0. and self.boxMin[2]> self.coordMax[2]):
nbPoints[self.North] = 0
boxDim = np.zeros(self.dim)
boxDim = np.zeros(self.dim, dtype=PARMES_INTEGER)
boxDim[0] = ind[self.Up] - ind[self.Down] + 1
boxDim[1] = ind[self.East] - ind[self.West] + 1
boxDim[2] = ind[self.North] - ind[self.South] + 1
nbPointsBox = boxDim[0] * boxDim[1] * boxDim[2]
tmp_box = np.zeros([self.dim,nbPointsBox])
tmp_box = np.zeros([self.dim,nbPointsBox], dtype=PARMES_INTEGER)
count_box=0
if(boxDim>0):
for k in xrange (ind[self.South],ind[self.North]): ## enlever les boucles !!
......@@ -134,12 +135,12 @@ class Compute_forces(object):
self.chi_box = tmp_box[:][0:count_box]
self.chi_up = np.zeros([self.dim,nbPoints[self.Up]])
self.chi_down = np.zeros([self.dim,nbPoints[self.Down]])
self.chi_east = np.zeros([self.dim,nbPoints[self.East]])
self.chi_west = np.zeros([self.dim,nbPoints[self.West]])
self.chi_north = np.zeros([self.dim,nbPoints[self.North]])
self.chi_south = np.zeros([self.dim,nbPoints[self.South]])
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)
if(nbPoints[self.South]!=0.): # South boundary is in the domain
count = 0
......@@ -197,7 +198,7 @@ class Compute_forces(object):
self.bufferForce = 0.
# Compute coef used to compute the drag coefficient
self.coef = 2./(uinf**2*np.pi*self.obstacle.radius**2)
self.coef = 2./(uinf**2*np.pi*self.obstacle.radius**2) # attention, uinf !!!!!!!!!
def apply(self, dt, velo, vort, Re):
......@@ -205,21 +206,28 @@ class Compute_forces(object):
Computation of the drag according to the "impulsion" formula presented by
Noca et. al in Journal of Fluids and Structures, 1999
"""
localForce = 0.
force = 0.
localForce = np.zeros(self.dim, dtype=PARMES_REAL)
force = np.zeros(self.dim, dtype=PARMES_REAL)
# computation of jacobian (du) and hessian (d2u) matrix
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.a = 2.*hessian1[0]+hessian1[4]+hessian2[1]+hessian1[8]+hessian3[2]
self.b = hessian2[0]+hessian1[1]+2.*hessian2[4]+hessian2[8]+hessian3[5]
self.c = hessian3[0]+hessian1[6]+hessian3[4]+hessian2[7]+2.*hessian3[8]
# Downstream and upstream surfaces
dsurf = step[1] * step[2]
localForce = integrateOnSurface(localForce,velo,vort,self.chi_down,self.normal[Down][:],X,Re,dsurf)
localForce = integrateOnSurface(localForce,velo,vort,self.chi_up,self.normal[Up][:],X,Re,dsurf)
localForce = integrateOnSurface(localForce,velo,vort,self.chi_down,self.normal[Down][:],Re,dsurf)
localForce = integrateOnSurface(localForce,velo,vort,self.chi_up,self.normal[Up][:],Re,dsurf)
# East and West surfaces
dsurf = step[0] * step[2]
localForce = integrateOnSurface(localForce,velo,vort,self.chi_east,self.normal[East][:],Y,Re,dsurf)
localForce = integrateOnSurface(localForce,velo,vort,self.chi_west,self.normal[West][:],Y,Re,dsurf)
localForce = integrateOnSurface(localForce,velo,vort,self.chi_east,self.normal[East][:],Re,dsurf)
localForce = integrateOnSurface(localForce,velo,vort,self.chi_west,self.normal[West][:],Re,dsurf)
# North and south surfaces
dsurf = step[0] * step[1]
localForce = integrateOnSurface(localForce,velo,vort,self.chi_south,self.normal[South][:],Z,Re,dsurf)
localForce = integrateOnSurface(localForce,velo,vort,self.chi_north,self.normal[North][:],Z,Re,dsurf)
localForce = integrateOnSurface(localForce,velo,vort,self.chi_south,self.normal[South][:],Re,dsurf)
localForce = integrateOnSurface(localForce,velo,vort,self.chi_north,self.normal[North][:],Re,dsurf)
# over the volume
dvol = step[0] * step[1] * step[2]
localForce = integrateOnBox(localForce,vort,self.chi_box,dvol,dt)
......@@ -236,7 +244,7 @@ class Compute_forces(object):
fact = -dvol/((self.dim-1)*dt)
int1 = 0.
# For all points in the box
for ind in xrange (int(chi[1].shape)): # ?????
for ind in xrange (chi.shape[1]):
i = chi[0][ind]
j = chi[0][ind]
k = chi[0][ind]
......@@ -248,7 +256,7 @@ class Compute_forces(object):
force = force + fact * (int1 - self.bufferForce)
self.bufferForce = int1 # Save for next time step ...
return force #??????????
return force
def integrateOnSurface(self,force,velo,vort,chi,NormalVec,direction,Re,dsurf):
"""
......@@ -259,10 +267,11 @@ class Compute_forces(object):
fact = 1./(self.dim-1)
# For each point of the current plane
int1 = 0.
int2 = 0.
diff_dir = 0.
for ind in xrange (int(chi[1].shape)): # ?????
int1 = np.zeros(self.dim, dtype=PARMES_REAL)
int2 = np.zeros(self.dim, dtype=PARMES_REAL)
X_n_DivT = np.zeros(self.dim, dtype=PARMES_REAL)
n_T = np.zeros(self.dim, dtype=PARMES_REAL)
for ind in xrange (chi.shape[1]):
i = chi[0][ind]
j = chi[0][ind]
k = chi[0][ind]
......@@ -281,6 +290,26 @@ class Compute_forces(object):
- 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, the one concerning T = nu(nabla u + nabla uT)
if (chi == self.chi_down):
X_n_DivT[0] = -coords[1]*self.b-coords[2]*self.c
X_n_DivT[1] = -coords[0]*self.b
X_n_DivT[2] = -coords[0]*self.c
n_T[0] = -1./Re*2*self.jacobian[0]
n_T[1] = -1./Re*self.jacobian[3]+self.jacobian[1]
n_T[2] = -1./Re*self.jacobian[6]+self.jacobian[2]
if (chi == self.chi_up):
X_n_DivT[0] = coords[1]*self.b+coords[2]*self.c
X_n_DivT[1] = coords[0]*self.b
X_n_DivT[2] = coords[0]*self.c
n_T[0] = 1./Re*2*self.jacobian[0]
n_T[1] = 1./Re*self.jacobian[3]+self.jacobian[1]
n_T[2] = 1./Re*self.jacobian[6]+self.jacobian[2]
int2 = int2 + X_n_DivT+n_T
# Product with element of surface and sum to the total (input) force
force = force + (int1 + int2) * dsurf
def __str__(self):
......
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