diff --git a/Docs/manual.pdf b/Docs/manual.pdf
new file mode 100644
index 0000000000000000000000000000000000000000..3f58dc12655ad1aa90334c96014e7db543a26632
Binary files /dev/null and b/Docs/manual.pdf differ
diff --git a/Docs/manual.tex b/Docs/manual.tex
new file mode 100644
index 0000000000000000000000000000000000000000..4e12e287fdd7506169d853a48256db7966869320
--- /dev/null
+++ b/Docs/manual.tex
@@ -0,0 +1,77 @@
+\documentclass[a4paper,10pt]{article}
+\usepackage{ucs,amsmath,a4wide,natbib, amssymb}
+\usepackage[utf8x]{inputenc}
+\usepackage{fancyhdr}
+\usepackage{lastpage}
+\fancyhead[C]{}
+\fancyhead[R]{\texttt{\thepage/\pageref{LastPage}}}
+\fancyfoot[L]{}
+\fancyfoot[C]{}
+\fancyfoot[R]{\texttt{file main.tex -- \isodayandtime}}
+\usepackage{multirow}
+\setlength{\parskip}{12pt}
+\usepackage{graphicx}
+\graphicspath{{./images/}}
+\usepackage{pdfsync}
+ \usepackage{cancel}
+\usepackage{xcolor}
+\usepackage[left=3cm,right=3cm,top=2.5cm,bottom=2.5cm]{geometry}
+\overfullrule 5pt
+\usepackage[left=3cm,right=3cm,top=2.5cm,bottom=2.5cm]{geometry}
+\overfullrule 5pt %Pour d\'esigner les listes couples d'un *full carre noir la ou l on sort
+
+\usepackage[pdftex]{hyperref}
+\hypersetup{
+ pdftitle={ParMeS User Guide},
+ pdfauthor={ParMeS Devel},
+ pdfcreator={ParMeS},
+ pdfkeywords={particular methods, GPU, MPI},
+% pdfsubject={}
+% pagebackref=false,
+%colorlinks=false,citecolor=red,
+%linkcolor=blue,
+%plainpages=false,
+%breaklinks=true,
+% pdfpagelabels=true
+}
+
+\usepackage{mathtools}
+
+\DeclarePairedDelimiter\abs{\lvert}{\rvert}%
+\DeclarePairedDelimiter\norm{\lVert}{\rVert}%
+% Swap the definition of \abs* and \norm*, so that \abs
+% and \norm resizes the size of the brackets, and the
+% starred version does not.
+\makeatletter
+\let\oldabs\abs
+\def\abs{\@ifstar{\oldabs}{\oldabs*}}
+%
+\let\oldnorm\norm
+\def\norm{\@ifstar{\oldnorm}{\oldnorm*}}
+\makeatother
+
+\title{ParMeS User Guide}
+\author{ParMeS Development team}
+
+\begin{document}
+\maketitle
+
+
+\part{Introduction}
+\part{Operators}
+\section{About operators}
+\section{Computational operators}
+\input{velocity_correction}
+\section{Monitors}
+\input{monitors}
+\section{Redistribute operators}
+
+\end{document}
+
+%%% Local Variables:
+%%% mode: latex
+%%% TeX-master: t
+%%% End:
+
+
+
diff --git a/Docs/monitors.tex b/Docs/monitors.tex
new file mode 100644
index 0000000000000000000000000000000000000000..c3edc57dea6c0bed57955628342cedb04a4c51d9
--- /dev/null
+++ b/Docs/monitors.tex
@@ -0,0 +1,41 @@
+
+\subsection{Introduction}
+\subsection{Printer}
+\subsection{DragAndLift}
+\subsection{Reprojection$\_$criterion}
+Purpose : compute a projection criterion, depending on the current vorticity value on the grid.
+
+IN : vorticity field, reproj$\_$cst\\
+OUT : reprojRate, projCriterion
+with\\
+ reprojRate : how often (in terms of iterations) this criteria must be taken into account \\
+ and checkCriterion : true if reprojRate must be checked and updated (if necessary).
+
+\begin{itemize}
+\item Case 1: checkCriterion == True.\\
+If checkCriterion is true, the operator computes a criterion projection every reprojRate iteration.
+If this criterion is greater than reproj$\_$cst (user defined constant), then reprojRate is set to 1.
+reprojRate is reset to defaultRate (user-defined) at the beginning of every apply.
+\item Case 2 : checkCriterion == False.\\
+The operator will do nothing but give the (user-defined) value for reprojRate.
+\end{itemize}
+
+If io$\_$params is set, then a file is appended with some diagnostics values (see below) each reprojRate iterations.\\
+
+$\Omega$ being the whole domain, $\omega$ the vorticity field, let us denote :
+\begin{eqnarray}
+d1 &=& \max_{\Omega}\abs{\sum_j \frac{\partial \omega_j}{\partial x_j}}\\
+d2 &=& \max_{\Omega}\abs{\frac{\partial \omega_i}{\partial x_j}} \ \ \forall i,j \in [0,2]
+\end{eqnarray}
+
+And so, diagnostics = simulation.time, projCriterion, d1, d2, counter \\
+
+counter being the number of projection done since the beginning of the simulation.
+
+
+\subsection{Energy$\_$enstrophy}
+%%% Local Variables:
+%%% mode: latex
+%%% TeX-master: "manual"
+%%% End:
+
diff --git a/Docs/velocity_correction.pdf b/Docs/velocity_correction.pdf
deleted file mode 100644
index 832eb6c5416b2f09a32efef01b8ab5b90d6800cb..0000000000000000000000000000000000000000
Binary files a/Docs/velocity_correction.pdf and /dev/null differ
diff --git a/Docs/velocity_correction.tex b/Docs/velocity_correction.tex
index 229f69ed7d3613013722b82a55affa2d57918f0f..a19616e22cb6100a0db8835836d5f5c8b5d733a7 100644
--- a/Docs/velocity_correction.tex
+++ b/Docs/velocity_correction.tex
@@ -1,43 +1,4 @@
-\documentclass[a4paper,10pt]{article}
-\usepackage[utf8x]{inputenc}
-\usepackage{amsmath}
-% \usepackage{showkeys} % affiche les labels et references
-\usepackage[T1]{fontenc}
-\usepackage{multirow}
-\usepackage{amsthm,amssymb}
-\setlength{\parskip}{12pt} %taille du saut entre deux paragraphes
-\usepackage{graphicx}
-\graphicspath{{./images/}}
-\usepackage{pdfsync} % permet de cliquer sur le pdf pour avoir l endroit ou il y a le code latex
- \usepackage{cancel}
-\usepackage{xcolor}
-\usepackage[left=3cm,right=3cm,top=2.5cm,bottom=2.5cm]{geometry}
-\overfullrule 5pt %Pour d\'esigner les listes couples d'un *full carre noir la ou l on sort
-
-\usepackage[pdftex]{hyperref}
-\hypersetup{
- pdftitle={Flow past 3D sphere},
- pdfauthor={Mimeau Chloe},
- pdfcreator={PARMES},
- pdfkeywords={sphere, potential flow, velocity correction, flow rate},
- pdfsubject={Flow past 3D sphere}
-% pagebackref=false,
-%colorlinks=false,citecolor=red,
-%linkcolor=blue,
-%plainpages=false,
-%breaklinks=true,
-% pdfpagelabels=true
-}
-
-
-\theoremstyle{plain}
-\newtheorem{Theoreme}{Th\'eor\`eme}[section]
-%opening
-\title{Correction de vitesse pour la modélisation d'un écoulement incompressible dans une boîte}
-\author{}
-\date{}
-
-\begin{document}
+\subsection{Correction de vitesse pour la modélisation d'un écoulement incompressible dans une boîte}
\maketitle
Lors de la résolution de l'équation de Poisson $\Delta \mathbf{u} = - \nabla \boldsymbol{\omega} $, le mode 0 des transformées de Fourier rapides de la vitesse est imposé à 0, ce qui entraîne une modification des valeurs moyennes des composantes de la vitesse et donc des composantes de la vitesse elles-mêmes lorsque la transformée de Fourier inverse est appliquée. Une correction de la vitesse est donc nécessaire après la résolution de l'équation de Poisson, à chaque itération.
@@ -219,5 +180,8 @@ u_z &= -u_{\infty} \dfrac{3R^3 xz}{2 r^5}
% \begin{equation}
% \underbrace{\iint\limits_S u_x \ dydz}_{\text{débit souhaité}} = \iint\limits_S \ \dfrac{-u_{\infty} x}{r} \left ( 1 - \dfrac{R^3}{r^3} \right ) \ dydz = -u_{\infty} x_0 \iint\limits_S \ \dfrac{1}{r} \left ( 1 - \dfrac{R^3}{r^3} \right ) \ dydz
% \end{equation}
+%%% Local Variables:
+%%% mode: latex
+%%% TeX-master: "manual"
+%%% End:
-\end{document}
diff --git a/Examples/NSDebug.py b/Examples/NSDebug.py
index 1217d0ab50b5615923667e5b0e31cbd26271197c..bbc6095af0bc65910349bd5e384927ecee78689d 100755
--- a/Examples/NSDebug.py
+++ b/Examples/NSDebug.py
@@ -11,7 +11,6 @@ import parmepy as pp
from parmepy.f2py import fftw2py
import numpy as np
from parmepy.fields.continuous import Field
-from parmepy.variables.variables import Variables
from parmepy.mpi.topology import Cartesian
from parmepy.domain.obstacle.sphere import Sphere
from parmepy.operator.advection import Advection
@@ -22,19 +21,17 @@ 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.printer import Printer
-from parmepy.operator.monitors.energy_enstrophy import Energy_enstrophy
-from parmepy.operator.monitors.compute_forces import Forces
+from parmepy.operator.monitors import Printer, Energy_enstrophy, DragAndLift,\
+ Reprojection_criterion
from parmepy.problem.simulation import Simulation
from parmepy.constants import VTK, HDF5
from parmepy.domain.obstacle.planes import PlaneBoundaries
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.operator.monitors.compute_forces import DragAndLift
from parmepy.domain.obstacle.controlBox import ControlBox
-
-
+import parmepy.tools.numpywrappers as npw
+import parmepy.tools.io_utils as io
## ----------- A 3d problem -----------
print " ========= Start Navier-Stokes 3D (Flow past bluff bodies) ========="
@@ -44,8 +41,8 @@ cos = np.cos
sin = np.sin
## Domain
-boxlength = [2.0 * pi, pi, pi]
-boxorigin = [0., 0., 0.]
+boxlength = npw.realarray([2.0 * pi, pi, pi])
+boxorigin = npw.realarray([0., 0., 0.])
box = pp.Box(dim, length=boxlength, origin=boxorigin)
## Global resolution
@@ -105,76 +102,79 @@ op['stretching'] = Stretching(velo, vorti,
op['diffusion'] = Diffusion(vorti, resolution=nbElem, viscosity=VISCOSITY)
-op['poisson'] = Poisson(velo, vorti, resolutions={velo: nbElem, vorti: nbElem},
- projection=PROJ)
+op['poisson'] = Poisson(velo, vorti, resolutions={velo: nbElem, vorti: nbElem})
op['curl'] = Curl(velo, vorti, resolutions={velo: nbElem, vorti: nbElem},
method=CURL_METHOD)
## Discretisation of computational operators
-for ope in op:
+for ope in op.values():
ope.discretize()
# Get fft topology and use it for penalization and correction operators
opfft = op['poisson']
topofft = opfft.discreteFields[opfft.vorticity].topology
+topocurl = op['curl'].discreteFields[op['curl'].invar].topology
-op['penal'] = Penalization(velo, [sphere, bc], coeff=[1e8], topo=topofft,
- resolutions={velo: nbElem})
+op['penalization'] = Penalization(velo, [sphere, bc], coeff=[1e8],
+ topo=topofft, resolutions={velo: nbElem})
-op['correc'] = VelocityCorrection(velo, vorti,
- resolutions={velo: nbElem, vorti: nbElem},
- uinf=uinf, topo=topofft)
+ref_rate = uinf * box.length[1] * box.length[2]
+op['correction'] = VelocityCorrection(velo, vorti,
+ resolutions={velo: nbElem,
+ vorti: nbElem},
+ req_flowrate=ref_rate, topo=topofft)
-op['penal'].discretize()
-op['correc'].discretize()
+op['penalization'].discretize()
+op['correction'].discretize()
## === Bridges between the different topologies ===
distr = {}
+
+distr['fft2curl'] = Redistribute([velo], op['penalization'], op['curl'])
+distr['curl2fft'] = Redistribute([vorti], op['curl'], op['poisson'])
+
# 1 - Curl_fft to stretching (velocity only)
distr['Curl2Str'] = Redistribute([velo], op['curl'], op['stretching'])
# 2 - Advection to stretching (vorticity only)
distr['Advec2Str_v'] = Redistribute([velo], op['advection'], op['stretching'])
distr['Advec2Str_w'] = Redistribute([vorti], op['advection'], op['stretching'])
# 3 - Stretching to Diffusion
-distr['stretch2diff'] = Redistribute([vorti], op['stretching'], op['diffusion'])
-# + - Brigdes required if Curl op. is solved using a FD method
-distrPoissCurl = Redistribute([vorti, velo], poisson, curl)
-distrCurlAdv = Redistribute([vorti, velo], curl, advec)
+distr['Str2diff'] = Redistribute([vorti], op['stretching'], op['diffusion'])
+distr['Curl2Advec'] = Redistribute([vorti], op['curl'], op['advection'])
-distrAdvStr_velo.discretize()
-distrAdvStr_vorti.discretize()
-distrStrDiff.discretize()
-distrCurlStr_velo.discretize()
-distrCurlAdv.discretize()
-
-## === Fields initialization ===
-velo.initialize(topo=topofft)
+for ope in distr.values():
+ ope.discretize()
-ind = sphere.discretize(topofft)
-vdfft = velo.discreteFields[topofft].data
+# === Simulation setup ===
+simu = Simulation(tinit=0.0, tend=0.2, timeStep=0.005, iterMax=1000000)
## === Monitoring operators ===
+monitors = {}
+# Save velo and vorti values on topofft
+monitors['printerFFT'] = Printer(variables=[velo, vorti], topo=topofft,
+ formattype=HDF5)
+
+monitors['printerCurl'] = Printer(variables=[velo, vorti], topo=topocurl,
+ formattype=HDF5, prefix='curl_io_vorti')
+# Compute and save enstrophy/Energy, on topofft
+io_default = {}
+monitors['energy'] = Energy_enstrophy(velo, vorti, topo=topofft, io_params={},
+ viscosity=VISCOSITY, isNormalized=False)
+# Setup for reprojection
+reprojection_constant = 0.04
+reprojection_rate = 1
+monitors["reprojection"] = Reprojection_criterion(vorti, reprojection_constant,
+ reprojection_rate, topofft,
+ io_params={})
+
+# Add reprojection for poisson operator
+op["poisson"].activateProjection(monitors["reprojection"])
-outputdir = 'res_' + str(topofft.size) + 'procs'
-pref = outputdir + '/vwfft'
-
-printerFFT = Printer(variables=[velo, vorti],
- topo=topofft,
- frequency=1,
- prefix=pref,
- formattype=HDF5)
-
-printerFFT.setUp()
-
-prefE = outputdir + '/ener'
-energy = Energy_enstrophy(velo, vorti,
- topo=topofft,
- viscosity=VISCOSITY,
- isNormalized=False,
- frequency=1,
- prefix=prefE)
+ind = sphere.discretize(topofft)
+vdfft = velo.discreteFields[topofft].data
+# Compute and save drag and lift on topofft
# 3D Control box
ref_step = topofft.mesh.space_step
cbpos = boxorigin + boxlength + 10 * ref_step
@@ -183,93 +183,118 @@ cb = ControlBox(box, cbpos, cblength)
nu = 0.1
coeffForce = 1.
# Monitor for forces
-forces = DragAndLift(velo, vorti, nu, coeffForce, topofft, cb,
- obstacle=[sphere],
- filename=outputdir + 'forces.dat')
-
-## Simulation with fixed time step
-simu = Simulation(tinit=0.0,
- tend=5., timeStep=0.005,
- iterMax=1000000)
-
-
-curl.setUp()
-advec.setUp()
-stretch.setUp()
-diffusion.setUp()
-poisson.setUp()
-correction.setUp()
-penal.setUp()
-distrAdvStr_velo.setUp()
-distrAdvStr_vorti.setUp()
-distrStrDiff.setUp()
-distrCurlStr_velo.setUp()
-distrCurlAdv.setUp()
-
-energy.setUp()
-
-
-## Initialization of velocity on topofft
-def run():
- penal.apply(simu)
- ## end of time loop ##
-
- ## from parmepy.mpi import MPI
- ## print "Start computation ..."
- ## time = MPI.Wtime()
- if topofft.rank == 0:
- print "start curl"
- curl.apply(simu)
-
- #printerFFT.apply(simu)
- # From topo fft to topo advec
- if topofft.rank == 0:
- print "start d0"
- distrCurlAdv.apply(simu)
- if topofft.rank == 0:
- print "start advec"
- advec.apply(simu)
-
- # From topo adv to topo stretch
- if topofft.rank == 0:
- print "start d1"
- distrAdvStr_velo.apply(simu)
- if topofft.rank == 0:
- print "start d1"
- distrAdvStr_vorti.apply(simu)
- if topofft.rank == 0:
- print "start stret"
- stretch.apply(simu)
-
- if topofft.rank == 0:
- print "start d2"
- distrStrDiff.apply(simu)
- if topofft.rank == 0:
- print "start diff"
- diffusion.apply(simu)
- if topofft.rank == 0:
- print "start poiss"
- poisson.apply(simu)
- if topofft.rank == 0:
- print "start correc"
- correction.apply(simu)
- if topofft.rank == 0:
- print "start penal"
- penal.apply(simu)
- energy.apply(simu)
- if topofft.rank == 0:
- print "start printer"
- printerFFT.apply(simu)
-
-
-#while not simu.isOver:
-for i in xrange(200):
+monitors['forces'] = DragAndLift(velo, vorti, nu, coeffForce, topofft, cb,
+ obstacles=[sphere], io_params={})
+
+# === Setup for all declared operators ===
+for ope in op.values():
+ ope.setUp()
+for ope in distr.values():
+ ope.setUp()
+# Set up for monitors
+for monit in monitors.values():
+ monit.setUp()
+
+allops = dict(op.items() + distr.items() + monitors.items())
+
+
+## === Fields initialization ===
+# Initialisation, mode 1:
+# - compute velocity on topofft
+# - penalize
+# - compute vorticity with curl
+def initFields_mode1():
+ # The velocity is initialized on the fftw topology
+ # 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 on topofft
+ distr['curl2fft'].apply(simu)
+ # From this point both velocity and vorticity are initialized
+ # on topofft and on topocurl
+
+## Call initialization
+initFields_mode1()
+
+##### Check initialize process results #####
+norm_vel_fft = velo.norm(topofft)
+norm_vel_curl = velo.norm(topocurl)
+norm_vort_fft = vorti.norm(topofft)
+norm_vort_curl = vorti.norm(topocurl)
+
+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)
+
+
+## Note Franck : tests init ok, same values on both topologies, for v and w.
+
+# === Definition of the 'one-step' sequence of operators ===
+#
+# 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
+#
+fullseq = ['curl', 'printerCurl', 'advection', 'stretching', 'diffusion',
+ 'reprojection', 'poisson', 'correction', 'energy', 'forces',
+ 'printerFFT']
+
+
+def run(sequence):
+ for name in sequence:
+ assert allops.keys().count(name) > 0 or distr.keys().count(name) > 0\
+ or monitors.keys().count(name) > 0
+ allops[name].apply(simu)
+
+seq = fullseq
+simu.initialize()
+while not simu.isOver:
simu.printState()
- run()
+ run(seq)
simu.advance()
## print 'total time (rank):', MPI.Wtime() - time, '(', topo.rank, ')'
+
+print 'aaaa', velo.norm(topofft)
+op['correction'].apply(simu)
+print 'iiii', velo.norm(topofft)
+
+sref = op['correction'].discreteOperator.surfRef
+
+flowrate = velo.integrateOnSurface(sref, topofft)
+print flowrate
+assert (flowrate - ref_rate) < 1e-10
+
+flowratey = velo.integrateOnSurface(sref, topofft, component=1)
+print flowratey
+flowratez = velo.integrateOnSurface(sref, topofft, component=2)
+print flowratez
+
+
+# === Finalize for all declared operators ===
+
+# Note FP : bug in fft finalize. To be checked ...
+# --> clean_fftw must be called only once
+#for ope in op.values():
+# ope.finalize()
## Clean memory buffers
fftw2py.clean_fftw_solver(box.dimension)
+for ope in distr.values():
+ ope.finalize()
+for monit in monitors.values():
+ monit.finalize()
+
diff --git a/HySoP/hysop/__init__.py.in b/HySoP/hysop/__init__.py.in
index f680e04c9d5e08c39b522340178e6727917cf288..7d7e223902ffec4e7a71b69c633e3ea8c8a99258 100755
--- a/HySoP/hysop/__init__.py.in
+++ b/HySoP/hysop/__init__.py.in
@@ -15,13 +15,13 @@ __FFTW_ENABLED__ = "@WITH_FFTW@" is "ON"
__SCALES_ENABLED__ = "@WITH_SCALES@" is "ON"
__VERBOSE__ = "@DEBUG@" in ["1", "3"]
__DEBUG__ = "@DEBUG@" in ["2", "3"]
-__PROFILE__= "@PROFILE@" in ["0", "1"]
-__OPTIMIZE__= "@OPTIM@" is "ON"
+__PROFILE__ = "@PROFILE@" in ["0", "1"]
+__OPTIMIZE__ = "@OPTIM@" is "ON"
# MPI
if __MPI_ENABLED__:
#import parmepy.mpi as mpi
# from mpi import main_rank
- import mpi
+ import parmepy.mpi
if(mpi.main_rank == 0):
print "Starting @PACKAGE_NAME@ version " + str(__version__) + ".\n"
else:
@@ -30,13 +30,20 @@ else:
version = "1.0.0"
## Box-type physical domain
-import domain.box
-Box = domain.box.Box
+import parmepy.domain.box
+Box = parmepy.domain.box.Box
## Fields
-import fields.continuous
-Field = fields.continuous.Field
+import parmepy.fields.continuous
+Field = parmepy.fields.continuous.Field
+## Variable parameters
+import parmepy.variable_parameter
+VariableParameter = parmepy.variable_parameter.VariableParameter
+
+## Simulation parameters
+import parmepy.problem.simulation
+Simulation = parmepy.problem.simulation.Simulation
# ## ## Problem
# import problem.problem
diff --git a/HySoP/hysop/domain/obstacle/planes.py b/HySoP/hysop/domain/obstacle/planes.py
index 6791bb23ff5d245e64fe7ca0a1756096ebf9a609..961ce968bb862d673929490262b6a37e6daf6fa6 100644
--- a/HySoP/hysop/domain/obstacle/planes.py
+++ b/HySoP/hysop/domain/obstacle/planes.py
@@ -115,6 +115,7 @@ class SubSpace(HalfSpace):
self.dist = dist
self.origin = npw.realarray(point)
self.max = self.origin + npw.realarray(lengths)
+ self.lengths = npw.realarray(lengths)
ndir = np.where(self.normal != 0)[0][0]
if normal[ndir] > 0:
self.max[ndir] = self.origin[ndir]
diff --git a/HySoP/hysop/domain/tests/test_obstacle.py b/HySoP/hysop/domain/tests/test_obstacle.py
index c40c2a37eaa1850ed85ac7e7633bac52e8ac9953..f7810b7ff11a0b5459aef5814d12111f459dec4d 100644
--- a/HySoP/hysop/domain/tests/test_obstacle.py
+++ b/HySoP/hysop/domain/tests/test_obstacle.py
@@ -33,7 +33,7 @@ dvol = np.prod(h3d)
ds = np.prod(h2d)
import math
pi = math.pi
-tol = 1e-3
+tol = 1e-6
lengths = np.asarray([20 * h3d[0], 22 * h3d[1], 31 * h3d[2]])
rlengths = lengths + h3d
rlengths2d = lengths[:2] + h2d
diff --git a/HySoP/hysop/fields/continuous.py b/HySoP/hysop/fields/continuous.py
index a60a3d6611f784786e4ef5d9fbc93065a3811dc5..1e5007e8a6fd06a35c36794ef18c44a348ff55fd 100644
--- a/HySoP/hysop/fields/continuous.py
+++ b/HySoP/hysop/fields/continuous.py
@@ -6,6 +6,8 @@ Continuous variable description.
from parmepy.constants import debug
from parmepy.fields.discrete import DiscreteField
from parmepy.mpi import main_rank
+import parmepy.tools.numpywrappers as npw
+import numpy as np
class Field(object):
@@ -256,3 +258,36 @@ class Field(object):
def finalize(self):
pass
+
+ def integrate(self, box, topo, useSlice=True,
+ component=0, root=0, mpiall=True):
+ """
+ integrate the field on a control box, on the current processus
+ @param box : a parmepy.domain.obstacles.controlBox.ControlBox
+ @param topo : discretization used for integration
+ @param useSlice : true if integrate with slices else integrate with
+ boolean numpy array
+ @param component : component number of the field to integrate
+ @param root : root process for mpi reduce
+ @param mpiall : true for allreduce, else reduce on root
+ """
+ return self.discreteFields[topo].integrate(box, useSlice, component,
+ root, mpiall)
+
+ def integrateOnSurface(self, surf, topo, useSlice=True, component=0,
+ root=0, mpiall=True):
+ """
+ integrate the field on a surface, on the current processus
+ @param surf : a parmepy.domain.obstacles.controlBox.planes.SubPlane
+ @param topo : discretization used for integration
+ @param useSlice : true if integrate with slices else integrate with
+ boolean numpy array
+ @param component : component number of the field to integrate
+ @param root : root process for mpi reduce
+ @param mpiall : true for allreduce, else reduce on root
+ Warning : surf slices or cond must have been computed for integration
+ purpose, that is last point in each dir must not be included.
+ """
+ return self.discreteFields[topo].integrateOnSurface(surf, useSlice,
+ component, root,
+ mpiall)
diff --git a/HySoP/hysop/fields/discrete.py b/HySoP/hysop/fields/discrete.py
index 60704c3ef653b582905006540d813a7f77429249..cdf9abfa679d1ef49a17acc7460ec9332dd12a1b 100644
--- a/HySoP/hysop/fields/discrete.py
+++ b/HySoP/hysop/fields/discrete.py
@@ -15,7 +15,8 @@ from parmepy.numerics.updateGhosts import UpdateGhosts
class DiscreteField(object):
- """Interface for discrete fields (scalar or vector).
+ """
+ Interface for discrete fields (scalar or vector).
A discrete field represents the data layout on each mpi process for
a given global resolution.
@@ -187,10 +188,10 @@ class DiscreteField(object):
# No formula, set all components to zero"
for d in xrange(self.nbComponents):
self.data[d][...] = 0.0
- if not all([(s == self.resolution).all()
- for s in [dat.shape for dat in self.data]]):
- print "WARNING: The shape of " + self.name + " has changed during"\
- " field initialisation."
+ assert all([(s == self.resolution).all()
+ for s in [dat.shape for dat in self.data]]),\
+ "WARNING: The shape of " + self.name + " has changed during"\
+ " field initialisation."
def norm(self):
"""
@@ -209,9 +210,6 @@ class DiscreteField(object):
self.topology.comm.Allreduce(result, gResult)
return gResult ** (1. / 2)
- def get_data_method(self):
- pass
-
@debug
@timed_function
def dump(self, filename, mode=None):
@@ -253,3 +251,79 @@ class DiscreteField(object):
""" set all components to zero"""
for dim in xrange(self.nbComponents):
self.data[dim][...] = 0.0
+
+ def integrate_on_proc(self, box, useSlice=True, component=0):
+ """
+ integrate the field on a control box, on the current processus
+ @param box : a parmepy.domain.obstacles.controlBox.ControlBox
+ @param useSlice : true if integrate with slices else integrate with
+ boolean numpy array
+ @param component : component number of the field to integrate
+ """
+ if useSlice:
+ cond = box.slices[self.topology]
+ else:
+ iC = self.topology.mesh.iCompute
+ cond = box.ind[self.topology][0][iC]
+ dvol = npw.prod(self.topology.mesh.space_step)
+ result = npw.sum(self.data[component][cond])
+ result *= dvol
+ return result
+
+ def integrate(self, box, useSlice=True, component=0,
+ root=0, mpiall=True):
+ """
+ integrate the field on a control box, on the current processus
+ @param box : a parmepy.domain.obstacles.controlBox.ControlBox
+ @param useSlice : true if integrate with slices else integrate with
+ boolean numpy array
+ @param component : component number of the field to integrate
+ @param root : root process for mpi reduce
+ @param mpiall : true for allreduce, else reduce on root
+ """
+ res = self.integrate_on_proc(box, useSlice, component)
+ if mpiall:
+ return self.topology.comm.allreduce(res)
+ else:
+ return self.topology.comm.reduce(res, root=root)
+
+ def integrateOnSurface(self, surf, useSlice=True, component=0,
+ root=0, mpiall=True):
+ """
+ integrate the field on a surface, on the current processus
+ @param surf : a parmepy.domain.obstacles.controlBox.planes.SubPlane
+ @param useSlice : true if integrate with slices else integrate with
+ boolean numpy array
+ @param component : component number of the field to integrate
+ @param root : root process for mpi reduce
+ @param mpiall : true for allreduce, else reduce on root
+ Warning : surf slices or cond must have been computed for integration
+ purpose, that is last point in each dir must not be included.
+ """
+ res = self.integrateOnSurf_proc(surf, useSlice, component)
+ if mpiall:
+ return self.topology.comm.allreduce(res)
+ else:
+ return self.topology.comm.reduce(res, root=root)
+
+ def integrateOnSurf_proc(self, surf, useSlice=True, component=0):
+ """
+ integrate the field on a surface, on the current processus
+ @param surf : a parmepy.domain.obstacles.controlBox.planes.Plane
+ or SubPlane
+ @param useSlice : true if integrate with slices else integrate with
+ boolean numpy array
+ @param component : component number of the field to integrate
+ Warning : surf slices or cond must have been computed for integration
+ purpose, that is last point in each dir must not be included.
+ """
+ if useSlice:
+ cond = surf.slices[self.topology]
+ else:
+ iC = self.topology.mesh.iCompute
+ cond = surf.ind[self.topology][0][iC]
+ dirs = np.where(surf.normal == 0)[0]
+ dS = npw.prod(self.topology.mesh.space_step[dirs])
+ result = npw.sum(self.data[component][cond])
+ result *= dS
+ return result
diff --git a/HySoP/hysop/fields/tests/test_field.py b/HySoP/hysop/fields/tests/test_field.py
index 631a81994415a5c59be108bbd276ae6c654a5e5c..020e8df0fd6ed393d8eae24d1f8ce805cadd6ff8 100644
--- a/HySoP/hysop/fields/tests/test_field.py
+++ b/HySoP/hysop/fields/tests/test_field.py
@@ -4,6 +4,7 @@ Testing parmepy.field.continuous.Field
from parmepy.fields.continuous import Field
from parmepy.domain.box import Box
from parmepy.mpi.topology import Cartesian
+from parmepy.domain.obstacle.planes import SubPlane
import sys
import numpy as np
@@ -63,3 +64,22 @@ def test_discretization():
cvf.discreteFields[topo].data[1].shape).all()
assert np.equal(cvf.discreteFields[topo].resolution,
cvf.discreteFields[topo].data[2].shape).all()
+
+
+def test_integrate_onSurf():
+ dom = Box()
+ myf = Field(dom, isVector=True)
+ resolTopo = [33, 33, 17]
+ topo = Cartesian(dom, 3, resolTopo)
+ fdiscr = myf.discretize(topo)
+ fdiscr[0][...] = 1.0
+ normal = [1, 0, 0]
+ hh = topo.mesh.space_step
+ surf = SubPlane(dom, normal=normal, point=dom.origin,
+ lengths=dom.length - hh,
+ epsilon=topo.mesh.space_step[0]/2.)
+
+ surf.discretize(topo)
+ res = myf.integrateOnSurface(surf, topo)
+ sref = dom.length[1] * dom.length[2]
+ assert(abs(res - sref) < 1e-6)
diff --git a/HySoP/hysop/gpu/QtRendering.py b/HySoP/hysop/gpu/QtRendering.py
index 81d7a02a05aa17a6fabbd3a37f44b011fd29e1b8..4a3ccfdbc702a636be66c2305e35ce1d66882843 100644
--- a/HySoP/hysop/gpu/QtRendering.py
+++ b/HySoP/hysop/gpu/QtRendering.py
@@ -14,6 +14,7 @@ from parmepy.gpu import cl
from parmepy.gpu.gpu_discrete import GPUDiscreteField
from parmepy.gpu.gpu_kernel import KernelLauncher
from parmepy.tools.timers import timed_function
+from parmepy.mpi import main_rank
class QtOpenGLRendering(Monitoring):
@@ -164,7 +165,8 @@ class QtOpenGLRendering(Monitoring):
"""
t = simulation.time
dt = simulation.timeStep
- simulation.printState()
+ if main_rank == 0:
+ simulation.printState()
# OpenCL update
self.numMethod(self.gpu_field.gpu_data[self.component],
self.color)
diff --git a/HySoP/hysop/gpu/tests/test_advection_nullVelocity.py b/HySoP/hysop/gpu/tests/test_advection_nullVelocity.py
index 1613faaf644845446e0bcf79fd54837bd6854f48..7dabdfc4e430fda423ab6b5d25703ac5714807bb 100644
--- a/HySoP/hysop/gpu/tests/test_advection_nullVelocity.py
+++ b/HySoP/hysop/gpu/tests/test_advection_nullVelocity.py
@@ -44,7 +44,7 @@ def assertion_2D(scal, velo, advec):
scal_d.toDevice()
velo_d.toDevice()
- advec.apply(Simulation(0., 0.01, 1))
+ advec.apply(Simulation(0., 1., 0.01))
scal_d.toHost()
@@ -69,8 +69,8 @@ def assertion_2D_withPython(scal, velo, advec, advec_py):
scal_d.toDevice()
velo_d.toDevice()
- advec_py.apply(Simulation(0., 0.01, 1))
- advec.apply(Simulation(0., 0.01, 1))
+ advec_py.apply(Simulation(0., 1., 0.01))
+ advec.apply(Simulation(0., 1., 0.01))
py_res = scal_d.data[0].copy()
scal_d.toHost()
@@ -99,7 +99,7 @@ def assertion_3D(scal, velo, advec):
scal_d.toDevice()
velo_d.toDevice()
- advec.apply(Simulation(0., 0.01, 1))
+ advec.apply(Simulation(0., 1., 0.01))
scal_d.toHost()
@@ -126,8 +126,8 @@ def assertion_3D_withPython(scal, velo, advec, advec_py):
scal_d.toDevice()
velo_d.toDevice()
- advec_py.apply(Simulation(0., 0.01, 1))
- advec.apply(Simulation(0., 0.01, 1))
+ advec_py.apply(Simulation(0., 1., 0.01))
+ advec.apply(Simulation(0., 1., 0.01))
py_res = scal_d.data[0].copy()
scal_d.toHost()
@@ -954,8 +954,8 @@ def test_rectangular_domain2D():
scal_d.toDevice()
velo_d.toDevice()
- advec.apply(Simulation(0., 0.01, 1))
- advec_py.apply(Simulation(0., 0.01, 1))
+ advec.apply(Simulation(0., 1., 0.01))
+ advec_py.apply(Simulation(0., 1., 0.01))
scal_py_res = scal_d.data[0].copy()
scal_d.toHost()
@@ -1008,8 +1008,8 @@ def test_rectangular_domain3D():
scal_d.toDevice()
velo_d.toDevice()
- advec.apply(Simulation(0., 0.01, 1))
- advec_py.apply(Simulation(0., 0.01, 1))
+ advec.apply(Simulation(0., 1., 0.01))
+ advec_py.apply(Simulation(0., 1., 0.01))
scal_py_res = scal_d.data[0].copy()
scal_d.toHost()
@@ -1063,8 +1063,8 @@ def test_2D_vector():
scal_d.toDevice()
velo_d.toDevice()
- advec.apply(Simulation(0., 0.01, 1))
- advec_py.apply(Simulation(0., 0.01, 1))
+ advec.apply(Simulation(0., 1., 0.01))
+ advec_py.apply(Simulation(0., 1., 0.01))
scal_py_res_X = scal_d.data[0].copy()
scal_py_res_Y = scal_d.data[1].copy()
@@ -1126,8 +1126,8 @@ def test_3D_vector():
scal_d.toDevice()
velo_d.toDevice()
- advec.apply(Simulation(0., 0.01, 1))
- advec_py.apply(Simulation(0., 0.01, 1))
+ advec.apply(Simulation(0., 1., 0.01))
+ advec_py.apply(Simulation(0., 1., 0.01))
scal_py_res_X = scal_d.data[0].copy()
scal_py_res_Y = scal_d.data[1].copy()
scal_py_res_Z = scal_d.data[2].copy()
diff --git a/HySoP/hysop/gpu/tests/test_advection_randomVelocity.py b/HySoP/hysop/gpu/tests/test_advection_randomVelocity.py
index 7e1d32982cc5872a3235fa2eda0820a4c6b9095f..e2165bfb0a40e38488bd8ea276203d24d8c9d578 100644
--- a/HySoP/hysop/gpu/tests/test_advection_randomVelocity.py
+++ b/HySoP/hysop/gpu/tests/test_advection_randomVelocity.py
@@ -2,7 +2,6 @@
@file parmepy.gpu.tests.test_advection_randomVelocity
Testing advection kernels with a random velocity field.
"""
-import parmepy
from parmepy.domain.box import Box
from parmepy.fields.continuous import Field
from parmepy.operator.advection import Advection
@@ -49,8 +48,8 @@ def assertion_2D_withPython(scal, velo, advec, advec_py):
scal_d.toDevice()
velo_d.toDevice()
- advec_py.apply(Simulation(0., 0.001, 1))
- advec.apply(Simulation(0., 0.001, 1))
+ advec_py.apply(Simulation(0., 1., 0.001))
+ advec.apply(Simulation(0., 1., 0.001))
py_res = scal_d.data[0].copy()
scal_d.toHost()
@@ -82,8 +81,8 @@ def assertion_3D_withPython(scal, velo, advec, advec_py):
scal_d.toDevice()
velo_d.toDevice()
- advec_py.apply(Simulation(0., 0.1, 1))
- advec.apply(Simulation(0., 0.1, 1))
+ advec_py.apply(Simulation(0., 1., 0.1))
+ advec.apply(Simulation(0., 1., 0.1))
py_res = scal_d.data[0].copy()
scal_d.toHost()
@@ -834,8 +833,8 @@ def test_rectangular_domain2D():
scal_d.toDevice()
velo_d.toDevice()
- advec_py.apply(Simulation(0., 0.1, 1))
- advec.apply(Simulation(0., 0.1, 1))
+ advec_py.apply(Simulation(0., 1., 0.1))
+ advec.apply(Simulation(0., 1., 0.1))
py_res = scal_d.data[0].copy()
scal_d.toHost()
@@ -889,8 +888,8 @@ def test_rectangular_domain3D():
scal_d.toDevice()
velo_d.toDevice()
- advec_py.apply(Simulation(0., 0.01, 1))
- advec.apply(Simulation(0., 0.01, 1))
+ advec_py.apply(Simulation(0., 1., 0.01))
+ advec.apply(Simulation(0., 1., 0.01))
py_res = scal_d.data[0].copy()
scal_d.toHost()
@@ -945,8 +944,8 @@ def test_vector_2D():
velo_d.toDevice()
print np.max(scal_d.data[0] - scal_d.data[1])
- advec_py.apply(Simulation(0., 0.01, 1))
- advec.apply(Simulation(0., 0.01, 1))
+ advec_py.apply(Simulation(0., 1., 0.01))
+ advec.apply(Simulation(0., 1., 0.01))
py_res_X = scal_d.data[0].copy()
py_res_Y = scal_d.data[1].copy()
@@ -1006,8 +1005,8 @@ def test_vector_3D():
scal_d.toDevice()
velo_d.toDevice()
- advec_py.apply(Simulation(0., 0.01, 1))
- advec.apply(Simulation(0., 0.01, 1))
+ advec_py.apply(Simulation(0., 1., 0.01))
+ advec.apply(Simulation(0., 1., 0.01))
py_res_X = scal_d.data[0].copy()
py_res_Y = scal_d.data[1].copy()
@@ -1018,4 +1017,3 @@ def test_vector_3D():
assert np.allclose(py_res_Y, scal_d.data[1], rtol=5e-02, atol=5e-05)
assert np.allclose(py_res_Z, scal_d.data[2], rtol=5e-02, atol=5e-05)
advec.finalize()
-
diff --git a/HySoP/hysop/gpu/tests/test_opencl_environment.py b/HySoP/hysop/gpu/tests/test_opencl_environment.py
index 1643ba3bda76d3a22382621053606afb50f2ca0a..46931e9af3231460586d6c9a87e2987fd1eff1b7 100644
--- a/HySoP/hysop/gpu/tests/test_opencl_environment.py
+++ b/HySoP/hysop/gpu/tests/test_opencl_environment.py
@@ -1,6 +1,5 @@
from parmepy.constants import np
-from parmepy.gpu.tools import get_opencl_environment, \
- get_opengl_shared_environment
+from parmepy.gpu.tools import get_opencl_environment
FLOAT_GPU = np.float32
diff --git a/HySoP/hysop/gpu/tests/test_transposition.py b/HySoP/hysop/gpu/tests/test_transposition.py
index 28873c2d34f584ac9193999e467d989569b5c4a3..fd245a3a5d973b7a64ffde136a22b2e6cd85bfe3 100644
--- a/HySoP/hysop/gpu/tests/test_transposition.py
+++ b/HySoP/hysop/gpu/tests/test_transposition.py
@@ -256,7 +256,8 @@ def test_transposition_xz3D():
src_transpose_xz = 'kernels/transpose_xz_noVec.cl'
build_options = ""
build_options += " -D NB_I=32 -D NB_II=32 -D NB_III=32"
- build_options += " -D TILE_DIM_XZ=16 -D BLOCK_ROWS_XZ=4 -D BLOCK_DEPH_XZ=4 -D PADDING_XZ=1"
+ build_options += " -D TILE_DIM_XZ=16 -D BLOCK_ROWS_XZ=4"
+ build_options += " -D BLOCK_DEPH_XZ=4 -D PADDING_XZ=1"
gwi = (int(resolution[0]),
int(resolution[1] / 4),
int(resolution[2] / 4))
@@ -302,7 +303,6 @@ def test_transposition_xz3D():
def test_transposition_xz3D_rect():
resolution = (32, 32, 64)
- resolutionT = (64, 32, 32)
cl_env = get_opencl_environment(0, 0, 'gpu', PARMES_REAL)
vec = 1
src_transpose_xz = 'kernels/transpose_xz_noVec.cl'
@@ -310,7 +310,8 @@ def test_transposition_xz3D_rect():
# Settings are taken from destination layout as current layout.
# gwi is computed form input layout (appears as transposed layout)
build_options += " -D NB_I=64 -D NB_II=32 -D NB_III=32"
- build_options += " -D TILE_DIM_XZ=16 -D BLOCK_ROWS_XZ=4 -D BLOCK_DEPH_XZ=4 -D PADDING_XZ=1"
+ build_options += " -D TILE_DIM_XZ=16 -D BLOCK_ROWS_XZ=4"
+ build_options += " -D BLOCK_DEPH_XZ=4 -D PADDING_XZ=1"
gwi = (int(resolution[0]),
int(resolution[1] / 4),
int(resolution[2] / 4))
@@ -321,7 +322,8 @@ def test_transposition_xz3D_rect():
build_options = ""
build_options += " -D NB_I=32 -D NB_II=32 -D NB_III=64"
- build_options += " -D TILE_DIM_XZ=16 -D BLOCK_ROWS_XZ=4 -D BLOCK_DEPH_XZ=4 -D PADDING_XZ=1"
+ build_options += " -D TILE_DIM_XZ=16 -D BLOCK_ROWS_XZ=4"
+ build_options += " -D BLOCK_DEPH_XZ=4 -D PADDING_XZ=1"
gwi = (int(resolution[2]),
int(resolution[1] / 4),
int(resolution[0] / 4))
@@ -373,7 +375,8 @@ def test_transposition_xz3Dslice():
src_transpose_xz = 'kernels/transpose_xz_slice_noVec.cl'
build_options = ""
build_options += " -D NB_I=32 -D NB_II=32 -D NB_III=32"
- build_options += " -D TILE_DIM_XZ=16 -D BLOCK_ROWS_XZ=1 -D BLOCK_DEPH_XZ=4 -D PADDING_XZ=1"
+ build_options += " -D TILE_DIM_XZ=16 -D BLOCK_ROWS_XZ=1"
+ build_options += " -D BLOCK_DEPH_XZ=4 -D PADDING_XZ=1"
gwi = (int(resolution[0]),
int(resolution[1]),
int(resolution[2] / 4))
@@ -419,7 +422,6 @@ def test_transposition_xz3Dslice():
def test_transposition_xz3Dslice_rect():
resolution = (32, 32, 64)
- resolutionT = (64, 32, 32)
cl_env = get_opencl_environment(0, 0, 'gpu', PARMES_REAL)
vec = 1
src_transpose_xz = 'kernels/transpose_xz_slice_noVec.cl'
@@ -427,7 +429,8 @@ def test_transposition_xz3Dslice_rect():
# Settings are taken from destination layout as current layout.
# gwi is computed form input layout (appears as transposed layout)
build_options += " -D NB_I=64 -D NB_II=32 -D NB_III=32"
- build_options += " -D TILE_DIM_XZ=16 -D BLOCK_ROWS_XZ=1 -D BLOCK_DEPH_XZ=4 -D PADDING_XZ=1"
+ build_options += " -D TILE_DIM_XZ=16 -D BLOCK_ROWS_XZ=1"
+ build_options += " -D BLOCK_DEPH_XZ=4 -D PADDING_XZ=1"
gwi = (int(resolution[0]),
int(resolution[1]),
int(resolution[2] / 4))
@@ -438,7 +441,8 @@ def test_transposition_xz3Dslice_rect():
build_options = ""
build_options += " -D NB_I=32 -D NB_II=32 -D NB_III=64"
- build_options += " -D TILE_DIM_XZ=16 -D BLOCK_ROWS_XZ=1 -D BLOCK_DEPH_XZ=4 -D PADDING_XZ=1"
+ build_options += " -D TILE_DIM_XZ=16 -D BLOCK_ROWS_XZ=1"
+ build_options += " -D BLOCK_DEPH_XZ=4 -D PADDING_XZ=1"
gwi = (int(resolution[2]),
int(resolution[1]),
int(resolution[0] / 4))
diff --git a/HySoP/hysop/gpu/tools.py b/HySoP/hysop/gpu/tools.py
index b33a41a7385c9203c333183cd1f007890e23db17..e6ec65a2942f622f38c485f4f63d8df35d002e6c 100644
--- a/HySoP/hysop/gpu/tools.py
+++ b/HySoP/hysop/gpu/tools.py
@@ -72,7 +72,6 @@ class OpenCLEnvironment(object):
self.kernels_config = kernel_cfg
self._locMem_Buffers = {}
-
def modify(self, platform_id, device_id, device_type,
precision, gl_sharing=False):
"""
diff --git a/HySoP/hysop/mpi/topology.py b/HySoP/hysop/mpi/topology.py
index 3f8fc6b37d953c4b366687fe9d1af4582e8ac471..4ddc7e3f62b03ffd3ed97be6511ca7f06a909b40 100644
--- a/HySoP/hysop/mpi/topology.py
+++ b/HySoP/hysop/mpi/topology.py
@@ -75,7 +75,7 @@ class Cartesian(object):
# to 0 if shape[:] = 1
## (Source) Communicator used to build the topology
if comm is None:
- from parmepy.mpi.main_var import main_comm as comm
+ from parmepy.mpi.main_var import main_comm as comm
self._comm_origin = comm
## number of process in comm_origin
self._origin_size = self._comm_origin.Get_size()
@@ -215,10 +215,6 @@ class Cartesian(object):
newId = self.domain.register(self)
self.isNew = True
if newId >= 0:
- s = "A similar topology already exist in the domain's list."
- s += "You'd rather destroy this one and use topo of id"
- s += str(newId)
- print s
self.isNew = False
def getParent(self):
@@ -227,12 +223,6 @@ class Cartesian(object):
"""
return self._comm_origin
- def setParent(self, comm):
- """
- returns the communicator used to build this topology
- """
- self._comm_origin = comm
-
@debug
def setUp(self):
"""Topology set up.
@@ -697,8 +687,6 @@ class Bridge(object):
return s
-
-
class Bridge_intercomm(object):
"""
This bridge is defined as the Bridge class.
diff --git a/HySoP/hysop/numerics/differential_operations.py b/HySoP/hysop/numerics/differential_operations.py
index 9773a8e19a7439f434a1bcfaef3a6ac645925154..1575827da586daae58298878362d80a76d17871c 100755
--- a/HySoP/hysop/numerics/differential_operations.py
+++ b/HySoP/hysop/numerics/differential_operations.py
@@ -6,7 +6,7 @@ Library of functions used to perform classical vector calculus
"""
from parmepy.constants import debug, XDIR, YDIR, ZDIR
from abc import ABCMeta, abstractmethod
-from parmepy.numerics.finite_differences import FD_C_4, FD2_C_2, FD_C_2
+from parmepy.numerics.finite_differences import FD_C_4, FD2_C_2
import numpy as np
@@ -326,107 +326,6 @@ class Laplacian(DifferentialOperation):
return result
-class DivStressTensor3D(DifferentialOperation):
- """
- Computes the 3D stress tensor T defined by
- \f$ T=div(\nabla.(U) + \nabla^{T}.(U)) \f$,
- \f$U\f$ being a velocity vector field.
-
- Methods :
- 1 - FD_C_4 (default) : 4th order, centered finite differences, with
- local workspace and based on fd.compute_and_add.
- init : func = GradS(topo, method)
- call : result = func(var, result)
-
- var stands for the vector field U.
- result must be a vector field (same number of components as var),
- each component with the same size as var[i].
- work must be a three component vector field
- each component with the same size as var[i].
-
- """
- def __init__(self, topo, method=FD_C_2):
- if method is FD_C_4:
- raise ValueError("4th order scheme Not yet implemented")
- else :
- # - 2nd ordered FD
- # - with work
- # - fd.compute_and_add method
-
- # connect call function
- self.fcall = self.FDCentral2
- self.fd_scheme = FD_C_2((topo.mesh.space_step))
-
- @staticmethod
- def getWorkLengths(nb_components=None, domain_dim=None, fd_method=None):
- return 0
-
- def __call__(self, var, result, ind):
- assert result[0].shape == var[0].shape
- return self.fcall(var, result, ind)
-
- def FDCentral2(self, var, result, ind):
- """
- """
- # compute indices according to "ind" list
- self._indices = ind
- self.fd_scheme.computeIndices(self._indices)
- self.fd_scheme.computeIndices_crossed(self._indices)
-
- # 1st component of the stress tensor
- # ----- 2 * d2ux/dx2
- self.fd_scheme.compute_2nd_deriv(var[XDIR], XDIR, result[XDIR])
- result[XDIR][...] *= 2.
- # ----- += d2ux/dy2
- self.fd_scheme.compute_and_add_2nd_deriv(var[XDIR], YDIR,
- result[XDIR])
- # ----- += d2ux/dz2
- self.fd_scheme.compute_and_add_2nd_deriv(var[XDIR], ZDIR,
- result[XDIR])
- # ----- += d2uy/dxdy
- self.fd_scheme.compute_and_add_2nd_crossed_deriv(var[YDIR], XDIR,
- YDIR, result[XDIR])
- # ----- += d2uz/dxdz
- self.fd_scheme.compute_and_add_2nd_crossed_deriv(var[ZDIR], XDIR,
- ZDIR, result[XDIR])
-
- # 2nd component of the stress tensor
- # ----- 2 * d2uy/dy2
- self.fd_scheme.compute_2nd_deriv(var[YDIR], YDIR, result[YDIR])
- result[YDIR][...] *= 2.
- # ----- += d2uy/dx2
- self.fd_scheme.compute_and_add_2nd_deriv(var[YDIR], XDIR,
- result[YDIR])
- # ----- += d2uy/dz2
- self.fd_scheme.compute_and_add_2nd_deriv(var[YDIR], ZDIR,
- result[YDIR])
- # ----- += d2ux/dydx
- self.fd_scheme.compute_and_add_2nd_crossed_deriv(var[XDIR], YDIR,
- XDIR, result[YDIR])
- # ----- += d2uz/dydz
- self.fd_scheme.compute_and_add_2nd_crossed_deriv(var[ZDIR], YDIR,
- ZDIR, result[YDIR])
-
- # 3rd component of the stress tensor
- # ----- 2 * d2uz/dz2
- self.fd_scheme.compute_2nd_deriv(var[ZDIR], ZDIR, result[ZDIR])
- result[ZDIR][...] *= 2.
- # ----- += d2uz/dx2
- self.fd_scheme.compute_and_add_2nd_deriv(var[ZDIR], XDIR,
- result[ZDIR])
- # ----- += d2uz/dy2
- self.fd_scheme.compute_and_add_2nd_deriv(var[ZDIR], YDIR,
- result[ZDIR])
- # ----- += d2ux/dzdx
- self.fd_scheme.compute_and_add_2nd_crossed_deriv(var[XDIR], ZDIR,
- XDIR, result[ZDIR])
- # ----- += d2uy/dzdy
- self.fd_scheme.compute_and_add_2nd_crossed_deriv(var[YDIR], ZDIR,
- YDIR, result[ZDIR])
-
- return result
-
-
class GradS(DifferentialOperation):
"""
Computes \f$ \nabla\rho \f$, \f$ \rho\f$ a scalar field
diff --git a/HySoP/hysop/numerics/finite_differences.py b/HySoP/hysop/numerics/finite_differences.py
index a260da267294f76d43900c6e447169cf78737993..bc05d712a6ce3ab699c29d48e5ed014bb1a53007 100644
--- a/HySoP/hysop/numerics/finite_differences.py
+++ b/HySoP/hysop/numerics/finite_differences.py
@@ -19,14 +19,14 @@ class FiniteDifference(object):
and initialize it with a step from domain discretisation description
\code
- >>> step = topo.mesh_space_step
- >>> scheme = FD_C_4(step)
+ >> step = topo.mesh_space_step
+ >> scheme = FD_C_4(step)
\endcode
- Compute scheme indices, based on topology grid points :
\code
- >>> scheme.computeIndices(topo.mesh.iCompute)
+ >> scheme.computeIndices(topo.mesh.iCompute)
\endcode
2 - Computation
@@ -36,14 +36,14 @@ class FiniteDifference(object):
\f$ result = \frac{\partial tab}{\partial dir}\f$ :
\code
- >>> scheme.compute(tab, dir, result)
+ >> scheme.compute(tab, dir, result)
\endcode
or :
\f$ result = result + \frac{\partial tab}{\partial dir}\f$
\code
- >>> scheme.compute_and_add(tab, dir, result)
+ >> scheme.compute_and_add(tab, dir, result)
\endcode
Notes FP :
@@ -270,4 +270,3 @@ class FD_C_4(FiniteDifference):
tab[self._m1[cdir]]) +
tab[self._m2[cdir]] -
tab[self._a2[cdir]])
-
diff --git a/HySoP/hysop/numerics/tests/test_diffOp.py b/HySoP/hysop/numerics/tests/test_diffOp.py
index ecfa4000777dfd0375c6158ed31bcbc83500c6cb..c1b5b768c4e50d8745df5166dc539b1e5469b67b 100755
--- a/HySoP/hysop/numerics/tests/test_diffOp.py
+++ b/HySoP/hysop/numerics/tests/test_diffOp.py
@@ -4,7 +4,6 @@ import numpy as np
from parmepy.fields.continuous import Field
from parmepy.mpi.topology import Cartesian
import parmepy.numerics.differential_operations as diffop
-from parmepy.numerics.updateGhosts import UpdateGhosts
import parmepy.tools.numpywrappers as npw
import math as m
pi = m.pi
@@ -61,11 +60,10 @@ velo = Field(domain=box, formula=computeVel,
name='Velocity', isVector=True)
vorti = Field(domain=box, formula=computeVort,
name='Vorticity', isVector=True)
-vorti.setTopoInit(topo)
velo.discretize(topo)
vorti.discretize(topo)
-velo.initialize()
-vorti.initialize()
+velo.initialize(topo=topo)
+vorti.initialize(topo=topo)
wd = vorti.discreteFields[topo]
vd = velo.discreteFields[topo]
ind = topo.mesh.iCompute
@@ -88,13 +86,13 @@ def testDivWV():
ref = Field(domain=box, formula=analyticalDivWV,
name='Analytical', isVector=True)
ref.discretize(topo)
- ref.initialize()
+ ref.initialize(topo=topo)
refd = ref.discreteFields[topo]
memshape = vd.data[0].shape
# Div operator
- lwork = diffop.DivV.getWorkLengths()
+ lwork = diffop.DivWV.getWorkLengths()
work = [npw.zeros(memshape) for i in xrange(lwork)]
- divOp = diffop.DivV(topo, work)
+ divOp = diffop.DivWV(topo, work)
result = [npw.zeros(memshape) for i in xrange(3)]
result = divOp(vd.data, wd.data, result)
@@ -110,7 +108,7 @@ def testGradVxW():
ref = Field(domain=box, formula=analyticalGradVxW,
name='Analytical', isVector=True)
ref.discretize(topo)
- ref.initialize()
+ ref.initialize(topo=topo)
refd = ref.discreteFields[topo]
memshape = vd.data[0].shape
lwork = diffop.GradVxW.getWorkLengths()
@@ -125,39 +123,3 @@ def testGradVxW():
errX = (Lx / (nb - 1)) ** 4
for i in xrange(3):
assert np.allclose(refd[i][ind], result[i][ind], rtol=errX)
-
-
-def testDivStressTensor():
- # Topologies
-# topo1G = Cartesian(box, box.dimension, nbElem,
-# ghosts=[2, 2, 2])
- topoNoG = Cartesian(box, box.dimension, nbElem)
- ind = topo.mesh.iCompute
- # Reference field
- ref = Field(domain=box, formula=analyticalDivStressTensor,
- name='Analytical', isVector=True)
- ref.discretize(topoNoG)
- ref.initialize()
- refd = ref.discreteFields[topoNoG]
- memshape = vd.data[0].shape
-
- # prepare ghost points synchro for velocity
- synchronize = UpdateGhosts(topo, velo.nbComponents)
- synchronize(vd.data)
-
- divTOp = diffop.DivStressTensor3D(topo)
- result = [npw.zeros(memshape) for i in xrange(3)]
- result = divTOp(vd.data, result, ind)
-
- # Numerical VS analytical
- Lx = box.length[0]
- errX = (Lx / (nb - 1)) ** 2
-
- for i in xrange(3):
- assert np.allclose(refd[i], result[i][ind], rtol=errX)
-
-if __name__ == "__main__":
- testDivWV()
- testGradVxW()
- testCurl()
- testDivStressTensor()
diff --git a/HySoP/hysop/operator/__init__.py b/HySoP/hysop/operator/__init__.py
index 09b8f2632344d26a2e52c155fcee2aab99c523e0..8cebe57a80a8ba046ec546043d544d3a651adec8 100644
--- a/HySoP/hysop/operator/__init__.py
+++ b/HySoP/hysop/operator/__init__.py
@@ -24,3 +24,4 @@
# - method["Support"] = GPU, CPU
# - method["Formulation"] = Conservative, GradUW (may depends on the operator)
# - ...
+
diff --git a/HySoP/hysop/operator/adapt_timestep.py b/HySoP/hysop/operator/adapt_timestep.py
index 6cb98714501470e5c259ae10efb572e2b79a48cb..70662d439d2d4f9830c97a33ff96ac5dd463d4ac 100755
--- a/HySoP/hysop/operator/adapt_timestep.py
+++ b/HySoP/hysop/operator/adapt_timestep.py
@@ -6,26 +6,27 @@ Definition of the adaptative time step according to the flow fields.
"""
from parmepy.constants import debug
-from parmepy.methods_keys import TimeIntegrator, SpaceDiscretisation, dtAdvecCrit
+from parmepy.methods_keys import TimeIntegrator, SpaceDiscretisation,\
+ dtAdvecCrit
from parmepy.numerics.integrators.runge_kutta3 import RK3
from parmepy.numerics.finite_differences import FD_C_4
from parmepy.operator.discrete.adapt_timestep import AdaptTimeStep_D
from parmepy.operator.continuous import Operator
+from parmepy.variable_parameter import VariableParameter
import numpy as np
class AdaptTimeStep(Operator):
"""
- The adaptative Time Step is computed according to the following expression :
+ The adaptative Time Step is computed according
+ to the following expression :
dt_adapt = min (dt_advection, dt_stretching, dt_cfl)
"""
@debug
- def __init__(self, velocity, vorticity, resolutions,
- dt_adapt=None,
- method=None,
- topo=None, ghosts=None, prefix='./dt.dat',
- task_id=None, comm=None, **other_config):
+ def __init__(self, velocity, vorticity, resolutions, dt_adapt,
+ method=None, topo=None, ghosts=None, filename='./dt.dat',
+ task_id=None, comm=None, time_range=None):
"""
Create a timeStep-evaluation operator from given
velocity and vorticity variables.
@@ -46,6 +47,11 @@ class AdaptTimeStep(Operator):
@param topo : a predefined topology to discretize velocity/vorticity
@param ghosts : number of ghosts points. Default depends on the method.
Automatically computed if not set.
+ @param time_range : [start, end] use to define a 'window' in which
+ the current operator is applied. Outside start-end, this operator
+ has no effect. Start/end are iteration numbers.
+ Default = [2, endofsimu]
+
"""
if method is None:
method = {TimeIntegrator: RK3,
@@ -53,7 +59,6 @@ class AdaptTimeStep(Operator):
dtAdvecCrit: 'vort'}
Operator.__init__(self, [velocity, vorticity], method, topo=topo,
ghosts=ghosts, task_id=task_id, comm=comm)
- self.config = other_config
## velocity variable (vector)
self.velocity = velocity
## vorticity variable (vector)
@@ -62,9 +67,12 @@ class AdaptTimeStep(Operator):
self.resolutions = resolutions
## adaptative time step variable ("variable" object)
self.dt_adapt = dt_adapt
+ assert isinstance(self.dt_adapt, VariableParameter)
+ # Check if 'dt' key is present in dt_adapt dict
+ assert 'dt' in self.dt_adapt.data
## Numerical methods for time and space discretization
self.method = method
- self.prefix = prefix
+ self.filename = filename
assert SpaceDiscretisation in self.method.keys()
assert TimeIntegrator in self.method.keys()
if not dtAdvecCrit in self.method.keys():
@@ -72,6 +80,7 @@ class AdaptTimeStep(Operator):
self.input = [self.velocity, self.vorticity]
self.output = [self.vorticity]
+ self.time_range = time_range
def discretize(self):
if self.method[SpaceDiscretisation] is FD_C_4:
@@ -105,9 +114,8 @@ class AdaptTimeStep(Operator):
self.discreteOperator =\
AdaptTimeStep_D(self.discreteFields[self.velocity],
self.discreteFields[self.vorticity],
- self.dt_adapt,
- method=self.method, prefix=self.prefix,
- **self.config)
+ self.dt_adapt, method=self.method,
+ filename=self.filename, time_range=self.time_range)
self.discreteOperator.setUp()
self._isUpToDate = True
diff --git a/HySoP/hysop/operator/analytic.py b/HySoP/hysop/operator/analytic.py
index a71394ea924b92e70d5e6eb87abf07f6855dc502..78a176c42ea9eae4db6cfe8b2c958b1df0aaf8fb 100644
--- a/HySoP/hysop/operator/analytic.py
+++ b/HySoP/hysop/operator/analytic.py
@@ -77,13 +77,6 @@ class Analytic(Operator):
for v in self.variables:
v.initialize(simulation.time)
- @debug
- def finalize(self):
- """
- Memory cleaning.
- """
- pass
-
def __str__(self):
s = super(Analytic, self).__str__()
if self.discreteOperator is None:
diff --git a/HySoP/hysop/operator/continuous.py b/HySoP/hysop/operator/continuous.py
index 641629c7f91cb6209b2cf73663df083204937c0d..7a8c0b5467c3196806704620259d0d0544d1e7f9 100644
--- a/HySoP/hysop/operator/continuous.py
+++ b/HySoP/hysop/operator/continuous.py
@@ -99,6 +99,8 @@ class Operator(object):
self.timer = Timer(self)
## Redistribute operator list that we must wait for.
self.requirements = []
+ ## time monitoring
+ self.time_info = None
@staticmethod
def getWorkLengths():
@@ -133,7 +135,7 @@ class Operator(object):
@abstractmethod
def setUp(self):
"""
- Last step of initializaton. After this, the operator must be
+ Last step of initialization. After this, the operator must be
ready for apply call.
Main step : setup for discrete operators.
@@ -144,8 +146,9 @@ class Operator(object):
"""
Memory cleaning.
"""
- self.discreteOperator.finalize()
- self.timer = self.timer + self.discreteOperator.timer
+ if self.discreteOperator is not None:
+ self.discreteOperator.finalize()
+ self.timer = self.timer + self.discreteOperator.timer
@debug
def apply(self, simulation=None):
@@ -169,6 +172,7 @@ class Operator(object):
shortName = str(self.__class__).rpartition('.')[-1][0:-2]
s = '[' + str(main_rank) + '] ' + shortName
s += " : operator not discretized --> no computation, time = 0."
+ print s
def isUp(self):
"""
diff --git a/HySoP/hysop/operator/discrete/adapt_timestep.py b/HySoP/hysop/operator/discrete/adapt_timestep.py
index e49dfe8355838653e9c81d05379d938ea340e874..c36d5d26606fd541ab36f1617b96adbd7bd66344 100755
--- a/HySoP/hysop/operator/discrete/adapt_timestep.py
+++ b/HySoP/hysop/operator/discrete/adapt_timestep.py
@@ -6,8 +6,9 @@ Evaluation of the adaptative time step according to the flow fields.
"""
from parmepy.constants import debug
-from parmepy.methods_keys import TimeIntegrator, SpaceDiscretisation, dtAdvecCrit
-from parmepy.numerics.finite_differences import FD_C_4, FD_C_2
+from parmepy.methods_keys import TimeIntegrator, SpaceDiscretisation,\
+ dtAdvecCrit
+from parmepy.numerics.finite_differences import FD_C_4
from parmepy.operator.discrete.discrete import DiscreteOperator
from parmepy.numerics.integrators.euler import Euler
from parmepy.numerics.integrators.runge_kutta2 import RK2
@@ -17,40 +18,48 @@ from parmepy.tools.timers import timed_function
from parmepy.numerics.differential_operations import GradV
import parmepy.tools.numpywrappers as npw
from parmepy.numerics.updateGhosts import UpdateGhosts
-from parmepy.mpi import main_rank
from parmepy.mpi import MPI
from parmepy.constants import np, PARMES_MPI_REAL
+import scitools.filetable as ft
import math
+import os
ceil = math.ceil
class AdaptTimeStep_D(DiscreteOperator):
"""
- The adaptative Time Step is computed according to the following expression :
+ The adaptative Time Step is computed according
+ to the following expression :
dt_adapt = min (dt_advection, dt_stretching, dt_cfl)
"""
@debug
- def __init__(self, velocity, vorticity,
- dt_adapt=None,
- method={TimeIntegrator: RK3, SpaceDiscretisation: FD_C_4,
- dtAdvecCrit: 'vort'},
- lcfl=0.125, cfl=0.5, prefix='./dt.dat'):
+ def __init__(self, velocity, vorticity, dt_adapt, method=None,
+ lcfl=0.125, cfl=0.5, filename='./dt.dat', time_range=None):
"""
@param velocity : discrete field
@param vorticity : discrete field
- @param dt_adapt : adaptative timestep variable
+ @param dt_adapt : adaptative timestep
+ (a parmepy.variable_parameter.VariableParameter)
@param method : numerical method for space/time discretizations
- @param lcfl : the lagrangian CFL coefficient used for advection stability
+ @param lcfl : the lagrangian CFL coefficient used
+ for advection stability
@param cfl : the CFL coefficient.
+ @param filename : output file name
+ @param time_range : [start, end] use to define a 'window' in which
+ the current operator is applied. Outside start-end, this operator
+ has no effect. Start/end are iteration numbers.
+ Default = [2, endofsimu]
"""
## velocity discrete field
self.velocity = velocity
## vorticity discrete field
self.vorticity = vorticity
- ## adaptative time step variable ("variable" object)
+ ## adaptative time step variable
self.dt_adapt = dt_adapt
-
+ if method is None:
+ method = {TimeIntegrator: RK3, SpaceDiscretisation: FD_C_4,
+ dtAdvecCrit: 'vort'}
DiscreteOperator.__init__(self, [self.velocity, self.vorticity],
method=method)
@@ -76,10 +85,17 @@ class AdaptTimeStep_D(DiscreteOperator):
# Definition of criterion for dt_advec computation
self.dtAdvecCrit = self.method[dtAdvecCrit]
-
- self.prefix = prefix
- if (main_rank == 0):
- self.f = open(self.prefix, 'w')
+ ## Output file name
+ self.filename = filename
+ ## Time range
+ if time_range is None:
+ time_range = [2, np.infty]
+ self.time_range = time_range
+
+ # local buffer :
+ # [time, dt, d1, d2, d3, d4, d5]
+ # for d1...d5 see computation details in apply.
+ self.diagnostics = npw.zeros((7))
def setUp(self):
@@ -87,26 +103,32 @@ class AdaptTimeStep_D(DiscreteOperator):
self._synchronize = UpdateGhosts(self.velocity.topology,
self.velocity.nbComponents)
+ if self.velocity.topology.rank == 0:
+ # Create output dir if required
+ d = os.path.dirname(self.filename)
+ if not os.path.exists(d):
+ os.makedirs(d)
+ self._file = open(self.filename, 'w')
# gradU function
self._function = GradV(self.velocity.topology,
method=self.method[SpaceDiscretisation])
memshape = self.velocity.data[0].shape
worklength = self.velocity.nbComponents ** 2
- # gradU result array
- self.grad = [npw.zeros(memshape)
- for i in xrange(worklength)]
+ # gradU result array.
+ self.grad = [npw.zeros(memshape) for i in xrange(worklength)]
self._isUpToDate = True
@timed_function
- def apply(self, simulation):
+ def apply(self, simulation=None):
if simulation is None:
raise ValueError("Missing simulation value for computation.")
# current time
time = simulation.time
iteration = simulation.currentIteration
-
- if iteration >= 2:
+ Nmax = min(simulation.iterMax, self.time_range[1])
+ self.diagnostics[0] = time
+ if iteration >= self.time_range[0] and iteration <= Nmax:
# Calling for requirements completion
DiscreteOperator.apply(self, simulation)
@@ -114,25 +136,30 @@ class AdaptTimeStep_D(DiscreteOperator):
self._synchronize(self.velocity.data)
# gradU computation
self.grad = self._function(self.velocity.data, self.grad)
-
- diagnostics = npw.zeros((5))
+ self.diagnostics[2:] = 0.0
nbComponents = self.velocity.nbComponents
for direc in xrange(nbComponents):
- # maxima of partial derivatives of velocity : needed for advection stability condition (1st option)
- diagnostics[0] = max(diagnostics[0],
- np.max(abs(self.grad[(nbComponents + 1) * direc])))
-
- # maxima of partial derivatives of velocity: needed for stretching stability condition
- diagnostics[1] = max(diagnostics[1],
- np.max(sum([abs(self.grad[i])
- for i in xrange(nbComponents * direc,
- nbComponents * (direc +1 ))])))
+ # maxima of partial derivatives of velocity :
+ # needed for advection stability condition (1st option)
+ tmp = np.max(abs(self.grad[(nbComponents + 1) * direc]))
+ self.diagnostics[2] = max(self.diagnostics[2], tmp)
+
+ # maxima of partial derivatives of velocity:
+ # needed for stretching stability condition
+ tmp = np.max(sum([abs(self.grad[i])
+ for i in xrange(nbComponents * direc,
+ nbComponents * (direc + 1))
+ ]))
+ self.diagnostics[3] = max(self.diagnostics[3], tmp)
+
# maxima of velocity : needed for CFL based time step
- diagnostics[2] = max(diagnostics[2],
- np.max(abs(self.velocity.data[direc])))
- # maxima of vorticity : needed for advection stability condition (2nd option)
- diagnostics[3] = max(diagnostics[3],
- np.max(abs(self.vorticity.data[direc])))
+ tmp = np.max(abs(self.velocity.data[direc]))
+ self.diagnostics[4] = max(self.diagnostics[4], tmp)
+
+ # maxima of vorticity :
+ # needed for advection stability condition (2nd option)
+ tmp = np.max(abs(self.vorticity.data[direc]))
+ self.diagnostics[5] = max(self.diagnostics[5], tmp)
# 1/2(gradU + gradU^T) computation
self.grad[1] += self.grad[3]
@@ -144,26 +171,28 @@ class AdaptTimeStep_D(DiscreteOperator):
self.grad[3][...] = self.grad[1][...]
self.grad[6][...] = self.grad[2][...]
self.grad[7][...] = self.grad[5][...]
+ # maxima of deformation tensor:
+ # needed for advection stability condition (3rd option)
for direc in xrange(nbComponents):
- # maxima of deformation tensor: needed for advection stability condition (3rd option)
- diagnostics[4] = max(diagnostics[4],
- np.max(sum([abs(self.grad[i])
- for i in xrange(nbComponents * direc,
- nbComponents * (direc +1 ))])))
+ tmp = np.max(sum([abs(self.grad[i])
+ for i in xrange(nbComponents * direc,
+ nbComponents * (direc + 1))
+ ]))
+ self.diagnostics[6] = max(self.diagnostics[6], tmp)
# definition of dt_advection
if self.dtAdvecCrit == 'gradU':
# => based on gradU
- dt_advec = self.lcfl / diagnostics[0]
+ dt_advec = self.lcfl / self.diagnostics[2]
elif self.dtAdvecCrit == 'deform':
# => based on the deformations
- dt_advec = self.lcfl / diagnostics[4]
+ dt_advec = self.lcfl / self.diagnostics[6]
else:
# => based on the vorticity
- dt_advec = self.lcfl / diagnostics[3]
+ dt_advec = self.lcfl / self.diagnostics[5]
# definition of dt_stretching
- dt_stretch = self.coef_stretch / diagnostics[1]
+ dt_stretch = self.coef_stretch / self.diagnostics[3]
# round dt_advec
# pow_adv = 0
@@ -183,24 +212,24 @@ class AdaptTimeStep_D(DiscreteOperator):
# redefine timestep
if self.cfl is None:
- self.dt_adapt.data[0] = min(dt_advec, dt_stretch)
- else :
+ self.dt_adapt['dt'] = min(dt_advec, dt_stretch)
+ else:
# definition of dt_cfl
- dt_cfl = (self.cfl * self.velocity.topology.mesh.space_step[0]) / (diagnostics[2])
- self.dt_adapt.data[0] = min(dt_advec, dt_stretch, dt_cfl)
+ h = self.velocity.topology.mesh.space_step[0]
+ dt_cfl = (self.cfl * h) / (self.diagnostics[4])
+ self.dt_adapt['dt'] = min(dt_advec, dt_stretch, dt_cfl)
- self.dt_adapt.data[0] = \
+ self.dt_adapt['dt'] = \
self.velocity.topology.comm.allreduce(self.dt_adapt.data[0],
- PARMES_MPI_REAL, op=MPI.MIN)
-
-
- if ((main_rank == 0)):
- self.f = open(self.prefix, 'a')
- self.f.write("%s %s %s %s %s %s %s\n" % (time,
- self.dt_adapt.data[0],
- diagnostics[0],
- diagnostics[1],
- diagnostics[2],
- diagnostics[3],
- diagnostics[4]))
- self.f.close()
+ PARMES_MPI_REAL,
+ op=MPI.MIN)
+ self.diagnostics[1] = self.dt_adapt['dt']
+ if self.velocity.topology.rank == 0:
+ self.writeToFile()
+ # Update simulation time step with the new dt
+ simulation.updateTimeStep(self.dt_adapt['dt'])
+
+ def writeToFile(self):
+ self._file = open(self.filename, 'a')
+ ft.write(self._file, self.diagnostics)
+ self._file.close()
diff --git a/HySoP/hysop/operator/discrete/discrete.py b/HySoP/hysop/operator/discrete/discrete.py
index 02d54ba6f81f7f6098385ce2e374c150110707b2..a752db902e728422419f6ccb33f1cb1fe4b6e96b 100644
--- a/HySoP/hysop/operator/discrete/discrete.py
+++ b/HySoP/hysop/operator/discrete/discrete.py
@@ -101,6 +101,7 @@ class DiscreteOperator(object):
return self._isUpToDate
@debug
+ @abstractmethod
def apply(self, simulation=None):
"""
Apply this operator to its variables.
diff --git a/HySoP/hysop/operator/discrete/poisson_fft.py b/HySoP/hysop/operator/discrete/poisson_fft.py
index ab2df6eebf2efc6bdaa004da09cd512e10de0f7c..9e24ff4f38e846060aa2f2ad2fe38388e50b26c8 100644
--- a/HySoP/hysop/operator/discrete/poisson_fft.py
+++ b/HySoP/hysop/operator/discrete/poisson_fft.py
@@ -5,10 +5,8 @@ Discrete operator for Poisson problem (fftw based)
"""
import parmepy.tools.numpywrappers as npw
from parmepy.f2py import fftw2py
-from discrete import DiscreteOperator
+from parmepy.operator.discrete.discrete import DiscreteOperator
from parmepy.constants import debug, prof
-from parmepy.tools.timers import timed_function
-from parmepy.variables.variables import Variables
from parmepy.mpi import MPI
@@ -19,111 +17,157 @@ class PoissonFFT(DiscreteOperator):
"""
@debug
- def __init__(self, velocity, vorticity, method=None, projection=None,
- multires=False, filterSize=None):
+ def __init__(self, velocity, vorticity, projection=None,
+ multires=False, filterSize=None, correction=None):
"""
Constructor.
@param[out] velocity : discretisation of the solution field
@param[in] vorticity : discretisation of the RHS (mind the minus rhs!)
+ @param projection :
+ @param[in] multires : true if velo/vorticity do not have the same
+ resolution
+ @param filterSize :
+ @param correction : operator used to shift velocity according
+ to a given input (fixed) flowrate.
+ See parmepy.operator.velocity_correction.
+ Default = None.
"""
## Solution field
self.velocity = velocity
## RHS field
self.vorticity = vorticity
## Solenoidal projection of vorticity ?
- if (isinstance(projection, Variables)):
- self.projection = projection.data
- else:
- self.projection = projection
- ## Multiresolution ?
- self.multires = multires
+ self.projection = projection
## Filter size array = domainLength/(CoarseRes-1)
self.filterSize = filterSize
# Base class initialisation
- DiscreteOperator.__init__(self, [velocity, vorticity], method)
+ DiscreteOperator.__init__(self, [velocity, vorticity])
# If 2D problem, vorticity must be a scalar
- if self.velocity.dimension == 2:
+ self.dim = self.velocity.domain.dimension
+ if self.dim == 2:
assert self.vorticity.nbComponents == 1
- ## fftw solver case (2 or 3D)
- self.case = '2D'
- elif self.velocity.dimension == 3:
- self.case = '3D'
- else:
+
+ if self.dim != 2 and self.dim != 3:
raise AttributeError("Wrong problem dimension: only 2D \
and 3D cases are implemented.")
self.input = [self.vorticity]
self.output = [self.velocity]
- # If there is a projection, vorticity is also an output
- if self.projection is not None and self.projection[0]:
- self.output.append(self.vorticity)
+ ## Multiresolution ?
+ self.multires = multires
- @debug
- @prof
- def apply(self, simulation):
- # Calling for requirements completion
- DiscreteOperator.apply(self, simulation)
- ctime = MPI.Wtime()
+ # connexion to the required apply function
+ # (to avoid 'if' calls during apply)
+ if self.dim == 2:
+ self._solve = self._solve2D
+ elif self.dim == 3:
+ # If there is a projection, vorticity is also an output
+ if self.projection is not None:
+ self.output.append(self.vorticity)
+ if self.multires:
+ self._solve = self._solve3D_proj_multires
+ else:
+ self._solve = self._solve3D_proj
+ else:
+ if self.multires:
+ self._solve = self._solve3D_multires
+ else:
+ self._solve = self._solve3D
+
+ # Operator to shift velocity according to an input required flowrate
+ if correction is not None:
+ self.correctionOp = correction
+ self.solve = self._solve_and_correct
+ else:
+ self.solve = self._solve
+
+ def _solve2D(self, simu=None):
+ """
+ Solve 2D poisson problem
+ """
+ self.velocity.data[0], self.velocity.data[1] =\
+ fftw2py.solve_poisson_2d(self.vorticity.data,
+ self.velocity.data[0],
+ self.velocity.data[1])
+
+ def _project(self):
+ """
+ apply projection onto vorticity
+ """
+ ghosts_w = self.vorticity.topology.ghosts
+ self.vorticity.data[0], self.vorticity.data[1], \
+ self.vorticity.data[2] = \
+ fftw2py.projection_om_3d(self.vorticity.data[0],
+ self.vorticity.data[1],
+ self.vorticity.data[2], ghosts_w)
- ite = simulation.currentIteration
+ def _solve3D_multires(self, simu=None):
+ """
+ 3D, multiresolution
+ """
+ # Projects vorticity values from fine to coarse grid
+ # in frequencies space by nullifying the smallest modes
+ vortFilter = npw.copy(self.vorticity.data)
+ vortFilter[0], vortFilter[1], vortFilter[2] = \
+ fftw2py.multires_om_3d(self.filterSize[0], self.filterSize[1],
+ self.filterSize[2], self.vorticity.data[0],
+ self.vorticity.data[1],
+ self.vorticity.data[2])
+
+ # Solves Poisson equation using filter vorticity
ghosts_v = self.velocity.topology.ghosts
ghosts_w = self.vorticity.topology.ghosts
- if self.case is '2D':
- self.velocity.data[0], self.velocity.data[1] =\
- fftw2py.solve_poisson_2d(self.vorticity.data,
- self.velocity.data[0],
- self.velocity.data[1])
-
- elif self.case is '3D':
- #=== SOLENOIDAL PROJECTION OF VORTICITY FIELD ===
- if (self.projection is not None and self.projection[0] \
- and ite % self.projection[1] == 0):
- self.vorticity.data[0], self.vorticity.data[1], \
- self.vorticity.data[2] = \
- fftw2py.projection_om_3d(self.vorticity.data[0],
- self.vorticity.data[1],
- self.vorticity.data[2],
- ghosts_w)
-
- if (self.multires):
- # Projects vorticity values from fine to coarse grid
- # in frequencies space by nullifying the smallest modes
- vortFilter = npw.copy(self.vorticity.data)
- vortFilter[0], vortFilter[1], \
- vortFilter[2] = \
- fftw2py.multires_om_3d(self.filterSize[0],
- self.filterSize[1],
- self.filterSize[2],
- self.vorticity.data[0],
- self.vorticity.data[1],
- self.vorticity.data[2])
-
- # Solves Poisson equation using filter vorticity
- self.velocity.data[0], self.velocity.data[1], \
- self.velocity.data[2] = \
- fftw2py.solve_poisson_3d(vortFilter[0],
- vortFilter[1],
- vortFilter[2],
- self.velocity.data[0],
- self.velocity.data[1],
- self.velocity.data[2],
- ghosts_w, ghosts_v)
+ self.velocity.data[0], self.velocity.data[1], self.velocity.data[2] = \
+ fftw2py.solve_poisson_3d(vortFilter[0], vortFilter[1],
+ vortFilter[2], self.velocity.data[0],
+ self.velocity.data[1],
+ self.velocity.data[2], ghosts_w, ghosts_v)
- else:
- # Solves Poisson equation using usual vorticity
- self.velocity.data[0], self.velocity.data[1], \
- self.velocity.data[2] = \
- fftw2py.solve_poisson_3d(self.vorticity.data[0],
- self.vorticity.data[1],
- self.vorticity.data[2],
- self.velocity.data[0],
- self.velocity.data[1],
- self.velocity.data[2],
- ghosts_w, ghosts_v)
+ def _solve3D_proj_multires(self, simu):
+ """
+ 3D, multiresolution, with projection
+ """
+ ite = simu.currentIteration
+ if self.projection.doProjection(ite):
+ self._project()
+ self._solve3D_multires()
- else:
- raise ValueError("invalid problem dimension")
+ def _solve3D_proj(self, simu):
+ """
+ 3D, with projection
+ """
+ ite = simu.currentIteration
+ if self.projection.doProjection(ite):
+ self._project()
+ self._solve3D()
+
+ def _solve3D(self,simu=None):
+ """
+ Basic solve
+ """
+ # Solves Poisson equation using usual vorticity
+ ghosts_v = self.velocity.topology.ghosts
+ ghosts_w = self.vorticity.topology.ghosts
+ self.velocity.data[0], self.velocity.data[1], self.velocity.data[2] =\
+ fftw2py.solve_poisson_3d(self.vorticity.data[0],
+ self.vorticity.data[1],
+ self.vorticity.data[2],
+ self.velocity.data[0],
+ self.velocity.data[1],
+ self.velocity.data[2], ghosts_w, ghosts_v)
+
+ def _solve_and_correct(self, simu):
+ self._solve(simu.currentIteration)
+ self.correctionOp.apply(simu)
+
+ @debug
+ @prof
+ def apply(self, simulation=None):
+ # Calling for requirements completion
+ DiscreteOperator.apply(self, simulation)
+ ctime = MPI.Wtime()
+ self.solve(simulation)
self._apply_timer.append_time(MPI.Wtime() - ctime)
def finalize(self):
diff --git a/HySoP/hysop/operator/discrete/velocity_correction.py b/HySoP/hysop/operator/discrete/velocity_correction.py
index f36f88ad470c7eb2afa5449169c1fd95aabba301..e574cf49dfa514389511d628d9b380db5b728cd5 100755
--- a/HySoP/hysop/operator/discrete/velocity_correction.py
+++ b/HySoP/hysop/operator/discrete/velocity_correction.py
@@ -15,116 +15,128 @@ from parmepy.mpi import MPI
class VelocityCorrection_D(DiscreteOperator):
"""
- The velocity field is corrected after solving the
- Poisson equation. For more details about calculations,
- see the "velocity_correction.pdf" explanations document
+ The velocity field is corrected after solving the
+ Poisson equation. For more details about calculations,
+ see the "velocity_correction.pdf" explanations document
in ParmeDoc directory.
"""
@debug
- def __init__(self, velocity, vorticity, uinf):
+ def __init__(self, velocity, vorticity, req_flowrate, cb=None):
"""
- @param velocity : discrete field
+ @param[in, out] velocity field to be corrected
+ @param[in] vorticity field used to compute correction
+ @param[in] req_flowrate : required value for the flowrate
+ @param[in] surf : surface (parmepy.domain.obstacle.planes.SubPlane)
+ used to compute reference flow rates. Default = surface at x_origin,
+ normal to x-dir.
"""
## velocity discrete field
self.velocity = velocity
## vorticity discrete field
self.vorticity = vorticity
- ## free stream velocity
- self.uinf = uinf
-
+ ## domain dimension
+ self.dim = self.velocity.domain.dimension
# If 2D problem, vorticity must be a scalar
- if self.velocity.dimension == 2:
+ if self.dim == 2:
assert self.vorticity.nbComponents == 1
- self.case = '2D'
- elif self.velocity.dimension == 3:
- self.case = '3D'
- else:
- raise AttributeError("Wrong problem dimension: only 2D \
- and 3D cases are implemented.")
+ assert (self.dim == 2 or self.dim == 3),\
+ "Wrong problem dimension: only 2D and 3D cases are implemented."
DiscreteOperator.__init__(self, [self.velocity, self.vorticity])
self.input = [self.velocity, self.vorticity]
self.output = [self.velocity]
+ ## Velocity from topology will be used as reference
+ self.topo = self.velocity.topology
+ hh = self.topo.mesh.space_step
+ self.cb = cb
+ self.cb.discretize(self.topo)
+ self.surfRef = cb.lowerS[XDIR]
+ normal = self.surfRef.normal
+ dirs = np.where(normal == 0)[0]
+ self.ds = np.prod(self.surfRef.lengths[dirs] + hh[dirs])
+ ## Expected value for the flow rate through self.surfRef
+ ## warning : divided with area of input surf.
+ self.req_flowrate = req_flowrate / self.ds
+ ## The correction that must be applied on each
+ ## component of the velocity.
+ self.velocity_shift = npw.zeros(self.dim)
+ nbf = self.velocity.nbComponents + self.vorticity.nbComponents
+ # temp buffer, used to save flow rates and mean
+ # values of vorticity
+ self.rates = npw.zeros(nbf)
def setUp(self):
- spaceStep = self.velocity.topology.mesh.space_step
- lengths = self.velocity.topology.domain.length
- self.coeff_mean = np.prod(spaceStep) / np.prod(lengths)
- self.x0_coord = self.velocity.topology.domain.origin[XDIR]
- self.x_coord = self.velocity.topology.mesh.coords[XDIR]
- if self.case is '2D':
- self.surf = lengths[YDIR]
- self.dvol = spaceStep[YDIR]
- elif self.case is '3D':
- self.surf = lengths[YDIR] * lengths[ZDIR]
- self.dvol = spaceStep[YDIR] * spaceStep[ZDIR]
- else:
- raise ValueError("invalid problem dimension")
- self.comput_flowrates = npw.zeros((3))
-
+ spaceStep = self.topo.mesh.space_step
+ lengths = self.topo.domain.length
+ self.coeff_mean = npw.prod(spaceStep) / npw.prod(lengths)
+ x0 = self.topo.domain.origin[XDIR]
+ self.x_coord = self.topo.mesh.coords[XDIR] - x0
self._isUpToDate = True
- @timed_function
- def apply(self, simulation):
- ## Calling for requirements completion
- DiscreteOperator.apply(self, simulation)
- ctime = MPI.Wtime()
-
- ## Computation of the flowrates evaluated from
+ def computeCorrection(self):
+ """
+ Compute the required correction for the current state
+ but do not apply it onto velocity.
+ """
+ ## Computation of the flowrates evaluated from
## current (ie non corrected) velocity
- local_comput_flowrates = npw.zeros((3))
+ ## local flow reduced only on proc 0
+ nbf = self.velocity.nbComponents + self.vorticity.nbComponents
+ localrates = npw.zeros((nbf))
for i in xrange(self.velocity.nbComponents):
- local_comput_flowrates[i] = np.sum(self.velocity[i][...])
-
- recvbuff = npw.zeros((3))
- self.velocity.topology.topo.Allreduce(local_comput_flowrates,
- recvbuff)
+ localrates[i] = self.velocity.integrateOnSurf_proc(self.surfRef,
+ component=i)
+ start = self.velocity.nbComponents
+ ## Integrate vorticity over the whole domain
+ for i in xrange(self.vorticity.nbComponents):
+ localrates[start + i] = \
+ self.vorticity.integrate_on_proc(self.cb, component=i)
+
+ # MPI reduction for rates
+ # rates = [flowrate[X], flowrate[Y], flowrate[Z],
+ # vort_mean[X], ..., vort_mean[Z]]
+ # or (in 2D) = [flowrate[X], flowrate[Y], vort_mean]
+ self.velocity.topology.comm.Allreduce(localrates, self.rates)
+
+ # Set velocity_shift == [Vx_shift, vort_mean[Y], vort_mean[Z]]
+ # or (in 2D) velocity_shift == [Vx_shift, vort_mean]
+ # Velocity shift for main dir component
+ self.velocity_shift[XDIR] = self.req_flowrate\
+ - self.rates[XDIR] / self.ds
+ # Shifts in other directions depend on x coord
+ # and will be computed during apply.
- for i in xrange(self.velocity.nbComponents):
- self.comput_flowrates[i] = recvbuff[i] * self.dvol
-
- ## Correction of the X-velocity component
- self.velocity[XDIR][...] += self.uinf - \
- self.comput_flowrates[XDIR] / self.surf
-
- ## Correction of the other velocity components
- if self.case is '2D':
- # Computation of vorticity mean (in space)
- local_vort_mean = np.sum(self.vorticity[0][...])
- recvbuff = 0.0
- self.velocity.topology.topo.Allreduce(local_vort_mean,
- recvbuff)
- vort_mean = recvbuff * self.coeff_mean
+ @timed_function
+ def apply(self, simulation=None):
+ # Calling for requirements completion
+ DiscreteOperator.apply(self, simulation)
+ ctime = MPI.Wtime()
+ # Computation of the required velocity shift
+ # for the current state
+ self.computeCorrection()
+ iCompute = self.topo.mesh.iCompute
+
+ # Apply shift to velocity
+ self.velocity[XDIR][iCompute] += self.velocity_shift[XDIR]
+ start = self.velocity.nbComponents
+ # reminder : self.rates =[vx_shift, flowrates[Y], flowrate[Z],
+ # vort_mean[X], vort_mean[Y], vort_mean[Z]]
+ # or (in 2D) [vx_shift, flowrates[Y], vort_mean]
+ vort_mean = self.rates[start:]
+ if self.dim == 2:
# Correction of the Y-velocity component
- self.velocity[YDIR][...] += vort_mean * \
- (self.x_coord - self.x0_coord) - \
- self.comput_flowrates[YDIR] \
- / self.surf
-
- elif self.case is '3D':
- # Computation of Y and Z-vorticity means (in space)
- local_vort_Y_mean = np.sum(self.vorticity[YDIR][...])
- local_vort_Z_mean = np.sum(self.vorticity[ZDIR][...])
-
- sendbuff = npw.zeros((2))
- recvbuff = npw.zeros((2))
- sendbuff[:] = [local_vort_Y_mean, local_vort_Z_mean]
- self.velocity.topology.topo.Allreduce(sendbuff, recvbuff)
- vort_mean = recvbuff * self.coeff_mean
+ self.velocity[YDIR][...] += vort_mean[XDIR] * self.x_coord \
+ - self.rates[YDIR] / self.ds
+ elif self.dim == 3:
# Correction of the Y and Z-velocity components
- self.velocity[YDIR][...] += vort_mean[1] * \
- (self.x_coord - self.x0_coord) - \
- self.comput_flowrates[YDIR] \
- / self.surf
- self.velocity[ZDIR][...] += vort_mean[0] * \
- (self.x_coord - self.x0_coord) - \
- self.comput_flowrates[ZDIR] \
- / self.surf
-
- else:
- raise ValueError("invalid problem dimension")
+ self.velocity[YDIR][...] += vort_mean[ZDIR] * self.x_coord - \
+ self.rates[YDIR] / self.ds
+ self.velocity[ZDIR][...] += -vort_mean[YDIR] * self.x_coord - \
+ self.rates[ZDIR] / self.ds
+
self._apply_timer.append_time(MPI.Wtime() - ctime)
+
diff --git a/HySoP/hysop/operator/energy_enstrophy.py b/HySoP/hysop/operator/energy_enstrophy.py
index cb68191782c23e9402079176812c2e7d4e4018b1..33644d0df331de7af3a880e5cb8aef46480d891b 100644
--- a/HySoP/hysop/operator/energy_enstrophy.py
+++ b/HySoP/hysop/operator/energy_enstrophy.py
@@ -8,8 +8,6 @@ from parmepy.constants import debug, XDIR
from parmepy.operator.monitors.monitoring import Monitoring
from parmepy.tools.timers import timed_function
import parmepy.tools.numpywrappers as npw
-import scitools.filetable as ft
-import os
class Energy_enstrophy(Monitoring):
@@ -25,8 +23,8 @@ class Energy_enstrophy(Monitoring):
"""
def __init__(self, velocity, vorticity,
- viscosity, isNormalized, topo, frequency, prefix=None,
- filename=None, safeOutput=True, task_id=None):
+ viscosity, isNormalized, topo,
+ io_params=None, task_id=None):
"""
Constructor.
@param velocity field
@@ -35,16 +33,20 @@ class Energy_enstrophy(Monitoring):
@param isNormalized : boolean indicating whether the enstrophy
and energy values have to be normalized by the domain lengths.
@param topo : the topology on which we want to monitor the fields
- @param frequency : output file producing frequency.
- @param filename : output file name. Default is None ==> no output.
- Full or relative path.
- @param safeOutput : boolean defining how often output file is written:
- True --> open/close file everytime data are written
- False --> open at init and close during finalize (cost less but if simu
- crashes, data are lost.)
+ @param io_params : parameters (dict) to set file output.
+ If None, no output. Set io_params = {} if you want output,
+ with default parameters values. Default file name = 'energy_enstrophy'
+ See parmepy.tools.io_utils.Writer for details
"""
- Monitoring.__init__(self, [velocity, vorticity], topo, frequency,
- task_id=task_id)
+ if io_params is not None:
+ if not "filename" in io_params:
+ io_params["filename"] = "energy_enstrophy"
+ # Set output buffer shape
+ io_params["writebuffshape"] = (1, 3)
+
+ io_params["writebuffshape"] = (1, 3)
+ Monitoring.__init__(self, [velocity, vorticity], topo,
+ io_params, task_id=task_id)
## velocity field
self.velocity = velocity
## vorticity field
@@ -63,70 +65,28 @@ class Energy_enstrophy(Monitoring):
# Note FP : for multiresolution case, it would probably be
# better to use two different operators for energy and enstrophy.
- # Output file
- if filename is None:
- self._writeOuput = False
- self.filename = None
- else:
- self._writeOuput = True
- self.filename = filename
- if (self._topo.rank == 0):
- # Create output dir if required
- d = os.path.dirname(self.filename)
- if not os.path.exists(d):
- os.makedirs(d)
- self.safeOutput = safeOutput
- if self.safeOutput:
- self._writeFile = self._fullwrite
- else:
- self._writeFile = self._partialwrite
-
- self._file = open(self.filename, 'w')
-
- def _fullwrite(self, data):
- self._file = open(self.filename, 'a')
- self._file.write("%s %s %s \n" % (data[0, 0],
- data[0, 1],
- data[0, 2]))
- self._file.close()
-
- def _partialwrite(self, data):
-# ft.write(self._file, data)
- self._file.write("%s %s %s \n" % (data[0, 0],
- data[0, 1],
- data[0, 2]))
-
def setUp(self):
- Monitoring.setUp(self)
- # Get discrete fields for velocity and vorticity.
- # Note FP : two options to get the discrete fields:
- # - 'field.discretization(topo)' that just try
- # to get the discrete field and return an error
- # if it does not exist.
- # - 'field.discretizet(topo) try to get the discrete field and create
- # a new discretization if it does not exist.
- # Current choice : no creation.
- self.discreteFields[self.velocity] = \
- self.velocity.discretization(self._topo)
- self.discreteFields[self.vorticity] =\
- self.vorticity.discretization(self._topo)
-
- self._isUpToDate = True
-
- spaceStep = self._topo.mesh.space_step
- length = self._topo.domain.length
- if self.isNormalized :
- self.coeffEnstrophy = (np.prod(spaceStep) /
- np.prod(length))
- else:
- self.coeffEnstrophy = np.prod(spaceStep)
- self.coeffEnergy = 0.5 * self.coeffEnstrophy
+ if not self._isUpToDate:
+ self.discreteFields[self.velocity] = \
+ self.velocity.discretization(self._topo)
+ self.discreteFields[self.vorticity] =\
+ self.vorticity.discretization(self._topo)
+
+ spaceStep = self._topo.mesh.space_step
+ length = self._topo.domain.length
+ if self.isNormalized:
+ self.coeffEnstrophy = (np.prod(spaceStep) /
+ np.prod(length))
+ else:
+ self.coeffEnstrophy = np.prod(spaceStep)
+ self.coeffEnergy = 0.5 * self.coeffEnstrophy
- # A work vector for local computations
- # Warning : we assume one topo for all variables
- self.ind = self._topo.mesh.iCompute
- shape = self.discreteFields[self.velocity].data[0][self.ind].shape
- self._work = npw.zeros(shape)
+ # A work vector for local computations
+ # Warning : we assume one topo for all variables
+ self.ind = self._topo.mesh.iCompute
+ shape = self.discreteFields[self.velocity].data[0][self.ind].shape
+ self._work = npw.zeros(shape)
+ self._isUpToDate = True
@debug
@timed_function
@@ -139,8 +99,6 @@ class Energy_enstrophy(Monitoring):
raise ValueError("Missing simulation value for computation.")
#time = MPI.Wtime()
- t = simulation.time
- ite = simulation.currentIteration
# --- Kinetic energy computation ---
vd = self.discreteFields[self.velocity]
@@ -167,6 +125,12 @@ class Energy_enstrophy(Monitoring):
recvbuff = npw.zeros((2))
sendbuff[:] = [local_energy, local_enstrophy]
self._topo.comm.Allreduce(sendbuff, recvbuff)
+ # the other way :
+ #energy = self._topovel.topo.allreduce(local_energy, PARMES_MPI_REAL,
+ # op=MPI.SUM)
+ #enstrophy = self._topovel.topo.allreduce(local_enstrophy,
+ # PARMES_MPI_REAL,
+ # op=MPI.SUM)
# Update global values
energy = recvbuff[0] * self.coeffEnergy
@@ -179,14 +143,9 @@ class Energy_enstrophy(Monitoring):
## self._buffer_1 = energy
# Print results, if required
- if self._writeOuput and self._topo.rank == 0 \
- and (ite % self.freq) == 0:
- dataout = npw.zeros((1, 3))
- dataout[0, 0] = simulation.time
- dataout[0, 1] = energy
- dataout[0, 2] = enstrophy
- self._writeFile(dataout)
-
- def finalize(self):
- if not self.safeOutput:
- self._file.close()
+ ite = simulation.currentIteration
+ if self._writer is not None and self._writer.doWrite(ite):
+ self._writer.buffer[0, 0] = simulation.time
+ self._writer.buffer[0, 1] = energy
+ self._writer.buffer[0, 2] = enstrophy
+ self._writer.write()
diff --git a/HySoP/hysop/operator/monitors/__init__.py b/HySoP/hysop/operator/monitors/__init__.py
index 905c2429f51593c82146ad38a9e3824e9c53003a..5728a8fe9c9feff86bc390bdcd5eac0b1daacb57 100644
--- a/HySoP/hysop/operator/monitors/__init__.py
+++ b/HySoP/hysop/operator/monitors/__init__.py
@@ -1,3 +1,19 @@
## @package parmepy.operator.monitors
# Parmes tools for data and fields monitoring.
#
+# import and alias so that monitors are
+# available with
+# from parmepy.operators.monitors import Printer, ...
+import printer
+Printer = printer.Printer
+import compute_forces
+DragAndLift = compute_forces.DragAndLift
+import reprojection_criterion
+Reprojection_criterion = reprojection_criterion.Reprojection_criterion
+import energy_enstrophy as energy_enstrophy
+Energy_enstrophy = energy_enstrophy.Energy_enstrophy
+
+# Set list for 'import *'
+__all__ = ['DragAndLift', 'Energy_enstrophy',
+ 'Reprojection_criterion', 'Printer']
+
diff --git a/HySoP/hysop/operator/monitors/compute_forces.py b/HySoP/hysop/operator/monitors/compute_forces.py
index 038600f8c245b0b7222bc9145be824f3b8549b2e..00b0c59f3a25d26b8d3d50c593c24c1fe7940cab 100644
--- a/HySoP/hysop/operator/monitors/compute_forces.py
+++ b/HySoP/hysop/operator/monitors/compute_forces.py
@@ -3,18 +3,12 @@
@file compute_forces.py
Compute forces
"""
-from parmepy.constants import debug, np, ORDER, \
- PARMES_INTEGER, PI, XDIR, YDIR, ZDIR, PARMES_MPI_REAL
+from parmepy.constants import np, XDIR, YDIR, ZDIR
from parmepy.numerics.updateGhosts import UpdateGhosts
-from parmepy.numerics.differential_operations import Laplacian,\
- DivStressTensor3D
+from parmepy.numerics.differential_operations import Laplacian
from parmepy.numerics.finite_differences import FD_C_2
from parmepy.operator.monitors.monitoring import Monitoring
-from parmepy.mpi import MPI
-from parmepy.tools.timers import timed_function
import parmepy.tools.numpywrappers as npw
-import scitools.filetable as ft
-import os
class DragAndLift(Monitoring):
@@ -25,8 +19,8 @@ class DragAndLift(Monitoring):
Integral inside the obstacle is not taken into account.
"""
def __init__(self, velocity, vorticity, nu, coefForce, topo,
- volumeOfControl, obstacles=None, frequency=1,
- filename=None, safeOutput=True, task_id=None):
+ volumeOfControl, obstacles=None, io_params=None,
+ task_id=None):
"""
@param velocity field
@param vorticity field
@@ -36,12 +30,19 @@ class DragAndLift(Monitoring):
(parmepy.domain.obstacle.controlBox.ControlBox object)
@param obstacles a list of parmepy.domain.obstacles inside
the box
- @param frequency : output rate
- @param filename output file name. Default is None ==> no output.
- Full or relative path.
+ @param io_params : parameters (dict) to set file output.
+ If None, no output. Set io_params = {} if you want output,
+ with default parameters values. Default file name = 'drag_and_lift'
+ See parmepy.tools.io_utils.Writer for details
"""
- Monitoring.__init__(self, [velocity, vorticity], topo, frequency,
- task_id=task_id)
+ if io_params is not None:
+ if not "filename" in io_params:
+ io_params["filename"] = "drag_and_lift"
+ # Set output buffer shape
+ io_params["writebuffshape"] = (1, 4)
+
+ Monitoring.__init__(self, [velocity, vorticity], topo,
+ io_params, task_id=task_id)
self.velocity = velocity
self.vorticity = vorticity
self._voc = volumeOfControl
@@ -85,40 +86,6 @@ class DragAndLift(Monitoring):
# Normalizing coefficient for forces
# (based on the physics of the flow)
self.coefForce = coefForce
- # Output file
- if filename is None:
- self._writeOuput = False
- self.filename = None
- else:
- self._writeOuput = True
- self.filename = filename
- if (self.topo.rank == 0):
- # Create output dir if required
- d = os.path.dirname(self.filename)
- if not os.path.exists(d):
- os.makedirs(d)
- ## Defines how often
- ## output file is written :
- ## True --> open/close file everytime
- ## data are written.
- ## False --> open at init and close
- ## during finalize. Cost less but if simu
- ## crashes, data are lost.
- self.safeOutput = safeOutput
- if self.safeOutput:
- self._writeFile = self._fullwrite
- else:
- self._writeFile = self._partialwrite
-
- self._file = open(self.filename, 'w')
-
- def _fullwrite(self, data):
- self._file = open(self.filename, 'a')
- ft.write(self._file, data)
- self._file.close()
-
- def _partialwrite(self, data):
- ft.write(self._file, data)
def _mpi_allsum(self):
"""
@@ -135,22 +102,24 @@ class DragAndLift(Monitoring):
"""
self.force = self.topo.comm.reduce(self.force, root=root)
- def discretize(self):
+ def setUp(self):
"""
"""
- for v in self.variables:
- # the discrete fields
- self.discreteFields[v] = v.discretize(self.topo)
- self.vd = self.discreteFields[self.velocity]
- self.wd = self.discreteFields[self.vorticity]
- self._voc.discretize(self.topo)
- for obs in self.obstacles:
- obs.discretize(self.topo)
- # prepare ghost points synchro for velocity and vorticity used
- # in fd schemes
- self._synchronize = UpdateGhosts(self.topo,
- self.vd.nbComponents
- + self.wd.nbComponents)
+ if not self._isUpToDate:
+ for v in self.variables:
+ # the discrete fields
+ self.discreteFields[v] = v.discretize(self.topo)
+ self.vd = self.discreteFields[self.velocity]
+ self.wd = self.discreteFields[self.vorticity]
+ self._voc.discretize(self.topo)
+ for obs in self.obstacles:
+ obs.discretize(self.topo)
+ # prepare ghost points synchro for velocity and vorticity used
+ # in fd schemes
+ self._synchronize = UpdateGhosts(self.topo,
+ self.vd.nbComponents
+ + self.wd.nbComponents)
+ self._isUpToDate = True
def apply(self, simulation=None):
"""
@@ -190,11 +159,10 @@ class DragAndLift(Monitoring):
self.force *= self.coefForce
# Print results, if required
- if self._writeOuput and self.topo.rank == 0 and (ite % self.freq) == 0:
- dataout = npw.zeros((1, 4))
- dataout[0, 0] = simulation.time
- dataout[0, 1:] = self.force
- self._writeFile(dataout)
+ if self._writer is not None and self._writer.doWrite(ite):
+ self._writer.buffer[0, 0] = simulation.time
+ self._writer.buffer[0, 1:] = self.force
+ self._writer.write()
return self.force
@@ -317,502 +285,3 @@ class DragAndLift(Monitoring):
res *= self._dvol
return res
-
- def finalize(self):
- if not self.safeOutput:
- self._file.close()
-
-class Forces(Monitoring):
- """
- Compute the forces according the Noca s formula
- """
-
- def __init__(self, velocity, vorticity, topo,
- obstacle, boxMin= None, boxMax=None, Reynolds=None,
- uinf=1.0, method=FD_C_2, frequency=None, prefix=None):
- """
- Constructor.
- @param velocity : Continuous velocity field
- @param vorticity : Continuous vorticity field
- @param obstacle : obstacle upon which the forces are exerced
- @param topo : the topology on which we want to monitor the fields
- @param boxMin : Global minimum coordinates of the control volume
- @param boxMax : Global maximum coordinates of the control volume
- @param Reynolds : Reynolds number value
- @param method : numerical method for Differential operators discretizations
- @param frequency : output file producing frequency.
- @param prefix : file name prefix, contains relative path.
- """
- Monitoring.__init__(self, [velocity, vorticity], topo, frequency)
- if prefix is None:
- self.prefix = './res/Noca_sphere.dat'
- else:
- self.prefix = prefix
- self.velocity = velocity
- self.vorticity = vorticity
- self.obstacle = obstacle
- self.boxMin = boxMin
- self.boxMax = boxMax
- if not (boxMin < boxMax):
- raise ValueError("Error, boxMin has to be lower than " +
- "boxMax to define control Noca box")
- self.Re = Reynolds
- self.method = method
- self.prefix = prefix
- self.input = [velocity, vorticity]
- self.output = []
- self.topo = self._predefinedTopo
-
- # Compute coef used to determine drag coefficient :
- # cD = coef.Fx = 2.Fx/rho.uinf**2.D
- self.bufferForce = 0.
- rho = 1.0
- self.coef = 2. / (rho * uinf ** 2 * PI * self.obstacle.radius ** 2)
-
- def setUp(self):
- Monitoring.setUp(self)
-
- # Variables Discretization
- for v in self.variables:
- self.discreteFields[v] = v.discretize(self.topo)
- self.velo = self.discreteFields[self.velocity]
- self.vort = self.discreteFields[self.vorticity]
-
- # Control box definition
- self.isForceComputationNeeded = True
- self.compute_control_box(self.obstacle)
-
- # prepare ghost points synchro for velocity
- self._synchronize = UpdateGhosts(self.topo,
- self.velocity.nbComponents)
-
- # Function for the computation of div(T) where
- # T is the stress tensor: T = Nabla(U) + Nabla(U)t
- self.divT = DivStressTensor3D(self.topo, method=self.method)
- self.fd_scheme = FD_C_2(self.topo.mesh.space_step)
-
- self._isUpToDate = True
-
- if (self.topo.comm.Get_rank() == 0):
- self.f = open(self.prefix, 'w')
-
- def compute_control_box(self, obst) :
- """
- Compute indicator functions for the control box
- (including the obstacle)
-
- @param obstacle : obstacle upon which the forces are exerced
-
- """
- # The control box definition is based on velocity topology
- self.dim = obst.dimension
- self.ghosts = self.velo.topology.ghosts
- self.res = self.velo.topology.mesh.resolution - 2 * self.ghosts
- self.step = self.velo.topology.mesh.space_step
- self.coords = self.velo.topology.mesh.coords
- self.coord_start = self.velo.topology.mesh.origin + \
- (self.ghosts * self.step)
- self.coord_end = self.velo.topology.mesh.end - \
- (self.ghosts * self.step)
- local_start = self.velo.topology.mesh.local_start
- local_end = self.velo.topology.mesh.local_end
-
- # normal vectors definitions
- self.normal = npw.zeros([self.dim * 2, self.dim])
- self.x_minus = 0
- self.x_plus = 1
- self.y_minus = 2
- self.y_plus = 3
- self.z_minus = 4
- self.z_plus = 5
- self.normal[self.x_minus][0] = - 1.0
- self.normal[self.x_plus][0] = 1.0
- self.normal[self.y_minus][1] = - 1.0
- self.normal[self.y_plus][1] = 1.0
- self.normal[self.z_minus][2] = - 1.0
- self.normal[self.z_plus][2] = 1.0
-
- # local indices and coordinates of the control box boundaries
- ind_boundaries = np.zeros([self.dim, 2], dtype=PARMES_INTEGER, order=ORDER)
-# coord_boundaries = npw.zeros([self.dim, 2])
-
- distMin = self.boxMin - self.coord_start
- distMax = self.boxMax - self.coord_end
-
- for i in xrange(self.dim):
- if (distMin[i]>=0. and distMax[i]<=0.):
- # the control volume is included inside the local domain
- ind_boundaries[i][0] = local_start[i] + distMin[i] // self.step[i] # interpolation to lower grid point
- ind_boundaries[i][1] = local_end[i] + distMax[i] // self.step[i] # interpolation to higher grid point
-# coord_boundaries[i][0] = self.coord_start[i] + distMin[i] // self.step[i] * self.step[i]
-# coord_boundaries[i][1] = self.coord_end[i] + distMax[i] // self.step[i] * self.step[i]
-
- elif (distMin[i]>=0. and self.boxMin[i]<= self.coord_end[i] and distMax[i]>0.):
- # only min corner of control volume is included inside the local domain
- ind_boundaries[i][0] = local_start[i] + distMin[i] // self.step[i]
- ind_boundaries[i][1] = local_end[i]
-# coord_boundaries[i][0] = self.coord_start[i] + distMin[i] // self.step[i] * self.step[i]
-# coord_boundaries[i][1] = self.coord_end[i]
-
- elif (distMin[i]<0. and self.boxMax[i]>= self.coord_start[i] and distMax[i]<=0.) :
- # only max corner of control volume is included inside the local domain
- ind_boundaries[i][0] = local_start[i]
- ind_boundaries[i][1] = local_end[i] + distMax[i] // self.step[i]
-# coord_boundaries[i][0] = self.coord_start[i]
-# coord_boundaries[i][1] = self.coord_end[i] + distMax[i] // self.step[i] * self.step[i]
-
- elif (distMin[i]<=0. and distMax[i]>=0.):
- # the local domain is included inside the control volume
- ind_boundaries[i][0] = local_start[i]
- ind_boundaries[i][1] = local_end[i]
-# coord_boundaries[i][0] = self.coord_start[i]
-# coord_boundaries[i][1] = self.coord_end[i]
- else:
- # there is no part of the volume control inside the local domain
- self.isForceComputationNeeded = False
-
-# print 'coord_start, coord_end', self.topo.comm.Get_rank(), self.coord_start, self.coord_end
-# print 'distMin, distMax, step, res', self.topo.comm.Get_rank(), distMin, distMax, self.step, self.res
-# print 'ind_bound_ChiX-', self.topo.comm.Get_rank(), ind_boundaries[0][0], coord_boundaries[0][0], self.coords[0][ind_boundaries[0][0], :, :]
-# print 'ind_bound_ChiX+', self.topo.comm.Get_rank(), ind_boundaries[0][1], coord_boundaries[0][1], self.coords[0][ind_boundaries[0][1], :, :]
-# print 'ind_bound_ChiY-', self.topo.comm.Get_rank(), ind_boundaries[1][0], coord_boundaries[1][0], self.coords[1][:, ind_boundaries[1][0], :]
-# print 'ind_bound_ChiY+', self.topo.comm.Get_rank(), ind_boundaries[1][1], coord_boundaries[1][1], self.coords[1][:, ind_boundaries[1][1], :]
-# print 'ind_bound_ChiZ-', self.topo.comm.Get_rank(), ind_boundaries[2][0], coord_boundaries[2][0], self.coords[2][:, :, ind_boundaries[2][0]]
-# print 'ind_bound_ChiZ+', self.topo.comm.Get_rank(), ind_boundaries[2][1], coord_boundaries[2][1], self.coords[2][:, :, ind_boundaries[2][1]]
-
- # definition of the indicator functions of box boundaries
- self.chi_x_minus = [slice(ind_boundaries[0][0], ind_boundaries[0][0] +1),
- slice(ind_boundaries[1][0], ind_boundaries[1][1] +1),
- slice(ind_boundaries[2][0], ind_boundaries[2][1] +1)]
-
- self.chi_x_plus = [slice(ind_boundaries[0][1], ind_boundaries[0][1] +1 ),
- slice(ind_boundaries[1][0], ind_boundaries[1][1] +1),
- slice(ind_boundaries[2][0], ind_boundaries[2][1] +1)]
-
- self.chi_y_minus = [slice(ind_boundaries[0][0], ind_boundaries[0][1] +1),
- slice(ind_boundaries[1][0], ind_boundaries[1][0] +1),
- slice(ind_boundaries[2][0], ind_boundaries[2][1] +1)]
-
- self.chi_y_plus = [slice(ind_boundaries[0][0], ind_boundaries[0][1] +1),
- slice(ind_boundaries[1][0], ind_boundaries[1][0] +1),
- slice(ind_boundaries[2][0], ind_boundaries[2][1] +1)]
-
- self.chi_z_minus = [slice(ind_boundaries[0][0], ind_boundaries[0][1] +1),
- slice(ind_boundaries[1][0], ind_boundaries[1][1] +1),
- slice(ind_boundaries[2][0], ind_boundaries[2][0] +1)]
-
- self.chi_z_plus = [slice(ind_boundaries[0][0], ind_boundaries[0][1] +1),
- slice(ind_boundaries[1][0], ind_boundaries[1][1] +1),
- slice(ind_boundaries[2][1], ind_boundaries[2][1] +1)]
-
- # definition of the indicator function of the whole box (includind the obstacle)
- # TODO : defining this indicator function without includind the obstacle
- self.chi_box = []
- for d in xrange(self.dim):
- self.chi_box.append(slice(ind_boundaries[d][0], ind_boundaries[d][1]))
-
- @debug
- @timed_function
- def apply(self, simulation):
- """
- Computation of the drag according to the "impulsion" formula presented by
- Noca et. al in Journal of Fluids and Structures, 1999.
- Integrals on the obstacle are neglected.
- \f$ F=-\frac{1}{N-1}\frac{d}{dt}\int_{V}\mathbf{x}\wedge\mathbf{\omega}\, dV
- +\int_{S}\mathbf{n}\cdot\gamma\, dS \f$
- """
- t = simulation.time
- dt = simulation.timeStep
- ite = simulation.currentIteration
-
- velo = self.velo.data
- vort = self.vort.data
-
- # Synchronize ghost points of velocity
- self._synchronize(velo)
-
- if not self.isForceComputationNeeded:
- pass
- else :
- localForce = npw.zeros(self.dim)
- force = npw.zeros(self.dim)
-
- # Surfaces normal to X axis
- # n=(-1, 0, 0)
- dsurf = self.step[YDIR] * self.step[ZDIR]
- localForce = self.integrateOnSurface(localForce, velo, vort, self.chi_x_minus,
- self.normal[self.x_minus][:],
- dsurf, direc=XDIR)
- # n=(1, 0, 0)
- localForce = self.integrateOnSurface(localForce, velo, vort, self.chi_x_plus,
- self.normal[self.x_plus][:],
- dsurf, direc=XDIR)
- # Surfaces normal to Y axis
- dsurf = self.step[XDIR] * self.step[ZDIR]
- # n=(0,-1, 0)
- localForce = self.integrateOnSurface(localForce, velo, vort, self.chi_y_minus,
- self.normal[self.y_minus][:],
- dsurf, direc=YDIR)
- # n=(0, 1, 0)
- localForce = self.integrateOnSurface(localForce, velo, vort, self.chi_y_plus,
- self.normal[self.y_plus][:],
- dsurf, direc=YDIR)
- # Surfaces normal to Z axis
- dsurf = self.step[XDIR] * self.step[YDIR]
- # n=(0, 0,-1)
- localForce = self.integrateOnSurface(localForce, velo, vort, self.chi_z_minus,
- self.normal[self.z_minus][:],
- dsurf, direc=ZDIR)
- # n=(0, 0, 1)
- localForce = self.integrateOnSurface(localForce, velo, vort, self.chi_z_plus,
- self.normal[self.z_plus][:],
- dsurf, direc=ZDIR)
- # Integration over the volume
- localForce = self.integrateOnBox(localForce, vort, self.chi_box, dt)
-
- # Reduction operation
- if self.topo.comm.Get_rank() == 0:
- comm = self.topo.comm
- comm.Reduce([localForce, PARMES_MPI_REAL],
- [force, PARMES_MPI_REAL],
- op=MPI.SUM, root=0)
- force = force * self.coef
- # Print drag, lift and side force coefficients in output file
- if ((ite % self.freq == 0) and (t > dt)):
- self.f = open(self.prefix, 'a')
- self.f.write("%s %s %s %s\n" % (t,
- force[XDIR],
- force[YDIR],
- force[ZDIR]))
- self.f.close()
-
- def integrateOnBox(self, force, vort, chi, dt):
- ## WARNING : this routine is available inly if velo and vort
- ## topologies have the same global mesh resolution !!!
- """
- Returns \f$-\frac{1}{2}\frac{d}{dt}
- \int_{V}\mathbf{x}\wedge\mathbf{\omega}\, dV \f$
- """
- dvol = self.step[XDIR] * self.step[YDIR] * self.step[ZDIR]
- fact = - dvol / ((self.dim - 1) * dt)
- intbox = npw.zeros(self.dim)
- work = [npw.zeros_like(vort[0]) for i in xrange(self.dim)]
-
- # (coord X vorticity)
- coords = np.asarray(self.coords)
- work[XDIR][...] = self.coords[YDIR] * vort[ZDIR][...] - \
- self.coords[ZDIR] * vort[YDIR][...]
- work[YDIR][...] = self.coords[ZDIR] * vort[XDIR][...] - \
- self.coords[XDIR] * vort[ZDIR][...]
- work[ZDIR][...] = self.coords[XDIR] * vort[YDIR][...] - \
- self.coords[YDIR] * vort[XDIR][...]
- for i in xrange(self.dim):
- intbox[i] = np.sum(work[i][chi])
- force += fact * (intbox - self.bufferForce) # 1st order time integration
- self.bufferForce = intbox # Save for next time step ...
-
- return force
-
-
- def integrateOnSurface(self, force, velo, vort, chi, NormalVec, dsurf, direc):
- ## WARNING : this routine is available inly if velo and vort
- ## topologies have the same global mesh resolution !!!
- """
- Compute integrals on surface normal to "direc" axis to calculate forces
- acting on the body. See (2.1) of Noca 1999 or (52) of Plouhmans 2002.
- Returns \f$ \int_{S}\mathbf{n}\cdot\gamma\, dS \f$, where
- \f$ \mathbf{n}\cdot\gamma = [(\frac{1}{2} (\mathbf{u}\cdot\mathbf{u})\mathbf{n} -
- (\mathbf{n}\cdot\mathbf{u})\mathbf{u}) -
- \frac{1}{N-1}(\mathbf{n}\cdot\mathbf{u})(\mathbf{x}\wedge\mathbf{\omega}) +
- \frac{1}{N-1}(\mathbf{n}\cdot\mathbf{\omega})(\mathbf{x}\wedge\mathbf{u})] +
- [\frac{1}{N-1}\mathbf{x}\wedge(\mathbf{n}\wedge\nabla\cdot T)] +
- [\frac{1}{N-1}\mathbf{n}\cdot T] \f$
- In the following, the computation of surface integrals is split in three parts :
- force = [int_part 1] + [int_part 2] + [int_part 3]
- """
- fact = 1. / (self.dim - 1)
- work = [npw.zeros_like(velo[0]) for i in xrange(self.dim)]
- divT = [npw.zeros_like(velo[0]) for i in xrange(self.dim)]
- # Divergence of stress tensor T, computed with 2nd order FD
- divT = self.divT(velo, divT, chi)
- # Functions to compute part 2 and part 3 of the surface integral, with:
- # part 2 = 1/(dim-1)(X^(n^Div(T))
- # part 3 = (1/Re) * n.T
- if direc==XDIR:
- self.X_n_DivT = self.part2_Xdir
- self.n_T = self.part3_Xdir
- elif direc==YDIR:
- self.X_n_DivT = self.part2_Ydir
- self.n_T = self.part3_Ydir
- else:
- assert direc==ZDIR
- self.X_n_DivT = self.part2_Zdir
- self.n_T = self.part3_Zdir
-
- # ====================== part 1 =======================
- # part 1 = 1/2(velocity.velocity)n -
- # (n.velocity)velocity -
- # 1/(dim-1)(n.velocity)(coord X vorticity) +
- # 1/(dim-1)(n.vorticity)(coord X velocity)
-
- # 1/2(velocity.velocity)n
- u_u = np.sum(velo[i][chi] * velo[i][chi] for i in xrange(self.dim))
- force[direc] += 0.5 * np.sum(u_u[...] * NormalVec[direc])
- # -(n.velocity)velocity
- n_u = velo[direc][chi] * NormalVec[direc]
- for i in xrange(self.dim):
- work[i][chi] = n_u[...] * velo[i][chi]
- for i in xrange(self.dim):
- force[i] -= np.sum(work[i][chi])
-
- # -1/(dim-1)(n.velocity)(coord X vorticity)
- work[XDIR][...] = self.coords[XDIR] * vort[ZDIR][...] - \
- self.coords[ZDIR] * vort[YDIR][...]
- work[YDIR][...] = self.coords[ZDIR] * vort[XDIR][...] - \
- self.coords[XDIR] * vort[ZDIR][...]
- work[ZDIR][...] = self.coords[XDIR] * vort[YDIR][...] - \
- self.coords[YDIR] * vort[XDIR][...]
- for i in xrange(self.dim):
- work[i][chi] *= n_u[...]
- for i in xrange(self.dim):
- force[i] -= fact * np.sum(work[i][chi])
-
- # 1/(dim-1)(n.vorticity)(coord X velocity)
- n_w = vort[direc][chi] * NormalVec[direc]
- work[XDIR][...] = self.coords[XDIR] * velo[ZDIR][...] - \
- self.coords[ZDIR] * velo[YDIR][...]
- work[YDIR][...] = self.coords[ZDIR] * velo[XDIR][...] - \
- self.coords[XDIR] * velo[ZDIR][...]
- work[ZDIR][...] = self.coords[XDIR] * velo[YDIR][...] - \
- self.coords[YDIR] * velo[XDIR][...]
- for i in xrange(self.dim):
- work[i][chi] *= n_w[...]
- for i in xrange(self.dim):
- force[i] += fact * np.sum(work[i][chi])
-
- # ====================== part 2 =======================
- # part 2 = 1/(dim-1)(X^(n^Div(T))
- work = self.X_n_DivT(chi, divT, work)
-
- if NormalVec[direc] == 1:
- for i in xrange(self.dim):
- force[i] += (fact / self.Re) * np.sum(work[i][chi])
- else:
- assert NormalVec[direc] == -1
- for i in xrange(self.dim):
- force[i] += (fact / self.Re) * np.sum(-work[i][chi])
-
- # ====================== part 3 =======================
- # part 3 = (1/Re) * n.T
- work = self.n_T(velo, work, chi)
-
- if NormalVec[direc] == 1:
- for i in xrange(self.dim):
- force[i] += (1.0 / self.Re) * np.sum(work[i][chi])
- else:
- assert NormalVec[direc] == -1
- for i in xrange(self.dim):
- force[i] += (1.0 / self.Re) * np.sum(-work[i][chi])
-
- force *= dsurf
-
- return force
-
- # Functions to compute part 2 of a surface integral
- def part2_Xdir(self, chi, divT, work):
- # X^(n^Div(T)) computation for normalVec = (1,0,0)
- work[XDIR][...] = self.coords[YDIR] * divT[YDIR][...] + \
- self.coords[ZDIR] * divT[ZDIR][...]
- work[YDIR][...] = - self.coords[XDIR] * divT[YDIR][...]
- work[ZDIR][...] = - self.coords[XDIR] * divT[ZDIR][...]
-
- return work
-
- def part2_Ydir(self, chi, divT, work):
- # X^(n^Div(T)) computation for normalVec = (0,1,0)
- work[XDIR][...] = - self.coords[YDIR] * divT[XDIR][...]
- work[YDIR][...] = self.coords[XDIR] * divT[XDIR][...] + \
- self.coords[ZDIR] * divT[ZDIR][...]
- work[ZDIR][...] = - self.coords[YDIR] * divT[ZDIR][...]
-
- return work
-
- def part2_Zdir(self, chi, divT, work):
- # X^(n^Div(T)) computation for normalVec = (0,0,1)
- work[XDIR][...] = - self.coords[ZDIR] * divT[XDIR][...]
- work[YDIR][...] = - self.coords[ZDIR] * divT[YDIR][...]
- work[ZDIR][...] = self.coords[XDIR] * divT[XDIR][...] + \
- self.coords[YDIR] * divT[YDIR][...]
-
- return work
-
- # Functions to compute part 3 of a surface integral
- def part3_Xdir(self, velo, work, chi):
- # n.T[XDIR] = 2 * dux/dx[i,j,k]
- # n.T[YDIR] = duy/dx[i,j,k] + dux/dy[i,j,k]
- # n.T[ZDIR] = duz/dx[i,j,k] + dux/dz[i,j,k]
-
- self.fd_scheme.computeIndices(chi)
- workJacobian = npw.zeros_like(velo[0])
-
- # n.T computation for normalVec = (1,0,0)
- self.fd_scheme.compute_1st_deriv(velo[XDIR], XDIR, workJacobian)
- work[XDIR][chi] = 2 * workJacobian[chi]
-
- self.fd_scheme.compute_1st_deriv(velo[YDIR], XDIR, workJacobian)
- self.fd_scheme.compute_and_add_1st_deriv(velo[XDIR], YDIR, workJacobian)
- work[YDIR][chi] = workJacobian[chi]
-
- self.fd_scheme.compute_1st_deriv(velo[ZDIR], XDIR, workJacobian)
- self.fd_scheme.compute_and_add_1st_deriv(velo[XDIR], ZDIR, workJacobian)
- work[ZDIR][chi] = workJacobian[chi]
-
- return work
-
- def part3_Ydir(self, velo, work, chi):
- # n.T[XDIR] = dux/dy[i,j,k] + duy/dx[i,j,k]
- # n.T[YDIR] = 2.0 * duy/dy[i,j,k]
- # n.T[ZDIR] = duz/dy[i,j,k] + duy/dz[i,j,k]
-
- self.fd_scheme.computeIndices(chi)
- workJacobian = npw.zeros_like(velo[0])
-
- # n.T computation for normalVec = (0,1,0)
- self.fd_scheme.compute_1st_deriv(velo[XDIR], YDIR, workJacobian)
- self.fd_scheme.compute_and_add_1st_deriv(velo[YDIR], XDIR, workJacobian)
- work[XDIR][chi] = workJacobian[chi]
-
- self.fd_scheme.compute_1st_deriv(velo[YDIR], YDIR, workJacobian)
- work[YDIR][chi] = 2.* workJacobian[chi]
-
- self.fd_scheme.compute_1st_deriv(velo[ZDIR], YDIR, workJacobian)
- self.fd_scheme.compute_and_add_1st_deriv(velo[YDIR], ZDIR, workJacobian)
- work[ZDIR][chi] = workJacobian[chi]
-
- return work
-
- def part3_Zdir(self, velo, work, chi):
- # n.T[XDIR] = dux/dz[i,j,k] + duz/dx[i,j,k]
- # n.T[YDIR] = duy/dz[i,j,k] + duz/dy[i,j,k]
- # n.T[ZDIR] = 2.0 * duz/dz[i,j,k]
-
- self.fd_scheme.computeIndices(chi)
- workJacobian = npw.zeros_like(velo[0])
-
- # n.T computation for normalVec = (0,0,1)
- self.fd_scheme.compute_1st_deriv(velo[XDIR], ZDIR, workJacobian)
- self.fd_scheme.compute_and_add_1st_deriv(velo[ZDIR], XDIR, workJacobian)
- work[XDIR][chi] = workJacobian[chi]
-
- self.fd_scheme.compute_1st_deriv(velo[YDIR], ZDIR, workJacobian)
- self.fd_scheme.compute_and_add_1st_deriv(velo[ZDIR], YDIR, workJacobian)
- work[YDIR][chi] = workJacobian[chi]
-
- self.fd_scheme.compute_1st_deriv(velo[ZDIR], ZDIR, workJacobian)
- work[ZDIR][chi] = 2.* workJacobian[chi]
-
- return work
-
- def __str__(self):
- s = "Compute_forces. "
- return s
diff --git a/HySoP/hysop/operator/monitors/monitoring.py b/HySoP/hysop/operator/monitors/monitoring.py
index c8703fbb466a5a3635e5e3f45a7e8bd8af29b12f..c96932fddfab2a3195bf21e01e8c99500adcccf6 100644
--- a/HySoP/hysop/operator/monitors/monitoring.py
+++ b/HySoP/hysop/operator/monitors/monitoring.py
@@ -5,6 +5,7 @@ Global interface for monitors.
from abc import ABCMeta, abstractmethod
from parmepy.operator.continuous import Operator
from parmepy.tools.timers import Timer
+import parmepy.tools.io_utils as io
class Monitoring(Operator):
@@ -13,27 +14,45 @@ class Monitoring(Operator):
__metaclass__ = ABCMeta
@abstractmethod
- def __init__(self, variables, topo, frequency, task_id=None):
+ def __init__(self, variables, topo, io_params=None, task_id=None):
""" Constructor
@param variables : list of fields to monitor.
@param topo : topo on which fields are to be monitored
- @param frequency : output rate.
- @param name : Monitor name.
+ @param filename output file name. Default is None ==> no output.
+ Full or relative path.
+ @param io_params : parameters (dict) to set file output.
+ If None, no output. Set io_params = {} if you want output,
+ with default parameters values.
+ See parmepy.tools.io_utils.Writer for details
"""
Operator.__init__(self, variables, topo=topo, task_id=task_id)
- ## Monitor frequency
- self.freq = frequency
## Object to store computational times of lower level functions
self.timer = Timer(self)
self.requirements = []
+ self._isUpToDate = False
+ # Output file.
+ if io_params is None:
+ self._writer = None
+ else:
+ self._writer = io.Writer(io_params, topo.comm)
def setUp(self):
- pass
+ """
+ Class-method required
+ to fit with operator base-class interface.
+ """
+ self._isUpToDate = True
def finalize(self):
- pass
+ if self._writer is not None:
+ self._writer.finalize()
def discretize(self):
+ """
+ Does nothing.
+ No discretization in monitors, but class-method required
+ to fit with operator base-class interface.
+ """
pass
def addRedistributeRequirement(self, red):
diff --git a/HySoP/hysop/operator/monitors/printer.py b/HySoP/hysop/operator/monitors/printer.py
index 7d9bbc0e7141ae270e23a98b8fc62ceaa9151f64..35022c94cf2fb84b1ca67a125f54ef3bd2782e5a 100644
--- a/HySoP/hysop/operator/monitors/printer.py
+++ b/HySoP/hysop/operator/monitors/printer.py
@@ -1,12 +1,14 @@
"""
@file printer.py
-Classes for handling ouputs.
+File output for field(s) value on a grid.
"""
from parmepy.constants import np, S_DIR, debug, VTK, HDF5, DATA
from parmepy.operator.monitors.monitoring import Monitoring
import parmepy.tools.numpywrappers as npw
+import parmepy.tools.io_utils as io
import os
+
try:
import evtk.hl as evtk
except ImportError as evtk_error:
@@ -20,29 +22,23 @@ from parmepy.tools.timers import timed_function
class Printer(Monitoring):
"""
- Basic printer.
-
- Performs outputs in VTK images.
+ Print field(s) values on a given topo, in HDF5 or VTK format.
"""
- def __init__(self, variables, topo, frequency=0,
- prefix=None, formattype=None, task_id=None):
+ def __init__(self, variables, topo, prefix=None,
+ frequency=1, formattype=None, task_id=None):
"""
- Create a results printer for given fields, filename
+ Create a results printer for given fields, prefix
prefix (relative path) and an output frequency.
- @param fields : list of variables to export.
+ @param variables : list of variables to export.
@param frequency : output rate (output every freq iteration)
- @param prefix : output file name prefix.
- @param ext : output file extension, default=.vtk.
+ @param prefix output file name. Default is None ==> no output.
+ Full or relative path.
+ @param formattype : output file format, default=vtk.
"""
- Monitoring.__init__(self, variables, topo, frequency,
- task_id=task_id)
- ## output file name prefix
- if prefix is None:
- self.prefix = './out_'
- else:
- self.prefix = prefix
-
+ Monitoring.__init__(self, variables, topo, task_id=task_id)
+ assert frequency > 0
+ self.frequency = frequency
## Default output type
if formattype is None:
self.formattype = VTK
@@ -61,41 +57,44 @@ class Printer(Monitoring):
self.input = self.variables
self.output = []
- ## Method to collect data in case of distributed data
- self.get_data_method = None
- if self.freq > 0:
- ## Printer core method
- self.step = self._printStep
+ # If no prefix is given, set it to
+ # the concatenation of variables'names.
+ if prefix is None:
+ prefix = io.io.defaultPath()
+ name = ''
+ for var in self.input:
+ name += var.name
+ name += '_'
+ prefix = os.path.join(prefix, name)
else:
- self.step = self._passStep
- # Internal counter
- self._call_number = 0
+ if not os.path.isabs(prefix):
+ prefix = os.path.join(io.io.defaultPath(), prefix)
## shortcut to topo name
self.topo = self._predefinedTopo
+ io.io.checkDir(prefix, 0, self.topo.comm)
+ self.prefix = prefix
self._xmf_data_files = []
-
-
- def setUp(self):
- """
- Print initial state
- """
- # Create output dir if required
- if self.topo.rank == 0:
- d = os.path.dirname(self.prefix)
- if not os.path.exists(d):
- os.makedirs(d)
- # Force synchro to be sure that all output dirs
- # have been created.
- self.topo.comm.barrier()
if self.formattype == VTK:
+ # filename = prefix_rk_N, rk = current process rank,
+ # N = counter value
self._get_filename = lambda i: self.prefix + \
- "{0}_iter_{1:03d}".format(self.topo.rank, i)
+ "_{0}_{1:05d}".format(self.topo.rank, i)
+ self.step = self._step_VTK
elif self.formattype == HDF5:
+ # filename = prefix_N, N = counter value
self._get_filename = lambda i: self.prefix + \
- "{0}procs_iter_{1:03d}".format(self.topo.size, i) + '.h5'
+ "_{0:05d}".format(i) + '.h5'
+ self.step = self._step_HDF5
elif self.formattype == DATA:
self._get_filename = lambda i: self.prefix + \
- "{0}_iter_{1:03d}.dat".format(self.topo.rank, i)
+ "_{0}_{1:05d}.dat".format(self.topo.rank, i)
+ self.step = self._step_DATA
+ self._isUpToDate = True
+ # count the number of calls
+ self._count = 0
+ ## Rank of 'leader' mpi process for output
+ ## Forced to 0. \todo : add a param for this?
+ self.io_rank = 0
@debug
@timed_function
@@ -103,21 +102,29 @@ class Printer(Monitoring):
if simulation is None:
raise ValueError("Missing simulation value for monitoring.")
- self._call_number += 1
- if self.freq > 0 and (simulation.currentIteration % self.freq) == 0:
+ if simulation.currentIteration % self.frequency == 0:
+ # Transfer from GPU to CPU if required
+ for f in self.variables:
+ df = f.discreteFields[self.topo]
+ try:
+ if not df.isBatch:
+ # To host only if data fit in the device memory
+ df.toHost()
+ df.wait()
+ except AttributeError:
+ pass
self.step(simulation)
+ self._count += 1
def finalize(self):
- if self.freq == -1:
- self._printStep(self._call_number)
- if self.formattype == HDF5 and self.topo.rank == 0:
+ if self.formattype == HDF5 and self.topo.rank == self.io_rank:
# Write the xmf file driving all h5 files.
# Writing only one file
# We have a temporal list of Grid => a single xmf file
# Manual writing of the xmf file because "XDMF library very
# difficult to compile and to use"
# [Advanced HDF5 & XDMF - Groupe Calcul]
- f = open(self.prefix + str(self.topo.size) + 'procs.xmf', 'w')
+ f = open(self.prefix + '.xmf', 'w')
f.write("<?xml version=\"1.0\" ?>\n")
f.write("<!DOCTYPE Xdmf SYSTEM \"Xdmf.dtd\">\n")
f.write("<Xdmf Version=\"2.0\">\n")
@@ -151,151 +158,119 @@ class Printer(Monitoring):
res[name] = npw.realarray(df.data[d])
return res
- def set_get_data_method(self, method):
- """Set data collect method."""
- if callable(method):
- self.get_data_method = method
- else:
- raise ValueError("Cannot set non callable method to get data " +
- "to print. Given method : " + str(method))
-
- def _passStep(self, simu=None):
- """A code method that do nothing."""
- pass
-
- def _printStep(self, simu=None):
+ def _step_VTK(self, simu):
"""
- Try to write data into VTK files.
- If fails, turn to classical ascii output.
- @param iter : current iteration number
"""
- ite = simu.currentIteration
- t = ite * simu.timeStep
- # Transfer from GPU to CPU if required
- for f in self.variables:
- df = f.discreteFields[self.topo]
- try:
- if not df.isBatch:
- # To host only if data fit in the device memory
- df.toHost()
- df.wait()
- except AttributeError:
- pass
+ orig = [0.] * 3
+ dim = self.topo.mesh.dim
+ orig[:dim] = [self.topo.mesh.origin[i] for i in xrange(dim)]
+ #orig = tuple(orig)
+ #ind = self.topo.mesh.local_start
+ ## orig = tuple([self.topo.mesh.coords[i].flatten()[ind[i]]
+ ## for i in xrange(self.topo.mesh.dim)])
+ spacing = [0] * 3
+ spacing[:dim] = [self.topo.mesh.space_step[i] for i in xrange(dim)]
+ evtk.imageToVTK(self._get_filename(self._count),
+ origin=orig, spacing=spacing,
+ pointData=self._build_vtk_dict())
- ## VTK output \todo: Need fix in 2D, getting an IOError.
- if self.formattype == VTK:
- orig = [0.] * 3
- dim = self.topo.mesh.dim
- orig[:dim] = [self.topo.mesh.origin[i] for i in xrange(dim)]
- #orig = tuple(orig)
- coords = self.topo.mesh.coords
- #ind = self.topo.mesh.local_start
- ## orig = tuple([self.topo.mesh.coords[i].flatten()[ind[i]]
- ## for i in xrange(self.topo.mesh.dim)])
- spacing = [0] * 3
- spacing[:dim] = [self.topo.mesh.space_step[i] for i in xrange(dim)]
- evtk.imageToVTK(self._get_filename(ite),
- origin=orig, spacing=spacing,
- pointData=self._build_vtk_dict())
-
- elif self.formattype == HDF5:
- filename = self._get_filename(ite)
- # Write the h5 file
- # (force np.float64, ParaView seems to not be able to read float32)
- # Writing compressed hdf5 files (gzip compression seems the best)
- # In parallel, one file is written, thus filename no longe contains
- # mpi rank.
- # TODO: Why gzip compression not working in parallel ??
- # Remark: h5py must be build with --mpi option
- if self.topo.size == 1:
- f = h5py.File(filename, "w")
- compression = 'gzip'
- else:
- f = h5py.File(filename, 'w',
- driver='mpio', comm=self.topo.comm)
- compression = None
+ def _step_HDF5(self, simu):
+ t = simu.time
+ filename = self._get_filename(self._count)
+ # Write the h5 file
+ # (force np.float64, ParaView seems to not be able to read float32)
+ # Writing compressed hdf5 files (gzip compression seems the best)
+ # In parallel, one file is written, thus filename no longe contains
+ # mpi rank.
+ # TODO: Why gzip compression not working in parallel ??
+ # Remark: h5py must be build with --mpi option
+ if self.topo.size == 1:
+ f = h5py.File(filename, "w")
+ compression = 'gzip'
+ else:
+ f = h5py.File(filename, 'w', driver='mpio', comm=self.topo.comm)
+ compression = None
+ g_start = self.topo.mesh.global_start
+ g_end = self.topo.mesh.global_end + 1
+ sl = [slice(g_start[i], g_end[i])
+ for i in xrange(self.domain.dimension)]
+ sl.reverse()
+ sl = tuple(sl)
+ datasetNames = []
+ for field in self.variables:
+ df = field.discreteFields[self.topo]
+ for d in xrange(df.nbComponents):
+ # creating datasets for the vector field
+ currentName = df.name + S_DIR[d]
+ datasetNames.append(currentName)
+ res = list(self.topo.globalMeshResolution - 1)
+ res.reverse()
+ ds = f.create_dataset(currentName,
+ res,
+ dtype=np.float64,
+ compression=compression)
+ # In parallel, each proc must write in the proper part
+ # Of the dataset (of site global resolution)
+ ds[sl] = npw.realarray(df.data[d][self.topo.mesh.iCompute].T)
+ self._xmf_data_files.append((datasetNames, self._count, t))
- g_start = self.topo.mesh.global_start
- g_end = self.topo.mesh.global_end + 1
- sl = [slice(g_start[i], g_end[i])
- for i in xrange(self.domain.dimension)]
- sl.reverse()
- sl = tuple(sl)
- datasetNames = []
- for field in self.variables:
- df = field.discreteFields[self.topo]
- for d in xrange(df.nbComponents):
- # creating datasets for the vector field
- currentName = df.name + S_DIR[d]
- datasetNames.append(currentName)
- res = list(self.topo.globalMeshResolution - 1)
- res.reverse()
- ds = f.create_dataset(currentName,
- res,
- dtype=np.float64,
- compression=compression)
- # In parallel, each proc must write in the proper part
- # Of the dataset (of site global resolution)
- ds[sl] = \
- npw.realarray(df.data[d][self.topo.mesh.iCompute].T)
- self._xmf_data_files.append((datasetNames, ite, t))
+ f.close()
- f.close()
- elif self.formattype == DATA:
- f = open(self._get_filename(ite), 'w')
- shape = self.topo.mesh.resolution
- coords = self.topo.mesh.coords
- pbDimension = self.domain.dimension
- if pbDimension == 2:
- if len(shape) == 2:
- for i in xrange(shape[0] - 1):
- for j in xrange(shape[1] - 1):
- f.write("{0:8.12} {1:8.12} ".format(
- coords[0][i, 0], coords[1][0, j]))
- for field in self.variables:
- df = field.discreteFields[self.topo]
- if field.isVector:
- f.write("{0:8.12} {1:8.12} ".format(
- df[0][i, j],
- df[1][i, j]))
- else:
- f.write("{0:8.12} ".format(
- df[0][i, j]))
- f.write("\n")
- elif pbDimension == 3:
+ def _step_DATA(self, simu):
+ f = open(self._get_filename(self._count), 'w')
+ shape = self.topo.mesh.resolution
+ coords = self.topo.mesh.coords
+ pbDimension = self.domain.dimension
+ if pbDimension == 2:
+ if len(shape) == 2:
for i in xrange(shape[0] - 1):
for j in xrange(shape[1] - 1):
- for k in xrange(shape[2] - 1):
- f.write(
- "{0:8.12} {1:8.12} {2:8.12} ".format(
- coords[0][i, 0, 0],
- coords[1][0, j, 0],
- coords[2][0, 0, k]))
- for field in self.variables:
- df = field.discreteFields[self.topo]
- if field.isVector:
- f.write(
- "{0:8.12} {1:8.12} " +
- "{2:8.12} ".format(
- df[0][i, j, k],
- df[1][i, j, k],
- df[2][i, j, k]))
- else:
- f.write("{0:8.12} ".format(
- df[0][i, j, k]))
- f.write("\n")
- else:
- for i in xrange(shape[0] - 1):
- f.write("{0:8.12} ".format(coords[0][i]))
- for field in self.variables:
- df = field.discreteFields[self.topo]
- if field.isVector:
- f.write(
- "{0:8.12} ".format(df[0][i]))
- else:
- f.write("{0:8.12} ".format(df[i]))
- f.write("\n")
- f.close()
+ f.write("{0:8.12} {1:8.12} ".format(
+ coords[0][i, 0], coords[1][0, j]))
+ for field in self.variables:
+ df = field.discreteFields[self.topo]
+ if field.isVector:
+ f.write("{0:8.12} {1:8.12} ".format(
+ df[0][i, j],
+ df[1][i, j]))
+ else:
+ f.write("{0:8.12} ".format(
+ df[0][i, j]))
+ f.write("\n")
+ elif pbDimension == 3:
+ for i in xrange(shape[0] - 1):
+ for j in xrange(shape[1] - 1):
+ for k in xrange(shape[2] - 1):
+ f.write(
+ "{0:8.12} {1:8.12} {2:8.12} ".format(
+ coords[0][i, 0, 0],
+ coords[1][0, j, 0],
+ coords[2][0, 0, k]))
+ for field in self.variables:
+ df = field.discreteFields[self.topo]
+ if field.isVector:
+ f.write(
+ "{0:8.12} {1:8.12} " +
+ "{2:8.12} ".format(
+ df[0][i, j, k],
+ df[1][i, j, k],
+ df[2][i, j, k]))
+ else:
+ f.write("{0:8.12} ".format(
+ df[0][i, j, k]))
+ f.write("\n")
+ else:
+ for i in xrange(shape[0] - 1):
+ f.write("{0:8.12} ".format(coords[0][i]))
+ for field in self.variables:
+ df = field.discreteFields[self.topo]
+ if field.isVector:
+ f.write(
+ "{0:8.12} ".format(df[0][i]))
+ else:
+ f.write("{0:8.12} ".format(df[i]))
+ f.write("\n")
+ f.close()
def _listFormat(l):
@@ -332,7 +307,7 @@ def _TemporalGridXMF(topo, datasetNames, ite, t, filename):
g += " " + _listFormat(ori) + "\n"
g += " </DataItem>\n"
g += " <DataItem Dimensions=\"" + str(topo.mesh.dim) + " \""
- g += " NumberType=\"Float\" Precision=\"4\" Format=\"XML\">\n"
+ g += " NumberType=\"Float\" Precision=\"8\" Format=\"XML\">\n"
step = list(topo.mesh.space_step)
step.reverse()
g += " " + _listFormat(step) + "\n"
diff --git a/HySoP/hysop/operator/poisson.py b/HySoP/hysop/operator/poisson.py
index 768d56466b919de611923599d182111de0c2c7d8..a9c7a454f8eb8709a0ee6b28dbf8c8ca433118de 100644
--- a/HySoP/hysop/operator/poisson.py
+++ b/HySoP/hysop/operator/poisson.py
@@ -7,9 +7,10 @@ Poisson problem.
from parmepy.operator.continuous import Operator
from parmepy.f2py import fftw2py
from parmepy.operator.discrete.poisson_fft import PoissonFFT
-from parmepy.variables.variables import Variables
from parmepy.constants import debug
import numpy as np
+from parmepy.mpi.main_var import main_size
+from parmepy.operator.velocity_correction import VelocityCorrection
class Poisson(Operator):
@@ -26,14 +27,19 @@ class Poisson(Operator):
@debug
def __init__(self, velocity, vorticity, resolutions, ghosts=None,
- method=None, topo=None, task_id=None, comm=None,
- **other_config):
+ method=None, flowrate=None, topo=None,
+ task_id=None, comm=None):
"""
Constructor for the Poisson problem.
@param[out] velocity : solution field
@param[in] vorticity : rhs field
@param[in] resolutions : list of resolutions (one per variable)
+ @param[in] ghosts : ghost layer (optional), default=[0, 0, ...]
+ @param[in] method : numerical method (default = fft)
+ @param[in] flowrate : a flow rate value (through input surf,
+ normal to xdir) used to compute a correction of the solution field.
+ Default = 0 (no correction). See parmepy.operator.velocity_correction.
"""
## solution of the problem
self.velocity = velocity
@@ -45,7 +51,6 @@ class Poisson(Operator):
Operator.__init__(self, [velocity, vorticity], method, ghosts=ghosts,
topo=topo, task_id=task_id, comm=comm)
- self.config = other_config
if method is None:
self.method = 'fftw'
if self.method is 'fftw':
@@ -53,25 +58,20 @@ class Poisson(Operator):
assert self.resolution == self.resolutions[self.vorticity],\
'Poisson error : for fftw, all variables must have\
the same global resolution.'
+ else:
+ raise AttributeError("Method not yet implemented.")
self.input = [self.vorticity]
self.output = [self.velocity]
-
- # If there is a projection, vorticity is also an output
- try:
- if (isinstance(other_config['projection'], Variables)):
- proj = other_config['projection'].data
- else:
- proj = other_config['projection']
- except KeyError:
- proj = None
- if proj is not None and proj[0]:
- self.output.append(self.vorticity)
+ if flowrate is not None:
+ self.withCorrection = True
+ self._flowrate = flowrate
+ else:
+ self.withCorrection = False
+ self.correction = None
+ self.projection = None
def discretize(self):
# Poisson solver based on fftw
- if self.method is not 'fftw':
- print self.method
- raise AttributeError("Method not yet implemented.")
if self.ghosts is not None:
raise AttributeError("Ghosts points not yet\
implemented for poisson operator.")
@@ -88,13 +88,13 @@ class Poisson(Operator):
self.resolution, self.domain.length, comm=comm.py2f())
if self._predefinedTopo is not None:
- main_size = self._predefinedTopo.size
+ commsize = self._predefinedTopo.size
elif self._comm is not None:
- main_size = self._comm.Get_size()
+ commsize = self._comm.Get_size()
else:
- from parmepy.mpi.main_var import main_size
+ commsize = main_size
topodims = np.ones((self.domain.dimension))
- topodims[-1] = main_size
+ topodims[-1] = commsize
#variables discretization
if self._predefinedTopo is not None:
topo = self._predefinedTopo
@@ -113,15 +113,32 @@ class Poisson(Operator):
ghosts=self.ghosts,
comm=self._comm)
self.discreteFields[v] = v.discretize(topo)
+ # Build and discretize correction operator, if necessary
+ if self.withCorrection:
+ self.correction = VelocityCorrection(self.velocity, self.vorticity,
+ resolutions=self.resolutions,
+ flowrate=self._flowrate,
+ topo=self.velocity.topology)
+ self.correction.discretize()
@debug
def setUp(self):
"""
"""
# Build and setup for the discrete operator
+ if self.withCorrection:
+ self.correction.setUp()
+ cd = self.correction.discreteOperator
+ else:
+ cd = None
self.discreteOperator = PoissonFFT(self.discreteFields[self.velocity],
self.discreteFields[self.vorticity],
- method=self.method,
- **self.config)
+ correction=cd,
+ projection=self.projection)
self.discreteOperator.setUp()
self._isUpToDate = True
+
+ def activateProjection(self, projection):
+ # If there is a projection, vorticity is also an output
+ self.output.append(self.vorticity)
+ self.projection = projection
diff --git a/HySoP/hysop/operator/redistribute.py b/HySoP/hysop/operator/redistribute.py
index d1f2a8c731cf3e6acf65ea530d8be93b697b0433..d805e01eaae9915a081a0d6ac59ca2a91fbf35d6 100644
--- a/HySoP/hysop/operator/redistribute.py
+++ b/HySoP/hysop/operator/redistribute.py
@@ -23,6 +23,7 @@ from parmepy import __VERBOSE__
from parmepy.constants import debug, PARMES_MPI_REAL, ORDERMPI, np, S_DIR
from parmepy.operator.continuous import Operator
from parmepy.mpi.topology import Bridge
+from parmepy.mpi import main_rank, main_comm
from parmepy.methods_keys import Support
@@ -365,9 +366,6 @@ class Redistribute(Operator):
return res
return res
- def finalize(self):
- pass
-
def addRedistributeRequirement(self, red):
raise ValueError(
"Cannot add a requirement to a Redistribute operator.")
diff --git a/HySoP/hysop/operator/redistribute_intercomm.py b/HySoP/hysop/operator/redistribute_intercomm.py
index 9b68569a044772aa8ec494acd30906b7bb1c76a9..959e4d1506b54427ec257275290691e356984892 100644
--- a/HySoP/hysop/operator/redistribute_intercomm.py
+++ b/HySoP/hysop/operator/redistribute_intercomm.py
@@ -346,5 +346,3 @@ class RedistributeIntercomm(Operator):
return res
return res
- def finalize(self):
- pass
diff --git a/HySoP/hysop/operator/reprojection.py b/HySoP/hysop/operator/reprojection.py
index eb3c89f9531af0aa06bdaa73c81e3f4e3d91d89a..a35bb739cf41f004624744fa136b1c80891d6a7e 100644
--- a/HySoP/hysop/operator/reprojection.py
+++ b/HySoP/hysop/operator/reprojection.py
@@ -4,10 +4,10 @@
Compute reprojection criterion and divergence maximum
"""
import numpy as np
-from parmepy.constants import np, debug, PARMES_REAL, PARMES_MPI_REAL
+from parmepy.constants import debug, PARMES_MPI_REAL
from parmepy.methods_keys import SpaceDiscretisation
from parmepy.operator.monitors.monitoring import Monitoring
-from parmepy.numerics.finite_differences import FD_C_4, FD_C_2
+from parmepy.numerics.finite_differences import FD_C_4
from parmepy.numerics.differential_operations import GradV
import parmepy.tools.numpywrappers as npw
from parmepy.numerics.updateGhosts import UpdateGhosts
@@ -21,150 +21,162 @@ class Reprojection_criterion(Monitoring):
See the related PDF called "vorticity_solenoidal_projection.pdf"
in ParmesDoc for more details.
"""
- def __init__(self, vorticity,
- need_reproj,
- reproj_cst,
- topo, frequency,
- method={SpaceDiscretisation: FD_C_4},
- prefix=None, task_id=None):
+ def __init__(self, vorticity, reproj_cst, reprojRate,
+ topo, checkCriterion=False, io_params=None,
+ method=None, task_id=None):
"""
Constructor.
@param vorticity field
@param method : finite difference scheme
@param topo : the topology on which we want to monitor the fields
- @param frequency : output file producing frequency.
@param prefix : output file name.
+ @param checkCriterion : set behavior of this operator :
+ if True, compute some criterion and force frequency
+ to 1 if reprojection is needed.
+ else (false) performs reprojection every reprojRate iterations.
+ @param reprojRate : set rate of execution of the reprojection
+ @param io_params : parameters (dict) to set file output.
+ If None, no output. Set io_params = {} if you want output,
+ with default parameters values. Default file name = 'reproj'
+ See parmepy.tools.io_utils.Writer for details
"""
- Monitoring.__init__(self, [vorticity], topo, frequency,
+ if io_params is not None:
+ if not "filename" in io_params:
+ io_params["filename"] = "reproj"
+ # Set output buffer shape
+ io_params["writebuffshape"] = (1, 5)
+
+ Monitoring.__init__(self, [vorticity], topo, io_params,
task_id=task_id)
- if prefix is None:
- self.prefix = './res/reproj_crit.dat'
- else:
- self.prefix = prefix
- ## vorticity field
- self.vorticity = vorticity
- ## Numerical methods for space discretization
- assert SpaceDiscretisation in method.keys()
- self.method = method[SpaceDiscretisation]
-
- self.input = [vorticity]
- self.output = []
# alias for local topo
self._topo = self._predefinedTopo
# \todo : rewrite for multiresolution case.
# Note FP : for multiresolution case, it would probably be
# better to use two different operators for energy and enstrophy.
-
- # Boolean : is the reprojection of vorticity needed for
- # current timestep ?
- self.need_reproj = need_reproj
-
+ ## Frequency for reprojection
+ self.reprojRate = reprojRate
+ ## The initial value will be used as default during
+ # simulation
+ self.default_rate = reprojRate
+ ## Set behavior of this operator :
+ ## if True, compute some criterion and force frequency
+ ## to 1 if reprojection is needed.
+ ## else (false) performs reprojection every frequency iterations.
+ self.checkCriterion = checkCriterion
# constant defining the reprojection criterion :
- # if the latter is greater than this constant, then a reprojection is needed
+ # if the latter is greater than this constant, then a reprojection
+ # is needed
self.reproj_cst = reproj_cst
- self.default_reproj_frequency = need_reproj.data[1]
self.proj_counter = 0
+ ## local buffer
+ ## diag = [time, projCriterion, d1, d2, proj_counter]
+ ## See apply for details about d1, d2.
+ self.diagnostics = npw.zeros(5)
+ # Connect writer buffer to diagnostics, if required
+ if self._writer is not None:
+ self._writer.buffer = self.diagnostics.reshape(1,5)
+
+ if self.checkCriterion:
+ ## vorticity field
+ self.vorticity = vorticity
+ if method is None:
+ method = {SpaceDiscretisation: FD_C_4}
+ ## Numerical methods for space discretization
+ assert SpaceDiscretisation in method.keys()
+ self.method = method[SpaceDiscretisation]
+
+ self.input = [vorticity]
+ else:
+ self.input = []
+
+ self.output = []
+
def setUp(self):
- Monitoring.setUp(self)
- # Get discrete fields for vorticity.
- # Note FP : two options to get the discrete fields:
- # - 'field.discretization(topo)' that just try
- # to get the discrete field and return an error
- # if it does not exist.
- # - 'field.discretizet(topo) try to get the discrete field and create
- # a new discretization if it does not exist.
- # Current choice : no creation.
- self.discreteFields[self.vorticity] =\
- self.vorticity.discretization(self._topo)
-
- # prepare ghost points synchro for vorticity
- self._synchronize = UpdateGhosts(self._topo,
- self.vorticity.nbComponents)
-
- # grad function
- self._function = GradV(self._topo,
- method=self.method)
-
- self.vort_d = self.discreteFields[self.vorticity]
-
- memshape = self.vort_d.data[0].shape
- worklength = self.vort_d.nbComponents ** 2
- # gradient result array
- self.grad = [npw.zeros(memshape)
- for i in xrange(worklength)]
+ if not self._isUpToDate and self.checkCriterion:
+ # Get discrete fields for vorticity.
+ # Note FP : two options to get the discrete fields:
+ # - 'field.discretization(topo)' that just try
+ # to get the discrete field and return an error
+ # if it does not exist.
+ # - 'field.discretizet(topo)
+ # try to get the discrete field and create
+ # a new discretization if it does not exist.
+ # Current choice : no creation.
+ self.discreteFields[self.vorticity] =\
+ self.vorticity.discretization(self._topo)
+
+ # prepare ghost points synchro for vorticity
+ self._synchronize = UpdateGhosts(self._topo,
+ self.vorticity.nbComponents)
+
+ # grad function
+ self._function = GradV(self._topo, method=self.method)
+
+ self.vort_d = self.discreteFields[self.vorticity]
+
+ memshape = self.vort_d.data[0].shape
+ worklength = self.vort_d.nbComponents ** 2
+ # gradient result array
+ self.grad = [npw.zeros(memshape) for i in xrange(worklength)]
self._isUpToDate = True
- # Open file for output
- if (self._topo.rank == 0):
- import os
- path = os.path.dirname(self.prefix)
- if not os.path.exists(path):
- os.mkdir(path)
-
- self.f = open(self.prefix, 'w')
-
@debug
@timed_function
def apply(self, simulation=None):
"""
Computes and prints reprojection criterion and divergence maximum
"""
- t = simulation.time
+ self.diagnostics[0] = simulation.time
ite = simulation.currentIteration
- self.need_reproj.data[1] = self.default_reproj_frequency
-
- # Synchronize ghost points of vorticity
- self._synchronize(self.vort_d.data)
- # gradU computation
- self.grad = self._function(self.vort_d.data, self.grad)
-
- diagnostics = npw.zeros((2))
- nbComponents = self.vorticity.nbComponents
+ # Reset reprojection frequency to default
+ self.reprojRate = self.default_rate
+
+ # If required, computation of a criterion
+ # and reset frequency
+ if self.checkCriterion:
+ # Synchronize ghost points of vorticity
+ self._synchronize(self.vort_d.data)
+ # gradU computation
+ self.grad = self._function(self.vort_d.data, self.grad)
+ nbComponents = self.vorticity.nbComponents
+ # maxima of vorticity divergence (abs)
+ self.diagnostics[2] = \
+ np.max(abs(sum([(self.grad[(nbComponents + 1) * i])
+ for i in xrange(nbComponents)])))
+
+ # maxima of partial derivatives of vorticity
+ for grad_n in self.grad:
+ self.diagnostics[3] = max(self.diagnostics[3],
+ np.max(abs(grad_n)))
+
+ # computation of the reprojection criterion and reduction
+ # among all proc.
+ projCriterion = self.diagnostics[2] / self.diagnostics[3]
+ projCriterion = \
+ self._topo.comm.allreduce(projCriterion, PARMES_MPI_REAL,
+ op=MPI.MAX)
+ self.diagnostics[1] = projCriterion
- # maxima of vorticity divergence (abs)
- diagnostics[0] = np.max(abs(sum([(self.grad[(nbComponents + 1) * i])
- for i in xrange(nbComponents)])))
+ # is reprojection of vorticity needed for the future time step ?
+ if (projCriterion > self.reproj_cst):
+ self.reprojRate = 1
- # maxima of partial derivatives of vorticity
- for direc1 in xrange(nbComponents):
- for direc2 in xrange(nbComponents):
- diagnostics[1] = max(diagnostics[1],
- np.max(abs(self.grad[nbComponents * \
- direc1 + direc2])))
+ # Note FP : is counter really useful? Maybe it should be increased
+ # by the operator that use projection (poisson)?
+ if (ite % self.reprojRate == 0):
+ self.proj_counter += 1
+ self.diagnostics[4] = self.proj_counter
- # computation of the reprojection criterion and reduction among all proc.
- projCriterion = diagnostics[0] / diagnostics[1]
- projCriterion = \
- self._topo.comm.allreduce(projCriterion,
- PARMES_MPI_REAL,
- op=MPI.MAX)
+ # Print results, if required
+ # Remark : writer buffer is (pointer) connected to diagnostics
+ if self._writer is not None and self._writer.doWrite(ite):
+ self._writer.write()
- if self.need_reproj.data[0]:
- # is reprojection of vorticity needed for the future time step ?
- if ((projCriterion > self.reproj_cst) and (self.need_reproj.data[2])):
- self.need_reproj.data[1] = 1
-
- if (ite % self.need_reproj.data[1] == 0):
- self.proj_counter += 1
-
- # printer
- if ((self._topo.comm.Get_rank() == 0) and (ite % self.freq == 0)):
- self.f = open(self.prefix, 'a')
- self.f.write("%s %s %s %s %s\n" % (t, projCriterion,
- diagnostics[0],
- diagnostics[1],
- self.proj_counter))
- self.f.close()
-
-
- def __str__(self):
- s = "Reprojection_criterion. "
- return s
-
-if __name__ == "__main__":
- print __doc__
- print "- Provided class : Reprojection_criterion"
- print Reprojection_criterion.__doc__
+ def doProjection(self, ite):
+ """
+ True if projection must be done
+ """
+ return ite % self.reprojRate == 0
diff --git a/HySoP/hysop/operator/tests/test_Stretching.py b/HySoP/hysop/operator/tests/test_Stretching.py
index 5af3da16e2d0b21d711ce940d88ca53abfa6abd6..1be61b32017611c2081f13a6eca6ee71e3f668b4 100755
--- a/HySoP/hysop/operator/tests/test_Stretching.py
+++ b/HySoP/hysop/operator/tests/test_Stretching.py
@@ -84,5 +84,5 @@ def test_Stretching():
velo.initialize()
vorti.initialize()
stretch.setUp()
- simulation = Simulation(0, timeStep, 20)
+ simulation = Simulation(0, 20, timeStep)
stretch.apply(simulation)
diff --git a/HySoP/hysop/operator/tests/test_advec_scales.py b/HySoP/hysop/operator/tests/test_advec_scales.py
index 87087a47206d9c31ea029d26f690b481a4743e2c..3001b05297a0d7b75dad5280b761f1847f54504c 100755
--- a/HySoP/hysop/operator/tests/test_advec_scales.py
+++ b/HySoP/hysop/operator/tests/test_advec_scales.py
@@ -48,8 +48,8 @@ def test_nullVelocity_m4():
dtype=PARMES_REAL, order=ORDER)
scal_ref_d.data[0][...] = scal_d.data[0][...]
- advec.apply(Simulation(0., 0.1, 1))
- advec_py.apply(Simulation(0., 0.1, 1))
+ advec.apply(Simulation(0., 1., 0.1))
+ advec_py.apply(Simulation(0., 1, 0.1))
assert np.allclose(scal_ref_d.data[0], scal_d.data[0])
@@ -97,8 +97,8 @@ def test_nullVelocity_vec_m4():
scal_ref_d.data[1][...] = scal_d.data[1][...]
scal_ref_d.data[2][...] = scal_d.data[2][...]
- advec.apply(Simulation(0., 0.1, 1))
- advec_py.apply(Simulation(0., 0.1, 1))
+ advec.apply(Simulation(0., 1., 0.1))
+ advec_py.apply(Simulation(0., 1., 0.1))
assert np.allclose(scal_ref_d.data[0], scal_d.data[0])
assert np.allclose(scal_ref_d.data[1], scal_d.data[1])
@@ -125,7 +125,7 @@ def test_nullVelocity_m6():
dtype=PARMES_REAL, order=ORDER)
scal_init = npw.copy(scal_d.data[0])
- advec.apply(Simulation(0., 0.1, 1))
+ advec.apply(Simulation(0., 1., 0.1))
assert np.allclose(scal_init, scal_d.data[0])
@@ -159,7 +159,7 @@ def test_nullVelocity_vec_m6():
scal_init1 = npw.copy(scal_d.data[1])
scal_init2 = npw.copy(scal_d.data[2])
- advec.apply(Simulation(0., 0.1, 1))
+ advec.apply(Simulation(0., 1., 0.1))
assert np.allclose(scal_init0, scal_d.data[0])
assert np.allclose(scal_init1, scal_d.data[1])
@@ -201,8 +201,8 @@ def test_nullVelocity_m8():
dtype=PARMES_REAL, order=ORDER)
scal_ref_d.data[0][...] = scal_d.data[0]
- advec.apply(Simulation(0., 0.1, 1))
- advec_py.apply(Simulation(0., 0.1, 1))
+ advec.apply(Simulation(0., 1., 0.1))
+ advec_py.apply(Simulation(0., 1., 0.1))
assert np.allclose(scal_ref_d.data[0], scal_d.data[0])
@@ -250,8 +250,8 @@ def test_nullVelocity_vec_m8():
scal_ref_d.data[1][...] = scal_d.data[1][...]
scal_ref_d.data[2][...] = scal_d.data[2][...]
- advec.apply(Simulation(0., 0.075, 1))
- advec_py.apply(Simulation(0., 0.075, 1))
+ advec.apply(Simulation(0., 1., 0.075))
+ advec_py.apply(Simulation(0., 1., 0.075))
assert np.allclose(scal_ref_d.data[0], scal_d.data[0])
assert np.allclose(scal_ref_d.data[1], scal_d.data[1])
@@ -301,8 +301,8 @@ def _randomVelocity_m4():
np.random.random(scal_d.data[0].shape),
dtype=PARMES_REAL, order=ORDER) / (2. * scal_d.resolution[1])
- advec.apply(Simulation(0., 0.1, 1))
- advec_py.apply(Simulation(0., 0.1, 1))
+ advec.apply(Simulation(0., 1., 0.1,))
+ advec_py.apply(Simulation(0., 1., 0.1))
assert np.allclose(scal_ref_d.data[0], scal_d.data[0])
@@ -358,8 +358,8 @@ def _randomVelocity_vec_m4():
np.random.random(scal_d.data[2].shape),
dtype=PARMES_REAL, order=ORDER) / (2. * scal_d.resolution[1])
- advec.apply(Simulation(0., 0.1, 1))
- advec_py.apply(Simulation(0., 0.1, 1))
+ advec.apply(Simulation(0., 1., 0.1))
+ advec_py.apply(Simulation(0., 1., 0.1))
assert np.allclose(scal_ref_d.data[0], scal_d.data[0])
assert np.allclose(scal_ref_d.data[1], scal_d.data[1])
@@ -409,8 +409,8 @@ def test_randomVelocity_m6():
np.random.random(scal_d.data[0].shape),
dtype=PARMES_REAL, order=ORDER) / (2. * scal_d.resolution[1])
- advec.apply(Simulation(0., 0.1, 1))
- advec_py.apply(Simulation(0., 0.1, 1))
+ advec.apply(Simulation(0., 1., 0.1))
+ advec_py.apply(Simulation(0., 1., 0.1))
assert np.allclose(scal_ref_d.data[0], scal_d.data[0])
@@ -466,8 +466,8 @@ def test_randomVelocity_vec_m6():
np.random.random(scal_d.data[2].shape),
dtype=PARMES_REAL, order=ORDER) / (2. * scal_d.resolution[1])
- advec.apply(Simulation(0., 0.1, 1))
- advec_py.apply(Simulation(0., 0.1, 1))
+ advec.apply(Simulation(0., 1., 0.1))
+ advec_py.apply(Simulation(0., 1., 0.1))
assert np.allclose(scal_ref_d.data[0], scal_d.data[0])
assert np.allclose(scal_ref_d.data[1], scal_d.data[1])
@@ -517,8 +517,8 @@ def test_randomVelocity_m8():
np.random.random(scal_d.data[0].shape),
dtype=PARMES_REAL, order=ORDER) / (2. * scal_d.resolution[1])
- advec.apply(Simulation(0., 0.1, 1))
- advec_py.apply(Simulation(0., 0.1, 1))
+ advec.apply(Simulation(0., 1., 0.1))
+ advec_py.apply(Simulation(0., 1., 0.1))
assert np.allclose(scal_ref_d.data[0], scal_d.data[0], atol=1e-07)
@@ -574,8 +574,8 @@ def test_randomVelocity_vec_m8():
np.random.random(scal_d.data[2].shape),
dtype=PARMES_REAL, order=ORDER) / (2. * scal_d.resolution[1])
- advec.apply(Simulation(0., 0.01, 1))
- advec_py.apply(Simulation(0., 0.01, 1))
+ advec.apply(Simulation(0., 1., 0.01))
+ advec_py.apply(Simulation(0., 1., 0.01))
assert np.allclose(scal_ref_d.data[0], scal_d.data[0], atol=1e-07)
assert np.allclose(scal_ref_d.data[1], scal_d.data[1], atol=1e-07)
diff --git a/HySoP/hysop/operator/tests/test_analytic.py b/HySoP/hysop/operator/tests/test_analytic.py
index b71cacc92009eab9468959d900690b4c59d4ef96..2d6240b713302099b197ad6d311b8127d0dc78c1 100644
--- a/HySoP/hysop/operator/tests/test_analytic.py
+++ b/HySoP/hysop/operator/tests/test_analytic.py
@@ -70,7 +70,7 @@ res2D = [33, 33]
L2D = [1., 1.]
origin2D = [0., 0.]
nbc = 4
-simu = Simulation(10., 0.1, 1)
+simu = Simulation(0., 1., 0.1)
## I - Tests with direct calls of field.initialize
diff --git a/HySoP/hysop/operator/tests/test_particle_advection.py b/HySoP/hysop/operator/tests/test_particle_advection.py
index 2de9116329aea1305bfedeb6c3d870e218b99b00..d4b604ca883d2a780b69691f4388f6c692ebe019 100644
--- a/HySoP/hysop/operator/tests/test_particle_advection.py
+++ b/HySoP/hysop/operator/tests/test_particle_advection.py
@@ -70,7 +70,7 @@ def assertion(scal, advec):
dtype=PARMES_REAL, order=ORDER)
scal_init = npw.copy(scal_d.data[0])
- advec.apply(Simulation(0., 0.01, 1))
+ advec.apply(Simulation(0., 1., 0.01))
return np.allclose(scal_init, scal_d.data[0])
@@ -86,7 +86,7 @@ def assertion_vector2D(scal, advec):
scal1_init = npw.copy(scal_d.data[0])
scal2_init = npw.copy(scal_d.data[1])
- advec.apply(Simulation(0., 0.01, 1))
+ advec.apply(Simulation(0., 1., 0.01))
return np.allclose(scal1_init, scal_d.data[0]) and \
np.allclose(scal2_init, scal_d.data[1])
@@ -106,7 +106,7 @@ def assertion_vector3D(scal, advec):
scal2_init = npw.copy(scal_d.data[1])
scal3_init = npw.copy(scal_d.data[2])
- advec.apply(Simulation(0., 0.01, 1))
+ advec.apply(Simulation(0., 1., 0.01))
return np.allclose(scal1_init, scal_d.data[0]) and \
np.allclose(scal2_init, scal_d.data[1]) and \
np.allclose(scal3_init, scal_d.data[2])
@@ -125,7 +125,7 @@ def assertion_list(scal, advec):
scal1_init = npw.copy(scal1_d.data[0])
scal2_init = npw.copy(scal2_d.data[0])
- advec.apply(Simulation(0., 0.01, 1))
+ advec.apply(Simulation(0., 1., 0.01))
print scal1_init, scal1_d.data[0]
print scal2_init, scal2_d.data[0]
return np.allclose(scal1_init, scal1_d.data[0]) and \
diff --git a/HySoP/hysop/operator/tests/test_velocity_correction.py b/HySoP/hysop/operator/tests/test_velocity_correction.py
new file mode 100755
index 0000000000000000000000000000000000000000..06c2ff23d91ae497816d208ae9a53d6264b16ca7
--- /dev/null
+++ b/HySoP/hysop/operator/tests/test_velocity_correction.py
@@ -0,0 +1,148 @@
+# -*- coding: utf-8 -*-
+
+import parmepy as pp
+from parmepy.operator.velocity_correction import VelocityCorrection
+from parmepy.problem.simulation import Simulation
+import numpy as np
+from parmepy.mpi.topology import Cartesian
+import math
+pi = math.pi
+sin = np.sin
+cos = np.cos
+import parmepy.tools.numpywrappers as npw
+
+## Physical Domain description
+print " ========= Start Navier-Stokes 3D (Flow past bluff bodies) ========="
+
+## pi constant
+pi = np.pi
+cos = np.cos
+sin = np.sin
+# Upstream flow velocity
+uinf = 1.0
+tol = 1e-12
+
+
+## Function to compute TG velocity
+def computeVel(res, x, y, z, t):
+ res[0][...] = sin(x) * cos(y) * cos(z)
+ res[1][...] = - cos(x) * sin(y) * cos(z)
+ res[2][...] = 0.
+ return res
+
+
+## Function to compute reference vorticity
+def computeVort(res, x, y, z, t):
+ res[0][...] = - cos(x) * sin(y) * sin(z)
+ res[1][...] = - sin(x) * cos(y) * sin(z)
+ res[2][...] = 2. * sin(x) * sin(y) * cos(z)
+ return res
+
+
+## Function to compute TG velocity
+def computeVel2D(res, x, y, t):
+ res[0][...] = sin(x) * cos(y)
+ res[1][...] = - cos(x) * sin(y)
+ return res
+
+
+## Function to compute reference vorticity
+def computeVort2D(res, x, y, t):
+ res[0][...] = - cos(x) * sin(y)
+ return res
+
+
+def test_velocity_correction_3D():
+ dim = 3
+ ## Domain
+ boxlength = npw.realarray([2.0 * pi, pi, pi])
+ boxorigin = npw.realarray([0., 0., 0.])
+ box = pp.Box(dim, length=boxlength, origin=boxorigin)
+
+ ## Global resolution
+ nbElem = [33, 33, 33]
+
+ ## Vector Fields
+ velo = pp.Field(domain=box, formula=computeVel,
+ name='Velocity', isVector=True)
+ vorti = pp.Field(domain=box, formula=computeVort,
+ name='Vorticity', isVector=True)
+
+ ## Usual Cartesian topology definition
+ NBGHOSTS = 2
+ ghosts = np.ones((box.dimension)) * NBGHOSTS
+ topo = Cartesian(box, box.dimension, nbElem,
+ ghosts=ghosts)
+
+ op = {}
+ ref_rate = uinf * box.length[1] * box.length[2]
+ op['correction'] = VelocityCorrection(velo, vorti,
+ resolutions={velo: nbElem,
+ vorti: nbElem},
+ req_flowrate=ref_rate, topo=topo)
+ op['correction'].discretize()
+ op['correction'].setUp()
+ # === Simulation setup ===
+ simu = Simulation(tinit=0.0, tend=5., timeStep=0.005, iterMax=1000000)
+ # init fields
+ velo.initialize(topo=topo)
+ vorti.initialize(topo=topo)
+
+ # Apply correction
+ op['correction'].apply(simu)
+ # check new flowrate values
+ sref = op['correction'].discreteOperator.surfRef
+ flowrate = velo.integrateOnSurface(sref, topo)
+ print flowrate
+ assert (flowrate - ref_rate) < tol
+ for i in xrange(1, dim):
+ flowrate2 = velo.integrateOnSurface(sref, topo, component=i)
+ assert flowrate2 < tol
+
+
+def test_velocity_correction_2D():
+ dim = 2
+ ## Domain
+ boxlength = npw.realarray([2.0 * pi, pi])
+ boxorigin = npw.realarray([0., 0.])
+ box = pp.Box(dim, length=boxlength, origin=boxorigin)
+
+ ## Global resolution
+ nbElem = [33, 33]
+
+ ## Vector Fields
+ velo = pp.Field(domain=box, formula=computeVel2D,
+ name='Velocity', isVector=True)
+ vorti = pp.Field(domain=box, formula=computeVort2D,
+ name='Vorticity', isVector=False)
+
+ ## Usual Cartesian topology definition
+ NBGHOSTS = 2
+ ghosts = np.ones((box.dimension)) * NBGHOSTS
+ topo = Cartesian(box, box.dimension, nbElem,
+ ghosts=ghosts)
+
+ op = {}
+ ref_rate = uinf * box.length[1]
+ op['correction'] = VelocityCorrection(velo, vorti,
+ resolutions={velo: nbElem,
+ vorti: nbElem},
+ req_flowrate=ref_rate, topo=topo)
+ op['correction'].discretize()
+ op['correction'].setUp()
+ # === Simulation setup ===
+ simu = Simulation(tinit=0.0, tend=5., timeStep=0.005, iterMax=1000000)
+ # init fields
+ velo.initialize(topo=topo)
+ vorti.initialize(topo=topo)
+
+ # Apply correction
+ op['correction'].apply(simu)
+ # check new flowrate values
+ sref = op['correction'].discreteOperator.surfRef
+ flowrate = velo.integrateOnSurface(sref, topo)
+ print flowrate
+ assert (flowrate - ref_rate) < tol
+ for i in xrange(1, dim):
+ flowrate2 = velo.integrateOnSurface(sref, topo, component=i)
+ assert flowrate2 < tol
diff --git a/HySoP/hysop/operator/velocity_correction.py b/HySoP/hysop/operator/velocity_correction.py
index 72de985cc96d0aa4d0bc7a884edd098f1c62a8aa..0735875e5bbfdcb848ce44e7abffbb283f324530 100755
--- a/HySoP/hysop/operator/velocity_correction.py
+++ b/HySoP/hysop/operator/velocity_correction.py
@@ -2,13 +2,13 @@
"""
@file operator/velocity_correction.py
-Correction of the velocity field according to the desired flowrate.
+Operator to shift velocity to fit with a required input flowrate.
"""
from parmepy.constants import debug
from parmepy.operator.discrete.velocity_correction import VelocityCorrection_D
from parmepy.operator.continuous import Operator
-import numpy as np
+from parmepy.domain.obstacle.controlBox import ControlBox
class VelocityCorrection(Operator):
@@ -20,20 +20,24 @@ class VelocityCorrection(Operator):
"""
@debug
- def __init__(self, velocity, vorticity, resolutions,
- topo=None, task_id=None, comm=None, **other_config):
+ def __init__(self, velocity, vorticity, resolutions, req_flowrate,
+ surf=None, topo=None, comm=None, task_id=None):
"""
Corrects the values of the velocity field after
solving Poisson equation in order to prescribe proper
mean flow and ensure the desired inlet flowrate.
- @param velocity field
- @param resolutions : grid resolution of velocity
- @param topo : a predefined topology to discretize velocity
+ @param[in, out] velocity field to be corrected
+ @param[in] vorticity field used to compute correction
+ @param resolutions : grid resolutions of velocity and vorticity
+ @param[in] req_flowrate : required value for the flowrate
+ @param[in] surf : surface (parmepy.domain.obstacle.planes.SubPlane)
+ used to compute reference flow rates. Default = surface at x_origin,
+ normal to x-dir.
+ @param topo : a predefined topology to discretize velocity/vorticity
"""
Operator.__init__(self, [velocity, vorticity], topo=topo,
task_id=task_id, comm=comm)
- self.config = other_config
## velocity variable (vector)
self.velocity = velocity
## velocity variable
@@ -43,6 +47,13 @@ class VelocityCorrection(Operator):
self.input = [self.velocity, self.vorticity]
self.output = [self.velocity]
+ ## A 'reference' surface on which flow rates will be computed
+ ## (usually, surface normal to the flow at origin)
+ self.surfRef = surf
+ ## Expected value for the flow rate through self.surfRef
+ self.req_flowrate = req_flowrate
+ dom = self.velocity.domain
+ self.cb = ControlBox(dom, dom.origin, dom.length)
def discretize(self):
if self._predefinedTopo is not None:
@@ -62,7 +73,20 @@ class VelocityCorrection(Operator):
self.discreteOperator =\
VelocityCorrection_D(self.discreteFields[self.velocity],
self.discreteFields[self.vorticity],
- **self.config)
+ self.req_flowrate, self.cb)
self.discreteOperator.setUp()
self._isUpToDate = True
+
+ def apply(self, simulation=None):
+ """
+ Compute correction and add it to current velocoty
+ """
+ self.discreteOperator.apply(simulation)
+
+ def computeCorrection(self):
+ """
+ Compute the required correction for the current state
+ but do not apply it onto velocity.
+ """
+ self.discreteOperator.computeCorrection()
diff --git a/HySoP/hysop/problem/problem.py b/HySoP/hysop/problem/problem.py
index 37227c62914bdc5349843e5238928b83af4f01d3..5483c7dd1b3e21a3f63c2cf598fea72ddfce7e37 100644
--- a/HySoP/hysop/problem/problem.py
+++ b/HySoP/hysop/problem/problem.py
@@ -181,7 +181,8 @@ class Problem(object):
if main_rank == 0:
print "\n\n Start solving ..."
while not self.simulation.isOver:
- self.simulation.printState()
+ if main_rank == 0:
+ self.simulation.printState()
for op in self.operators:
if __VERBOSE__:
diff --git a/HySoP/hysop/problem/problem_tasks.py b/HySoP/hysop/problem/problem_tasks.py
index 9749aa4e178c610078e60e3d7a1e0edc46310943..3ec307560d1e3f284811ba7e49b6d017f4a579ab 100644
--- a/HySoP/hysop/problem/problem_tasks.py
+++ b/HySoP/hysop/problem/problem_tasks.py
@@ -178,7 +178,8 @@ class ProblemTasks(Problem):
if self._main_rank == 0:
print "\n\n Start solving ..."
while not self.simulation.isOver:
- self.simulation.printState()
+ if self._main_rank == 0:
+ self.simulation.printState()
for op in self.operators:
op.apply(self.simulation)
testdump = \
diff --git a/HySoP/hysop/problem/simulation.py b/HySoP/hysop/problem/simulation.py
index df6b3850fc7f318318216277ca2812e49b625dc8..fb16b807648bede2ff8fdd1171c8c090dfa6cbe0 100644
--- a/HySoP/hysop/problem/simulation.py
+++ b/HySoP/hysop/problem/simulation.py
@@ -3,8 +3,6 @@
Description of the simulation parameters (time, iteration ...)
"""
-from parmepy.mpi import main_rank
-from parmepy.variables.variables import Variables
class Simulation(object):
@@ -29,10 +27,13 @@ class Simulation(object):
@param bridgeOp : Redistribute operator
@param computeTimeStepOp : Operator which evaluates
the adaptative time step according to the flow fields
+
+ Notation:
+ iteration number 'currentIteration'
+ between tk and tkp1 = tk + timeStep
"""
## Simulation final time
self.end = tend
- print 'final time:', tend, self.end, iterMax
## Starting time
self.start = tinit
## Simulation current time
@@ -41,25 +42,21 @@ class Simulation(object):
self.isOver = False
## Iteration counter
self.currentIteration = 0
- ## Simulation time step variable
- if (isinstance(timeStep, Variables)):
- ## Set isAdaptative to True
- self.isAdaptative = True
- ## Simulation time step variable
- self.timeStep_variable = timeStep
- ## Simulation time step adaptative value
- self.timeStep = timeStep.data[0]
- ## Save the initial timestep value
- self.timeStepInit = timeStep.data[0]
- else:
- ## Set isAdaptative to False
- self.isAdaptative = False
- ## Simulation time step fixed value
- self.timeStep = timeStep
-
+ ## Simulation time step
+ self.timeStep = timeStep
+ ## Starting time for the current iteration
+ self.tk = tinit
+ ## tk + dt
+ self.tkp1 = tinit + self.timeStep
## Maximum number of iterations
self.iterMax = iterMax
self._lastStep = False
+ print self.start, self.timeStep, self.end
+ print self.start + self.timeStep
+ assert self.end > self.start, \
+ 'Final time must be greater than initial time'
+ assert (self.start + self.timeStep) <= self.end,\
+ 'start + step is bigger than end.'
def init(self):
self.currentIteration = 1
@@ -71,6 +68,8 @@ class Simulation(object):
Compute a timestep for the incoming iteration (from a vairable and
to reach the end of simulation)
"""
+ # Increment iteration counter
+ self.currentIteration += 1
if self._lastStep:
# The timestep was adjusted to reach end in the previous call
# So now the simulation is over
@@ -78,42 +77,33 @@ class Simulation(object):
else:
if self.currentIteration < self.iterMax:
# Advance time for the iteration just ended
- self.time += self.timeStep
- self.currentIteration += 1
-
- # Prepare the timestep for the following itertation
- if self.isAdaptative and self.currentIteration >= 2:
- self.timeStep = self.timeStep_variable.data[0]
+ self.tk = self.tkp1
+ self.tkp1 = self.tk + self.timeStep
- # Axdjust last timestep to reach self.end
- if self.end - self.time < self.timeStep:
- self.timeStep = self.end - self.time
+ # Adjust last timestep to reach self.end
+ if self.tkp1 > self.end:
+ self.timeStep = self.end - self.tk
+ self.tkp1 = self.end
self._lastStep = True
else:
# iteration number is reached
self.isOver = True
- # Time is close enough to the end to consider that the simulation
- # is over.
- if abs(self.time - self.end) < 1e-10:
- self.isOver = True
- def printState(self):
+ self.time = self.tkp1
+
+ def updateTimeStep(self, newDt):
"""
- print current state
+ Update current time step.
+ This function is usually called from Adapt_timestep operator.
"""
- if main_rank == 0:
- print "=== Iteration : {0:3d}, start from t ={1:6.3f} \
-to t ={2:6.3f} ===".format(
- self.currentIteration, self.time, self.time + self.timeStep)
+ self.timeStep = newDt
- def reset(self):
+ def printState(self):
"""
- set initial time to current and reset iteration counter.
- Used to initialize simulation after a restart.
+ print current state
"""
- self.start = self.time
- self.isOver = False
- self.currentIteration = 1
+ msg = "== Iteration : {0:3d}, from t = {1:6.5} to t = {2:6.5f} =="
+ print msg.format(self.currentIteration, self.tk, self.time)
def __str__(self):
s = "Simulation parameters : "
@@ -123,3 +113,8 @@ to t ={2:6.3f} ===".format(
s += str(self.currentIteration) + ', max number of iterations : '
s += str(self.iterMax)
return s
+
+ def initialize(self):
+ self.time = self.tkp1
+ self.isOver = False
+ self.currentIteration = 0
diff --git a/HySoP/hysop/tools/io_utils.py b/HySoP/hysop/tools/io_utils.py
new file mode 100644
index 0000000000000000000000000000000000000000..1cced51835f5256a3f195bac2f4e7bab0017e891
--- /dev/null
+++ b/HySoP/hysop/tools/io_utils.py
@@ -0,0 +1,181 @@
+"""
+@file io_utils.py
+Tools related to i/o in parmes.
+"""
+import os
+import scitools.filetable as ft
+import parmepy.tools.numpywrappers as npw
+import inspect
+import parmepy.mpi as mpi
+
+
+class io(object):
+ """
+ Static class with utilities for file i/o
+ """
+ @staticmethod
+ def defaultPath():
+ """
+ Get the default path for io.
+ Used name if set. If not, get basename of the
+ top exec. file.
+ Example : if you run python MyExe.py for your parmes simulation,
+ results will be saved in ./MyExe/p1/, ./MyExe/p4/
+ the number after 'p' being the number of mpi processus
+ set for the simulation.
+ """
+ # Memo for getouterframe usage:
+ # frame, filename, line_number, function_name, lines, index =\
+ # inspect.getouterframes(inspect.currentframe())[-1]
+ a = inspect.getouterframes(inspect.currentframe())[-1]
+ return os.path.join(a[1].split('.')[0], 'p' + str(mpi.main_size))
+
+ @staticmethod
+ def checkDir(filename, io_rank=0, comm=None):
+ """
+ Check if the directory of 'filename' exists.
+ Creates it if not.
+ @param filename : file name with full or relative path
+ @param io_rank : processus rank that does the check.
+ @param comm : the mpi communicator that does the check
+ """
+ # Create output dir if required
+ if comm is None:
+ comm = mpi.main_comm
+ if comm.Get_rank() == io_rank:
+ d = os.path.dirname(filename)
+ if not os.path.exists(d):
+ os.makedirs(d)
+
+
+class Writer(object):
+ """
+ @param params : dictionnary to set parameters
+ for output file (name, path, ...)
+ available keys are :
+ - filename : file name with full or relative path.
+ Default = data.out, in filepath.
+ If filename contains absolute path, filepath is ignored.
+ - filepath : path to the i/o file.
+ Default = parmepy.tools.io_utils.io.defaultPath.
+ - frequency : how often output file is written. Default = 1
+ (i.e. every iteration). N means writes every N iteration and
+ 0 means write only after last iteration.
+ - io_rank : mpi process rank (in comm) that performs writting. Default = 0
+ - writebuffshape : shape (numpy array) of the output buffer,
+ written in filename. Default = 1.
+ - safeOutput : boolean. Default = True.
+ True --> open/close file everytime data are written.
+ False --> open at init and close during finalize. Cost less but if simu
+ crashes, data are lost.
+ @param comm : mpi communicator that handles this
+ writer. Default = parmepy.mpi.main_comm.
+
+ Usage :
+ \code
+ >>> import parmepy.tools.io_utils as io
+ >>> defaultPath = io.io.defaultPath()
+ >>> params = {"filename": defaultPath + 'r.dat', "writebuffshape": (1, 2)}
+ >>> wr = Writer(params)
+ >>> ite = 3 # current iteration number
+ >>> if wr.doWrite(ite):
+ ... wr.buffer[...] = 3.
+ ... wr.write()
+ >>> wr.finalize()
+
+ \endcode
+ result : buffer is written into r.dat.
+ """
+ def __init__(self, params, comm=None):
+ """
+ """
+ # Directory for output
+ if "filepath" in params:
+ filepath = params["filepath"]
+ else:
+ filepath = io.defaultPath()
+ # file name
+ if "filename" in params:
+ filename = params["filename"]
+ else:
+ filename = 'data.out'
+ ## Absolute path + name for i/o file
+ ## Note that if filename contains absolute path
+ ## filepath is ignored
+ self.filename = os.path.join(filepath, filename)
+
+ ## How often i/o must occur
+ if "frequency" in params:
+ self.frequency = params["frequency"]
+ else:
+ self.frequency = 1
+ ## Rank of the mpi process that performs i/o (if any)
+ if "io_rank" in params:
+ self.io_rank = params["io_rank"]
+ else:
+ self.io_rank = 0
+
+ if comm is None:
+ comm = mpi.main_comm
+ ## A reference communicator, just to identify a
+ ## process rank for io.
+ self.comm = comm
+
+ # check if output dir exists, create it if not.
+ io.checkDir(self.filename, self.io_rank, self.comm)
+
+ ## Shape of the output buffer (must be a 2D numpy array)
+ if "writebuffshape" in params:
+ self.buffshape = params["writebuffshape"]
+ else:
+ self.buffshape = (1, 1)
+ assert len(self.buffshape) == 2,\
+ '2D shape required : set writebuffshape: (x,y)'
+ ## The buffer (numpy array) that will be printed to a file
+ self.buffer = npw.zeros(self.buffshape)
+
+ ## Defines how often
+ ## output file is written :
+ ## True --> open/close file everytime
+ ## data are written.
+ ## False --> open at init and close
+ ## during finalize. Cost less but if simu
+ ## crashes, data are lost.
+ if "safeOutput" in params:
+ self.safeOutput = params["safeOutput"]
+ else:
+ self.safeOutput = True
+
+ if self.safeOutput:
+ self.write = self._fullwrite
+ else:
+ self.write = self._partialwrite
+ # Force synchro to be sure that all output dirs
+ # have been created.
+ self.comm.barrier()
+ if self.comm.Get_rank() == self.io_rank:
+ self._file = open(self.filename, 'w')
+
+ def doWrite(self, ite):
+ """
+ True if output is required
+ for iteration ite
+ @param ite : current iteration number
+ """
+ rank = self.comm.Get_rank()
+ return rank == self.io_rank and (ite % self.frequency) == 0
+
+ def _fullwrite(self):
+ if self.comm.Get_rank() == self.io_rank:
+ self._file = open(self.filename, 'a')
+ ft.write(self._file, self.buffer)
+ self._file.close()
+
+ def _partialwrite(self):
+ if self.comm.Get_rank() == self.io_rank:
+ ft.write(self._file, self.buffer)
+
+ def finalize(self):
+ if self.comm.Get_rank() == self.io_rank:
+ if not self._file.closed:
+ self._file.close()
diff --git a/HySoP/hysop/tools/numpywrappers.py b/HySoP/hysop/tools/numpywrappers.py
index 11d8011ba74ed46a2fea3d0465e6739b7d8ccca2..333d6527dcca5e41ad1c6889b80f56c3488c7adc 100644
--- a/HySoP/hysop/tools/numpywrappers.py
+++ b/HySoP/hysop/tools/numpywrappers.py
@@ -6,6 +6,7 @@ Tools to build numpy arrays based on parmepy setup (float type ...)
"""
from parmepy.constants import PARMES_REAL, ORDER, PARMES_INTEGER, PARMES_INDEX
import numpy as np
+import scitools.filetable as ft
def zeros(shape, dtype=PARMES_REAL):
@@ -100,3 +101,11 @@ def prod(tab, dtype=PARMES_REAL):
"""
return np.prod(tab, dtype=dtype)
+
+def writeToFile(fname, data, mode='a'):
+ """
+ write data (numpy array) to file fname
+ """
+ fout = open(fname, mode)
+ ft.write(fout, data)
+ fout.close()
diff --git a/HySoP/hysop/variable_parameter.py b/HySoP/hysop/variable_parameter.py
new file mode 100644
index 0000000000000000000000000000000000000000..3cd0bea10fb771cc6585d87807265dbff3dcce84
--- /dev/null
+++ b/HySoP/hysop/variable_parameter.py
@@ -0,0 +1,75 @@
+#
+# Definition of the adaptative time step (for example) on the 3D domain:
+# \code
+# dt_adapt = Variable(dom3D, name="adaptative_timestep, data=[0.1]")
+# ...
+# \endcode
+#
+#
+"""
+@file variables/variables.py
+
+Tools to define parameters that do not depend on space but may
+be used as input-output parameter during simulation.
+
+Example : the time step.
+
+\code
+data = {'timestep': 0.1 }
+dt = VariableParameter(data)
+
+## An operator that can change dt
+op = Adaptative_timestep(..., dt)
+simu = Simulation(..., dt, ...)
+
+op.apply(simu) ===> change simu.time_step
+
+\endcode
+
+"""
+
+
+class VariableParameter(object):
+ """
+ Class to embed a user-defined parameter (a dictionnary indeed).
+
+ VariableParameter has a member data which is mutable.
+ """
+
+ def __init__(self, data=None, name=None):
+ """
+ Creates a dictionnary with data
+ @param data: the data used as parameter
+ @param name : optional name, used if data is not a dict.
+
+ If data is a dictionnary, self.data = data, that's it.
+ Else data will be added into self.data dict, with name as a key.
+ """
+ ## name of the variable
+ self.name = name
+ ## data values
+ if isinstance(data, dict):
+ self.data = data
+ else:
+ if name is None:
+ msg = "Name arg is required when data is not a dict."
+ raise AttributeError(msg)
+ self.data = {name: data}
+
+ def __getitem__(self, key):
+ """ Access to the content of the data[key]
+ @param key : requested key.
+ @return component 'key' of the dict
+ """
+ return self.data[key]
+
+ def __setitem__(self, key, value):
+ """
+ Set the data[key] = value
+ """
+ self.data[key] = value
+
+ def __str__(self):
+ s = "Variable parameter, which contains: "
+ s += str(self.data)
+ return s
diff --git a/HySoP/hysop/variables/__init__.py b/HySoP/hysop/variables/__init__.py
deleted file mode 100644
index c876f935ed708d186e1e8652c1583183e777da0c..0000000000000000000000000000000000000000
--- a/HySoP/hysop/variables/__init__.py
+++ /dev/null
@@ -1,27 +0,0 @@
-## @package parmepy.variables
-# Tools to build and use numpy arrays of 0_d_variables : all the types are allowed
-# to define a variable (real, integer, boolean ...)
-#
-# To define a variable, one must provide :
-# - the domain on which it is defined
-# and optionnaly:
-# - a user defined initialization of the numpy array containing the 0_d_subvariables
-# - a name
-#
-# First of all, one must have defined a domain :
-# \code
-# from parmepy.domain.box import Box
-# from parmepy.fields.continuous import Field
-#
-# # First, the domain : 3D and 2D boxes
-# dom3D = Box()
-# dom2D = Box(dimension=2, length=[1.0, 1.0], origin=[0., 0.])
-#\endcode
-#
-# Definition of the adaptative time step (for example) on the 3D domain:
-# \code
-# dt_adapt = Variable(dom3D, name="adaptative_timestep, data=[0.1]")
-# ...
-# \endcode
-#
-#
diff --git a/HySoP/hysop/variables/variables.py b/HySoP/hysop/variables/variables.py
deleted file mode 100644
index 5163d63e746cf10b1ddb7f51a0dd3e5b78456497..0000000000000000000000000000000000000000
--- a/HySoP/hysop/variables/variables.py
+++ /dev/null
@@ -1,47 +0,0 @@
-"""
-@file variables/variables.py
-
-0_d_Variables description.
-"""
-from parmepy.constants import debug
-#from parmepy.fields.vector import VectorField
-from parmepy.mpi import main_rank
-from parmepy.constants import np
-
-
-class Variables(object):
- """
- 0_d variable defined on a physical domain.
- """
-
- @debug
- def __new__(cls, *args, **kw):
- return object.__new__(cls, *args, **kw)
-
- @debug
- def __init__(self, domain, name=None, data=None):
- """
- Create a scalar of vector field which dimension is determined
- by its domain of definition.
-
- @param domain : physical domain where this field is defined.
- @param name : name of the variable.
- @param data : list of data subvalues of the variable.
-
- """
- ## Physical domain of definition of the variable.
- self.domain = domain
- ## name of the variable
- if(name is not None):
- self.name = name
- else:
- self.name = 'unamed'
- ## data values
- self.data = np.asarray(data)
-
- def __str__(self):
- """ Variable info display """
-
- s = '[' + str(main_rank) + '] " ' + self.name + ' ", '
- s += ' 0D variable with ' + str(len(self.data)) + 'sub_values '
- return s
diff --git a/HySoP/setup.py.in b/HySoP/setup.py.in
index eef3c0bc9bea2c24d2afd24ec7e2822029e87a13..3c75c97226a0edb1984fd3a8d8261a6be509fcda 100644
--- a/HySoP/setup.py.in
+++ b/HySoP/setup.py.in
@@ -16,7 +16,6 @@ packages = ['parmepy',
'parmepy.domain',
'parmepy.domain.obstacle',
'parmepy.fields',
- 'parmepy.variables',
'parmepy.operator',
'parmepy.operator.discrete',
'parmepy.operator.monitors',
diff --git a/HySoP/src/CMakeLists.txt b/HySoP/src/CMakeLists.txt
index 88df76a15bd81f3bab2bd5004b17e1dbcaab1dd1..c03b48f1d21157d920cf0b3356bb1ad10bee29a4 100644
--- a/HySoP/src/CMakeLists.txt
+++ b/HySoP/src/CMakeLists.txt
@@ -23,15 +23,18 @@ if(WITH_FFTW)
)
endif()
+#set(SCALES_DIR scalesReduced)
+set(SCALES_DIR scalesInterface)
+
# --- scales ---
if(WITH_SCALES)
set(${PARMES_LIBRARY_NAME}_SRCDIRS
${${PARMES_LIBRARY_NAME}_SRCDIRS}
- scalesInterface/
- scalesInterface/layout
- scalesInterface/particles
- scalesInterface/particles/advec_line
- scalesInterface/output
+ ${SCALES_DIR}/
+ ${SCALES_DIR}/layout
+ ${SCALES_DIR}/particles
+ ${SCALES_DIR}/particles/advec_line
+ ${SCALES_DIR}/output
)
endif()