From 06761b0d9e17a8e4d191b8487e7f164b9a0d8872 Mon Sep 17 00:00:00 2001
From: Chloe Mimeau <Chloe.Mimeau@imag.fr>
Date: Mon, 7 Apr 2014 14:54:28 +0000
Subject: [PATCH] update example for flow past sphere
---
Examples/NSDebug.py | 274 +++++++++++++++++++++++++-------------------
1 file changed, 158 insertions(+), 116 deletions(-)
diff --git a/Examples/NSDebug.py b/Examples/NSDebug.py
index be6101baa..43328a9bb 100755
--- a/Examples/NSDebug.py
+++ b/Examples/NSDebug.py
@@ -12,24 +12,26 @@ from parmepy.f2py import fftw2py
import numpy as np
from parmepy.fields.continuous import Field
from parmepy.mpi.topology import Cartesian
+from parmepy.fields.variable_parameter import VariableParameter
from parmepy.domain.obstacle.sphere import Sphere
from parmepy.operator.advection import Advection
from parmepy.operator.stretching import Stretching
from parmepy.operator.poisson import Poisson
from parmepy.operator.velocity_correction import VelocityCorrection
+#from parmepy.operator.vorticity_correction import VorticityCorrection
from parmepy.operator.diffusion import Diffusion
from parmepy.operator.penalization import Penalization
from parmepy.operator.differential import Curl
from parmepy.operator.redistribute import Redistribute
from parmepy.operator.monitors import Printer, Energy_enstrophy, DragAndLift,\
- Reprojection_criterion
+ Reprojection_criterion, ProfilesOutput
from parmepy.problem.simulation import Simulation
from parmepy.constants import VTK, HDF5
-from parmepy.domain.obstacle.planes import PlaneBoundaries
+from parmepy.domain.obstacle.planes import PlaneBoundaries, SubPlane
+from parmepy.domain.obstacle.controlBox import ControlBox
from dataNS_bb import dim, nb, NBGHOSTS, ADVECTION_METHOD, VISCOSITY, \
OUTPUT_FREQ, FILENAME, PROJ, LCFL, CFL, CURL_METHOD, \
TIMESTEP_METHOD, OBST_Ox, OBST_Oy, OBST_Oz, RADIUS
-from parmepy.domain.obstacle.controlBox import ControlBox
import parmepy.tools.numpywrappers as npw
import parmepy.tools.io_utils as io
## ----------- A 3d problem -----------
@@ -41,17 +43,35 @@ cos = np.cos
sin = np.sin
## Domain
-boxlength = npw.realarray([2.0 * pi, 2 * pi, 2 * pi])
-boxorigin = npw.realarray([0., 0., 0.])
+#boxlength = npw.realarray([8.0 * pi, 2.0 * pi, 2.0 * pi])
+#boxorigin = npw.realarray([-2.0, -pi, -pi])
+boxlength = npw.realarray([10.24, 5.12, 5.12])
+boxorigin = npw.realarray([-2.0, -2.56, -2.56])
box = pp.Box(dim, length=boxlength, origin=boxorigin)
## Global resolutionsimu
-nbElem = [129]*3
+#nbElem = [1025, 513, 513]
+nbElem = [129, 65, 65]
-# Upstream flow velocity
+
+# Initial upstream flow velocity
uinf = 1.0
-# Function to compute potential flow past a sphere
+def computeFlowrate(simu):
+ ##Time-dependant flow rate
+ t = simu.tk
+ Tstart = 3.0
+ flowrate = np.zeros(3)
+ flowrate[0] = uinf * box.length[1] * box.length[2]
+ if t >= Tstart and t <= Tstart + 1.0:
+ flowrate[1] = sin(pi * (t - Tstart)) * \
+ box.length[1] * box.length[2]
+ ##constant flow rate
+# flowrate = np.zeros(3)
+# flowrate[0] = uinf * box.length[1] * box.length[2]
+ return flowrate
+
+# Function to compute initial velocity
def computeVel(res, x, y, z, t):
res[0][...] = uinf
res[1][...] = 0.
@@ -72,16 +92,23 @@ velo = Field(domain=box, formula=computeVel,
vorti = Field(domain=box, formula=computeVort,
name='Vorticity', isVector=True)
-# Topo for operators that required ghosts points (stretching)
+## Usual Cartesian topology definition
+# At the moment we use two (or three?) topologies :
+# - "topo" for Stretching and all operators based on finite differences.
+# --> ghost layer = 2
+# - topo from Advection operator for all the other operators.
+# --> no ghost layer
+# - topo from fftw for Poisson and Diffusion.
+# Todo : check compat between scales and fft operators topologies.
ghosts = np.ones((box.dimension)) * NBGHOSTS
-topostr = Cartesian(box, box.dimension, nbElem, ghosts=ghosts)
+topo = Cartesian(box, box.dimension, nbElem,
+ ghosts=ghosts)
## === Obstacles (sphere + up and down plates) ===
-sphere_pos = (boxlength - boxorigin) * 0.5
+#sphere_pos = (boxlength - boxorigin) * 0.5
+sphere_pos = [0., 0., 0.]
sphere = Sphere(box, position=sphere_pos, radius=RADIUS)
-bc = PlaneBoundaries(box, 2, thickness=0.1)
-
## === Computational Operators ===
op = {}
op['advection'] = Advection(velo, vorti,
@@ -90,7 +117,7 @@ op['advection'] = Advection(velo, vorti,
op['stretching'] = Stretching(velo, vorti,
resolutions={velo: nbElem, vorti: nbElem},
- topo=topostr)
+ topo=topo)
op['diffusion'] = Diffusion(vorti, resolution=nbElem, viscosity=VISCOSITY)
@@ -108,15 +135,20 @@ opfft = op['poisson']
topofft = opfft.discreteFields[opfft.vorticity].topology
topocurl = op['curl'].discreteFields[op['curl'].invar].topology
topoadvec = op['advection'].discreteFields[op['advection'].velocity].topology
+# space step
+ref_step = topofft.mesh.space_step
+step_dir = ref_step[0]
-op['penalization'] = Penalization(velo, [sphere, bc], coeff=[1e8],
+## Penalization velocity
+op['penalization'] = Penalization(velo, [sphere], coeff=[1e8],
topo=topofft, resolutions={velo: nbElem})
-ref_rate = uinf * box.length[1] * box.length[2]
+## Time-dependant required-flowrate (Variable Parameter)
+req_flowrate = VariableParameter(formula=computeFlowrate)
op['correction'] = VelocityCorrection(velo, vorti,
resolutions={velo: nbElem,
vorti: nbElem},
- req_flowrate=ref_rate, topo=topofft)
+ req_flowrate=req_flowrate, topo=topofft)
op['penalization'].discretize()
op['correction'].discretize()
@@ -126,37 +158,73 @@ distr = {}
distr['fft2curl'] = Redistribute([velo], op['penalization'], op['curl'])
distr['curl2fft'] = Redistribute([vorti], op['curl'], op['poisson'])
+
+# 1 - Curl to stretching (velocity only)
+#distr['curl2str'] = Redistribute([velo], op['curl'], op['stretching'])
+# 2 - Advection to stretching
+# velocity only
+#distr['adv2str_v'] = Redistribute([velo], op['advection'], op['stretching'])
+# vorticity only
+#distr['adv2str_w'] = Redistribute([vorti], op['advection'], op['stretching'])
+# Both
distr['adv2str'] = Redistribute([velo, vorti],
op['advection'], op['stretching'])
+# 3 - Stretching to Diffusion
distr['str2diff'] = Redistribute([vorti], op['stretching'], op['diffusion'])
distr['curl2adv'] = Redistribute([velo, vorti], op['curl'], op['advection'])
distr['curl2str'] = Redistribute([velo, vorti], op['curl'], op['stretching'])
+#distr['adv2corr'] = Redistribute([velo, vorti],
+# op['advection'], op['correc_vorti'])
+
for ope in distr.values():
ope.discretize()
# === Simulation setup ===
-simu = Simulation(tinit=0.0, tend=10., nbIter=2000, iterMax=1000000)
+simu = Simulation(tinit=0.0, tend=75., timeStep=0.0125, iterMax=1000000)
## === Monitoring operators ===
monitors = {}
# Save velo and vorti values on topofft
monitors['printerFFT'] = Printer(variables=[velo, vorti], topo=topofft,
- formattype=HDF5, prefix='fft_io')
+ formattype=HDF5)
monitors['printerCurl'] = Printer(variables=[velo, vorti], topo=topocurl,
- formattype=HDF5, prefix='curl_io')
-monitors['printerAdvec'] = Printer(variables=[velo, vorti], topo=topoadvec,
- formattype=HDF5, prefix='advec_io',
+ formattype=HDF5, prefix='curl_io_vorti')
+
+monitors['printerAdvec'] = Printer(variables=[velo, vorti], topo=topofft,
+ frequency=1, formattype=HDF5, prefix='advec_io',
xmfalways=True)
-monitors['printerStr'] = Printer(variables=[velo, vorti], topo=topostr,
- formattype=HDF5, prefix='str_io',
- xmfalways=True)
+
+thickSliceXY = ControlBox(box, origin=[-2.0, -2.56, -2.0 * step_dir],
+ lengths=[10.24, 5.12, 4.0 * step_dir])
+#thickSliceXY = ControlBox(box, origin=[-2.0, -pi, -2.0 * step_dir],
+# lengths=[8.0*pi, 2.0*pi, 4.0 * step_dir])
+monitors['printerSliceXY'] = Printer(variables=[velo, vorti], topo=topofft,
+ frequency=100, formattype=HDF5, prefix='sliceXY',
+ subset=thickSliceXY, xmfalways=True)
+
+thickSliceXZ = ControlBox(box, origin=[-2.0, -2.0 * step_dir, -2.56],
+ lengths=[10.24, 4.0 * step_dir, 5.12])
+monitors['printerSliceXZ'] = Printer(variables=[velo, vorti], topo=topofft,
+ frequency=100, formattype=HDF5, prefix='sliceXZ',
+ subset=thickSliceXZ, xmfalways=True)
+
+subBox = ControlBox(box, origin=[-0.7, -1.25, -1.25], lengths=[7.0, 2.5, 2.5])
+monitors['printerSubBox'] = Printer(variables=[velo, vorti], topo=topofft,
+ frequency=100, formattype=HDF5, prefix='subBox',
+ subset=subBox, xmfalways=True)
# Compute and save enstrophy/Energy, on topofft
io_default = {}
-monitors['energy'] = Energy_enstrophy(velo, vorti, topo=topoadvec,
+monitors['energy'] = Energy_enstrophy(velo, vorti, topo=topofft,
io_params={},
viscosity=VISCOSITY, isNormalized=False)
+
+# Compute velocity and vorticity magnitude on topofft and print it in data file
+io_default = {}
+monitors['profiles'] = ProfilesOutput(velo, vorti, obstacles=[sphere],
+ topo=topofft, io_params={"frequency":1000})
+
# Setup for reprojection
reprojection_constant = 0.04
reprojection_rate = 1
@@ -172,19 +240,17 @@ vdfft = velo.discreteFields[topofft].data
# Compute and save drag and lift on topofft
# 3D Control box
-ref_step = topofft.mesh.space_step
cbpos = boxorigin + 10 * ref_step
-cblength = boxlength - 20 * ref_step
+cblength = boxlength - 25 * ref_step
cb = ControlBox(box, cbpos, cblength)
coeffForce = 1. / (0.5 * uinf ** 2 * pi * RADIUS ** 2)
-print coeffForce
+#print coeffForce
# Monitor for forces
-## monitors['forces'] = DragAndLift(velo, vorti, VISCOSITY, coeffForce,
-## topo, cb,
-## obstacles=[sphere], io_params={})
monitors['forces'] = DragAndLift(velo, vorti, VISCOSITY, coeffForce,
- topostr, cb, io_params={})
+ topo, cb, obstacles=[sphere], io_params={})
+#monitors['forces'] = DragAndLift(velo, vorti, VISCOSITY, coeffForce,
+# topo, cb, io_params={})
# === Setup for all declared operators ===
for ope in op.values():
@@ -205,21 +271,18 @@ allops = dict(op.items() + distr.items() + monitors.items())
# - compute vorticity with curl
def initFields_mode1():
# The velocity is initialized on the fftw topology
- # penalized and distributed to curl topo
+ # penalized and distributed on curl topo
velo.initialize(topo=topofft)
op['penalization'].apply(simu)
distr['fft2curl'].apply(simu)
# Vorticity is initialized after penalization of velocity
op['curl'].apply(simu)
- # Distribute vorticity and velocity on topostr
- distr['curl2str'].apply(simu)
- # Warning : distribute operators do not update ghost points. This is done
- # by computational operators. So a synchronize call may be required
- # to monitor proper values.
- #op['stretching'].discreteOperator._synchronize(op['stretching'].discreteOperator.velocity.data + op['stretching'].discreteOperator.vorticity.data)
- op['stretching'].updateGhosts()
+ op['advection'].apply(simu)
+ # Distribute vorticity on topofft
+ distr['adv2str'].apply(simu)
+# distr['curl2adv'].apply(simu)
# From this point both velocity and vorticity are initialized
- # on topostr and topocurl. velocity is also init. on topofft.
+ # on topocurl and topoadvec. velocity is also init. on topofft.
## Call initialization
initFields_mode1()
@@ -227,73 +290,70 @@ initFields_mode1()
##### Check initialize process results #####
norm_vel_fft = velo.norm(topofft)
norm_vel_curl = velo.norm(topocurl)
-norm_vel_str = velo.norm(topostr)
-
+norm_vort_fft = vorti.norm(topofft)
norm_vort_curl = vorti.norm(topocurl)
-norm_vort_str = vorti.norm(topostr)
-
-print norm_vort_curl, norm_vort_str
-print norm_vel_curl, norm_vel_fft
-assert np.allclose(norm_vort_curl, norm_vort_str)
-assert np.allclose(norm_vel_curl, norm_vel_fft)
-assert np.allclose(norm_vel_curl, norm_vel_str)
-
-wref = np.asarray([0., 4354.98713193, 612.8061093], dtype=np.float64)
-vref = np.asarray([1418.04337028, 0., 0.], dtype=np.float64)
-assert np.allclose(norm_vel_curl, vref)
-assert np.allclose(norm_vort_str, wref)
+
+#print norm_vort_curl, norm_vort_fft
+#print norm_vel_curl, norm_vel_fft
+#assert np.allclose(norm_vort_curl, norm_vort_fft)
+#assert np.allclose(norm_vel_curl, norm_vel_fft)
+
# Print initial state
-for mon in monitors.values():
- mon.apply(simu)
+#for mon in monitors.values():
+# mon.apply(simu)
## Note Franck : tests init ok, same values on both topologies, for v and w.
# === After init ===
#
-# v,w init on topostr
+# v init on topocurl (which may be equal to topofft)
#
# === Definition of the 'one-step' sequence of operators ===
#
-# A - w = stretching(v,w), topostr
-# B
-# 1 - w = diffusion(w), topofft
-# 2 - w = reprojection(w), topofft
-# 3 - v = poisson(w), topofft
-# 4 - v = correction(v, w), topofft
-# 5 - v = penalization(v), topofft
-# C - w = curl(v), topocurl
-# D - w = advection(v,w), topoadvec
+# 1 - w = curl(v), topocurl
+# 2 - w = advection(v,w), topo-advec
+# 3 - w = stretching(v,w), topostr
+# 4 - w = diffusion(w), topofft
+# 5 - v = poisson(w), topofft
+# 6 - v = correction(v, w), topofft
+# 7 - v = penalization(v), topofft
#
-
-# Required bridges and in/out var:
-# A, IN : v, w on topostr
-# OUT : w on topostr
-# A - B bridge topostr-topofft for w
-# B, IN : w
-# OUT : w, v on topofft
-# B - C : bridge topofft - topocurl for v
-# C, IN : v on topocurl
-# OUT : w on topocurl
-# C - D : bridge topocurl - topoadvec for v and w
-# D, IN : v,w on topoadvec
-# OUT : w on topoadvec
-# D - A : bridge topoadvec - topostr for v and w
-
-
-# Monitors :
-# - energy(v,w), works for any topo
-# - reprojection(w), needs ghost points
-# - printer(v,w), works for any topo
-# - forces(v,w), needs ghost points
-
-fullseq = ['stretching',
- 'str2diff', 'diffusion', 'reprojection', 'poisson', 'correction',
+# Required bridges :
+# 1 --> 2 : (v,w) on topo-advec
+# 2 --> 3 : (v,w) on topostr
+# 3 --> 4 : w on topofft
+# 7 --> 1 : v on topocurl
+
+#fullseq = ['advection', # in: v,w out: w, on topoadvec
+# 'adv2str', 'stretching', # in: v,w out: w on topostr
+# # in: w, out: v, w on topofft
+# 'str2diff', 'diffusion', 'reprojection', 'poisson', 'correction',
+# 'penalization',
+# 'fft2curl', 'curl', 'forces',
+# 'curl2adv', 'energy', 'profiles', 'printerAdvec',
+# ]
+
+fullseq = ['stretching',
+ 'str2diff', 'diffusion', 'poisson', 'correction',
'penalization',
- 'fft2curl', 'curl',
- 'curl2adv',
+ 'fft2curl', 'curl',
+ 'curl2fft', 'poisson', 'correction',
'advection',
- 'energy', 'printerAdvec',
+# 'printerSliceXY', 'printerSliceXZ', 'printerSubBox',
+ 'energy',
'adv2str', 'forces']
+#fullseq = ['advection',
+# 'diffusion',
+# 'poisson',
+# 'correction',
+# 'penalization',
+# 'fft2curl',
+# 'curl',
+# 'curl2fft',
+# 'poisson',
+# 'correction',
+# 'printerSlice']
+
def run(sequence):
for name in sequence:
@@ -302,12 +362,12 @@ def run(sequence):
allops[name].apply(simu)
seq = fullseq
-seq = []
+
simu.initialize()
while not simu.isOver:
- #if topocurl.rank == 0:
- # simu.printState()
+ if topocurl.rank == 0:
+ simu.printState()
run(seq)
simu.advance()
@@ -343,21 +403,3 @@ for ope in distr.values():
ope.finalize()
for monit in monitors.values():
monit.finalize()
-
-from abc import ABCMeta, abstractmethod
-
-
-class A(object):
-
- def toto(self):
- print "toto A "
-
- def titi(self):
- self.toto()
-
-
-class B(A):
-
- def toto(self):
- print "toto B"
-
--
GitLab