From 9578a1e3168c32fe072ec977eb0d3f99a07ba700 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Franck=20P=C3=A9rignon?= <franck.perignon@imag.fr>
Date: Fri, 31 Jan 2014 13:33:56 +0000
Subject: [PATCH] Fix hdf5 bug + file output in compute_forces
---
.../hysop/operator/monitors/compute_forces.py | 54 ++++++++++++++++++-
HySoP/hysop/operator/monitors/printer.py | 37 ++++++++-----
2 files changed, 76 insertions(+), 15 deletions(-)
diff --git a/HySoP/hysop/operator/monitors/compute_forces.py b/HySoP/hysop/operator/monitors/compute_forces.py
index 28fd33b9a..da84c6436 100644
--- a/HySoP/hysop/operator/monitors/compute_forces.py
+++ b/HySoP/hysop/operator/monitors/compute_forces.py
@@ -13,6 +13,8 @@ from parmepy.operator.monitors.monitoring import Monitoring
from parmepy.mpi import main_comm, main_rank, MPI
from parmepy.tools.timers import timed_function
import parmepy.tools.numpywrappers as npw
+import scitools.filetable as ft
+import os
class DragAndLift(Monitoring):
@@ -23,14 +25,18 @@ class DragAndLift(Monitoring):
Integral inside the obstacle is not taken into account.
"""
def __init__(self, velocity, vorticity, topo, volumeOfControl,
- obstacles=None, frequency=1):
+ obstacles=None, frequency=1, filename=None, safeOutput=None):
"""
@param velocity field
@param vorticity field@
@param the topology on which forces will be computed
@param a volume of control
(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.
"""
Monitoring.__init__(self, [velocity, vorticity], topo, frequency)
self.velocity = velocity
@@ -71,6 +77,41 @@ class DragAndLift(Monitoring):
self.vd = None
# discrete vorticity field
self.wd = None
+ # 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)
+ if safeOutput is None:
+ ## 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 = True
+ self._writeFile = self._fullwrite
+ else:
+ self.safeOutput = False
+ 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):
"""
@@ -113,6 +154,7 @@ class DragAndLift(Monitoring):
assert simulation is not None,\
"Simulation parameter is required for DragAndLift apply."
dt = simulation.timeStep
+ ite = simulation.currentIteration
# Synchro of ghost points is required for fd schemes
self._synchronize(self.vd.data + self.wd.data)
@@ -138,6 +180,13 @@ class DragAndLift(Monitoring):
# Reduce results over all MPI process in topo
self.mpi_sum()
+ # 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)
+
return self.force
def _integrateOnSurface(self, surf, res):
@@ -260,6 +309,9 @@ class DragAndLift(Monitoring):
res *= self._dvol
return res
+ def finalize(self):
+ if not self.safeOutput:
+ self._file.close()
class Forces(Monitoring):
"""
diff --git a/HySoP/hysop/operator/monitors/printer.py b/HySoP/hysop/operator/monitors/printer.py
index 70b2cbc83..0549e6811 100644
--- a/HySoP/hysop/operator/monitors/printer.py
+++ b/HySoP/hysop/operator/monitors/printer.py
@@ -3,7 +3,7 @@
Classes for handling ouputs.
"""
-from parmepy.constants import np, S_DIR, PARMES_REAL, debug, VTK, HDF5, DATA
+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 os
@@ -116,11 +116,9 @@ class Printer(Monitoring):
else:
name = df.name
if len(df.data[d].shape) == 2:
- res[name] = np.expand_dims(
- df.data[d], axis=2).astype(PARMES_REAL)
+ res[name] = npw.realarray(df.data[d][:, :, np.newaxis])
else:
- res[name] = \
- df.data[d].astype(PARMES_REAL)
+ res[name] = npw.realarray(df.data[d])
return res
def set_get_data_method(self, method):
@@ -174,6 +172,7 @@ class Printer(Monitoring):
spacing[:dim] = [self.topo.mesh.space_step[i] for i in xrange(dim)]
evtk.imageToVTK(filename, origin=orig, spacing=spacing,
pointData=self._build_vtk_dict())
+
elif self.formattype == HDF5:
filename = self.prefix + str(self.topo.size)
filename += "procs_iter_{0:03d}".format(ite) + '.h5'
@@ -194,8 +193,10 @@ class Printer(Monitoring):
g_start = self.topo.mesh.global_start
g_end = self.topo.mesh.global_end + 1
- sl = tuple([slice(g_start[i], g_end[i])
- for i in xrange(self.domain.dimension)])
+ 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]
@@ -203,14 +204,16 @@ class Printer(Monitoring):
# 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,
- self.topo.globalMeshResolution - 1,
+ 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] = np.asarray(df.data[d][self.topo.mesh.iCompute],
- dtype=np.float64)
+ ds[sl] = \
+ npw.realarray(df.data[d][self.topo.mesh.iCompute].T)
f.close()
if self.topo.rank == 0:
@@ -316,15 +319,21 @@ def _TemporalGridXMF(topo, datasetNames, ite, t, filename):
g += " <Time Value=\"{0}\" />\n".format(t)
g += " <Topology TopologyType=\"" + str(topoType) + "\""
g += " NumberOfElements=\""
- g += _listFormat(topo.globalMeshResolution - 1) + " \"/>\n"
+ res = list(topo.globalMeshResolution - 1)
+ res.reverse()
+ g += _listFormat(res) + " \"/>\n"
g += " <Geometry GeometryType=\"" + geoType + "\">\n"
g += " <DataItem Dimensions=\"" + str(topo.mesh.dim) + " \""
g += " NumberType=\"Float\" Precision=\"4\" Format=\"XML\">\n"
- g += " " + _listFormat(topo.domain.origin) + "\n"
+ ori = list(topo.domain.origin)
+ ori.reverse()
+ g += " " + _listFormat(ori) + "\n"
g += " </DataItem>\n"
g += " <DataItem Dimensions=\"" + str(topo.mesh.dim) + " \""
g += " NumberType=\"Float\" Precision=\"4\" Format=\"XML\">\n"
- g += " " + _listFormat(topo.mesh.space_step) + "\n"
+ step = list(topo.mesh.space_step)
+ step.reverse()
+ g += " " + _listFormat(step) + "\n"
g += " </DataItem>\n"
g += " </Geometry>\n"
for name in datasetNames:
@@ -332,7 +341,7 @@ def _TemporalGridXMF(topo, datasetNames, ite, t, filename):
g += name + "\""
g += " AttributeType=\"Scalar\" Center=\"Node\">\n"
g += " <DataItem Dimensions=\""
- g += _listFormat(topo.globalMeshResolution - 1) + " \""
+ g += _listFormat(res) + " \""
g += " NumberType=\"Float\" Precision=\"8\" Format=\"HDF\""
g += " Compression=\"Raw\">\n" #
g += " " + filename.split('/')[-1]
--
GitLab