Merge branch 'feature/FC-streamfunction' into 'master'

Move StreamFunction utility to a FieldConvert module

See merge request !809

CHANGELOG.md

@@ -28,6 +28,7 @@ v5.0.0
 **FieldConvert**:
 - Add input module for Semtex field files (!777)
 - Fixed interppoints module (!760)
+- Move StreamFunction utility to a FieldConvert module (!809)
 
 **Documentation**:
 - Added the developer-guide repository as a submodule (!751)

docs/user-guide/utilities/fieldconvert.tex

@@ -178,6 +178,7 @@ possibly also Reynolds stresses) into single file;
 \item \inltt{scalargrad}: Computes scalar gradient field;
 \item \inltt{scaleinputfld}: Rescale input field by a constant factor;
 \item \inltt{shear}: Computes time-averaged shear stress metrics: TAWSS, OSI, transWSS, TAAFI, TACFI, WSSG;
+\item \inltt{streamfunction}: Calculates stream function of a 2D incompressible flow.
 \item \inltt{surfdistance}: Computes height of a prismatic boundary layer mesh and projects onto the surface (for e.g. $y^+$ calculation).
 \item \inltt{vorticity}: Computes the vorticity field.
 \item \inltt{wss}: Computes wall shear stress field.
@@ -873,6 +874,26 @@ The argument \inltt{N} and \inltt{fromfld} are compulsory arguments that respect
 
 The input \inltt{.fld} files are the outputs of the \textit{wss} module. If they do not contain the surface normals (an optional output of the \textit{wss} modle), then the \textit{shear} module will not compute the last metric, |WSSG|.
 
+%
+%
+%
+\subsection{Stream function of a 2D incompressible flow: \textit{streamfunction} module}
+
+The streamfunction module calculates the stream function of a 2D incompressible flow, by
+solving the Poisson equation
+\[
+\nabla^2 \psi = -\omega
+\]
+where $\omega$ is the vorticity. Note that this module applies the same boundary conditions
+specified for the y-direction velocity component \inltt{v} to the stream function,
+what may not be the most appropriate choice.
+
+To use this module, the user can run
+\begin{lstlisting}[style=BashInputStyle]
+FieldConvert -m streamfunction test.xml test.fld test-streamfunc.fld
+\end{lstlisting}
+where the file \inltt{test-streamfunc.fld} can be processed in a similar
+way as described in section \ref{s:utilities:fieldconvert:sub:convert}.
 
 %
 %

library/FieldUtils/CMakeLists.txt

@@ -41,6 +41,7 @@ SET(FieldUtilsHeaders
     ProcessModules/ProcessPointDataToFld.h
     ProcessModules/ProcessPrintFldNorms.h
     ProcessModules/ProcessScaleInFld.h
+    ProcessModules/ProcessStreamFunction.h
     ProcessModules/ProcessSurfDistance.h
     ProcessModules/ProcessVorticity.h
     ProcessModules/ProcessScalGrad.h
@@ -95,6 +96,7 @@ SET(FieldUtilsSources
     ProcessModules/ProcessScaleInFld.cpp
     ProcessModules/ProcessVorticity.cpp
     ProcessModules/ProcessScalGrad.cpp
+    ProcessModules/ProcessStreamFunction.cpp
     ProcessModules/ProcessSurfDistance.cpp
     ProcessModules/ProcessMultiShear.cpp
     ProcessModules/ProcessWSS.cpp

library/FieldUtils/ProcessModules/ProcessStreamFunction.cpp

+////////////////////////////////////////////////////////////////////////////////
+//
+//  File: ProcessStreamFunction.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: Computes stream-function for 2D field.
+//
+////////////////////////////////////////////////////////////////////////////////
+
+#include "ProcessStreamFunction.h"
+
+using namespace std;
+
+namespace Nektar
+{
+namespace FieldUtils
+{
+
+ModuleKey ProcessStreamFunction::className =
+    GetModuleFactory().RegisterCreatorFunction(
+        ModuleKey(eProcessModule, "streamfunction"),
+        ProcessStreamFunction::create,
+        "Computes stream-function for 2D field.");
+
+ProcessStreamFunction::ProcessStreamFunction(FieldSharedPtr f)
+    : ProcessModule(f)
+{
+    m_f->m_declareExpansionAsContField = true;
+}
+
+ProcessStreamFunction::~ProcessStreamFunction()
+{
+}
+
+void ProcessStreamFunction::Process(po::variables_map &vm)
+{
+    int nfields  = m_f->m_variables.size();
+    // Append field name
+    m_f->m_variables.push_back("StreamFunc");
+
+    // Skip in case of empty partition
+    if (m_f->m_exp[0]->GetNumElmts() == 0)
+    {
+        return;
+    }
+
+    // Check if dimension is correct
+    int expdim   = m_f->m_graph->GetMeshDimension();
+    int spacedim = expdim + m_f->m_numHomogeneousDir;
+    ASSERTL0(spacedim == 2, "Stream function can only be obtained in 2D.");
+
+    // Allocate arrays
+    int npoints = m_f->m_exp[0]->GetNpoints();
+    Array<OneD, NekDouble> vx(npoints);
+    Array<OneD, NekDouble> uy(npoints);
+    Array<OneD, NekDouble> vortNeg(npoints);
+
+    // Resize expansion
+    m_f->m_exp.resize(nfields + 1);
+    m_f->m_exp[nfields] = m_f->AppendExpList(m_f->m_numHomogeneousDir);
+
+    // Calculate vorticity: -W_z = Uy - Vx
+    m_f->m_exp[0]->PhysDeriv(MultiRegions::DirCartesianMap[1],
+                             m_f->m_exp[0]->GetPhys(),uy);
+    m_f->m_exp[1]->PhysDeriv(MultiRegions::DirCartesianMap[0],
+                             m_f->m_exp[1]->GetPhys(),vx);
+
+    Vmath::Vsub(npoints, uy, 1, vx, 1, vortNeg, 1);
+
+    // Calculate the Stream Function as the solution of the
+    //   Poisson equation: \nabla^2 StreamFunction = -Vorticity
+    //   Use same boundary conditions as v
+    StdRegions::ConstFactorMap factor;
+    factor[StdRegions::eFactorLambda] = 0.0;
+
+    m_f->m_exp[1]->HelmSolve(vortNeg,
+                            m_f->m_exp[nfields]->UpdateCoeffs(),
+                            NullFlagList, factor);
+    m_f->m_exp[1]->BwdTrans(m_f->m_exp[nfields]->GetCoeffs(),
+                            m_f->m_exp[nfields]->UpdatePhys());
+}
+
+}
+}

library/FieldUtils/ProcessModules/ProcessStreamFunction.h

+////////////////////////////////////////////////////////////////////////////////
+//
+//  File: ProcessStreamFunction.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: Computes stream-function for 2D field.
+//
+////////////////////////////////////////////////////////////////////////////////
+
+#ifndef FIELDUTILS_PROCESSSTREAMFUNCTION
+#define FIELDUTILS_PROCESSSTREAMFUNCTION
+
+#include "../Module.h"
+
+namespace Nektar
+{
+namespace FieldUtils
+{
+/**
+ * @brief This processing module calculates the stream function of a 2D field
+ * and adds it as an extra-field to the output file.
+ */
+class ProcessStreamFunction : public ProcessModule
+{
+public:
+    /// Creates an instance of this class
+    static boost::shared_ptr<Module> create(FieldSharedPtr f)
+    {
+        return MemoryManager<ProcessStreamFunction>::AllocateSharedPtr(f);
+    }
+    static ModuleKey className;
+
+    ProcessStreamFunction(FieldSharedPtr f);
+    virtual ~ProcessStreamFunction();
+
+    /// Write mesh to output file.
+    virtual void Process(po::variables_map &vm);
+
+    virtual std::string GetModuleName()
+    {
+        return "ProcessStreamFunction";
+    }
+
+    virtual std::string GetModuleDescription()
+    {
+        return "Calculating stream-function";
+    }
+
+    virtual ModulePriority GetModulePriority()
+    {
+        return eModifyExp;
+    }
+
+};
+}
+}
+
+#endif

utilities/Extras/CMakeLists.txt

 ADD_UTILITIES_EXECUTABLE(CheckXmlFile COMPONENT utilities-extra SOURCES CheckXmlFile.cpp)
-ADD_UTILITIES_EXECUTABLE(StreamFunction2D COMPONENT utilities-extra SOURCES StreamFunction2D.cpp)
 
 IF (NEKTAR_USE_VTK)
     ADD_UTILITIES_EXECUTABLE(VtkToFld COMPONENT utilities-extra SOURCES VtkToFld.cpp)

utilities/Extras/StreamFunction2D.cpp

-/**
- * This function calculate the vorticity vector starting from an .fld file.
- * It is meant to be used with solutions produced by the incompressible Navier-Stokes solver.
- * To use it with solutions coming form another solver further generalisations are required.
- */
-#include <cstdio>
-#include <cstdlib>
-
-#include <MultiRegions/ExpList.h>
-#include <MultiRegions/ExpList1D.h>
-#include <MultiRegions/ExpList2D.h>
-#include <MultiRegions/ExpList3D.h>
-#include <MultiRegions/ExpList2DHomogeneous1D.h>
-#include <MultiRegions/ExpList3DHomogeneous1D.h>
-#include <MultiRegions/ExpList3DHomogeneous2D.h>
-
-
-#include <SolverUtils/EquationSystem.h>
-
-#include <LocalRegions/MatrixKey.h>
-#include <LibUtilities/BasicUtils/Equation.h>
-#include <MultiRegions/ContField1D.h>
-#include <MultiRegions/ContField2D.h>
-#include <MultiRegions/ContField3D.h>
-#include <MultiRegions/ContField3DHomogeneous1D.h>
-#include <MultiRegions/ContField3DHomogeneous2D.h>
-
-using namespace std;
-using namespace Nektar;
-
-int main(int argc, char *argv[])
-{
-    int i,j;
-    
-    if(argc != 3)
-    {
-        fprintf(stderr,"Usage: ./StreamFunction2D file.xml file.fld\n");
-        exit(1);
-    }
-    
-    LibUtilities::SessionReaderSharedPtr vSession
-    = LibUtilities::SessionReader::CreateInstance(argc, argv);
-    
-    
-    //----------------------------------------------
-    // Read in mesh from input file
-    string meshfile(argv[argc-2]);
-    SpatialDomains::MeshGraphSharedPtr graphShPt = SpatialDomains::MeshGraph::Read(vSession);//meshfile);
-    //----------------------------------------------
-    
-    //----------------------------------------------
-    // Import field file.
-    string fieldfile(argv[argc-1]);
-    vector<LibUtilities::FieldDefinitionsSharedPtr> fielddef;
-    vector<vector<NekDouble> > fielddata;
-    LibUtilities::Import(fieldfile,fielddef,fielddata);
-    
-    //----------------------------------------------
-    // Define Expansion
-    int expdim  = graphShPt->GetMeshDimension();
-    int nfields = fielddef[0]->m_fields.size();
-    int vorticitydim;
-
-    vorticitydim = 2;
-    
-    
-    
-    Array<OneD, MultiRegions::ExpListSharedPtr> Exp(nfields + vorticitydim);
-    
-    Array<OneD, MultiRegions::ExpListSharedPtr> FieldVar(1);
-    
-    switch(expdim)
-    {
-
-        case 2:
-        {
-
-            {
-                MultiRegions::ExpList2DSharedPtr Exp2D;
-                Exp2D = MemoryManager<MultiRegions::ExpList2D>
-                ::AllocateSharedPtr(vSession,graphShPt);
-                Exp[0] =  Exp2D;
-                
-                FieldVar[0] = MemoryManager<MultiRegions::ContField2D>::AllocateSharedPtr(vSession, graphShPt, "v");
-                                
-                for(i = 1; i < nfields + vorticitydim; ++i)
-                {
-                    Exp[i] = MemoryManager<MultiRegions::ExpList2D>
-                    ::AllocateSharedPtr(*Exp2D);
-                }
-            }
-        }
-        break;
-
-        default:
-            ASSERTL0(false,"The input file must be two-dimensional");
-            break;
-    }
-    
-    //----------------------------------------------
-    // Copy data from field file
-    for(j = 0; j < nfields; ++j)
-    {
-        for(unsigned int i = 0; i < fielddata.size(); ++i)
-        {
-            Exp[j]->ExtractDataToCoeffs(fielddef [i],
-                                        fielddata[i],
-                                        fielddef [i]->m_fields[j],
-                                        Exp[j]->UpdateCoeffs());
-        }
-        Exp[j]->BwdTrans(Exp[j]->GetCoeffs(),Exp[j]->UpdatePhys());
-    }
-    //----------------------------------------------
-    
-    int nq = Exp[0]->GetNpoints();
-    
-    Array<OneD, NekDouble> Uy(nq);
-    Array<OneD, NekDouble> Vx(nq);
-
-    Array<OneD, NekDouble> StreamFunc(nq);
-    Array<OneD, NekDouble> Qz(nq);
-    
-    switch(expdim)
-    {
-
-        case 2:
-        {
-
-            {
-                Exp[0]->PhysDeriv(MultiRegions::DirCartesianMap[0],Exp[1]->GetPhys(),Vx);
-                Exp[0]->PhysDeriv(MultiRegions::DirCartesianMap[1],Exp[0]->GetPhys(),Uy);
-                
-                Vmath::Vsub(nq,Vx,1,Uy,1,Qz,1);
-                
-                //The Vorticity is stored
-                Exp[3]->FwdTrans(Qz, Exp[3]->UpdateCoeffs());
-                
-                //We now calculate the Stream Function as the solution of the
-                //Poisson equation: Vorticity = - \nabla^2 StreamFunction
-                StdRegions::ConstFactorMap factor;
-                factor[StdRegions::eFactorLambda] = 0.0;
-                
-                Vmath::Smul(nq,-1.0,Qz,1,Qz,1);
-                Exp[4]->SetPhys(Qz);
-                
-                FieldVar[0]->HelmSolve(Qz, Exp[4]->UpdateCoeffs(), NullFlagList, factor);
-            }
-        }
-        break;
-
-        default:
-        {
-            ASSERTL0(false,"The input file must be two-dimensional");
-        }
-        break;
-    }
-    
-    //-----------------------------------------------
-    // Write solution to file with additional computed fields
-    string   fldfilename(argv[2]);
-    string   out = fldfilename.substr(0, fldfilename.find_last_of("."));
-    string   endfile("_with_2DstremFunction.fld");
-    out += endfile;
-    std::vector<LibUtilities::FieldDefinitionsSharedPtr> FieldDef
-    = Exp[0]->GetFieldDefinitions();
-    std::vector<std::vector<NekDouble> > FieldData(FieldDef.size());
-    
-    for(j = 0; j < nfields + vorticitydim; ++j)
-    {
-        for(i = 0; i < FieldDef.size(); ++i)
-        {
-            if (j >= nfields)
-            {
-                if(j == 4)
-                {
-                    FieldDef[i]->m_fields.push_back("StreamFunc");
-                }
-                else
-                {
-                    FieldDef[i]->m_fields.push_back("Qz");
-                }
-            }
-            else
-            {
-                FieldDef[i]->m_fields.push_back(fielddef[i]->m_fields[j]);
-            }
-            Exp[j]->AppendFieldData(FieldDef[i], FieldData[i]);
-        }
-    }
-    LibUtilities::Write(out, FieldDef, FieldData);
-    //-----------------------------------------------
-    
-    return 0;
-}
-