Updates to introduce part-only options and make inputXml more efficient

docs/user-guide/utilities/fieldconvert.tex

@@ -662,12 +662,16 @@ field files placed within the directory.
 %
 %
 %
-\section{Processing large files}
+\section{Processing large files in serial}
 When processing large files it is not always convenient to run in parallel 
 but process each parallel partition in serial, for example when interpolating 
-one field to another. To do this we can use the \inltt{--nprocs} and 
-\inltt{--procid} options. For example the following option will interpolate 
-partition 2 of a decomposition into 10 partitions of \inltt{fiile2.xml} from 
+one field to another.
+
+\subsection{Using the nprocs and procid options}
+
+One option is to use the \inltt{--nprocs} and \inltt{--procid}
+options. For example the following option will interpolate partition 2
+of a decomposition into 10 partitions of \inltt{fiile2.xml} from
 \inltt{file1.fld}
 \begin{lstlisting}[style=BashInputStyle] 
 FieldConvert --nprocs 10 --procid 2 \
@@ -700,6 +704,34 @@ tell FieldConvert that you wish to generate an info file but the extension
 to this file is .xml. This syntax allows the routine not to get confused with 
 the input/output xml files.
 
+\subsection{ Using the --part-only and --part-only-overlapping options}
+
+Another approach to serially proessing a large file is to intiially process the file into multiple partitions. This can be done with the \inltt{part-only} option. So the command
+\begin{lstlisting}[style=BashInputStyle] 
+FieldConvert --part-only 10 file.xml file.fld
+\end{lstlisting}
+will call the mesh partitioner for 10 paritions and write the
+partitions into a director called file\_xml. If you enter this
+directory you will find parititoned xml files \inltt{P0000000.xml},
+\inltt{P0000001.xml},...,\inltt{P0000009.xml} which can then be
+processed individually.
+
+There is also a \inltt{part-only-overlapping} option
+\begin{lstlisting}[style=BashInputStyle] 
+FieldConvert --part-only-overlapping 10 file.xml file.fld
+\end{lstlisting}
+which will partition the mesh into 10 overlapping partitions which can
+then be processed individually. This is sometime helpful when for
+example trying to produce a global isocontour which has been
+smoothed. Applying the smoothed isocontour extraction routine with the
+\inltt{part-only} option will produce a series of isocontour where
+there will be a gap between partitions as the smoother tends to shrink
+the isocontour within a partition. using the
+\inltt{part-only-overlapping} will still have shrinking isocnotour but
+the overlapping nature of these isocontour can be helpful when trying
+to mask the partiiton boundaries.
+
+
 %%% Local Variables:
 %%% mode: latex
 %%% TeX-master: "../user-guide"

library/LibUtilities/BasicUtils/MeshPartition.cpp

@@ -93,7 +93,8 @@ namespace Nektar
 
         }
 
-        void MeshPartition::PartitionMesh(int nParts, bool shared)
+        void MeshPartition::PartitionMesh(int nParts, bool shared,
+                                          bool overlapping)
         {
             ASSERTL0(m_meshElements.size() >= nParts,
                      "Too few elements for this many processes.");
@@ -104,7 +105,7 @@ namespace Nektar
                 WeightElements();
             }
             CreateGraph(m_mesh);
-            PartitionGraph(m_mesh, nParts, m_localPartition);
+            PartitionGraph(m_mesh, nParts, m_localPartition, overlapping);
         }
 
         void MeshPartition::WriteLocalPartition(LibUtilities::SessionReaderSharedPtr& pSession)
@@ -1143,7 +1144,9 @@ namespace Nektar
 
         void MeshPartition::PartitionGraph(BoostSubGraph& pGraph,
                                            int nParts,
-                                           std::vector<BoostSubGraph>& pLocalPartition)
+                                           std::vector<BoostSubGraph>& pLocalPartition,
+                                           bool overlapping)
+
         {
             int i;
             int nGraphVerts = boost::num_vertices(pGraph);
@@ -1151,6 +1154,7 @@ namespace Nektar
 
             // Convert boost graph into CSR format
             BoostVertexIterator    vertit, vertit_end;
+            BoostAdjacencyIterator adjvertit, adjvertit_end;
             Array<OneD, int> part(nGraphVerts,0);
 
             if (m_comm->GetRowComm()->TreatAsRankZero())
@@ -1256,6 +1260,25 @@ namespace Nektar
                 pGraph[*vertit].partition = part[i];
                 boost::add_vertex(i, pLocalPartition[part[i]]);
             }
+
+            if(overlapping)
+            {
+                // 
+                for ( boost::tie(vertit, vertit_end) = boost::vertices(pGraph);
+                      vertit != vertit_end;
+                      ++vertit)
+                {
+                    for (boost::tie(adjvertit, adjvertit_end) = boost::adjacent_vertices(*vertit,pGraph);
+                         adjvertit != adjvertit_end; ++adjvertit)
+                    {
+                        if(part[*adjvertit] != part[*vertit])
+                        {
+                            boost::add_vertex(*adjvertit, pLocalPartition[part[*vertit]]);
+                            
+                        }
+                    }
+                }
+            }
         }
 
 

library/LibUtilities/BasicUtils/MeshPartition.h

@@ -66,7 +66,7 @@ namespace Nektar
             LIB_UTILITIES_EXPORT MeshPartition(const SessionReaderSharedPtr& pSession);
             LIB_UTILITIES_EXPORT virtual ~MeshPartition();
 
-            LIB_UTILITIES_EXPORT void PartitionMesh(int nParts, bool shared = false);
+            LIB_UTILITIES_EXPORT void PartitionMesh(int nParts, bool shared = false, bool overlapping = false);
             LIB_UTILITIES_EXPORT void WriteLocalPartition(
                     SessionReaderSharedPtr& pSession);
             LIB_UTILITIES_EXPORT void WriteAllPartitions(
@@ -213,7 +213,8 @@ namespace Nektar
             void CreateGraph(BoostSubGraph& pGraph);
             void PartitionGraph(BoostSubGraph& pGraph,
                                 int nParts,
-                                std::vector<BoostSubGraph>& pLocalPartition);
+                                std::vector<BoostSubGraph>& pLocalPartition,
+                                bool overlapping = false);
 
             virtual void PartitionGraphImpl(
                     int&                              nVerts,

library/LibUtilities/BasicUtils/SessionReader.cpp

@@ -338,6 +338,8 @@ namespace Nektar
                                  "number of slices in Z-dir")
                 ("part-only",    po::value<int>(),
                                  "only partition mesh into N partitions.")
+                ("part-only-overlapping",    po::value<int>(),
+                                 "only partition mesh into N overlapping partitions.")
                 ("part-info",    "Output partition information")
             ;
 
@@ -1600,7 +1602,8 @@ namespace Nektar
             // Mesh has not been partitioned so do partitioning if required.
             // Note in the serial case nothing is done as we have already loaded
             // the mesh.
-            if (DefinesCmdLineArgument("part-only"))
+            if (DefinesCmdLineArgument("part-only")||
+                DefinesCmdLineArgument("part-only-overlapping"))
             {
                 // Perform partitioning of the mesh only. For this we insist
                 // the code is run in serial (parallel execution is pointless).
@@ -1608,12 +1611,21 @@ namespace Nektar
                         "The 'part-only' option should be used in serial.");
 
                 // Number of partitions is specified by the parameter.
-                int nParts = GetCmdLineArgument<int>("part-only");
+                int nParts; 
                 SessionReaderSharedPtr vSession     = GetSharedThisPtr();
                 MeshPartitionSharedPtr vPartitioner =
                         GetMeshPartitionFactory().CreateInstance(
                                             vPartitionerName, vSession);
-                vPartitioner->PartitionMesh(nParts, true);
+                if(DefinesCmdLineArgument("part-only"))
+                {
+                    nParts = GetCmdLineArgument<int>("part-only");
+                    vPartitioner->PartitionMesh(nParts, true);
+                }
+                else  
+                {
+                    nParts = GetCmdLineArgument<int>("part-only-overlapping");
+                    vPartitioner->PartitionMesh(nParts, true, true);
+                }
                 vPartitioner->WriteAllPartitions(vSession);
                 vPartitioner->GetCompositeOrdering(m_compOrder);
                 vPartitioner->GetBndRegionOrdering(m_bndRegOrder);

utilities/FieldConvert/FieldConvert.cpp

@@ -67,6 +67,10 @@ int main(int argc, char* argv[])
         ("onlyshape", po::value<string>(),
                  "Only use element with defined shape type i.e. -onlyshape "
                  " Tetrahedron")
+        ("part-only", po::value<int>(),
+                "Partition into specfiied npart partitions and exit")
+        ("part-only-overlapping", po::value<int>(),
+                "Partition into specfiied npart overlapping partitions and exit")
         ("procid", po::value<int>(),
                 "Process as single procid of a partition of size nproc "
                 "(-nproc must be specified).")
@@ -235,6 +239,15 @@ int main(int argc, char* argv[])
     modcmds.push_back(*(inout.end()-1));
     int nInput = inout.size()-1;
 
+    // For special case of part-only or part-only-overlapping options
+    // only require a single input file and so reset the nInputs to be
+    // of size inout.size(). Since the code will exit before reaching
+    // any output module this seems to work as expected
+    if(vm.count("part-only")||vm.count("part-only-overlapping"))
+    {
+        nInput = inout.size();
+    }
+
     InputModuleSharedPtr inputModule;
 
     for (int i = 0; i < modcmds.size(); ++i)

utilities/FieldConvert/InputModules/InputXml.cpp

@@ -84,11 +84,28 @@ InputXml::~InputXml()
  */
 void InputXml::Process(po::variables_map &vm)
 {
+    Timer timer, timerpart;
+
+    //check for multiple calls to inputXml due to split xml
+    //files. If so just return
+    int expsize = m_f->m_exp.size();
+    m_f->m_comm->AllReduce(expsize,LibUtilities::ReduceMax);
+    
+    if(expsize != 0)
+    {
+        return; 
+    }
 
     if(m_f->m_verbose)
     {
-        cout << "Processing input xml file" << endl;
+        if(m_f->m_comm->GetRank() == 0)
+        {
+            cout << "Processing input xml file" << endl;
+            timer.Start();
+            timerpart.Start();
+        }
     }
+
     // check to see if fld file defined so can use in
     // expansion defintion if required
     string fldending;
@@ -207,31 +224,95 @@ void InputXml::Process(po::variables_map &vm)
     }
 
 
+    // set up command lines options for session read to pass through
+    string firstarg = "FieldConvert";
+    int argc = 1;
+    char **argv;
+    argv = (char**)malloc(6*sizeof(char*));
+
+    argv[0] = (char *)malloc(firstarg.size()*sizeof(char));
+    sprintf(argv[0],"%s",firstarg.c_str());
+    
     if(m_f->m_verbose)
     {
-        string firstarg = "FieldConvert";
-        string verbose = "-v";
-        char **argv;
-        argv = (char**)malloc(2*sizeof(char*));
-        argv[0] = (char *)malloc(firstarg.size()*sizeof(char));
-        argv[1] = (char *)malloc(verbose.size()*sizeof(char));
-
-        sprintf(argv[0],"%s",firstarg.c_str());
-        sprintf(argv[1],"%s",verbose.c_str());
-
-        m_f->m_session = LibUtilities::SessionReader::
-            CreateInstance(2, (char **)argv, files, m_f->m_comm);
+        string verbose = "--verbose";
+        argv[argc] = (char *)malloc(verbose.size()*sizeof(char));
+
+        sprintf(argv[argc],"%s",verbose.c_str());
+        ++argc;
     }
-    else
+
+    if(vm.count("shared-filesystem"))
+    {
+        string sharedfs = "--shared-filesystem";
+        argv[argc] = (char *)malloc(sharedfs.size()*sizeof(char));
+        sprintf(argv[argc],"%s",sharedfs.c_str());
+        argc ++;
+    }
+
+    if(vm.count("part-only"))
+    {
+        string argstr = "--part-only";
+        argv[argc] = (char *)malloc(argstr.size()*sizeof(char));
+        sprintf(argv[argc],"%s",argstr.c_str());
+        argc ++;
+        
+        argstr = boost::lexical_cast<string>(vm["part-only"].as<int>());
+        argv[argc] = (char *)malloc(argstr.size()*sizeof(char));
+        sprintf(argv[argc],"%s",argstr.c_str());
+        argc ++;
+    }
+    
+    if(vm.count("part-only-overlapping"))
     {
-        m_f->m_session = LibUtilities::SessionReader::
-            CreateInstance(0, 0, files, m_f->m_comm);
+        string argstr = "--part-only-overlapping";
+        argv[argc] = (char *)malloc(argstr.size()*sizeof(char));
+        sprintf(argv[argc],"%s",argstr.c_str());
+        argc ++;
+
+        argstr = boost::lexical_cast<string>(vm["part-only-overlapping"].as<int>());
+        argv[argc] = (char *)malloc(argstr.size()*sizeof(char));
+        sprintf(argv[argc],"%s",argstr.c_str());
+        argc ++;
     }
 
+    m_f->m_session = LibUtilities::SessionReader::
+        CreateInstance(argc, (char **) argv, files, m_f->m_comm);
+
+    if(m_f->m_verbose)
+    {
+        if(m_f->m_comm->GetRank() == 0)
+        {
+            timerpart.Stop();
+            NekDouble cpuTime = timerpart.TimePerTest(1);
+            
+            stringstream ss;
+            ss << cpuTime << "s";
+            cout << "\t InputXml session reader CPU Time: " << setw(8) << left
+                 << ss.str() << endl;
+            timerpart.Start();
+        }
+    }
+    
     m_f->m_graph = SpatialDomains::MeshGraph::Read(m_f->m_session,rng);
     m_f->m_fld = MemoryManager<LibUtilities::FieldIO>
                     ::AllocateSharedPtr(m_f->m_session->GetComm());
 
+    if(m_f->m_verbose)
+    {
+        if(m_f->m_comm->GetRank() == 0)
+        {
+            timerpart.Stop();
+            NekDouble cpuTime = timerpart.TimePerTest(1);
+            
+            stringstream ss;
+            ss << cpuTime << "s";
+            cout << "\t InputXml mesh graph setup  CPU Time: " << setw(8) << left
+                 << ss.str() << endl;
+            timerpart.Start();
+        }
+    }
+
     // currently load all field (possibly could read data from
     // expansion list but it is re-arranged in expansion)
     const SpatialDomains::ExpansionMap &expansions = m_f->m_graph->GetExpansions();
@@ -251,7 +332,22 @@ void InputXml::Process(po::variables_map &vm)
     int NumHomogeneousDir = 0;
     if(fldfilegiven)
     {
-        m_f->m_fld->Import(m_f->m_inputfiles[fldending][0],m_f->m_fielddef);
+        // use original expansion to identify which elements are in
+        // this partition/subrange
+
+        Array<OneD,int> ElementGIDs(expansions.size());
+        SpatialDomains::ExpansionMap::const_iterator expIt;
+
+        int i = 0; 
+        for (expIt = expansions.begin(); expIt != expansions.end(); ++expIt)
+        {
+            ElementGIDs[i++] = expIt->second->m_geomShPtr->GetGlobalID();
+        }
+        
+        m_f->m_fld->Import(m_f->m_inputfiles[fldending][0],m_f->m_fielddef,
+                           LibUtilities::NullVectorNekDoubleVector,
+                           LibUtilities::NullFieldMetaDataMap,
+                           ElementGIDs);
         NumHomogeneousDir = m_f->m_fielddef[0]->m_numHomogeneousDir;
 
         //----------------------------------------------
@@ -291,6 +387,21 @@ void InputXml::Process(po::variables_map &vm)
         m_f->m_graph->SetExpansionsToEvenlySpacedPoints(nPointsNew);
     }
 
+    if(m_f->m_verbose)
+    {
+        if(m_f->m_comm->GetRank() == 0)
+        {
+            timerpart.Stop();
+            NekDouble cpuTime = timerpart.TimePerTest(1);
+            
+            stringstream ss;
+            ss << cpuTime << "s";
+            cout << "\t InputXml setexpansion CPU Time: " << setw(8) << left
+                 << ss.str() << endl;
+            timerpart.Start();
+        }
+    }
+
     // Override number of planes with value from cmd line
     if(NumHomogeneousDir == 1 && vm.count("output-points-hom-z"))
     {
@@ -299,7 +410,31 @@ void InputXml::Process(po::variables_map &vm)
     }
 
     m_f->m_exp[0] = m_f->SetUpFirstExpList(NumHomogeneousDir,fldfilegiven);
-}
+    
+    if(m_f->m_verbose)
+    {
+        if(m_f->m_comm->GetRank() == 0)
+        {
+            timerpart.Stop();
+            NekDouble cpuTime = timerpart.TimePerTest(1);
+            
+            stringstream ss1;
+
+            ss1 << cpuTime << "s";
+            cout << "\t InputXml set first exp CPU Time: " << setw(8) << left
+                 << ss1.str() << endl;
+
+            
+            timer.Stop();
+            cpuTime = timer.TimePerTest(1);
+            
+            stringstream ss;
+            ss << cpuTime << "s";
+            cout << "InputXml  CPU Time: " << setw(8) << left
+                 << ss.str() << endl;
 
+        }
+    }
+}
 }
 }

utilities/FieldConvert/Module.h

@@ -39,6 +39,7 @@
 #include <string>
 #include <map>
 #include <iostream>
+#include <iomanip>
 #include <fstream>
 #include <vector>
 #include <set>
@@ -46,6 +47,7 @@
 #include <LibUtilities/BasicUtils/NekFactory.hpp>
 #include <StdRegions/StdNodalTriExp.h>
 #include <LibUtilities/Communication/CommSerial.h>
+#include <LibUtilities/BasicUtils/Timer.h>
 
 #include "Field.hpp"