From 66191dab6c51510a4dd2f7b828452226c9961f8d Mon Sep 17 00:00:00 2001
From: Chloe Mimeau <Chloe.Mimeau@imag.fr>
Date: Wed, 28 Aug 2013 08:20:18 +0000
Subject: [PATCH] general updates
---
.../hysop/numerics/tests/test_integrators.py | 2 --
.../hysop/operator/discrete/diffusion_fft.py | 1 +
.../operator/discrete/scales_advection.py | 30 ++++++++++++++++++-
HySoP/hysop/operator/energy_enstrophy.py | 11 +++----
HySoP/hysop/operator/monitors/printer.py | 1 +
HySoP/hysop/problem/navier_stokes.py | 6 ++++
HySoP/hysop/problem/problem.py | 3 ++
HySoP/hysop/problem/simulation.py | 1 +
HySoP/hysop/tools/cpu_data_transfer.py | 3 +-
HySoP/hysop/tools/cpu_data_transfer_S.py | 9 +++---
10 files changed, 50 insertions(+), 17 deletions(-)
diff --git a/HySoP/hysop/numerics/tests/test_integrators.py b/HySoP/hysop/numerics/tests/test_integrators.py
index ad53ca3ee..3137fb9ac 100644
--- a/HySoP/hysop/numerics/tests/test_integrators.py
+++ b/HySoP/hysop/numerics/tests/test_integrators.py
@@ -110,5 +110,3 @@ def test_RK4_1D():
nbSteps = 100
dt, err = integrate(RK4(1, func1D), nbSteps)
assert err < dt ** 4
-
-
diff --git a/HySoP/hysop/operator/discrete/diffusion_fft.py b/HySoP/hysop/operator/discrete/diffusion_fft.py
index 105564a2c..6f3bfb23f 100644
--- a/HySoP/hysop/operator/discrete/diffusion_fft.py
+++ b/HySoP/hysop/operator/discrete/diffusion_fft.py
@@ -41,6 +41,7 @@ class DiffusionFFT(DiscreteOperator):
def apply(self, simulation):
if simulation is None:
raise ValueError("Missing dt value for diffusion computation.")
+
dt = simulation.timeStep
if (self.vorticity.dimension == 2):
self.vorticity.data = fftw2py.solve_diffusion_2d(
diff --git a/HySoP/hysop/operator/discrete/scales_advection.py b/HySoP/hysop/operator/discrete/scales_advection.py
index 038e30618..2d770b59d 100644
--- a/HySoP/hysop/operator/discrete/scales_advection.py
+++ b/HySoP/hysop/operator/discrete/scales_advection.py
@@ -7,6 +7,7 @@ Discrete Advection operator based on scales library (Jean-Baptiste)
from parmepy.f2py import scales2py
from discrete import DiscreteOperator
from parmepy.constants import debug, np
+from parmepy.numerics.differential_operations import GradVxW
from parmepy.tools.timers import timed_function
@@ -46,7 +47,19 @@ class ScalesAdvection(DiscreteOperator):
raise ValueError("Missing simulation value for computation.")
dt = simulation.timeStep
-
+# # Compute the number of required subcycles
+# ndt, subdt = self.checkStability(dt)
+# # Call scales advection
+# print "time steps ...", dt, subdt
+# assert sum(subdt) == dt
+# for i in xrange(ndt):
+# for adF in self.advectedFields:
+# for d in xrange(adF.nbComponents):
+# adF[d][...] = scales2py.solve_advection(
+# subdt[i], self.velocity.data[0],
+# self.velocity.data[1],
+# self.velocity.data[2],
+# adF[d][...])
for adF in self.advectedFields:
for d in xrange(adF.nbComponents):
adF[d][...] = scales2py.solve_advection(
@@ -55,6 +68,21 @@ class ScalesAdvection(DiscreteOperator):
self.velocity.data[2],
adF[d][...])
+ def checkStability(self, dt):
+ """
+ Computes a stability condition depending on velocity gradient
+ @param current time step
+ @return int,float : the number of required subcycles and
+ the subcycle time-step.
+ """
+ dt_stab = min(dt, self.cststretch / self.diagnostics[0])
+ nb_cycles = int(ceil(dt / dt_stab))
+ subdt = np.zeros((nb_cycles))
+ subdt[:] = dt_stab
+ subdt[-1] = dt - (nb_cycles - 1) * dt_stab
+ print 'dt_advec', dt_stab
+ return nb_cycles, subdt
+
def finalize(self):
"""
\todo check memory deallocation in scales???
diff --git a/HySoP/hysop/operator/energy_enstrophy.py b/HySoP/hysop/operator/energy_enstrophy.py
index 5f91121f3..22208d0cc 100644
--- a/HySoP/hysop/operator/energy_enstrophy.py
+++ b/HySoP/hysop/operator/energy_enstrophy.py
@@ -102,15 +102,12 @@ class Energy_enstrophy(Monitoring):
enstrophy = self.coeffEnstrophy * np.sum([vortd[i][ind] ** 2
for i in xrange(nbc)])
# Collective communications
-
- energy = self.topovel.comm.reduce(energy, PARMES_MPI_REAL,
- op=MPI.SUM, root=0)
- energyBuff1 = self.topovel.comm.reduce(self._buffer_1, PARMES_MPI_REAL,
- op=MPI.SUM, root=0)
- energyBuff2 = self.topovel.comm.reduce(self._buffer_2, PARMES_MPI_REAL,
- op=MPI.SUM, root=0)
+ energy = self.topovel.comm.allreduce(energy, PARMES_MPI_REAL,
+ op=MPI.SUM)
enstrophy = self.topovort.comm.reduce(enstrophy, PARMES_MPI_REAL,
op=MPI.SUM, root=0)
+ energyBuff1 = self._buffer_1
+ energyBuff2 = self._buffer_2
# Update buffers
diff --git a/HySoP/hysop/operator/monitors/printer.py b/HySoP/hysop/operator/monitors/printer.py
index 374bcd932..e9c121495 100644
--- a/HySoP/hysop/operator/monitors/printer.py
+++ b/HySoP/hysop/operator/monitors/printer.py
@@ -130,6 +130,7 @@ class Printer(Monitoring):
orig = tuple([topo.mesh.coords[i].flatten()[ind[i]]
for i in xrange(topo.mesh.dim)])
spacing = tuple(topo.mesh.space_step)
+ print 'printing output vtk file ...'
evtk.imageToVTK(filename, origin=orig, spacing=spacing,
pointData=self._build_vtk_dict())
elif self.ext == '.dat':
diff --git a/HySoP/hysop/problem/navier_stokes.py b/HySoP/hysop/problem/navier_stokes.py
index f40c5284e..4ec4fb45c 100644
--- a/HySoP/hysop/problem/navier_stokes.py
+++ b/HySoP/hysop/problem/navier_stokes.py
@@ -4,9 +4,11 @@
from parmepy.problem.problem import Problem
from parmepy.operator.analytic import Analytic
from parmepy.operator.advection import Advection
+from parmepy.operator.density import DensityVisco
from parmepy.operator.stretching import Stretching
from parmepy.operator.poisson import Poisson
from parmepy.operator.diffusion import Diffusion
+from parmepy.operator.multiphase import Baroclinic
from parmepy.operator.penalization import Penalization
@@ -22,12 +24,16 @@ class NSProblem(Problem):
for op in operators:
if isinstance(op, Advection):
self.advection = op
+ if isinstance(op, DensityVisco):
+ self.density = op
if isinstance(op, Stretching):
self.stretch = op
if isinstance(op, Diffusion):
self.diffusion = op
if isinstance(op, Poisson):
self.poisson = op
+ if isinstance(op, Baroclinic):
+ self.multiphase = op
if isinstance(op, Penalization):
self.penal = op
if isinstance(op, Analytic):
diff --git a/HySoP/hysop/problem/problem.py b/HySoP/hysop/problem/problem.py
index 35a3a80f3..db8783c6d 100644
--- a/HySoP/hysop/problem/problem.py
+++ b/HySoP/hysop/problem/problem.py
@@ -6,6 +6,7 @@ Complete problem description.
from parmepy.constants import debug
from parmepy import __VERBOSE__
from parmepy.operator.monitors.monitoring import Monitoring
+from parmepy.operator.continuous import Operator
from parmepy.operator.redistribute import Redistribute
from parmepy.tools.timers import Timer, timed_function
from parmepy.mpi import main_rank
@@ -114,6 +115,8 @@ class Problem(object):
for op in self.operators:
if isinstance(op, Monitoring):
op.setUp()
+
+ for op in self.operators:
if isinstance(op, Redistribute):
op.setUp()
diff --git a/HySoP/hysop/problem/simulation.py b/HySoP/hysop/problem/simulation.py
index 35cdb9f42..7584b4bc7 100644
--- a/HySoP/hysop/problem/simulation.py
+++ b/HySoP/hysop/problem/simulation.py
@@ -22,6 +22,7 @@ class Simulation(object):
"""
## Simulation final time
self.end = tend
+ print 'final time:', tend, self.end, iterMax
## Starting time
self.start = tinit
## Simulation current time
diff --git a/HySoP/hysop/tools/cpu_data_transfer.py b/HySoP/hysop/tools/cpu_data_transfer.py
index 49d6ced8c..576ae0ab7 100644
--- a/HySoP/hysop/tools/cpu_data_transfer.py
+++ b/HySoP/hysop/tools/cpu_data_transfer.py
@@ -3,8 +3,7 @@
@file cpu_data_transfer.py
Communicator MPI to synchronize ghost
"""
-from parmepy.constants import ORDER
-import numpy as np
+from parmepy.constants import np, ORDER
from parmepy.mpi import main_comm
import time
diff --git a/HySoP/hysop/tools/cpu_data_transfer_S.py b/HySoP/hysop/tools/cpu_data_transfer_S.py
index 96c8c1877..90b908a45 100644
--- a/HySoP/hysop/tools/cpu_data_transfer_S.py
+++ b/HySoP/hysop/tools/cpu_data_transfer_S.py
@@ -3,9 +3,8 @@
@file cpu_data_transfer_S.py
Communicator MPI to synchronize ghosts of scalar fields
"""
-from ..constants import *
-import numpy as np
-import mpi4py.MPI as MPI
+from parmepy.constants import np, ORDER
+from parmepy.mpi import main_comm
import time
@@ -67,8 +66,8 @@ class SynchronizeS(object):
# NORTH SOUTH
if(self.topo.tabSort[i,1] == 4):
- f[0][:,:,self.topo.mesh.resolution[2]-self.topo.ghosts[2]:self.topo.mesh.resolution[2]] =\
- f[0][:,:,self.topo.ghosts[2]:2*self.topo.ghosts[2]]
+ f[:,:,self.topo.mesh.resolution[2]-self.topo.ghosts[2]:self.topo.mesh.resolution[2]] =\
+ f[:,:,self.topo.ghosts[2]:2*self.topo.ghosts[2]]
if(self.topo.tabSort[i,1] == 5):
f[:,:,0:self.topo.ghosts[2]] =\
--
GitLab