From ca2b3d4bc9f4f30bb5186736620d45b7b2e50991 Mon Sep 17 00:00:00 2001
From: Chloe Mimeau <Chloe.Mimeau@imag.fr>
Date: Wed, 27 Feb 2013 12:24:29 +0000
Subject: [PATCH]
---
Examples/NavierStokes3d_sphere.py | 245 ++++++++++++++++++++++++++++
Examples/input.dat | 38 +++--
HySoP/hysop/operator/diffusion_d.py | 25 +--
3 files changed, 284 insertions(+), 24 deletions(-)
create mode 100755 Examples/NavierStokes3d_sphere.py
diff --git a/Examples/NavierStokes3d_sphere.py b/Examples/NavierStokes3d_sphere.py
new file mode 100755
index 000000000..d0ffc052a
--- /dev/null
+++ b/Examples/NavierStokes3d_sphere.py
@@ -0,0 +1,245 @@
+#!/usr/bin/python
+import parmepy as pp
+import parmepy.f2py
+import numpy as np
+import mpi4py.MPI as MPI
+import math as m
+
+PARMES_REAL = pp.constants.PARMES_REAL
+ORDER = pp.constants.ORDER
+
+#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 ========="
+
+## Opening/reading input data file
+inputData = open("input.dat", "r")
+
+Ox = np.float64(inputData.readline().rstrip('\n\r'))
+Oy = np.float64(inputData.readline().rstrip('\n\r'))
+Oz = np.float64(inputData.readline().rstrip('\n\r'))
+Lx = np.float64(inputData.readline().rstrip('\n\r'))
+Ly = np.float64(inputData.readline().rstrip('\n\r'))
+Lz = np.float64(inputData.readline().rstrip('\n\r'))
+resX = np.uint32(inputData.readline().rstrip('\n\r'))
+resY = np.uint32(inputData.readline().rstrip('\n\r'))
+resZ = np.uint32(inputData.readline().rstrip('\n\r'))
+timeStep = np.float64(inputData.readline().rstrip('\n\r'))
+finalTime = np.float64(inputData.readline().rstrip('\n\r'))
+outputModulo = np.uint32(inputData.readline().rstrip('\n\r'))
+outputFilePrefix = inputData.readline().rstrip('\n\r')
+nbGhosts = np.uint32(inputData.readline().rstrip('\n\r'))
+visco = np.float64(inputData.readline().rstrip('\n\r'))
+uinf = np.float64(inputData.readline().rstrip('\n\r'))
+obstName = inputData.readline().rstrip('\n\r')
+lambdFluid = np.float64(inputData.readline().rstrip('\n\r'))
+lambdPorous = np.float64(inputData.readline().rstrip('\n\r'))
+lambdSolid = np.float64(inputData.readline().rstrip('\n\r'))
+obstRadius = np.float64(inputData.readline().rstrip('\n\r'))
+obstOx = np.float64(inputData.readline().rstrip('\n\r'))
+obstOy = np.float64(inputData.readline().rstrip('\n\r'))
+obstOz = np.float64(inputData.readline().rstrip('\n\r'))
+obstOrientation = inputData.readline().rstrip('\n\r')
+thicknPorous = np.float64(inputData.readline().rstrip('\n\r'))
+thicknZbound = np.float64(inputData.readline().rstrip('\n\r'))
+
+inputData.close()
+
+## Physical Domain description
+myDomain3d = pp.Box(dimension=3, length=[Lx, Ly, Lz], origin=[Ox, Oy, Oz])
+resolution3d = np.asarray((resX, resY, resZ))
+ncells = resolution3d - 1
+hx = Lx / ncells[0]
+hy = Ly / ncells[1]
+hz = Lz / ncells[2]
+
+## Obstacle
+lambd = np.array([lambdFluid, lambdPorous, lambdSolid],
+ dtype=PARMES_REAL, order=ORDER)
+sphere = pp.Obstacle(myDomain3d,
+ name=obstName,
+ zlayer=thicknZbound,
+ radius=obstRadius,
+ center=[obstOx, obstOy, obstOz],
+ orientation=obstOrientation,
+ porousLayer=thicknPorous)
+
+### Post
+#outputFilePrefix = './res/NS_'
+#outputModulo = 50
+
+### Simulation parameters
+#timeStep = 1.#1e-1
+#finalTime = 10000.#20.*timeStep
+
+
+## Topologies definitions
+# 1---- FFT topology + diffusion/poisson solvers initialization -----------
+localres, localoffset = ppfft.init_fftw_solver(resolution3d,
+ myDomain3d.length)
+print "FFT solver local resolution/offset: ", localres, localoffset
+##topofft = poisson.getTopology()
+
+# ---- Cartesian topology -----
+topoCart = pp.CartesianTopology(domain=myDomain3d,
+ resolution=resolution3d, dim=3, ghosts=[nbGhosts,nbGhosts,nbGhosts])
+
+# ---- JB's topology ----------
+nbProcs = MPI.COMM_WORLD.Get_size()
+topodims = [1, 1, nbProcs] # fits with fftw topology
+
+
+## Fields declaration
+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.
+#------------------
+ vx = 0.
+ module = (x-obstOx)**2 + (y-obstOy)**2 + (z-obstOz)**2
+ if (module >= obstRadius**2):
+ vx = uinf * (1-(obstRadius**2/module))
+ vy = 0.
+ vz = 0.
+#-----------------
+# vx = 0.
+# vy = 0.
+# 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)
+#-----------------
+ wx = 0.
+ wy = 0.
+ wz = 0.
+#-----------------
+# xc = 6./2.
+# yc = 6./2.
+# zc = 6./6.
+# R = 1.5
+# sigma = R / 3.
+# Gamma = 0.0075
+# dist = m.sqrt((x-xc)**2 + (y-yc)**2)
+# s2 = (z - zc)**2 + (dist - R)**2
+# wx = 0.
+# wy = 0.
+# wz = 0.
+# if (dist != 0.):
+# cosTheta = (x-xc) / dist
+# sinTheta = (y-yc) / dist
+# wTheta = Gamma / (pi * sigma**2) * m.exp(-(s2 / sigma**2))
+# wx = - wTheta * sinTheta
+# wy = wTheta * cosTheta
+# wz = 0.
+ return wx, wy, wz
+
+def computeVortNorm(x, y, z):
+ norm = 0.
+#----------------
+# xc = 6./2.
+# yc = 6./2.
+# zc = 6./6.
+# R = 1.5
+# sigma = R / 3.
+# Gamma = 0.0075
+# dist = m.sqrt((x-xc)**2 + (y-yc)**2)
+# s2 = (z - zc)**2 + (dist - R)**2
+# norm = 0.
+# if (dist != 0.):
+# cosTheta = (x-xc) / dist
+# sinTheta = (y-yc) / dist
+# wTheta = Gamma / (pi * sigma**2) * m.exp(-(s2 / sigma**2))
+# wx = - wTheta * sinTheta
+# wy = wTheta * cosTheta
+# wz = 0.
+# norm = np.sqrt(wx ** 2 + wy ** 2 + wz ** 2)
+ return norm
+
+#def computeDensity(x, y, z):
+# if (x>=0):
+# rho = 1000.
+# else :
+# rho = 1.
+# return rho
+
+velocity = pp.AnalyticalField(domain=myDomain3d, formula=computeVel,
+ name='Velocity', vector=True)
+vorticity = pp.AnalyticalField(domain=myDomain3d, formula=computeVort,
+ name='Vorticity', vector=True)
+vortNorm = pp.AnalyticalField(domain=myDomain3d, formula=computeVortNorm,
+ name='VortNorm', vector=False)
+#density = pp.AnalyticalField(domain=myDomain3d, formula=computeDensity,
+# name='Density', vector=False)
+
+
+
+## 2 - Advection solver initialization. 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,
+ topodims, order='p_O2')
+
+print 'advection stability coeff', stab_coeff
+advecVort = pp.Advec_scales(stab_coeff, velocity, vorticity, vortNorm)
+#advecDensity = pp.Advec_scales(stab_coeff, velocity, density)
+
+# ----> Change TOPO from Scales to Cartesian
+
+## 3 - Stretching
+stretch = pp.Stretching(velocity, vorticity)
+
+# ----> Change TOPO from Cartesian to FFT
+
+## 4 - Diffusion
+diffusion = pp.Diffusion(velocity, vorticity, viscosity=visco)
+
+## 4 - Poisson
+poisson = pp.Poisson(velocity, vorticity)
+
+# ----> Change TOPO from FFT to Cartesian
+
+## 6 - Penalization
+penal = pp.Penalization(velocity, vorticity, sphere, lambd)
+
+# ----> Change TOPO from Cartesian to Scales
+
+## Define the problem to solve
+navierStokes = pp.Problem(topoCart, [advecVort, stretch, diffusion, poisson, penal])
+#navierStokes = pp.Problem(topoCart, [advecVort, stretch, poisson])
+
+printer = pp.Printer(fields=[vorticity, velocity, vortNorm], frequency=outputModulo,
+ outputPrefix=outputFilePrefix)
+
+navierStokes.setSolver(finalTime, timeStep, solver_type='basic', io=printer)
+
+## Problem => ParticularSover = basic.initialize()
+navierStokes.initSolver()
+
+penal.apply(0., timeStep)
+#poisson.apply(0., timeStep)
+
+## solve stretching
+navierStokes.solve()
+
+## end of time loop ##
+
+# Clean memory buffers
+ppfft.clean_fftw_solver(myDomain3d.dimension)
diff --git a/Examples/input.dat b/Examples/input.dat
index c7ac525a0..defc7d24f 100644
--- a/Examples/input.dat
+++ b/Examples/input.dat
@@ -1,15 +1,27 @@
-0.0
-0.0
-0.0
-6.0
-6.0
-6.0
-257
-257
-257
-1.
-10000.
-50
+-6.0
+-6.0
+-6.0
+12.0
+12.0
+12.0
+65
+65
+65
+0.02
+20000.
+10
./res/NS_
2
-0.000001
+0.005
+1.
+sphere
+0.
+10.
+10E8
+1.0
+0.
+0.
+0.
+West
+0.
+0.05
diff --git a/HySoP/hysop/operator/diffusion_d.py b/HySoP/hysop/operator/diffusion_d.py
index 7abd019c7..f3af5ed3c 100644
--- a/HySoP/hysop/operator/diffusion_d.py
+++ b/HySoP/hysop/operator/diffusion_d.py
@@ -8,7 +8,6 @@ from discrete import DiscreteOperator
from parmepy.constants import ORDER, PARMES_REAL
import numpy as np
import time
-import numpy.linalg as la
class Diffusion_d(DiscreteOperator):
@@ -39,7 +38,6 @@ class Diffusion_d(DiscreteOperator):
Apply operator.
"""
self.compute_time = time.time()
- print "Norm de vz pre ... ", la.norm(self.velocity[2])
ind0a = self.ghosts[0]
ind0b = self.resolution[0] - self.ghosts[0]
@@ -60,22 +58,27 @@ class Diffusion_d(DiscreteOperator):
velocityNoG[i][...] = self.velocity[i][ind0a:ind0b,
ind1a:ind1b, ind2a:ind2b]
- # Vorticity diffusion
+ # Curl + Vorticity diffusion
+# vorticityNoG[0][...], vorticityNoG[1][...], vorticityNoG[2][...] = \
+# fftw2py.solve_curl_diffusion_3d(self.viscosity * dt,
+# velocityNoG[0][...],
+# velocityNoG[1][...],
+# velocityNoG[2][...],
+# vorticityNoG[0][...],
+# vorticityNoG[1][...],
+# vorticityNoG[2][...])
+
+ # Pure Vorticity diffusion
vorticityNoG[0][...], vorticityNoG[1][...], vorticityNoG[2][...] = \
fftw2py.solve_diffusion_3d(self.viscosity * dt,
- velocityNoG[0][...],
- velocityNoG[1][...],
- velocityNoG[2][...],
- vorticityNoG[0][...],
- vorticityNoG[1][...],
- vorticityNoG[2][...])
+ vorticityNoG[0][...],
+ vorticityNoG[1][...],
+ vorticityNoG[2][...])
for i in xrange(self.dim):
self.vorticity[i][ind0a:ind0b, ind1a:ind1b, ind2a:ind2b] = \
vorticityNoG[i][...]
- print "Norm de vz post ... ", la.norm(self.velocity[2])
-
self.compute_time = time.time() - self.compute_time
self.total_time += self.compute_time
return self.compute_time
--
GitLab