From 477dccac2cc7cafd96094af9552a08726ff59980 Mon Sep 17 00:00:00 2001
From: Chloe Mimeau <Chloe.Mimeau@imag.fr>
Date: Fri, 8 Feb 2013 14:19:50 +0000
Subject: [PATCH]
---
HySoP/hysop/obstacle/obstacle_d.py | 15 ++++++++-------
HySoP/hysop/operator/differentialOperator_d.py | 4 ++--
HySoP/hysop/operator/penalization_d.py | 16 ++++++++++++++--
HySoP/hysop/particular_solvers/basic.py | 18 ++++++++++++++++++
.../hysop/test/test_obstacle/test_obstacle.py | 5 +++--
5 files changed, 45 insertions(+), 13 deletions(-)
diff --git a/HySoP/hysop/obstacle/obstacle_d.py b/HySoP/hysop/obstacle/obstacle_d.py
index a21bcd13f..f56899d70 100644
--- a/HySoP/hysop/obstacle/obstacle_d.py
+++ b/HySoP/hysop/obstacle/obstacle_d.py
@@ -68,17 +68,18 @@ class Obstacle_d(object):
chiPorousTmp = []
step = self.topology.mesh.size
+ ghosts = self.topology.ghosts
coord_start = self.topology.mesh.origin
- coord_end = self.topology.mesh.end * step + self.topology.domain.origin
- print 'start, end', coord_start, coord_end
- layerMin = coord_start[2] + self.zlayer
- layerMax = coord_end[2] - self.zlayer
+ coord_end = self.topology.mesh.end * step + self.topology.domain.origin - (ghosts * step)
+ layerMin = coord_start[2] + ghosts[2] * step[2] + self.zlayer
+ layerMax = coord_end[2] - ghosts[2] * step[2] - self.zlayer
radiusMinuslayer = self.radius- self.porousLayer
- for k in xrange (self.resolution[2]):
+ print 'start, end, layerMin, layerMax', coord_start, coord_end, layerMin, layerMax
+ for k in xrange (ghosts[2], self.resolution[2] - ghosts[2]):
rz = coord_start[2] + step[2] * k
- for j in xrange (self.resolution[1]):
+ for j in xrange (ghosts[1], self.resolution[1] - ghosts[1]):
ry = coord_start[1] + step[1] * j
- for i in xrange (self.resolution[0]):
+ for i in xrange (ghosts[0], self.resolution[0] - ghosts[0]):
if (rz > layerMax or rz < layerMin):
# we are in the z-layer (North or South):
chiBoundaryTmp.append([i,j,k])
diff --git a/HySoP/hysop/operator/differentialOperator_d.py b/HySoP/hysop/operator/differentialOperator_d.py
index c900aed22..0d2f90135 100755
--- a/HySoP/hysop/operator/differentialOperator_d.py
+++ b/HySoP/hysop/operator/differentialOperator_d.py
@@ -178,8 +178,8 @@ class DifferentialOperator_d(DiscreteOperator):
elif self.choice == 'gradV':
# Ghosts synchronization
-# OpSynchronize = Synchronize(self.topology)
-# OpSynchronize.apply(self.field2)
+ OpSynchronize = Synchronize(self.topology)
+ OpSynchronize.apply(self.field2)
maxgersh = np.zeros(2, dtype=PARMES_REAL, order=ORDER)
diff --git a/HySoP/hysop/operator/penalization_d.py b/HySoP/hysop/operator/penalization_d.py
index 7cbd1d0b0..24654fdda 100644
--- a/HySoP/hysop/operator/penalization_d.py
+++ b/HySoP/hysop/operator/penalization_d.py
@@ -6,7 +6,7 @@ Discrete penalization representation
from ..constants import *
from .discrete import DiscreteOperator
from .differentialOperator_d import DifferentialOperator_d
-#from ..obstacle.obstacle import Obstacle
+from ..physics.compute_forces import Compute_forces
import pyopencl as cl
import numpy as np
import time
@@ -42,12 +42,15 @@ class Penalization_d(DiscreteOperator):
# ## output field
# self.output = [self.res_velocity]
self.obstacle.chiFunctions()
+ self.noca = Compute_forces(self.topology, self.obstacle, boxMin= [0.2, 0.2, 0.2], boxMax=[0.8, 0.8, 0.8])
+ if (self.topology.rank == 0):
+ self.f = open('./res/NocaForces.dat', 'w')
def apply(self, t, dt):
self.compute_time = time.time()
-# self.obstacle.chiFunctions()
+ # self.obstacle.chiFunctions()
coef = 1.0 / (1.0 + dt * self.penal.lambd[:])
for k in self.obstacle.chiBoundary[:]:
@@ -66,7 +69,16 @@ class Penalization_d(DiscreteOperator):
self.velocity[2][k[0],k[1],k[2]] = coef[2] * self.velocity[2][k[0],k[1],k[2]]
self.vorticity.data = DifferentialOperator_d(self.velocity, self.vorticity, choice='curl', topology=self.topology).apply()
+
+ # Computation of forces exerced on the obstacle (according to Noca's formula)
+ Re = 200.
+ nocares = self.noca.apply(t, dt, self.velocity, self.vorticity, Re)
+ if (self.topology.rank == 0):
+ # print time and forces values in the following order : time, cX, cY, cZ
+ self.f.write("%s %s %s %s\n" % (t, nocares[0], nocares[1], nocares[2]))
+
self.compute_time = time.time() - self.compute_time
+ self.total_time += self.compute_time
return self.compute_time
diff --git a/HySoP/hysop/particular_solvers/basic.py b/HySoP/hysop/particular_solvers/basic.py
index de31f57ad..1fdbf3e70 100644
--- a/HySoP/hysop/particular_solvers/basic.py
+++ b/HySoP/hysop/particular_solvers/basic.py
@@ -14,6 +14,9 @@ from ..operator.transport import Transport
from ..operator.stretching import Stretching
from ..operator.remeshing import Remeshing
from ..operator.penalization import Penalization
+from ..operator.advec_scales import Advec_scales
+from ..operator.diffusion import Diffusion
+from ..operator.poisson import Poisson
from ..operator.splitting import Splitting
from ..tools.timer import Timer
from ..tools.printer import Printer
@@ -52,6 +55,9 @@ class ParticleSolver(Solver):
self.advection = None
self.stretch = None
self.penal = None
+ self.advec_scales = None
+ self.diffusion = None
+ self.poisson = None
## ODE Solver method.
self.ODESolver = ODESolver
## Interpolation method.
@@ -77,6 +83,12 @@ class ParticleSolver(Solver):
self.stretch = op
if isinstance(op, Penalization):
self.penal = op
+ if isinstance(op, Advec_scales):
+ self.advec_scales = op
+ if isinstance(op, Diffusion):
+ self.diffusion = op
+ if isinstance(op, Poisson):
+ self.poisson = op
#if isinstance(op, Velocity)
#self.velocity = op
#if self.advection is None:
@@ -109,6 +121,12 @@ class ParticleSolver(Solver):
if op is self.penal:
op.obstacle.discretize(self.problem.topology)
op.discretize()
+ if op is self.advec_scales:
+ op.discretize(topology=self.problem.topology)
+ if op is self.diffusion:
+ op.discretize(topology=self.problem.topology)
+ if op is self.poisson:
+ op.discretize(topology=self.problem.topology)
## Velocity is only program on gpu for instance
#elif op is self.velocity:
# op.discretize(result_position=self.p_position, result_scalar=self.p_scalar, method=self.ODESolver)
diff --git a/HySoP/hysop/test/test_obstacle/test_obstacle.py b/HySoP/hysop/test/test_obstacle/test_obstacle.py
index 94905b67e..7154dac26 100644
--- a/HySoP/hysop/test/test_obstacle/test_obstacle.py
+++ b/HySoP/hysop/test/test_obstacle/test_obstacle.py
@@ -11,10 +11,11 @@ from math import *
def run():
# Parameters
- nb = 11
+ nb = 65
timeStep = 0.02
finalTime = 1.
- outputFilePrefix = './parmepy/test/test_obstacle/Domain_'
+# outputFilePrefix = './parmepy/test/test_obstacle/Domain_'
+ outputFilePrefix = './res/Domain_'
outputModulo = 1
t0 = time.time()
--
GitLab