Merge branch 'master' into fix/hdf5-deadlock

CHANGELOG.md

@@ -11,6 +11,9 @@ v4.5.0
 - Added in sum factorisation version for pyramid expansions and orthogonal
   expansion in pyramids (!750)
 
+**FieldConvert**:
+- Add input module for Semtex field files (!777)
+
 **Documentation**:
 - Added the developer-guide repository as a submodule (!751)
 
@@ -22,9 +25,15 @@ v4.4.1
 - Remove the duplicate output of errorutil (!756)
 - Fix BLAS CMake dependencies (!763)
 - Fix interpolation issue with Lagrange basis functions (!768)
+- Fix issue with average fields not working with different polynomial order
+  fields (!776)
+- Fix Hdf5 output in FilterFieldConvert (!781)
+- Fixed extreme memory consumption of Interpolator when interpolating from pts
+  to fld or between different meshes (!783)
 - Fix deadlock with HDF5 input (!786)
 
-**FieldConvert**:
+**FieldConvert:**
+- Fix issue with field ordering in the interppointdatatofld module (!754)
 - Fix issue with FieldConvert when range flag used (!761)
 
 **NekMesh**:
@@ -34,8 +43,8 @@ v4.4.1
 - Add manifold meshing option (!756)
 - Fix issue with older rea input files (!765)
 
-**FieldConvert:**
-- Fix issue with field ordering in the interppointdatatofld module (!754)
+**IncNavierStokesSolver**
+- Fix an initialisation issue when using an additional advective field (!779)
 
 v4.4.0
 ------

docs/user-guide/utilities/fieldconvert.tex

 \chapter{FieldConvert}
 \label{s:utilities:fieldconvert}
-FieldConvert is a utility embedded in \nekpp with the primary aim of
-allowing the user to convert the \nekpp output binary files (.chk and
-.fld) into a format which can be read by two common visualisation
-softwares: Paraview/VisIt (.vtu format) or Tecplot/VisIt (.dat or .plt
-formats). FieldConvert also allows the user to manipulate the \nekpp
-output binary files by using some additional modules which can be
-called with the option \inltt{-m} which stands for
-\inltt{m}odule. Note that another flag, \inltt{-r} (which stand for
-\inltt{r}ange) allows the user to specify a sub-range of the domain on
-which the conversion or manipulation of the \nekpp output binary files
-will be performed.
+FieldConvert is a utility embedded in \nekpp with the primary aim of allowing
+the user to convert the \nekpp output binary files (\inltt{.chk} and
+\inltt{.fld}) into formats which can be read by common visualisation and
+post-processing software, primarily Paraview/VisIt (in unstructured VTK
+\inltt{.vtu} format) or Tecplot/VisIt (in ASCII \inltt{.dat} or binary
+\inltt{.plt} formats). FieldConvert also allows the user to manipulate the
+\nekpp output binary files by using some additional modules which can be called
+with the option \inltt{-m} which stands for \inltt{m}odule. Note that another
+flag, \inltt{-r} (which stand for \inltt{r}ange) allows the user to specify a
+sub-range of the domain on which the conversion or manipulation of the \nekpp
+output binary files will be performed.
 
 Almost all of the FieldConvert functionalities can be run in parallel if \nekpp
 is compiled using MPI (see the installation documentation for additional info on
@@ -19,6 +19,49 @@ how to implement \nekpp using MPI). \footnote{Modules that do not have parallel
 %
 %
 %
+\section{Basic usage}
+FieldConvert expects at least one input specification (such as a session file
+and its corresponding field file) and one output specification. These are
+specified on the command line as
+%
+\begin{lstlisting}[style=BashInputStyle]
+  FieldConvert in1.xml in2.fld out.dat
+\end{lstlisting}
+%
+These can be combined with a processing module by adding the \inltt{-m} command
+line option. There can be more than one module specified, and they can appear
+anywhere in the command line arguments, although the order of execution is
+inferred from their order in the command line. For example, the command
+%
+\begin{lstlisting}[style=BashInputStyle]
+  FieldConvert in1.xml -m module1 in2.fld -m module2 out.dat
+\end{lstlisting}
+%
+causes \inltt{in1.xml} and \inltt{in2.fld} to be read, followed by the
+\inltt{module1} processing module, the \inltt{module2} processing module, and
+finally output to the \inltt{out.dat} Tecplot file.
+
+\subsection{Input formats}
+
+FieldConvert supports XML and FLD-format files as produced by \nekpp. It also
+supports the reading of data files from two external spectral element codes:
+\emph{Semtex}\footnote{http://users.monash.edu.au/~bburn/semtex.html} and
+\emph{Nek5000}\footnote{https://nek5000.mcs.anl.gov}. These files can be
+directly converted to \nekpp format files by using the command
+%
+\begin{lstlisting}[style=BashInputStyle]
+  FieldConvert input.fld output.fld
+\end{lstlisting}
+%
+Note that even though the \inltt{.fld} extension is typically associated with
+\nekpp files, FieldConvert can automatically identify \emph{Semtex} and
+\emph{Nek5000} input field files.
+
+To use these files in a simulation, or to post-process the results of a
+simulation, an appropriate mesh must also be defined in the \nekpp XML format.
+\nm can be used to convert these input files to XML, as outlined in
+section~\ref{s:utilities:nekmesh}.
+
 \section{Convert .fld / .chk files into Paraview, VisIt or Tecplot format}
 \label{s:utilities:fieldconvert:sub:convert}
 To convert the \nekpp output binary files (.chk and .fld) into a

library/FieldUtils/CMakeLists.txt

@@ -7,6 +7,7 @@ SET(FieldUtilsHeaders
     InputModules/InputXml.h
     InputModules/InputPts.h
     InputModules/InputNek5000.h
+    InputModules/InputSemtex.h
     OutputModules/OutputInfo.h
     OutputModules/OutputTecplot.h
     OutputModules/OutputVtk.h
@@ -55,6 +56,7 @@ SET(FieldUtilsSources
     InputModules/InputXml.cpp
     InputModules/InputPts.cpp
     InputModules/InputNek5000.cpp
+    InputModules/InputSemtex.cpp
     OutputModules/OutputInfo.cpp
     OutputModules/OutputTecplot.cpp
     OutputModules/OutputVtk.cpp

library/FieldUtils/InputModules/InputNek5000.cpp

@@ -54,28 +54,6 @@ ModuleKey InputNek5000::m_className[1] = {
         "Reads Nek5000 field file.")
 };
 
-/**
- * @brief Swap endian ordering of the input variable.
- */
-template <typename T>
-void swap_endian(T &u)
-{
-    union
-    {
-        T u;
-        unsigned char u8[sizeof(T)];
-    } source, dest;
-
-    source.u = u;
-
-    for (size_t k = 0; k < sizeof(T); k++)
-    {
-        dest.u8[k] = source.u8[sizeof(T) - k - 1];
-    }
-
-    u = dest.u;
-}
-
 /**
  * @brief Set up InputNek5000 object.
  *

library/FieldUtils/InputModules/InputSemtex.cpp

+////////////////////////////////////////////////////////////////////////////////
+//
+//  File: InputSemtex.cpp
+//
+//  For more information, please see: http://www.nektar.info/
+//
+//  The MIT License
+//
+//  Copyright (c) 2006 Division of Applied Mathematics, Brown University (USA),
+//  Department of Aeronautics, Imperial College London (UK), and Scientific
+//  Computing and Imaging Institute, University of Utah (USA).
+//
+//  License for the specific language governing rights and limitations under
+//  Permission is hereby granted, free of charge, to any person obtaining a
+//  copy of this software and associated documentation files (the "Software"),
+//  to deal in the Software without restriction, including without limitation
+//  the rights to use, copy, modify, merge, publish, distribute, sublicense,
+//  and/or sell copies of the Software, and to permit persons to whom the
+//  Software is furnished to do so, subject to the following conditions:
+//
+//  The above copyright notice and this permission notice shall be included
+//  in all copies or substantial portions of the Software.
+//
+//  THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
+//  OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+//  FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
+//  THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+//  LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
+//  FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
+//  DEALINGS IN THE SOFTWARE.
+//
+//  Description: Reads a Semtex checkpoint file.
+//
+////////////////////////////////////////////////////////////////////////////////
+
+#include <iostream>
+#include <fstream>
+#include <string>
+#include <vector>
+using namespace std;
+
+#include <LibUtilities/BasicUtils/CompressData.h>
+#include <boost/algorithm/string.hpp>
+
+#include "InputSemtex.h"
+
+namespace Nektar
+{
+namespace FieldUtils
+{
+
+ModuleKey InputSemtex::m_className[1] = {
+    GetModuleFactory().RegisterCreatorFunction(
+        ModuleKey(eInputModule, "fldsem"), InputSemtex::create,
+        "Reads Semtex field file.")
+};
+
+/**
+ * @brief Set up InputSemtex object.
+ *
+ */
+InputSemtex::InputSemtex(FieldSharedPtr f) : InputModule(f)
+{
+    m_allowedFiles.insert("fldsem");
+}
+
+/**
+ *
+ */
+InputSemtex::~InputSemtex()
+{
+}
+
+/**
+ * @brief Process Semtex input file.
+ *
+ * This routine reads a binary-format Semtex field file, loads the data into
+ * memory and populates the field definitions to match the data format. Semtex
+ * is a classic nodal-Lagrangian spectral element code at a single polynomial
+ * order, meaning that the field data are set up according to this structure.
+ */
+void InputSemtex::Process(po::variables_map &vm)
+{
+    if (m_f->m_verbose)
+    {
+        if (m_f->m_comm->TreatAsRankZero())
+        {
+            cout << "Processing Semtex field file" << endl;
+        }
+    }
+
+    // Variables to be read from session file
+    string sessionName, date, fields, endian;
+    int nr, ns, nz, nelmt, step;
+    NekDouble time, dt, kinvis, beta;
+
+    ifstream file(m_f->m_inputfiles["fldsem"][0].c_str(), ios::binary);
+
+    // -- Read header information.
+    char buf[25];
+    string line;
+
+    // Session name
+    file.read(buf, 25);
+    sessionName = string(buf, 25);
+    boost::trim(sessionName);
+    getline(file, line);
+    m_f->m_fieldMetaDataMap["SessionName0"] = sessionName;
+
+    // Date
+    file.read(buf, 25);
+    date = string(buf, 25);
+    boost::trim(date);
+    getline(file, line);
+
+    // nP, nZ, nElmt
+    file >> nr >> ns >> nz >> nelmt;
+    getline(file, line);
+
+    // Step
+    file >> step;
+    getline(file, line);
+
+    // Time
+    file >> time;
+    getline(file, line);
+    m_f->m_fieldMetaDataMap["Time"] = boost::lexical_cast<string>(time);
+
+    // Timestep
+    file >> dt;
+    getline(file, line);
+    m_f->m_fieldMetaDataMap["TimeStep"] = boost::lexical_cast<string>(dt);
+
+    // Viscosity
+    file >> kinvis;
+    getline(file, line);
+    m_f->m_fieldMetaDataMap["Kinvis"] = boost::lexical_cast<string>(kinvis);
+
+    // Beta
+    file >> beta;
+    getline(file, line);
+
+    // Fields
+    file.read(buf, 25);
+    fields = string(buf, 25);
+    boost::trim(fields);
+    getline(file, line);
+
+    // Endian-ness
+    LibUtilities::EndianType systemEndian =  LibUtilities::Endianness();
+    std::string endianSearch;
+    if (systemEndian == LibUtilities::eEndianBig)
+    {
+        endianSearch = "big";
+    }
+    else if (systemEndian == LibUtilities::eEndianLittle)
+    {
+        endianSearch = "little";
+    }
+    else
+    {
+        ASSERTL0(false, "Only little- or big-endian systems are supported");
+    }
+
+    file.read(buf, 25);
+    endian = string(buf, 25);
+    bool byteSwap = endian.find(endianSearch) == string::npos;
+    getline(file, line);
+
+    // Print some basic information for input if in verbose mode.
+    if (m_f->m_verbose)
+    {
+        cout << "Found header information:" << endl;
+        cout << " -- From session         : " << sessionName << endl;
+        cout << " -- File generated       : " << date        << endl;
+        cout << " -- Polynomial order     : " << nr-1        << endl;
+        cout << " -- Number of planes     : " << nz          << endl;
+        cout << " -- Number of elements   : " << nelmt       << endl;
+        cout << " -- Simulation time      : " << time        << endl;
+        cout << " -- Timestep             : " << dt          << endl;
+        cout << " -- Viscosity            : " << kinvis      << endl;
+        cout << " -- Fields               : " << fields
+             << " (" << fields.size() << " total)" << endl;
+
+        if (nz > 1)
+        {
+            cout << " -- Homogeneous length   : " << 2*M_PI/beta << endl;
+        }
+
+        cout << " -- " << (byteSwap ? "" : "do not ") << "need to swap endian"
+             << endl;
+    }
+
+    ASSERTL0(nr == ns, "Semtex reader assumes values of nr and ns are equal");
+
+    // Set up a field definition
+    LibUtilities::FieldDefinitionsSharedPtr fielddef = MemoryManager<
+        LibUtilities::FieldDefinitions>::AllocateSharedPtr();
+    fielddef->m_shapeType         = LibUtilities::eQuadrilateral;
+    fielddef->m_homoStrips        = false;
+    fielddef->m_pointsDef         = false;
+    fielddef->m_uniOrder          = true;
+    fielddef->m_numPointsDef      = false;
+
+    // Set up basis
+    fielddef->m_basis.push_back(LibUtilities::eGLL_Lagrange);
+    fielddef->m_basis.push_back(LibUtilities::eGLL_Lagrange);
+    fielddef->m_numModes.push_back(nr);
+    fielddef->m_numModes.push_back(nr);
+
+    // Set up elements
+    fielddef->m_elementIDs.resize(nelmt);
+    for (int i = 0; i < nelmt; ++i)
+    {
+        fielddef->m_elementIDs[i] = i;
+    }
+
+    // Deal with homogeneous direction.
+    if (nz > 1)
+    {
+        fielddef->m_numHomogeneousDir = 1;
+        fielddef->m_homogeneousLengths.push_back(2 * M_PI / beta);
+        fielddef->m_numModes.push_back(nz);
+        fielddef->m_basis.push_back(LibUtilities::eFourier);
+
+        for (int i = 0; i < nz; ++i)
+        {
+            fielddef->m_homogeneousZIDs.push_back(i);
+        }
+    }
+    else
+    {
+        fielddef->m_numHomogeneousDir = 0;
+    }
+
+    for (string::size_type i = 0; i < fields.size(); ++i)
+    {
+        fielddef->m_fields.push_back(string(&fields[i], 1));
+    }
+
+    // Size of data to read.
+    size_t elmtSize  = nr * ns;
+    size_t planeSize = elmtSize * nelmt;
+    size_t fieldSize = planeSize * nz;
+    size_t dataSize = fieldSize * fields.size();
+
+    // Allocate our storage.
+    m_f->m_data.resize(1);
+    m_f->m_data[0].resize(dataSize);
+
+    // Temporary storage for one plane of data.
+    vector<NekDouble> tmp(planeSize);
+    size_t offset = nz * nr * ns;
+
+    // Now reorder data; Semtex ordering traverses memory fastest over planes,
+    // whereas Nektar++ expects it over elements
+    for (int i = 0; i < fields.size(); ++i)
+    {
+        NekDouble *data = &m_f->m_data[0][i * fieldSize];
+        for (int j = 0; j < nz; ++j)
+        {
+            size_t elSizeJ = j * elmtSize;
+            file.read((char *)&tmp[0], planeSize * sizeof(NekDouble));
+
+            if (byteSwap)
+            {
+                swap_endian(tmp);
+            }
+
+            for (int k = 0; k < nelmt; ++k)
+            {
+                std::copy(&tmp[k * elmtSize], &tmp[(k+1) * elmtSize],
+                          data + k * offset + elSizeJ);
+
+            }
+        }
+    }
+
+    m_f->m_fielddef.push_back(fielddef);
+
+    // If we have defined expansion lists, then create additional expansion
+    // lists.
+    if (m_f->m_exp.size())
+    {
+        int nfields = m_f->m_fielddef[0]->m_fields.size();
+        m_f->m_exp.resize(nfields);
+
+        for (string::size_type i = 0; i < fields.size(); ++i)
+        {
+            if (!m_f->m_exp[i])
+            {
+                m_f->m_exp[i] = m_f->AppendExpList(
+                    m_f->m_fielddef[0]->m_numHomogeneousDir);
+            }
+
+            m_f->m_exp[i]->ExtractDataToCoeffs(
+                m_f->m_fielddef[0],
+                m_f->m_data[0],
+                m_f->m_fielddef[0]->m_fields[i],
+                m_f->m_exp[i]->UpdateCoeffs());
+
+            m_f->m_exp[i]->BwdTrans(
+                m_f->m_exp[i]->GetCoeffs(), m_f->m_exp[i]->UpdatePhys());
+        }
+    }
+}
+
+}
+}

library/FieldUtils/InputModules/InputSemtex.h

+////////////////////////////////////////////////////////////////////////////////
+//
+//  File: InputSemtex.h
+//
+//  For more information, please see: http://www.nektar.info/
+//
+//  The MIT License
+//
+//  Copyright (c) 2006 Division of Applied Mathematics, Brown University (USA),
+//  Department of Aeronautics, Imperial College London (UK), and Scientific
+//  Computing and Imaging Institute, University of Utah (USA).
+//
+//  License for the specific language governing rights and limitations under
+//  Permission is hereby granted, free of charge, to any person obtaining a
+//  copy of this software and associated documentation files (the "Software"),
+//  to deal in the Software without restriction, including without limitation
+//  the rights to use, copy, modify, merge, publish, distribute, sublicense,
+//  and/or sell copies of the Software, and to permit persons to whom the
+//  Software is furnished to do so, subject to the following conditions:
+//
+//  The above copyright notice and this permission notice shall be included
+//  in all copies or substantial portions of the Software.
+//
+//  THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
+//  OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+//  FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
+//  THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+//  LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
+//  FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
+//  DEALINGS IN THE SOFTWARE.
+//
+//  Description: Reads a Semtex checkpoint file.
+//
+////////////////////////////////////////////////////////////////////////////////
+
+#ifndef FIELDUTILS_INPUTSEMTEX
+#define FIELDUTILS_INPUTSEMTEX
+
+#include "../Module.h"
+
+namespace Nektar
+{
+namespace FieldUtils
+{
+
+/**
+ * Converter for Fld files.
+ */
+class InputSemtex : public InputModule
+{
+public:
+    InputSemtex(FieldSharedPtr f);
+    virtual ~InputSemtex();
+    virtual void Process(po::variables_map &vm);
+
+    /// Creates an instance of this class
+    static ModuleSharedPtr create(FieldSharedPtr f)
+    {
+        return MemoryManager<InputSemtex>::AllocateSharedPtr(f);
+    }
+    /// %ModuleKey for class.
+    static ModuleKey m_className[];
+
+    virtual std::string GetModuleName()
+    {
+        return "InputSemtex";
+    }
+
+private:
+};
+}
+}
+
+#endif

library/FieldUtils/Interpolator.cpp

@@ -64,9 +64,6 @@ void Interpolator::CalcWeights(const LibUtilities::PtsFieldSharedPtr ptsInField,
     int nOutPts  = m_ptsOutField->GetNpoints();
     int lastProg = 0;
 
-    m_weights   = Array<OneD, Array<OneD, float> >(nOutPts);
-    m_neighInds = Array<OneD, Array<OneD, unsigned int> >(nOutPts);
-
     // set a default method
     if (m_method == eNoMethod)
     {
@@ -89,6 +86,9 @@ void Interpolator::CalcWeights(const LibUtilities::PtsFieldSharedPtr ptsInField,
     {
         case eNearestNeighbour:
         {
+            m_weights   = Array<TwoD, float>(nOutPts, 1, 0.0);
+            m_neighInds = Array<TwoD, unsigned int>(nOutPts, 1, (unsigned int) 0);
+
             for (int i = 0; i < nOutPts; ++i)
             {
                 Array<OneD, NekDouble> tmp(m_dim, 0.0);
@@ -116,6 +116,9 @@ void Interpolator::CalcWeights(const LibUtilities::PtsFieldSharedPtr ptsInField,
             ASSERTL0(m_ptsInField->GetDim() == 1 || m_coordId >= 0,
                      "not implemented");
 
+            m_weights   = Array<TwoD, float>(nOutPts, 3, 0.0);
+            m_neighInds = Array<TwoD, unsigned int>(nOutPts, 3, (unsigned int) 0);
+
             if (m_ptsInField->GetDim() == 1)
             {
                 m_coordId = 0;
@@ -152,6 +155,12 @@ void Interpolator::CalcWeights(const LibUtilities::PtsFieldSharedPtr ptsInField,
 
         case eShepard:
         {
+            int numPts = pow(double(2), m_ptsInField->GetDim());
+            numPts     = min(numPts, int(m_ptsInField->GetNpoints() / 2));
+
+            m_weights   = Array<TwoD, float>(nOutPts, numPts, 0.0);
+            m_neighInds = Array<TwoD, unsigned int>(nOutPts, numPts, (unsigned int) 0);
+
             for (int i = 0; i < nOutPts; ++i)
             {
                 Array<OneD, NekDouble> tmp(m_dim, 0.0);
@@ -161,7 +170,7 @@ void Interpolator::CalcWeights(const LibUtilities::PtsFieldSharedPtr ptsInField,
                 }
                 PtsPoint searchPt(i, tmp, 1E30);
 
-                CalcW_Shepard(searchPt);
+                CalcW_Shepard(searchPt, numPts);
 
                 int progress = int(100 * i / nOutPts);
                 if (m_progressCallback && progress > lastProg)
@@ -181,6 +190,11 @@ void Interpolator::CalcWeights(const LibUtilities::PtsFieldSharedPtr ptsInField,
             // use m_filtWidth as FWHM
             NekDouble sigma = m_filtWidth * 0.4246609001;
 
+            m_maxPts = min(m_maxPts, int(m_ptsInField->GetNpoints() / 2));
+
+            m_weights   = Array<TwoD, float>(nOutPts, m_maxPts, 0.0);
+            m_neighInds = Array<TwoD, unsigned int>(nOutPts, m_maxPts, (unsigned int) 0);
+