From 5de52ec1a6552265ec0d4b8692c1d0ad9f58df13 Mon Sep 17 00:00:00 2001
From: Chloe Mimeau <Chloe.Mimeau@imag.fr>
Date: Tue, 4 Jun 2013 15:32:13 +0000
Subject: [PATCH] Updates concerning Navier-Stokes problems
---
Docs/force_computation.tex | 25 ++--
Docs/vorticity_solenoidal_projection.tex | 8 ++
Examples/Attic/NavierStokes3d_sphere.py | 26 ++--
Examples/Attic/NavierStokes3d_vortRing.py | 71 +++++++---
Examples/testDiffusion.py | 47 +++++--
HySoP/hysop/f2py/fftw2py.f90 | 16 +++
HySoP/hysop/operator/discrete/poisson_fft.py | 57 ++++++--
.../operator/discrete/scales_advection.py | 104 ++++++++++++--
HySoP/hysop/operator/energy_enstrophy.py | 81 ++++++-----
.../hysop/operator/monitors/compute_forces.py | 2 +
HySoP/hysop/operator/monitors/printer.py | 4 +-
HySoP/hysop/operator/reprojection.py | 129 ++++++++++++++++++
HySoP/src/fftw/fft3d.f90 | 31 ++++-
13 files changed, 488 insertions(+), 113 deletions(-)
create mode 100644 HySoP/hysop/operator/reprojection.py
diff --git a/Docs/force_computation.tex b/Docs/force_computation.tex
index 5ad09d22c..24ee46fe8 100644
--- a/Docs/force_computation.tex
+++ b/Docs/force_computation.tex
@@ -121,19 +121,19 @@ Ainsi, en remplaçant $n$ par cette expression puis en effectuant une intégrati
\int_{\partial B} p \cdot n \ d\text{s}= R \int_{\partial B} p \cdot \dfrac{\partial \Gamma}{\partial s} \ d\text{s} = -R \int_{\partial B} \dfrac{\partial p}{\partial s} \cdot \Gamma \ d\text{s}
\end{equation}
%
-En supposant $\dfrac{\partial p}{\partial s} = \dfrac{\partial p}{\partial \Gamma}$, on a alors :
+On a alors :
\begin{equation}
-\int_{\partial B} p \cdot n \ d\text{s}= -R \int_{\partial B} \dfrac{\partial p}{\partial \Gamma} \cdot \Gamma \ d\text{s}
+\int_{\partial B} p \cdot n \ d\text{s}= -R \int_{\partial B} \dfrac{\partial p}{\partial s} \cdot \Gamma \ d\text{s}
\end{equation}
%
D'après l'équation de la quantité de mouvement l'égalité suivante est vérifiée :
\begin{equation}
-\nu (\nabla \times \omega) \cdot \Gamma + \dfrac{\partial p}{\partial \Gamma} \cdot \Gamma = 0
+\nu (\nabla \times \omega) \cdot \Gamma + \dfrac{\partial p}{\partial s} \cdot \Gamma = 0
\end{equation}
%
Ainsi,
\begin{equation}
--R \nu \int_{\partial B} (\nabla \times \omega) \cdot \Gamma \ d\text{s} - R \int_{\partial B} \dfrac{\partial p}{\partial \Gamma} \cdot \Gamma \ d\text{s} = 0
+-R \nu \int_{\partial B} (\nabla \times \omega) \cdot \Gamma \ d\text{s} - R \int_{\partial B} \dfrac{\partial p}{\partial s} \cdot \Gamma \ d\text{s} = 0
\end{equation}
%
Et par conséquent,
@@ -170,7 +170,7 @@ Dans le cas du cylindre 2D on a donc :
%
Et dans ce cas, le coefficient de pression est donné par
\begin{equation}
-C_p = \dfrac{F_p \cdot e_x}{u_{\infty}^2 H R} = -\dfrac{\nu}{u_{\infty}^2 H}\int_{0}^{2\pi} \dfrac{\partial\omega}{\partial n} e_{\theta} \ d\theta \ \cdot e_x
+C_p = \dfrac{F_p \cdot e_x}{u_{\infty}^2 \rho R} = -\dfrac{\nu}{u_{\infty}^2 \rho}\int_{0}^{2\pi} \dfrac{\partial\omega}{\partial n} e_{\theta} \ d\theta \ \cdot e_x
\end{equation}
%
or
@@ -180,7 +180,7 @@ e_{\theta} \cdot e_x = -(e_x \cdot e_x) \text{sin} \theta + (e_y \cdot e_x)\text
%
donc
\begin{equation}
-\boxed{C_p = \dfrac{\nu}{u_{\infty}^2 H}\int_{0}^{2\pi} \dfrac{\partial\omega}{\partial n} \ \text{sin} \theta \ d\theta}
+\boxed{C_p = \dfrac{\nu}{u_{\infty}^2 \rho}\int_{0}^{2\pi} \dfrac{\partial\omega}{\partial n} \ \text{sin} \theta \ d\theta}
\end{equation}
\underline{Calcul de la traînée de friction $F_f=- \nu \int_{\partial B} \omega \times n \ d\text{s}$}
@@ -213,12 +213,21 @@ Dans le cas du cylindre 2D on a donc :
%
Et dans ce cas, le coefficient de friction est donné par
\begin{equation}
-C_f = \dfrac{F_f \cdot e_x}{u_{\infty}^2 H R} = \dfrac{\nu}{u_{\infty}^2 H}\int_{0}^{2\pi} \omega e_{\theta} \ d\theta \ \cdot e_x
+C_f = \dfrac{F_f \cdot e_x}{u_{\infty}^2 \rho R} = \dfrac{\nu}{u_{\infty}^2 \rho}\int_{0}^{2\pi} \omega e_{\theta} \ d\theta \ \cdot e_x
\end{equation}
%
donc
\begin{equation}
-\boxed{C_f = -\dfrac{\nu}{u_{\infty}^2 H}\int_{0}^{2\pi} \omega \ \text{sin} \theta \ d\theta}
+\boxed{C_f = -\dfrac{\nu}{u_{\infty}^2 \rho}\int_{0}^{2\pi} \omega \ \text{sin} \theta \ d\theta}
\end{equation}
+\vspace{2cm}
+
+\begin{figure}[ht]
+\centering
+ \includegraphics[width=0.5\linewidth]{images/pressure_friction_drag.png}
+\caption{Les forces de pression ($PdA$) et de friction ($\tau_w dA$) s'exerçant sur l'élément de surface $dA$ d'un corps bidimensionnel.}
+\label{tore}
+\end{figure}
+
\end{document}
diff --git a/Docs/vorticity_solenoidal_projection.tex b/Docs/vorticity_solenoidal_projection.tex
index b4165ffe6..e18de789f 100644
--- a/Docs/vorticity_solenoidal_projection.tex
+++ b/Docs/vorticity_solenoidal_projection.tex
@@ -99,6 +99,14 @@ or $\widehat{\nabla \boldsymbol{\pi}}_k = i \boldsymbol{\xi}_k \widehat{\boldsym
\boxed{\widehat{\boldsymbol{\omega}'}_k = \widehat{\boldsymbol{\omega}}_k - \boldsymbol{\xi}_k \dfrac{\sum_{j=1}^3 \boldsymbol{\xi}_j \widehat{\boldsymbol{\omega}}_j}{\mid \boldsymbol{\xi} \mid^2}}, \ 1\leq k \leq 3
\end{equation}
+Un critère de reprojection peut-être donné par :
+\begin{equation}
+\max\limits_{x} \left ( \dfrac{\rvert \text{div} (\omega) \rvert }{\max\limits_{x} \left \rvert\dfrac{\partial \omega_i (x)}{\partial x_j} \right \rvert } \right ) \geq c
+\end{equation}
+
+avec $c$ une constante fixée inférieure à 1, de l'ordre de 1-10$\%$.
+
+
\end{document}
diff --git a/Examples/Attic/NavierStokes3d_sphere.py b/Examples/Attic/NavierStokes3d_sphere.py
index 2206ab9f8..ca7ca59e8 100755
--- a/Examples/Attic/NavierStokes3d_sphere.py
+++ b/Examples/Attic/NavierStokes3d_sphere.py
@@ -78,7 +78,7 @@ hy = Ly / ncells[1]
hz = Lz / ncells[2]
## Obstacle
-lambd = np.array([lambdFluid, lambdPorous, lambdSolid],
+lambd = np.array([lambdSolid, lambdPorous],
dtype=PARMES_REAL, order=ORDER)
sphere = Obstacle(myDomain3d,
name=obstName,
@@ -105,8 +105,10 @@ print "FFT solver local resolution/offset: ", localres, localoffset
##topofft = poisson.getTopology()
# ---- Cartesian topology -----
-topoCart = pp.CartesianTopology(domain=myDomain3d,
- resolution=resolution3d, dim=3, ghosts=[nbGhosts,nbGhosts,nbGhosts])
+topoCart = pp.Cartesian(domain=myDomain3d,
+ dim=3,
+ globalMeshResolution=resolution3d,
+ ghosts=[nbGhosts,nbGhosts,nbGhosts])
# ---- JB's topology ----------
nbProcs = MPI.COMM_WORLD.Get_size()
@@ -243,21 +245,21 @@ navierStokes = pp.Problem(topoCart, [advecVort, stretch, diffusion, poisson, pen
printer = Printer(fields=[vorticity, velocity, vortNorm], frequency=outModuloPrinter,
outputPrefix=outputPrinter)
-forces = Compute_forces(velocity=velocity, vorticity=vorticity,
- topology=topoCart,
- obstacle=sphere,
- boxMin=[Ox + Lx / 10., Oy + Ly / 10., Oz + Lz / 10.],
- boxMax=[Lx - Lx / 10., Ly - Ly / 10., Lz - Lz / 10.],
- Reynolds=uinf * obstRadius/ visco,
- frequency=outModuloForces,
- outputPrefix=outputForces)
+#forces = Compute_forces(velocity=velocity, vorticity=vorticity,
+# topology=topoCart,
+# obstacle=sphere,
+# boxMin=[Ox + Lx / 10., Oy + Ly / 10., Oz + Lz / 10.],
+# boxMax=[Lx - Lx / 10., Ly - Ly / 10., Lz - Lz / 10.],
+# Reynolds=uinf * obstRadius/ visco,
+# frequency=outModuloForces,
+# outputPrefix=outputForces)
energy = Energy_enstrophy(velocity=velocity, vorticity=vorticity,
topology=topoCart,
frequency=outModuloEnergies,
outputPrefix=outputEnergies)
-navierStokes.addDiagnostics([forces, energy])
+navierStokes.addDiagnostics([energy])
navierStokes.setSolver(finalTime, timeStep, solver_type='basic', io=printer)
diff --git a/Examples/Attic/NavierStokes3d_vortRing.py b/Examples/Attic/NavierStokes3d_vortRing.py
index ab8803fcf..fd7053228 100755
--- a/Examples/Attic/NavierStokes3d_vortRing.py
+++ b/Examples/Attic/NavierStokes3d_vortRing.py
@@ -1,4 +1,6 @@
#!/usr/bin/python
+import sys
+sys.path.insert(0,'/scratch/mimeau/install-parmes3/')
import time
import parmepy as pp
from parmepy.f2py import fftw2py
@@ -6,7 +8,6 @@ import numpy as np
import mpi4py.MPI as MPI
import math as m
from parmepy.constants import PARMES_REAL, ORDER
-from parmepy.domain.obstacle.obstacle_3D import Obstacle3D
from parmepy.operator.advection import Advection
from parmepy.operator.stretching import Stretching
from parmepy.operator.poisson import Poisson
@@ -14,9 +15,12 @@ from parmepy.operator.diffusion import Diffusion
from parmepy.operator.penalization import Penalization
from parmepy.operator.redistribute import Redistribute
from parmepy.problem.navier_stokes import NSProblem
+from parmepy.problem.simulation import Simulation
from parmepy.operator.monitors.energy_enstrophy import Energy_enstrophy
+from parmepy.operator.monitors.reprojection_criterion import Reprojection_criterion
from parmepy.operator.monitors.printer import Printer
+
#from numpy import linalg as LA
pi = m.pi
@@ -39,8 +43,11 @@ Lz = np.float64(inputData.readline().rstrip('\n\r'))
resX = np.uint32(inputData.readline().rstrip('\n\r'))
resY = np.uint32(inputData.readline().rstrip('\n\r'))
resZ = np.uint32(inputData.readline().rstrip('\n\r'))
-timeStep = np.float64(inputData.readline().rstrip('\n\r'))
+dt = np.float64(inputData.readline().rstrip('\n\r'))
finalTime = np.float64(inputData.readline().rstrip('\n\r'))
+restart = np.uint32(inputData.readline().rstrip('\n\r'))
+dumpfreq = np.uint32(inputData.readline().rstrip('\n\r'))
+dumpfile = inputData.readline().rstrip('\n\r')
outModuloPrinter = np.uint32(inputData.readline().rstrip('\n\r'))
outputPrinter = inputData.readline().rstrip('\n\r')
outModuloForces = np.uint32(inputData.readline().rstrip('\n\r'))
@@ -51,12 +58,17 @@ visco = np.float64(inputData.readline().rstrip('\n\r'))
vortProj = np.uint32(inputData.readline().rstrip('\n\r'))
methodAdv = inputData.readline().rstrip('\n\r')
methodStretch = inputData.readline().rstrip('\n\r')
+multires = np.uint32(inputData.readline().rstrip('\n\r'))
+resFilterX = np.uint32(inputData.readline().rstrip('\n\r'))
+resFilterY = np.uint32(inputData.readline().rstrip('\n\r'))
+resFilterZ = np.uint32(inputData.readline().rstrip('\n\r'))
inputData.close()
# Parameters
dim = 3
nbElem = [resX, resY, resZ]
+nbElemFilter = [resFilterX, resFilterY, resFilterZ]
t0 = time.time()
@@ -153,36 +165,53 @@ poisson = Poisson(velo, vorti,
resolutions={velo: nbElem,
vorti: nbElem},
method='fftw',
- projection=vortProj
+ projection=vortProj,
+ multires=multires,
+ filterSize=box.length / \
+ (np.asarray(nbElemFilter,
+ dtype=PARMES_REAL) - 1)
)
+## Diagnostics related to the problem
+energy = Energy_enstrophy(velo, vorti,
+ resolutions={velo: nbElem,
+ vorti: nbElem},
+ viscosity=visco,
+ frequency=outModuloEnergies,
+ prefix=outputEnergies)
+
+reproj = Reprojection_criterion(velo, vorti,
+ resolutions={velo: nbElem,
+ vorti: nbElem},
+ method='FD_order4',
+ frequency=1,
+ prefix='./res/Reproj.dat')
+
+printer = Printer(fields=[vorti, velo],
+ frequency=outModuloPrinter,
+ prefix=outputPrinter)
+
distrAdvStr = Redistribute([vorti, velo], advec, stretch)
distrStrPoiss = Redistribute([vorti, velo], stretch, poisson)
-## Define the problem to solve
-#pb = NSProblem([advec, diffusion, poisson],
-pb = NSProblem([advec, distrAdvStr, stretch, distrStrPoiss, diffusion, poisson],
- monitors=[Energy_enstrophy(velo, vorti,
- resolutions={velo: nbElem,
- vorti: nbElem},
- viscosity=visco,
- frequency=outModuloEnergies,
- outputPrefix=outputEnergies),
- Printer(fields=[vorti, velo],
- resolutions={velo: nbElem,
- vorti: nbElem},
- frequency=outModuloPrinter,
- outputPrefix=outputPrinter)])
-
-
+## Definition of simulation parameters
+simu = Simulation(tinit=0.0, tend=finalTime, timeStep=dt, iterMax=1000000)
+## Define the problem to solve
+#pb = NSProblem([diffusion, poisson],
+pb = NSProblem(operators=[advec, distrAdvStr, stretch, distrStrPoiss, diffusion, poisson],
+ simulation=simu, monitors=[printer, energy, reproj], dumpFreq=int(dumpfreq), name=dumpfile)
## Setting solver to Problem
-pb.setUp(finalTime, timeStep)
+pb.setUp()
+
+## Restarting Problem from dumped field values
+if restart :
+ pb.restart(dumpfile)
t1 = time.time()
## Solve problem
-#poisson.apply(0., timeStep, 0)
+#poisson.apply(simu)
timings = pb.solve()
tf = time.time()
diff --git a/Examples/testDiffusion.py b/Examples/testDiffusion.py
index 3872fc894..6ef8bfe40 100755
--- a/Examples/testDiffusion.py
+++ b/Examples/testDiffusion.py
@@ -2,8 +2,8 @@
import parmepy as pp
from parmepy.operator.diffusion import Diffusion
from parmepy.operator.poisson import Poisson
-from parmepy.problem.navier_stokes import NSProblem
from parmepy.operator.monitors.energy_enstrophy import Energy_enstrophy
+from parmepy.problem.simulation import Simulation
import numpy as np
import math
pi = math.pi
@@ -48,36 +48,57 @@ vorticity = pp.Field(domain=dom, formula=computeVort,
## Definition of the Diffusion and Poisson operators
+output = 'energy.dat'
nu = 1e-3
+finalTime = 1.0
+dt = 0.01
+restartTime = 0.5
+filenameVel = 'dump_velo_T=_'
+filenameVel += str(restartTime)
+
+filenameVort = 'dump_vort_T=_'
+filenameVort += str(restartTime)
+
diffusion = Diffusion(vorticity, resolution=resol, viscosity=nu)
poisson = Poisson(velocity, vorticity, resolutions={velocity: resol,
vorticity: resol})
-output = 'energy.dat'
-
-timeStep = 1e-2
-
+## Definition of monitor
monitor = Energy_enstrophy(velocity, vorticity,
resolutions={velocity: resol,
vorticity: resol},
viscosity=nu,
frequency=1,
- outputPrefix=output)
+ prefix=output)
+
+## Definition of simulation parameters
+tini = restartTime
+niter = (finalTime - tini) / dt
+simu = Simulation(tinit=tini, tend=finalTime, nbiter=niter)
diffusion.setUp()
poisson.setUp()
monitor.setUp()
-t = 0.0
-niter = 1
+velocity.initialize()
vorticity.initialize()
print "pre ", vorticity.norm()
+velocity.load('dump_velo_T=_0.5')
+vorticity.load('dump_vort_T=_0.5')
+
+while simu.time <= finalTime:
+ print 'time, niter : ', simu.time, simu.currentIteration
+# if (simu.currentIteration == 50):
+# print 'dump fields ...'
+# velocity.dump(filenameVel)
+# vorticity.dump(filenameVort)
+# print "norm restart ", vorticity.norm()
+ diffusion.apply(simu)
+ poisson.apply(simu)
+ monitor.apply(simu)
+ simu.time += simu.timeStep
+ simu.currentIteration += 1
-while t <= 1.0:
- diffusion.apply(timeStep)
- poisson.apply()
- monitor.apply(t, timeStep, niter)
- t += timeStep
print "post ", vorticity.norm()
diff --git a/HySoP/hysop/f2py/fftw2py.f90 b/HySoP/hysop/f2py/fftw2py.f90
index 9225ff8b6..42764aca1 100755
--- a/HySoP/hysop/f2py/fftw2py.f90
+++ b/HySoP/hysop/f2py/fftw2py.f90
@@ -192,4 +192,20 @@ contains
end subroutine projection_om_3d
+ !> Projects vorticity values from fine to coarse grid :
+ !! @param[in] dxf, dyf, dzf: grid filter size = domainLength/(CoarseRes-1)
+ !! in the following, omega is the 3D vorticity vector field.
+ subroutine multires_om_3d(dxf, dyf, dzf, omega_x,omega_y,omega_z)
+ real(pk), intent(in) :: dxf, dyf, dzf
+ real(pk),dimension(:,:,:),intent(inout):: omega_x,omega_y,omega_z
+ !f2py intent(in,out) :: omega_x,omega_y,omega_z
+
+ call r2c_3d(omega_x,omega_y,omega_z)
+
+ call filter_multires_om_3d(dxf, dyf, dzf)
+
+ call c2r_3d(omega_x,omega_y,omega_z)
+
+ end subroutine multires_om_3d
+
end module fftw2py
diff --git a/HySoP/hysop/operator/discrete/poisson_fft.py b/HySoP/hysop/operator/discrete/poisson_fft.py
index 6620e3aee..f0616c3ef 100644
--- a/HySoP/hysop/operator/discrete/poisson_fft.py
+++ b/HySoP/hysop/operator/discrete/poisson_fft.py
@@ -3,8 +3,9 @@
@file poisson_fft.py
Discrete operator for Poisson problem (fftw based)
"""
-
+import numpy as np
from parmepy.f2py import fftw2py
+from parmepy.constants import PARMES_REAL, ORDER
from discrete import DiscreteOperator
from parmepy.constants import debug, prof
from parmepy.tools.timers import timed_function
@@ -17,7 +18,8 @@ class PoissonFFT(DiscreteOperator):
"""
@debug
- def __init__(self, velocity, vorticity, method=None, projection=False):
+ def __init__(self, velocity, vorticity, method=None, projection=False,
+ multires=False, filterSize=None):
"""
Constructor.
@param[out] velocity : discretisation of the solution field
@@ -29,6 +31,10 @@ class PoissonFFT(DiscreteOperator):
self.vorticity = vorticity
## Solenoidal projection on vorticity ?
self.projection = projection
+ ## Multiresolution ?
+ self.multires = multires
+ ## Filter size array = domainLength/(CoarseRes-1)
+ self.filterSize = filterSize
# Base class initialisation
DiscreteOperator.__init__(self, [velocity, vorticity], method,
name="Poisson FFT")
@@ -51,6 +57,8 @@ class PoissonFFT(DiscreteOperator):
@prof
@timed_function
def apply(self, simulation=None):
+ ite = simulation.currentIteration
+
if self.case is '2D':
self.velocity.data[0], self.velocity.data[1] =\
fftw2py.solve_poisson_2d(self.vorticity.data,
@@ -58,22 +66,47 @@ class PoissonFFT(DiscreteOperator):
self.velocity.data[1])
elif self.case is '3D':
- #=== SOLENOIDAL PROJECTION ON VORTICITY FIELD ===
- if (self.projection):
+ #=== SOLENOIDAL PROJECTION OF VORTICITY FIELD ===
+ if (self.projection and ite % 10 == 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])
- 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])
+ if (self.multires):
+ # Projects vorticity values from fine to coarse grid
+ # in frequencies space by nullifying the smallest modes
+ vortFilter = np.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])
+
+ 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])
else:
raise ValueError("invalid problem dimension")
diff --git a/HySoP/hysop/operator/discrete/scales_advection.py b/HySoP/hysop/operator/discrete/scales_advection.py
index 845a55d57..d35a2eb50 100644
--- a/HySoP/hysop/operator/discrete/scales_advection.py
+++ b/HySoP/hysop/operator/discrete/scales_advection.py
@@ -7,7 +7,7 @@ Discrete Advection operator based on scales library (Jean-Baptiste)
from parmepy.f2py import scales2py
from discrete import DiscreteOperator
from parmepy.fields.vector import VectorField
-from parmepy.constants import debug
+from parmepy.constants import debug, np
from parmepy.tools.timers import timed_function
@@ -47,20 +47,98 @@ class ScalesAdvection(DiscreteOperator):
raise ValueError("Missing simulation value for computation.")
dt = simulation.timeStep
- for adF in self.advectedFields:
- if isinstance(adF, VectorField):
- for d in xrange(adF.dimension):
- adF[d][...] = scales2py.solve_advection(
- dt, self.velocity.data[0],
+
+# for adF in self.advectedFields:
+# if isinstance(adF, VectorField):
+# for d in xrange(adF.dimension):
+# adF[d][...] = scales2py.solve_advection(
+# dt, self.velocity.data[0],
+# self.velocity.data[1],
+# self.velocity.data[2],
+# adF[d][...])
+# else:
+# adF[...] = scales2py.solve_advection(
+# dt, self.velocity.data[0],
+# self.velocity.data[1],
+# self.velocity.data[2],
+# adF.data)
+
+
+ ## Space discretization of grad(u)
+ mesh_size = self.velocity.topology.mesh.space_step
+ ind=np.arange(self.velocity.topology.mesh.resolution[0])
+
+### First order scheme
+### X components of temp and result
+ temp1 = (self.velocity[0][...] -
+ self.velocity[0][np.roll(ind,+1),...]) / (mesh_size[0])
+
+ temp2 = (self.velocity[0][...] -
+ self.velocity[0][:,np.roll(ind,+1),:]) / (mesh_size[1])
+
+ temp3 = (self.velocity[0][...] -
+ self.velocity[0][...,np.roll(ind,+1)]) / (mesh_size[2])
+
+ maxstr1= np.max(abs(temp1)+abs(temp2)+abs(temp3))
+ maxadv1= np.max(abs(temp1))
+
+# Y components of temp and result
+ temp4 = (self.velocity[1][...] -
+ self.velocity[1][np.roll(ind,+1),...]) / (mesh_size[0])
+
+ temp5 = (self.velocity[1][...] -
+ self.velocity[1][:,np.roll(ind,+1),:]) / (mesh_size[1])
+
+ temp6 = (self.velocity[1][...] -
+ self.velocity[1][...,np.roll(ind,+1)]) / (mesh_size[2])
+
+ maxstr2= np.max(abs(temp4)+abs(temp5)+abs(temp6))
+ maxadv2= np.max(abs(temp5))
+
+# Z components of temp and result
+ temp7 = (self.velocity[2][...] -
+ self.velocity[2][np.roll(ind,+1),...]) / (mesh_size[0])
+
+ temp8 = (self.velocity[2][...] -
+ self.velocity[2][:,np.roll(ind,+1),:]) / (mesh_size[1])
+
+ temp9 = (self.velocity[2][...] -
+ self.velocity[2][...,np.roll(ind,+1)]) / (mesh_size[2])
+
+ maxstr3= np.max(abs(temp7)+abs(temp8)+abs(temp9))
+ maxadv3= np.max(abs(temp9))
+ maxAdv = max(maxadv1,maxadv2,maxadv3)
+
+ stabcoeff = 1.0 / 2.0
+
+ print "dtAdvec", stabcoeff / maxAdv
+
+ ## Time Subcycles
+ dtadv = dt
+ dtadv = min(dtadv, stabcoeff / maxAdv)
+ if dt == dtadv:
+ print "dtadv, ndt, stab/max ", dtadv, int(dt / dtadv), stabcoeff / maxAdv
+ ndt = int(dt / dtadv)
+ else:
+ print "dtadv, ndt", dtadv, int(dt / dtadv) + 1, stabcoeff / maxAdv
+ ndt = int(dt / dtadv) + 1
+ dtadv = dt / float(ndt)
+
+ for i in range(ndt):
+ for adF in self.advectedFields:
+ if isinstance(adF, VectorField):
+ for d in xrange(adF.dimension):
+ adF[d][...] = scales2py.solve_advection(
+ dtadv, self.velocity.data[0],
+ self.velocity.data[1],
+ self.velocity.data[2],
+ adF[d][...])
+ else:
+ adF[...] = scales2py.solve_advection(
+ dtadv, self.velocity.data[0],
self.velocity.data[1],
self.velocity.data[2],
- adF[d][...])
- else:
- adF[...] = scales2py.solve_advection(
- dt, self.velocity.data[0],
- self.velocity.data[1],
- self.velocity.data[2],
- adF.data)
+ adF.data)
def finalize(self):
"""
diff --git a/HySoP/hysop/operator/energy_enstrophy.py b/HySoP/hysop/operator/energy_enstrophy.py
index 1486fc1c1..83725d448 100644
--- a/HySoP/hysop/operator/energy_enstrophy.py
+++ b/HySoP/hysop/operator/energy_enstrophy.py
@@ -3,9 +3,11 @@
@file energy_enstrophy.py
Compute Energy and Enstrophy
"""
-from parmepy.constants import np, debug, PARMES_REAL
+import numpy as np
+from numpy import linalg as la
+from parmepy.constants import np, debug, PARMES_REAL, PARMES_MPI_REAL
from parmepy.operator.monitors.monitoring import Monitoring
-from parmepy.mpi import main_rank, main_comm, MPI
+from parmepy.mpi import main_rank, main_comm, main_size, MPI
from parmepy.tools.timers import timed_function
@@ -50,12 +52,15 @@ class Energy_enstrophy(Monitoring):
Monitoring.setUp(self)
discreteFields = {}
self.spaceStep = {}
+ topodims = np.ones((self.domain.dimension))
+ topodims[-1] = main_size
# Variables discretization
for v in self.variables:
# the topology for v ...
topo = self.domain.getOrCreateTopology(self.domain.dimension,
- self.resolutions[v])
+ self.resolutions[v],
+ topodims)
# ... and the corresponding discrete field
discreteFields[v] = v.discretize(topo)
self.spaceStep[v] = topo.mesh.space_step
@@ -87,57 +92,71 @@ class Energy_enstrophy(Monitoring):
else:
raise ValueError("invalid problem dimension")
- energy = (0.5 * squareNormVel) / (2.0 * np.pi) ** 3
+ myEnergy = (0.5 * squareNormVel) / (2.0 * np.pi) ** 3
# Enstrophy computation
dvol = np.prod(self.spaceStep[self.vorticity])
if (self.vort.dimension == 2):
- enstrophy = np.sum(self.vort.data ** 2) * dvol
+ myEnstrophy = np.sum(self.vort.data ** 2) * dvol
elif (self.vort.dimension == 3):
- enstrophy = np.sum(self.vort.data[0] ** 2 +
+ myEnstrophy = np.sum(self.vort.data[0] ** 2 +
self.vort.data[1] ** 2 +
self.vort.data[2] ** 2) * dvol
else:
raise ValueError("invalid problem dimension")
- enstrophy = enstrophy / (2.0 * np.pi) ** 3
+ myEnstrophy = myEnstrophy / (2.0 * np.pi) ** 3
- # 1st order
-# effectiveViscosity = - (energy - self.buffer_1) / (dt * enstrophy)
-# timeDerivativeEnergy = - (energy - self.buffer_1) / dt
+ # Collective communications
- # 2nd order
- effectiveViscosity = - (3.0 * energy - 4.0 *
- self.buffer_1 + 1.0 * self.buffer_2) / \
- (2.0 * dt * enstrophy)
+ topovelo = self.velo.topology
+ energy = topovelo.comm.reduce(myEnergy, PARMES_MPI_REAL, op=MPI.SUM, root=0)
+ energyBuff1 = topovelo.comm.reduce(self.buffer_1, PARMES_MPI_REAL, op=MPI.SUM, root=0)
+ energyBuff2 = topovelo.comm.reduce(self.buffer_2, PARMES_MPI_REAL, op=MPI.SUM, root=0)
- timeDerivativeEnergy = - (3.0 * energy - 4.0 * self.buffer_1 + 1.0 *
- self.buffer_2) / (2.0 * dt)
+ topovort = self.vort.topology
+ enstrophy = topovort.comm.reduce(myEnstrophy, PARMES_MPI_REAL, op=MPI.SUM, root=0)
- ratio = effectiveViscosity / self.viscosity
+ # Update buffers
- nuEnstrophy = self.viscosity * enstrophy
+ self.buffer_2 = self.buffer_1
+ self.buffer_1 = myEnergy
- nu_effEnstrophy = effectiveViscosity * enstrophy
+ # Effective viscosity computation
- self.buffer_2 = self.buffer_1
- self.buffer_1 = energy
+ if main_rank == 0 :
+
+ # 1st order
+ # effectiveViscosity = - (energy - energyBuff1) / (dt * enstrophy)
+ # timeDerivativeEnergy = - (energy - energyBuff1) / dt
+
+ # 2nd order
+ effectiveViscosity = - (3.0 * energy - 4.0 *
+ energyBuff1 + 1.0 * energyBuff2) / \
+ (2.0 * dt * enstrophy)
+
+ ratio = effectiveViscosity / self.viscosity
+
+ # Energy dissipation computation
+
+ timeDerivativeEnergy = - (3.0 * energy - 4.0 * energyBuff1 + 1.0 *
+ energyBuff2) / (2.0 * dt)
+
+ nuEnstrophy = self.viscosity * enstrophy
-# dissipE = 0.
-# main_comm.Reduce(nuEnstrophy, dissipE, op=MPI.SUM, root=0)
- dissipE = main_comm.reduce(nuEnstrophy, op=MPI.SUM, root=0)
+ nu_effEnstrophy = effectiveViscosity * enstrophy
+ # dissipE = main_comm.reduce(nuEnstrophy, PARMES_MPI_REAL, op=MPI.SUM, root=0)
if ((main_rank == 0) and (ite % self.freq == 0) and (t > dt)):
self.f = open(self.prefix, 'a')
-# self.f.write("%s %s %s %s %s %s %s\n" % (t, energy,
-# enstrophy,
-# ratio,
-# timeDerivativeEnergy,
-# nuEnstrophy,
-# nu_effEnstrophy))
- self.f.write("%s %s\n" % (t, dissipE))
+ self.f.write("%s %s %s %s %s %s %s\n" % (t, energy,
+ enstrophy,
+ ratio,
+ timeDerivativeEnergy,
+ nuEnstrophy,
+ nu_effEnstrophy))
self.f.close()
def __str__(self):
diff --git a/HySoP/hysop/operator/monitors/compute_forces.py b/HySoP/hysop/operator/monitors/compute_forces.py
index 8a9acbed7..5e99259b0 100644
--- a/HySoP/hysop/operator/monitors/compute_forces.py
+++ b/HySoP/hysop/operator/monitors/compute_forces.py
@@ -399,6 +399,8 @@ class Compute_forces(Monitoring):
"""
t = simulation.time
dt = simulation.timeStep
+ ite = simulation.currentIteration
+
if self.nbPointsBox == 0:
pass
else :
diff --git a/HySoP/hysop/operator/monitors/printer.py b/HySoP/hysop/operator/monitors/printer.py
index dfa31245f..4f1656e4d 100644
--- a/HySoP/hysop/operator/monitors/printer.py
+++ b/HySoP/hysop/operator/monitors/printer.py
@@ -115,8 +115,8 @@ class Printer(Monitoring):
except AttributeError:
pass
# Set output file name
- filename = self.prefix + "iter_{0:06d}".format(ite)
- filename += '_' + str(main_rank) + self.ext
+ filename = self.prefix + str(main_rank)
+ filename += '_' + "iter_{0:02d}".format(ite) + self.ext
print "Print to file " + filename
## VTK output \todo: Need fix in 2D, getting an IOError.
if self.ext == '.vtk':
diff --git a/HySoP/hysop/operator/reprojection.py b/HySoP/hysop/operator/reprojection.py
new file mode 100644
index 000000000..87e3a2e99
--- /dev/null
+++ b/HySoP/hysop/operator/reprojection.py
@@ -0,0 +1,129 @@
+# -*- coding: utf-8 -*-
+"""
+@file reprojection_criterion.py
+Compute reprojection criterion and divergence maximum
+"""
+import numpy as np
+from numpy import linalg as la
+from parmepy.constants import np, debug, PARMES_REAL, PARMES_MPI_REAL
+from parmepy.operator.monitors.monitoring import Monitoring
+from parmepy.numerics.differentialOperator import DifferentialOperator
+from parmepy.mpi import main_rank, main_comm, main_size, MPI
+from parmepy.tools.timers import timed_function
+
+
+class Reprojection_criterion(Monitoring):
+ """
+ Computes and prints reprojection criterion and divergence maximum
+ """
+
+ def __init__(self, velocity, vorticity,
+ resolutions=None,
+ method=None,
+ frequency=None, prefix=None):
+ """
+ Constructor.
+ @param fields : List of fields
+ @param resolutions : list of resolutions (one per variable)
+ @param frequency : output file producing frequency.
+ @param prefix : file name prefix, contains relative path.
+ """
+ Monitoring.__init__(self, [velocity, vorticity], frequency,
+ name="Reprojection")
+ self.prefix = prefix
+ self.velocity = velocity
+ self.vorticity = vorticity
+ self.resolutions = resolutions
+ self.method = method
+ if (main_rank == 0):
+ self.f = open(self.prefix, 'w')
+ self.input = [velocity, vorticity]
+ self.output = []
+
+ def setUp(self):
+ Monitoring.setUp(self)
+ discreteFields = {}
+ topodims = np.ones((self.domain.dimension))
+ topodims[-1] = main_size
+ nbGhosts = 0
+ if self.method.find('FD_order2') >= 0:
+ nbGhosts = 1
+ elif self.method.find('FD_order4') >= 0:
+ nbGhosts = 2
+ else:
+ nbGhosts = 2
+
+ # Variables discretization
+ ghosts = np.zeros((self.domain.dimension))
+ ghosts[:] = nbGhosts
+ for v in self.variables:
+ # the topology for v ...
+ topo = self.domain.getOrCreateTopology(self.domain.dimension,
+ self.resolutions[v],
+ ghosts=ghosts)
+ # ... and the corresponding discrete field
+ discreteFields[v] = v.discretize(topo)
+ self.velo = discreteFields[self.velocity]
+ self.vort = discreteFields[self.vorticity]
+
+ self._isUpToDate = True
+
+ @debug
+ @timed_function
+ def apply(self, simulation=None):
+ """
+ Computes and prints reprojection criterion and divergence maximum
+ """
+ t = simulation.time
+ dt = simulation.timeStep
+ ite = simulation.currentIteration
+ gradVort, maxArrayVort, divVort = \
+ DifferentialOperator(self.vort, self.vort, self.method,
+ choice='gradV').__call__()
+
+ gradVelo, maxArrayVelo, divVelo = \
+ DifferentialOperator(self.velo, self.velo, self.method,
+ choice='gradV').__call__()
+
+ maxDiv = np.max(abs(divVort))
+ maxdwx_dx = np.max(abs(gradVort[0]))
+ maxdwx_dy = np.max(abs(gradVort[1]))
+ maxdwx_dz = np.max(abs(gradVort[2]))
+ maxdwy_dx = np.max(abs(gradVort[3]))
+ maxdwy_dy = np.max(abs(gradVort[4]))
+ maxdwy_dz = np.max(abs(gradVort[5]))
+ maxdwz_dx = np.max(abs(gradVort[6]))
+ maxdwz_dy = np.max(abs(gradVort[7]))
+ maxdwz_dz = np.max(abs(gradVort[8]))
+ maxmax = max(maxdwx_dx, maxdwx_dy, maxdwx_dz,
+ maxdwy_dx, maxdwy_dy, maxdwy_dz,
+ maxdwz_dx, maxdwz_dy, maxdwz_dz)
+ projCriterion = maxDiv / maxmax
+ dtstabAdv = 0.5 / maxArrayVelo[1]
+ dtstabStretch = 2.5127 / maxArrayVelo[0]
+
+ if ((main_rank == 0) and (ite % self.freq == 0) and (t > dt)):
+ self.f = open(self.prefix, 'a')
+ self.f.write("%s %s %s %s %s %s %s %s %s %s %s %s %s %s\n" % (t, projCriterion,
+ maxDiv,
+ maxdwx_dx,
+ maxdwx_dy,
+ maxdwx_dz,
+ maxdwy_dx,
+ maxdwy_dy,
+ maxdwy_dz,
+ maxdwz_dx,
+ maxdwz_dy,
+ maxdwz_dz,
+ dtstabAdv,
+ dtstabStretch))
+ 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__
diff --git a/HySoP/src/fftw/fft3d.f90 b/HySoP/src/fftw/fft3d.f90
index d06ecd737..5a4a8878a 100755
--- a/HySoP/src/fftw/fft3d.f90
+++ b/HySoP/src/fftw/fft3d.f90
@@ -25,7 +25,8 @@ module fft3d
public :: init_r2c_3d,init_c2c_3d,r2c_3d,c2c_3d,c2r_3d,cleanFFTW_3d,&
getParamatersTopologyFFTW3d,filter_poisson_3d,filter_curl_diffusion_3d, &
init_r2c_3d_many, r2c_3d_many, c2r_3d_many, filter_diffusion_3d_many,&
- filter_poisson_3d_many, filter_diffusion_3d, filter_projection_om_3d
+ filter_poisson_3d_many, filter_diffusion_3d, filter_projection_om_3d,&
+ filter_multires_om_3d
!> plan for fftw "c2c" forward or r2c transform
type(C_PTR) :: plan_forward1, plan_forward2, plan_forward3
@@ -682,6 +683,34 @@ contains
end do
end subroutine filter_projection_om_3d
+
+ !> Projects vorticity values from fine to coarse grid :
+ !> the smallest modes of vorticity are nullified
+ !! @param[in] dxf, dyf, dzf: grid filter size = domainLength/(CoarseRes-1)
+ subroutine filter_multires_om_3d(dxf, dyf, dzf)
+
+ real(C_DOUBLE), intent(in) :: dxf, dyf, dzf
+ integer(C_INTPTR_T) :: i,j,k
+ real(C_DOUBLE) :: kxc, kyc, kzc
+
+ kxc = pi / dxf
+ kyc = pi / dyf
+ kzc = pi / dzf
+
+ !! mind the transpose -> index inversion between y and z
+ do j = 2,local_resolution(c_Y)
+ do k = 2, fft_resolution(c_Z)
+ do i = 2, local_resolution(c_X)
+ if ((abs(kx(i))>kxc) .and. (abs(ky(j))>kyc) .and. (abs(kz(k))>kzc)) then
+ dataout1(i,k,j) = 0.
+ dataout2(i,k,j) = 0.
+ dataout3(i,k,j) = 0.
+ endif
+ end do
+ end do
+ end do
+
+ end subroutine filter_multires_om_3d
!> Solve Poisson problem in the Fourier space :
!! \f{eqnarray*} \Delta \psi &=& - \omega \\ v &=& \nabla\times\psi \f}
--
GitLab