diff --git a/HySoP/hysop/__init__.py b/HySoP/hysop/__init__.py
index 6d822a99102b8e2fa460599faef11c89d03b2218..63f9a5f7d676374f3b8bab7fe4ba83d375facadf 100644
--- a/HySoP/hysop/__init__.py
+++ b/HySoP/hysop/__init__.py
@@ -49,6 +49,8 @@ Obstacle_d = obstacle.obstacle_d.Obstacle_d
import operator.transport
import operator.velocity
import operator.stretching
+import operator.penalization
+Penalization = operator.penalization.Penalization
Transport = operator.transport.Transport
Velocity = operator.velocity.Velocity
Stretching = operator.stretching.Stretching
diff --git a/HySoP/hysop/__init__.py.in b/HySoP/hysop/__init__.py.in
index d952fd4ebb476a77deee9d746bb383603b1a7608..b5634f2c18d6a858dbc803f388b27e9e3643fd51 100755
--- a/HySoP/hysop/__init__.py.in
+++ b/HySoP/hysop/__init__.py.in
@@ -49,6 +49,8 @@ Obstacle_d = obstacle.obstacle_d.Obstacle_d
import operator.transport
import operator.velocity
import operator.stretching
+import operator.penalization
+Penalization = operator.penalization.Penalization
Transport = operator.transport.Transport
Velocity = operator.velocity.Velocity
Stretching = operator.stretching.Stretching
diff --git a/HySoP/hysop/obstacle/obstacle_d.py b/HySoP/hysop/obstacle/obstacle_d.py
index ef38319af3b1017307f8104ccdf0653c5f87413f..804beed23c0df178dfc5dc800b68d9ed7d2b299a 100644
--- a/HySoP/hysop/obstacle/obstacle_d.py
+++ b/HySoP/hysop/obstacle/obstacle_d.py
@@ -103,7 +103,6 @@ class Obstacle_d(object):
# chiSolidTmp1.append(j)
# chiSolidTmp2.append(k)
else :
-# print "je suis passe ici, oups "
if (radiusMinuslayer < dist and dist <= self.radius and rx <= self.center[0]):
# we are in the porous area of the hemisphere:
chiPorousTmp.append([i,j,k])
diff --git a/HySoP/hysop/operator/penalization.py b/HySoP/hysop/operator/penalization.py
index 678b23bbb800ee7de32769a4c3400eb463f89106..a33e19bea05cbce33fe3b10e82c3110f6e8f5f6e 100644
--- a/HySoP/hysop/operator/penalization.py
+++ b/HySoP/hysop/operator/penalization.py
@@ -13,7 +13,7 @@ class Penalization(ContinuousOperator):
Penalization operator representation
"""
- def __init__(self, velocity, obstacle):
+ def __init__(self, velocity, obstacle, lambd):
"""
Constructor.
Create a Penalization operator from given velocity variavle.
@@ -25,16 +25,18 @@ class Penalization(ContinuousOperator):
self.obstacle = obstacle
## velocity variable to penalize
self.velocity = velocity
+ ## lambda penalization function
+ self.lambd = lambd
self.addVariable([velocity])
- def discretize(self, idVelocityD=0, result_velocity=None):
+ def discretize(self, idVelocityD=0, idObstacleD=0):
"""
Penalization operator discretization method.
Create a PenalizationDOp.PenalizationDOp from given specifications.
"""
if self.discreteOperator is None:
- self.discreteOperator = Penalization_d(self, idVelocityD, result_velocity)
+ self.discreteOperator = Penalization_d(self, idVelocityD, idObstacleD)
def __str__(self):
"""ToString method"""
diff --git a/HySoP/hysop/operator/penalization_d.py b/HySoP/hysop/operator/penalization_d.py
index a9405fa3307e664107549f7e957597861044277f..b56d7b45e0caf88f097c7066575d89c98bbbbeb0 100644
--- a/HySoP/hysop/operator/penalization_d.py
+++ b/HySoP/hysop/operator/penalization_d.py
@@ -5,18 +5,19 @@ Discrete penalization representation
"""
from ..constants import *
from .discrete import DiscreteOperator
+#from ..obstacle.obstacle import Obstacle
import pyopencl as cl
import numpy as np
import time
-class PenalizationDOp(DiscreteOperator):
+class Penalization_d(DiscreteOperator):
"""
Penalization operator representation.
DiscreteOperator.DiscreteOperator specialization.
"""
- def __init__(self, penal, idVelocityD=0, result_velocity=None):
+ def __init__(self, penal, idVelocityD=0, idObstacleD=0):
"""
Constructor.
Create a Penalization operator on a given continuous domain.
@@ -25,42 +26,55 @@ class PenalizationDOp(DiscreteOperator):
@param penal : Penalization operator.
"""
DiscreteOperator.__init__(self)
+ ## Continous penalization operator
+ self.penal = penal
## Velocity.
self.velocity = penal.velocity.discreteField[idVelocityD]
- ## Result position
- self.res_velocity = result_velocity
+ ## Obstacle represented by characteristic function Chi
+ self.obstacle = self.penal.obstacle.discreteObstacle[idObstacleD]
+# ## Topology
+# self.topology = self.velocity.topology
## input field
- self.input = [self.velocity]
- ## output field
- self.output = [self.res_velocity]
+# self.input = [self.velocity]
+# ## output field
+# self.output = [self.res_velocity]
+ self.obstacle.chiFunctions()
- def apply(self, t, dt, lambd, myObstacle, topology): # faut-il vraiment mettre lambd, myObstacle, topology comme arguments ici ? NON
+ def apply(self, t, dt):
- myObstacle.discreteObstacles(topology)
- myObstacle.chiFunctions()
- coef = 1.0 / (1.0 + dt * lambd)
+ self.compute_time = time.time()
- for k in xrange (myObstacle.ptsInBoundary):
- self.velocity[0][myObstacle.chiBoundary[0][k],myObstacle.chiBoundary[1][k],myObstacle.chiBoundary[2][k]] =
- coef * self.velocity[0][myObstacle.chiBoundary[0][k],myObstacle.chiBoundary[1][k],myObstacle.chiBoundary[2][k]]
- self.velocity[1][myObstacle.chiBoundary[0][k],myObstacle.chiBoundary[1][k],myObstacle.chiBoundary[2][k]] =
- coef * self.velocity[1][myObstacle.chiBoundary[0][k],myObstacle.chiBoundary[1][k],myObstacle.chiBoundary[2][k]]
- self.velocity[2][myObstacle.chiBoundary[0][k],myObstacle.chiBoundary[1][k],myObstacle.chiBoundary[2][k]] =
- coef * self.velocity[2][myObstacle.chiBoundary[0][k],myObstacle.chiBoundary[1][k],myObstacle.chiBoundary[2][k]]
+# self.obstacle.chiFunctions()
+ coef = 1.0 / (1.0 + dt * self.penal.lambd[:])
- for k in xrange (myObstacle.ptsInSolid):
- self.velocity[0][myObstacle.chiSolid[0][k],myObstacle.chiSolid[1][k],myObstacle.chiSolid[2][k]] =
- coef * self.velocity[0][myObstacle.chiSolid[0][k],myObstacle.chiSolid[1][k],myObstacle.chiSolid[2][k]]
- self.velocity[1][myObstacle.chiSolid[0][k],myObstacle.chiSolid[1][k],myObstacle.chiSolid[2][k]] =
- coef * self.velocity[1][myObstacle.chiSolid[0][k],myObstacle.chiSolid[1][k],myObstacle.chiSolid[2][k]]
- self.velocity[2][myObstacle.chiSolid[0][k],myObstacle.chiSolid[1][k],myObstacle.chiSolid[2][k]] =
- coef * self.velocity[2][myObstacle.chiSolid[0][k],myObstacle.chiSolid[1][k],myObstacle.chiSolid[2][k]]
+ for k in self.obstacle.chiBoundary[:]:
+ self.velocity[0][k[0],k[1],k[2]] = coef[2] * self.velocity[0][k[0],k[1],k[2]]
+ self.velocity[1][k[0],k[1],k[2]] = coef[2] * self.velocity[1][k[0],k[1],k[2]]
+ self.velocity[2][k[0],k[1],k[2]] = coef[2] * self.velocity[2][k[0],k[1],k[2]]
+
+ for k in self.obstacle.chiPorous[:]:
+ self.velocity[0][k[0],k[1],k[2]] = coef[1] * self.velocity[0][k[0],k[1],k[2]]
+ self.velocity[1][k[0],k[1],k[2]] = coef[1] * self.velocity[1][k[0],k[1],k[2]]
+ self.velocity[2][k[0],k[1],k[2]] = coef[1] * self.velocity[2][k[0],k[1],k[2]]
+
+ for k in self.obstacle.chiSolid[:]:
+ self.velocity[0][k[0],k[1],k[2]] = coef[2] * self.velocity[0][k[0],k[1],k[2]]
+ self.velocity[1][k[0],k[1],k[2]] = coef[2] * self.velocity[1][k[0],k[1],k[2]]
+ self.velocity[2][k[0],k[1],k[2]] = coef[2] * self.velocity[2][k[0],k[1],k[2]]
+
+ return self.compute_time
+
+ def printComputeTime(self):
+ self.timings_info[0] = "\"Penalization total\""
+ self.timings_info[1] = str(self.total_time)
+ print "Penalization total time : ", self.total_time
+ print "Time of the last Penalization iteration :", self.compute_time
def __str__(self):
- s = "PenalizationDOp (DiscreteOperator). " + DiscreteOperator.__str__(self)
+ s = "Penalization_d (DiscreteOperator). " + DiscreteOperator.__str__(self)
return s
if __name__ == "__main__":
print __doc__
- print "- Provided class : PenalizationDOp"
- print PenalizationDOp.__doc__
+ print "- Provided class : Penalization_d"
+ print Penalization_d.__doc__
diff --git a/HySoP/hysop/particular_solvers/basic.py b/HySoP/hysop/particular_solvers/basic.py
index e694ecabc75f16c57bde3472abffccefed8053f3..fb62e489083d60e271ae2eca2442da3deb64ecd0 100644
--- a/HySoP/hysop/particular_solvers/basic.py
+++ b/HySoP/hysop/particular_solvers/basic.py
@@ -4,7 +4,6 @@
Particular solver description.
"""
from solver import Solver
-#from integrator.runge_kutta4 import RK4
from .integrator.integrator import ODESolver
from .integrator.euler import Euler
from .integrator.runge_kutta2 import RK2
@@ -14,6 +13,7 @@ from ..fields.continuous import ContinuousField
from ..operator.transport import Transport
from ..operator.stretching import Stretching
from ..operator.remeshing import Remeshing
+from ..operator.penalization import Penalization
from ..operator.splitting import Splitting
from ..tools.timer import Timer
from ..tools.printer import Printer
@@ -50,6 +50,8 @@ class ParticleSolver(Solver):
self.problem = problem
## Advection operator of the problem. Advection need special treatment in particle methods.
self.advection = None
+ self.stretch = None
+ self.penal = None
## ODE Solver method.
self.ODESolver = ODESolver
## Interpolation method.
@@ -73,6 +75,8 @@ class ParticleSolver(Solver):
self.advection = op
if isinstance(op, Stretching):
self.stretch = op
+ if isinstance(op, Penalization):
+ self.penal = op
#if isinstance(op, Velocity)
#self.velocity = op
#if self.advection is None:
@@ -100,8 +104,11 @@ class ParticleSolver(Solver):
for op in self.problem.operators:
if op is self.advection:
op.discretize(result_position=self.p_position, result_scalar=self.p_scalar, method=self.ODESolver)
- elif op is self.stretch:
+ if op is self.stretch:
op.discretize(method=self.ODESolver)
+ if op is self.penal:
+ op.obstacle.discretize(self.problem.topology)
+ op.discretize()
## 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 00e18d7ae2a5837cb986b8dc3d4c5c35d24b9e2f..1c60c2feba7d49d61d7e8e1459d87fe9d4ce7732 100644
--- a/HySoP/hysop/test/test_obstacle/test_obstacle.py
+++ b/HySoP/hysop/test/test_obstacle/test_obstacle.py
@@ -11,10 +11,10 @@ from math import *
def run():
# Parameters
- nb = 50
+ nb = 128
timeStep = 0.02
finalTime = 1.
- outputFilePrefix = './parmepy/test/test_domain/Domain_'
+ outputFilePrefix = './parmepy/test/test_obstacle/Domain_'
outputModulo = 1
t0 = time.time()
@@ -37,15 +37,21 @@ def run():
sphere.discretize(topo3D)
sphereD = sphere.discreteObstacle[0]
sphereD.chiFunctions()
- for k in xrange (topo3D.mesh.resolution[2]):
- for j in xrange (topo3D.mesh.resolution[2]):
- for i in xrange (topo3D.mesh.resolution[2]):
- if ([i,j,k] in sphereD.chiBoundary) :
- chiDomainD[i,j,k]=1.
- if ([i,j,k] in sphereD.chiSolid) :
- chiDomainD[i,j,k]=1.
- if ([i,j,k] in sphereD.chiPorous) :
- chiDomainD[i,j,k]=0.5
+ for x in sphereD.chiBoundary[:] :
+ chiDomainD[x[0], x[1], x[2]]=1.
+ for x in sphereD.chiSolid[:] :
+ chiDomainD[x[0], x[1], x[2]]=1.
+ for x in sphereD.chiPorous[:] :
+ chiDomainD[x[0], x[1], x[2]]=0.5
+# for k in xrange (topo3D.mesh.resolution[2]):
+# for j in xrange (topo3D.mesh.resolution[2]):
+# for i in xrange (topo3D.mesh.resolution[2]):
+# if ([i,j,k] in sphereD.chiBoundary) :
+# chiDomainD[i,j,k]=1.
+# if ([i,j,k] in sphereD.chiSolid) :
+# chiDomainD[i,j,k]=1.
+# if ([i,j,k] in sphereD.chiPorous) :
+# chiDomainD[i,j,k]=0.5
io=pp.Printer(fields=[chiDomain], frequency=outputModulo, outputPrefix=outputFilePrefix)
io.step()
diff --git a/HySoP/hysop/test/test_operator/test_Penalization.py b/HySoP/hysop/test/test_operator/test_Penalization.py
new file mode 100644
index 0000000000000000000000000000000000000000..62b552a730a8cc0f8aaf13a33b93a9f73afaa783
--- /dev/null
+++ b/HySoP/hysop/test/test_operator/test_Penalization.py
@@ -0,0 +1,67 @@
+# -*- coding: utf-8 -*-
+import unittest
+import time
+import parmepy as pp
+import numpy as np
+from parmepy.particular_solvers.integrator.runge_kutta4 import RK4
+import numpy.testing as npt
+from parmepy.constants import *
+from math import *
+
+def vitesse(x, y, z):
+ vx = 2.
+ vy = 2.
+ vz = 2.
+ return vx, vy, vz
+
+def run():
+ # Parameters
+ nb = 128
+ timeStep = 0.09
+ finalTime = 0.09
+ outputFilePrefix = './parmepy/test/test_operator/Penalization_'
+ outputModulo = 1
+
+ t0 = time.time()
+
+ ## Domain
+ box = pp.Box(3, length=[1., 1., 1.], origin=[0., 0., 0.])
+
+ ## Obstacle
+ lambd=np.array([0, 1, 10**8],dtype = PARMES_REAL, order=ORDER)
+ sphere = pp.Obstacle(box, name='sphere', zlayer=0.1, radius=0.1, center=[0.5,0.5,0.5], orientation='West', porousLayer=0.05)
+
+ ## ChiDomain
+ chiDomain = pp.ContinuousField(domain=box, name='ChiDomain', vector=False)
+
+ ## Fields
+ velo = pp.AnalyticalField(domain=box, formula=vitesse, name='Velocity', vector=True)
+
+ ## Operators
+ penal = pp.Penalization(velo, sphere, lambd)
+
+ ## Solver creation (discretisation of objects is done in solver initialisation)
+ topo3D = pp.CartesianTopology(domain=box, resolution=[nb, nb, nb], dim=3)
+
+ pb = pp.Problem(topo3D, [penal])
+
+ ## Setting solver to Problem
+ pb.setSolver(finalTime, timeStep, solver_type='basic', io=pp.Printer(fields=[velo], frequency=outputModulo, outputPrefix=outputFilePrefix))
+ pb.solver.ODESolver= RK4# RK3# RK2# Euler#
+ pb.initSolver()
+ t1 = time.time()
+
+ ## Solve problem
+ pb.solve()
+
+ tf = time.time()
+
+
+ print "\n"
+ print "Total time : ", tf - t0, "sec (CPU)"
+ print "Init time : ", t1 - t0, "sec (CPU)"
+ print "Solving time : ", tf - t1, "sec (CPU)"
+
+
+if __name__ == "__main__":
+ run()