From 2d793b78c8841432fd7fc28218d95884b9a5e86d Mon Sep 17 00:00:00 2001
From: Eric Dalissier <Eric.Dalissier@imag.fr>
Date: Mon, 17 Dec 2012 13:32:12 +0000
Subject: [PATCH] Add vorticity penalization, and penalization for NS
---
Examples/NavierStokes3dStretch.py | 142 ++++++++++++++++++
.../hysop/operator/differentialOperator_d.py | 1 +
HySoP/hysop/operator/penalization.py | 10 +-
HySoP/hysop/operator/penalization_d.py | 8 +-
.../test/test_operator/test_Penalization.py | 12 +-
5 files changed, 166 insertions(+), 7 deletions(-)
create mode 100644 Examples/NavierStokes3dStretch.py
diff --git a/Examples/NavierStokes3dStretch.py b/Examples/NavierStokes3dStretch.py
new file mode 100644
index 000000000..a8b300e86
--- /dev/null
+++ b/Examples/NavierStokes3dStretch.py
@@ -0,0 +1,142 @@
+#!/usr/bin/python
+import parmepy as pp
+import parmepy.f2py
+import numpy as np
+import mpi4py.MPI as MPI
+import math as m
+import time
+
+#from numpy import linalg as LA
+
+pi = m.pi
+
+# Import scales and fftw solvers
+ppfft = parmepy.f2py.fftw2py
+scales = parmepy.f2py.scales2py
+
+rank = MPI.COMM_WORLD.Get_rank()
+print "Mpi process number ", rank
+
+# ----------- A 3d problem -----------
+print " ========= Start test for Navier-Stokes 3D ========="
+
+# Physical Domain description
+Lx = Ly = Lz = 2 * pi
+myDomain3d = pp.Box(dimension=3, length=[Lx, Ly, Lz], origin=[0., 0., 0.])
+resolution3d = np.asarray((65, 65, 65))
+ncells = resolution3d - 1
+hx = Lx / ncells[0]
+hy = Ly / ncells[1]
+hz = Lz / ncells[2]
+
+# Simulation parameters
+finalTime = 0.05
+timeStep = 1e-3
+outputFilePrefix = './res/Stretching_'
+outputModulo = 1
+
+t0 = time.time()
+
+## 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)
+
+# Fields declaration
+
+# 1 - Poisson/diffusion solvers initialisation.
+# See poisson3d.py for a working example.
+# poisson = pp.Poisson(vorticity,velocity)
+#
+localres, localoffset = ppfft.init_fftw_solver(resolution3d,
+ myDomain3d.length)
+
+print "FFT solver local resolution/offset: ", localres, localoffset
+
+##topofft = poisson.getTopology()
+topofft = pp.CartesianTopology(domain=myDomain3d, resolution=resolution3d, dim=3)
+
+
+def computeVel(x, y, z):
+ vx = 2./np.sqrt(3) * np.sin(2.* pi/3.) * np.sin(x) * np.cos(y) * np.cos(z)
+ vy = 2./np.sqrt(3) * np.sin(-2.* pi/3.) * np.cos(x) * np.sin(y) * np.cos(z)
+ vz = 0.
+ return vx , vy , vz
+
+
+def computeVort(x, y, z):
+ wx = - np.cos(x) * np.sin(y) * np.sin(z)
+ wy = - np.sin(x) * np.cos(y) * np.sin(z)
+ wz = 2. * np.sin(x) * np.sin(y) * np.cos(z)
+ return wx , wy , wz
+
+
+velocity = pp.AnalyticalField(domain=myDomain3d, formula=computeVel,
+ name='Velocity', vector=True)
+vorticity = pp.AnalyticalField(domain=myDomain3d, formula=computeVort,
+ name='Vorticity', vector=True)
+
+## 2 - Advection solver initialisation. See testScales for a working example.
+## Based on scales JB solver
+## Warning : fields input for scales should be of size (ncells), not localres.
+#scalesres, scalesoffset, stab_coeff = \
+# scales.init_advection_solver(ncells, myDomain3d.length, order='p_O2')
+
+# 3 - Stretching
+stretch = pp.Stretching(vorticity, velocity)
+
+# 4 - Penalization
+penal = pp.Penalization(velo, vorti, sphere, lambd)
+
+navierStokes = pp.Problem(topofft, [stretch,penal])
+printer = pp.Printer(fields=[vorticity, velocity], frequency=outputModulo,
+ outputPrefix=outputFilePrefix)
+
+navierStokes.setSolver(finalTime, timeStep, solver_type='basic', io=printer)
+
+## Problem => ParticularSover = basic.initialize()
+navierStokes.initSolver()
+
+## end of init ##
+## Mind that scales works on arrays of size resolution - 1
+## --> pointers to subarrays of velocity/vorticity
+#svx = velocity[0][:-1, :-1, :-1]
+#svy = velocity[1][:-1, :-1, :-1]
+#svz = velocity[2][:-1, :-1, :-1]
+#somega_x = vorticity[0][:-1, :-1, :-1]
+#somega_y = vorticity[1][:-1, :-1, :-1]
+#somega_z = vorticity[2][:-1, :-1, :-1]
+
+## solve advection
+#scales.solve_advection(timeStep, svx, svy, svz, somega_x)
+#scales.solve_advection(timeStep, svx, svy, svz, somega_y)
+#scales.solve_advection(timeStep, svx, svy, svz, somega_z)
+
+t1 = time.time()
+
+## solve stretching
+navierStokes.solve()
+
+## solve diffusion
+#vx = velocity[0]
+#vy = velocity[1]
+#vz = velocity[2]
+#omega_x = vorticity[0]
+#omega_y = vorticity[1]
+#omega_z = vorticity[2]
+
+#omega_x, omega_y, omega_z = ppfft.solve_diffusion_3d(vx, vy, vz)
+
+## solve poisson
+#vx, vy, vz = ppfft.solve_poisson_3d(omega_x, omega_y, omega_z)
+
+## end of time loop ##
+
+# Clean memory buffers
+#ppfft.clean_fftw_solver(myDomain3d.dimension)
+
+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)"
diff --git a/HySoP/hysop/operator/differentialOperator_d.py b/HySoP/hysop/operator/differentialOperator_d.py
index 247ad5472..8955cbf12 100755
--- a/HySoP/hysop/operator/differentialOperator_d.py
+++ b/HySoP/hysop/operator/differentialOperator_d.py
@@ -238,6 +238,7 @@ class DifferentialOperator_d(DiscreteOperator):
1.0 * self.field1[0][:,np.roll(ind,-2),:] ) / (12. * self.meshSize[1])
tmp3 = temp1 - temp2
result = np.concatenate((np.array([tmp1]), np.array([tmp2]), np.array([tmp3])))
+ print 'curly' , result
return result
else:
diff --git a/HySoP/hysop/operator/penalization.py b/HySoP/hysop/operator/penalization.py
index a33e19bea..663463a8f 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, lambd):
+ def __init__(self, velocity, vorticity, obstacle, lambd):
"""
Constructor.
Create a Penalization operator from given velocity variavle.
@@ -25,18 +25,20 @@ class Penalization(ContinuousOperator):
self.obstacle = obstacle
## velocity variable to penalize
self.velocity = velocity
+ ## vorticity variable to penalize
+ self.vorticity = vorticity
## lambda penalization function
self.lambd = lambd
- self.addVariable([velocity])
+ self.addVariable([velocity, vorticity])
- def discretize(self, idVelocityD=0, idObstacleD=0):
+ def discretize(self, idVelocityD=0, idVorticityD=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, idObstacleD)
+ self.discreteOperator = Penalization_d(self, idVelocityD, idVorticityD, idObstacleD)
def __str__(self):
"""ToString method"""
diff --git a/HySoP/hysop/operator/penalization_d.py b/HySoP/hysop/operator/penalization_d.py
index b56d7b45e..ff1895940 100644
--- a/HySoP/hysop/operator/penalization_d.py
+++ b/HySoP/hysop/operator/penalization_d.py
@@ -5,6 +5,7 @@ Discrete penalization representation
"""
from ..constants import *
from .discrete import DiscreteOperator
+from .differentialOperator_d import DifferentialOperator_d
#from ..obstacle.obstacle import Obstacle
import pyopencl as cl
import numpy as np
@@ -17,7 +18,7 @@ class Penalization_d(DiscreteOperator):
DiscreteOperator.DiscreteOperator specialization.
"""
- def __init__(self, penal, idVelocityD=0, idObstacleD=0):
+ def __init__(self, penal, idVelocityD=0, idVorticityD=0, idObstacleD=0):
"""
Constructor.
Create a Penalization operator on a given continuous domain.
@@ -30,6 +31,8 @@ class Penalization_d(DiscreteOperator):
self.penal = penal
## Velocity.
self.velocity = penal.velocity.discreteField[idVelocityD]
+ ## Vorticity.
+ self.vorticity = penal.vorticity.discreteField[idVorticityD]
## Obstacle represented by characteristic function Chi
self.obstacle = self.penal.obstacle.discreteObstacle[idObstacleD]
# ## Topology
@@ -62,6 +65,9 @@ class Penalization_d(DiscreteOperator):
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]]
+ self.vorticity.data = DifferentialOperator_d(self.velocity, self.vorticity, choice='curl', resolution=self.vorticity.topology.resolution, meshSize=self.vorticity.topology.mesh.size).apply()
+
+
return self.compute_time
def printComputeTime(self):
diff --git a/HySoP/hysop/test/test_operator/test_Penalization.py b/HySoP/hysop/test/test_operator/test_Penalization.py
index 62b552a73..4ec5a0be3 100644
--- a/HySoP/hysop/test/test_operator/test_Penalization.py
+++ b/HySoP/hysop/test/test_operator/test_Penalization.py
@@ -14,6 +14,13 @@ def vitesse(x, y, z):
vz = 2.
return vx, vy, vz
+def vorticite(x, y, z):
+ wx = 3.
+ wy = 3.
+ wz = 3.
+ return wx, wy, wz
+
+
def run():
# Parameters
nb = 128
@@ -36,9 +43,10 @@ def run():
## Fields
velo = pp.AnalyticalField(domain=box, formula=vitesse, name='Velocity', vector=True)
+ vorti = pp.AnalyticalField(domain=box, formula=vorticite, name='Vorticity', vector=True)
## Operators
- penal = pp.Penalization(velo, sphere, lambd)
+ penal = pp.Penalization(velo, vorti, sphere, lambd)
## Solver creation (discretisation of objects is done in solver initialisation)
topo3D = pp.CartesianTopology(domain=box, resolution=[nb, nb, nb], dim=3)
@@ -46,7 +54,7 @@ def run():
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.setSolver(finalTime, timeStep, solver_type='basic', io=pp.Printer(fields=[velo, vorti], frequency=outputModulo, outputPrefix=outputFilePrefix))
pb.solver.ODESolver= RK4# RK3# RK2# Euler#
pb.initSolver()
t1 = time.time()
--
GitLab