Velocity penalization and Obstacle validation

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

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

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])

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"""

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__

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)

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()
 

HySoP/hysop/test/test_operator/test_Penalization.py

+# -*- 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()