Merge branch 'feature/ExpFromFile' into 'master'

Expansion from restart file and improved partitioning

This MR makes two changes which might be useful for simulations with variable polynomial order:
- Allow the expansions to be loaded directly from the field definitions in a restart file (issue #14), using:

```xml
<EXPANSIONS>
     <F FILE="restartfile.fld" />
</EXPANSIONS>
```
This is useful when using DriverAdaptive, since in that case the final distribution of polynomial orders is only available in the field file.

- Introduce new options for load balancing in the mesh partitioning, through the solver info `WEIGHTPARTITIONS`. Possible values are `DOF` (load balanced by number of degrees of freedom in each element), `BOUNDARY` (load balanced by number of boundary degrees of freedom in each element) or `BOTH`. By default, the current behaviour of using no weights is maintained. 


See merge request !617

CHANGELOG.md

@@ -21,6 +21,8 @@ v4.4.0
 - Add support for HDF5 as an alternative output to XML-based output, including
   refactoring of FieldIO, improvements to MPI interface and added communicators
   to boundary conditions (!615)
+- Allow expansions to be loaded directly from field file (!617)
+- New options for load balancing (DOF or BOUNDARY) in mesh partitioner (!617)
 
 **ADRSolver:**
 - Add a projection equation system for C^0 projections (!675)

docs/user-guide/solvers/incompressible-ns.tex

@@ -1107,6 +1107,20 @@ element is taken as the average of the function in the quadrature
 points of the element. If these functions are defined, the values set
 for the tolerances in the \texttt{PARAMETERS} section are ignored.
 
+\subsection{Restarting the simulation}
+
+The simulation can be restarted using the final distribution of
+polynomial orders obtained from the adaptive procedure by setting
+the expansions as
+\begin{lstlisting}[style=XMLStyle]
+<EXPANSIONS>
+  <F FILE="restartfile.fld" />
+</EXPANSIONS>
+\end{lstlisting}
+note that this will only affect the polynomial order. The initial condition
+still needs to be set correctly, and does not need to come from the same file
+used for the expansions.
+
 \section{Advecting extra passive scalar fields}
 
 In some cases, it might be useful to advect passive scalar fields with the velocity obtained from the

library/LibUtilities/BasicUtils/MeshPartition.h

@@ -131,6 +131,8 @@ namespace Nektar
                 int id;             ///< Universal ID of the vertex
                 int partition;      ///< Index of the partition to which it belongs
                 MultiWeight weight; ///< Weightings to this graph vertex
+                MultiWeight bndWeight;
+                MultiWeight edgeWeight;
             };
 
             // Face/Edge/Vertex between two adjacent elements
@@ -189,12 +191,17 @@ namespace Nektar
             std::vector<unsigned int>           m_domain;
             std::map<std::string, std::string>  m_vertexAttributes;
 
-            // hierarchial mapping: composite id -> field name -> integer list
+            // hierarchial mapping: elmt id -> field name -> integer list
             // of directional nummodes described by expansion type clause.
             std::map<int, NummodesPerField>     m_expansions;
 
+            // map of each elements shape
+            std::map<int, char>                 m_shape;
+
             std::map<std::string, int>          m_fieldNameToId;
             std::map<int, MultiWeight>          m_vertWeights;
+            std::map<int, MultiWeight>          m_vertBndWeights;
+            std::map<int, MultiWeight>          m_edgeWeights;
 
             BndRegionOrdering                   m_bndRegOrder;
 
@@ -204,6 +211,8 @@ namespace Nektar
             CommSharedPtr                       m_comm;
 
             bool                                m_weightingRequired;
+            bool                                m_weightBnd;
+            bool                                m_weightDofs;
             bool                                m_shared;
 
             void ReadExpansions(const SessionReaderSharedPtr& pSession);
@@ -223,6 +232,7 @@ namespace Nektar
                     Nektar::Array<Nektar::OneD, int>& adjcy,
                     Nektar::Array<Nektar::OneD, int>& vertWgt,
                     Nektar::Array<Nektar::OneD, int>& vertSize,
+                    Nektar::Array<Nektar::OneD, int>& edgeWgt,
                     int&                              nparts,
                     int&                              volume,
                     Nektar::Array<Nektar::OneD, int>& part) = 0;
@@ -230,6 +240,7 @@ namespace Nektar
             void OutputPartition(SessionReaderSharedPtr& pSession, BoostSubGraph& pGraph, TiXmlElement* pGeometry);
             void CheckPartitions(int nParts, Array<OneD, int> &pPart);
             int CalculateElementWeight(char elmtType, bool bndWeight, int na, int nb, int nc);
+            int CalculateEdgeWeight(char elmtType, int na, int nb, int nc);
         };
 
         typedef boost::shared_ptr<MeshPartition> MeshPartitionSharedPtr;

library/LibUtilities/BasicUtils/MeshPartitionMetis.cpp

@@ -70,11 +70,12 @@ namespace LibUtilities
             Nektar::Array<Nektar::OneD, int>& adjcy,
             Nektar::Array<Nektar::OneD, int>& vertWgt,
             Nektar::Array<Nektar::OneD, int>& vertSize,
+            Nektar::Array<Nektar::OneD, int>& edgeWgt,
             int&                              nparts,
             int&                              volume,
             Nektar::Array<Nektar::OneD, int>& part)
     {
-        Metis::PartGraphVKway(nVerts, nVertConds, xadj, adjcy, vertWgt, vertSize, nparts, volume, part);
+        Metis::PartGraphVKway(nVerts, nVertConds, xadj, adjcy, vertWgt, vertSize, edgeWgt, nparts, volume, part);
     }
 }
 }

library/LibUtilities/BasicUtils/MeshPartitionMetis.h

@@ -68,6 +68,7 @@ namespace LibUtilities
                     Nektar::Array<Nektar::OneD, int>& adjcy,
                     Nektar::Array<Nektar::OneD, int>& vertWgt,
                     Nektar::Array<Nektar::OneD, int>& vertSize,
+                    Nektar::Array<Nektar::OneD, int>& edgeWgt,
                     int&                              nparts,
                     int&                              volume,
                     Nektar::Array<Nektar::OneD, int>& part);

library/LibUtilities/BasicUtils/MeshPartitionScotch.cpp

@@ -70,6 +70,7 @@ namespace LibUtilities
             Nektar::Array<Nektar::OneD, int>& adjcy,
             Nektar::Array<Nektar::OneD, int>& vertWgt,
             Nektar::Array<Nektar::OneD, int>& vertSize,
+            Nektar::Array<Nektar::OneD, int>& edgeWgt,
             int&                              nparts,
             int&                              volume,
             Nektar::Array<Nektar::OneD, int>& part)

library/LibUtilities/BasicUtils/MeshPartitionScotch.h

@@ -72,6 +72,7 @@ namespace LibUtilities
                     Nektar::Array<Nektar::OneD, int>& adjcy,
                     Nektar::Array<Nektar::OneD, int>& vertWgt,
                     Nektar::Array<Nektar::OneD, int>& vertSize,
+                    Nektar::Array<Nektar::OneD, int>& edgeWgt,
                     int&                              nparts,
                     int&                              volume,
                     Nektar::Array<Nektar::OneD, int>& part);

library/LibUtilities/BasicUtils/Metis.hpp

@@ -75,12 +75,14 @@ namespace Metis
             Nektar::Array<Nektar::OneD, int>& adjcy,
             Nektar::Array<Nektar::OneD, int>& vertWgt,
             Nektar::Array<Nektar::OneD, int>& vertSize,
+            Nektar::Array<Nektar::OneD, int>& edgeWgt,
             int&                              nparts,
             int&                              volume,
             Nektar::Array<Nektar::OneD, int>& part)
     {
-        int *vwgt = 0;
-        int *vsize = 0;
+        int *vwgt   = 0;
+        int *vsize  = 0;
+        int *adjwgt = 0;
         if (vertWgt.num_elements() > 0)
         {
             vwgt = &vertWgt[0];
@@ -89,12 +91,16 @@ namespace Metis
         {
             vsize = &vertSize[0];
         }
+        if (edgeWgt.num_elements() > 0)
+        {
+            adjwgt = &edgeWgt[0];
+        }
         // number of balancing conditions (size of vertex multi-weight)
         int ncon = nVertConds;
         int options[METIS_NOPTIONS];
         METIS_SetDefaultOptions(options);
         METIS_PartGraphKway(&nVerts, &ncon, &xadj[0], &adjcy[0], vwgt, vsize,
-                            0, &nparts, 0, 0, options, &volume, &part[0]);
+                            adjwgt, &nparts, 0, 0, options, &volume, &part[0]);
     }
 }
 #endif //NEKTAR_LIB_UTILITIES_BASICUTILS_METIS_HPP

library/SpatialDomains/MeshGraph.cpp

@@ -1057,6 +1057,22 @@ namespace Nektar
                     cout << "    Number of elements: " << fielddefs.size() << endl;
                     SetExpansions(fielddefs);
                 }
+                else if(expType == "F")
+                {
+                    ASSERTL0(expansion->Attribute("FILE"),
+                                "Attribute FILE expected for type F expansion");
+                    std::string filenameStr = expansion->Attribute("FILE");
+                    ASSERTL0(!filenameStr.empty(),
+                                "A filename must be specified for the FILE "
+                                "attribute of expansion");
+
+                    std::vector<LibUtilities::FieldDefinitionsSharedPtr> fielddefs;
+                    LibUtilities::FieldIOSharedPtr f =
+                            LibUtilities::FieldIO::CreateForFile(
+                            m_session, filenameStr);
+                    f->Import(filenameStr, fielddefs);
+                    SetExpansions(fielddefs);
+                }
                 else
                 {
                     ASSERTL0(false,"Expansion type not defined");
@@ -2361,7 +2377,7 @@ namespace Nektar
                     std::string field = fielddef[i]->m_fields[j];
                     if(m_expansionMapShPtrMap.count(field) == 0)
                     {
-                        expansionMap = MemoryManager<ExpansionMap>::AllocateSharedPtr();
+                        expansionMap = SetUpExpansionMap();
                         m_expansionMapShPtrMap[field] = expansionMap;
 
                         // check to see if DefaultVar also not set and
@@ -2370,26 +2386,6 @@ namespace Nektar
                         {
                             m_expansionMapShPtrMap["DefaultVar"] = expansionMap;
                         }
-
-                        // loop over all elements and set expansion
-                        for(k = 0; k < fielddef.size(); ++k)
-                        {
-                            for(int h = 0; h < fielddef[k]->m_fields.size(); ++h)
-                            {
-                                if(fielddef[k]->m_fields[h] == field)
-                                {
-                                    expansionMap = m_expansionMapShPtrMap.find(field)->second;
-                                    LibUtilities::BasisKeyVector def;
-
-                                    for(int g = 0; g < fielddef[k]->m_elementIDs.size(); ++g)
-                                    {
-                                        ExpansionShPtr tmpexp =
-                                                MemoryManager<Expansion>::AllocateSharedPtr(geom, def);
-                                        (*expansionMap)[fielddef[k]->m_elementIDs[g]] = tmpexp;
-                                    }
-                                }
-                            }
-                        }
                     }
                 }
             }
@@ -2424,6 +2420,11 @@ namespace Nektar
                             if(m_segGeoms.count(fielddef[i]->m_elementIDs[j]) == 0)
                             {
                                 // skip element likely from parallel read
+                                if(!UniOrder)
+                                {
+                                    cnt++;
+                                    cnt += fielddef[i]->m_numHomogeneousDir;
+                                }
                                 continue;
                             }
                             geom = m_segGeoms[fielddef[i]->m_elementIDs[j]];
@@ -2461,6 +2462,11 @@ namespace Nektar
                             if(m_triGeoms.count(fielddef[i]->m_elementIDs[j]) == 0)
                             {
                                 // skip element likely from parallel read
+                                if(!UniOrder)
+                                {
+                                    cnt += 2;
+                                    cnt += fielddef[i]->m_numHomogeneousDir;
+                                }
                                 continue;
                             }
                             geom = m_triGeoms[fielddef[i]->m_elementIDs[j]];
@@ -2516,6 +2522,11 @@ namespace Nektar
                             if(m_quadGeoms.count(fielddef[i]->m_elementIDs[j]) == 0)
                             {
                                 // skip element likely from parallel read
+                                if(!UniOrder)
+                                {
+                                    cnt += 2;
+                                    cnt += fielddef[i]->m_numHomogeneousDir;
+                                }
                                 continue;
                             }
 
@@ -2561,6 +2572,10 @@ namespace Nektar
                             // parallel runs
                             if(m_tetGeoms.count(k) == 0)
                             {
+                                if(!UniOrder)
+                                {
+                                    cnt += 3;
+                                }
                                 continue;
                             }
                             geom = m_tetGeoms[k];
@@ -2647,6 +2662,10 @@ namespace Nektar
                             k = fielddef[i]->m_elementIDs[j];
                             if(m_prismGeoms.count(k) == 0)
                             {
+                                if(!UniOrder)
+                                {
+                                    cnt += 3;
+                                }
                                 continue;
                             }
                             geom = m_prismGeoms[k];
@@ -2770,6 +2789,10 @@ namespace Nektar
                             k = fielddef[i]->m_elementIDs[j];
                             if(m_hexGeoms.count(k) == 0)
                             {
+                                if(!UniOrder)
+                                {
+                                    cnt += 3;
+                                }
                                 continue;
                             }
 
@@ -2831,7 +2854,7 @@ namespace Nektar
                 std::vector<LibUtilities::FieldDefinitionsSharedPtr> &fielddef,
                 std::vector< std::vector<LibUtilities::PointsType> > &pointstype)
         {
-            int i,j,k,g,h,cnt,id;
+            int i,j,k,cnt,id;
             GeometrySharedPtr geom;
 
             ExpansionMapShPtr expansionMap;
@@ -2845,7 +2868,7 @@ namespace Nektar
                     std::string field = fielddef[i]->m_fields[j];
                     if(m_expansionMapShPtrMap.count(field) == 0)
                     {
-                        expansionMap = MemoryManager<ExpansionMap>::AllocateSharedPtr();
+                        expansionMap = SetUpExpansionMap();
                         m_expansionMapShPtrMap[field] = expansionMap;
 
                         // check to see if DefaultVar also not set and
@@ -2854,26 +2877,6 @@ namespace Nektar
                         {
                             m_expansionMapShPtrMap["DefaultVar"] = expansionMap;
                         }
-
-                        // loop over all elements and set expansion
-                        for(k = 0; k < fielddef.size(); ++k)
-                        {
-                            for(h = 0; h < fielddef[k]->m_fields.size(); ++h)
-                            {
-                                if(fielddef[k]->m_fields[h] == field)
-                                {
-                                    expansionMap = m_expansionMapShPtrMap.find(field)->second;
-                                    LibUtilities::BasisKeyVector def;
-
-                                    for(g = 0; g < fielddef[k]->m_elementIDs.size(); ++g)
-                                    {
-                                        ExpansionShPtr tmpexp =
-                                                MemoryManager<Expansion>::AllocateSharedPtr(geom, def);
-                                        (*expansionMap)[fielddef[k]->m_elementIDs[g]] = tmpexp;
-                                    }
-                                }
-                            }
-                        }
                     }
                 }
             }
@@ -2908,6 +2911,7 @@ namespace Nektar
                         if(!UniOrder)
                         {
                             cnt++;
+                            cnt += fielddef[i]->m_numHomogeneousDir;
                         }
                         bkeyvec.push_back(bkey);
                     }
@@ -2928,6 +2932,7 @@ namespace Nektar
                         if(!UniOrder)
                         {
                             cnt += 2;
+                            cnt += fielddef[i]->m_numHomogeneousDir;
                         }
                     }
                     break;
@@ -2948,6 +2953,7 @@ namespace Nektar
                         if(!UniOrder)
                         {
                             cnt += 2;
+                            cnt += fielddef[i]->m_numHomogeneousDir;
                         }
                     }
                     break;
@@ -2967,7 +2973,7 @@ namespace Nektar
 
                         if(!UniOrder)
                         {
-                            cnt += 2;
+                            cnt += 3;
                         }
                     }
                     break;
@@ -2987,7 +2993,7 @@ namespace Nektar
 
                         if(!UniOrder)
                         {
-                            cnt += 2;
+                            cnt += 3;
                         }
                     }
                     break;
@@ -3007,7 +3013,7 @@ namespace Nektar
 
                         if(!UniOrder)
                         {
-                            cnt += 2;
+                            cnt += 3;
                         }
                     }
                     break;
@@ -3027,7 +3033,7 @@ namespace Nektar
 
                         if(!UniOrder)
                         {
-                            cnt += 2;
+                            cnt += 3;
                         }
                     }
                     break;

solvers/IncNavierStokesSolver/CMakeLists.txt

@@ -123,6 +123,7 @@ IF( NEKTAR_SOLVER_INCNAVIERSTOKES )
         ADD_NEKTAR_TEST(ChanFlow_3DH1D_Parallel_mode1)
         ADD_NEKTAR_TEST(ChanFlow_3DH1D_Parallel_mode2)
         ADD_NEKTAR_TEST(ChanFlow_m3_par)
+        ADD_NEKTAR_TEST(KovaFlow_expFromFile_par)
         ADD_NEKTAR_TEST_LENGTHY(ChanFlow_m8_BodyForce_par)
         ADD_NEKTAR_TEST(KovaFlow_varP_per_par)
         ADD_NEKTAR_TEST_LENGTHY(Hex_channel_m8_par)

solvers/IncNavierStokesSolver/Tests/KovaFlow_expFromFile.xml

+<?xml version="1.0" encoding="utf-8" ?>
+
+<NEKTAR xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance"
+    xsi:noNamespaceSchemaLocation="http://www.nektar.info/schema/nektar.xsd">
+
+    <EXPANSIONS>
+        <F FILE="KovaFlow_m8.rst" />
+    </EXPANSIONS>
+
+    <CONDITIONS>
+        <SOLVERINFO>
+            <I PROPERTY="SolverType" VALUE="VelocityCorrectionScheme" />
+            <I PROPERTY="EQTYPE" VALUE="UnsteadyNavierStokes" />
+            <I PROPERTY="AdvectionForm" VALUE="Convective" />
+            <I PROPERTY="Projection" VALUE="Galerkin" />
+            <I PROPERTY="TimeIntegrationMethod" VALUE="IMEXOrder1" />
+            <I PROPERTY="GlobalSysSoln" VALUE="XxtMultiLevelStaticCond" />
+            <I PROPERTY="WEIGHTPARTITIONS" VALUE="DOF" />
+        </SOLVERINFO>
+
+        <PARAMETERS>
+            <P> TimeStep      = 0.001        </P>
+            <P> NumSteps      = 100       </P>
+            <P> IO_CheckSteps = 100       </P>
+            <P> IO_InfoSteps  = 100       </P>
+            <P> Kinvis        = 0.025        </P>
+        </PARAMETERS>
+
+        <VARIABLES>
+            <V ID="0"> u </V>
+            <V ID="1"> v </V>
+            <V ID="2"> p </V>
+        </VARIABLES>
+
+        <BOUNDARYREGIONS>
+            <B ID="0"> C[1] </B>
+            <B ID="1"> C[2] </B>
+            <B ID="2"> C[3] </B>
+        </BOUNDARYREGIONS>
+
+        <BOUNDARYCONDITIONS>
+            <REGION REF="0">
+                <D VAR="u" VALUE="1-1.619099729265964*cos(2*PI*y)" />
+                <D VAR="v" VALUE="-0.248344108585656*sin(2*PI*y)" />
+                <N VAR="p" USERDEFINEDTYPE="H" VALUE="0" />
+            </REGION>
+            <REGION REF="1">
+                <D VAR="u" VALUE="1-0.381463333531742*cos(2*PI*y)" />
+                <D VAR="v" VALUE="-0.058510399212408*sin(2*PI*y)" />
+                <D VAR="p" VALUE="0.427242862585425" />
+            </REGION>
+            <REGION REF="2">
+                <N VAR="u" VALUE="0" />
+                <D VAR="v" VALUE="0" />
+                <N VAR="p" VALUE="0" />
+            </REGION>
+        </BOUNDARYCONDITIONS>
+
+        <FUNCTION NAME="InitialConditions">
+            <F VAR="u,v,p" FILE="KovaFlow_m8.rst" />
+        </FUNCTION>
+
+        <FUNCTION NAME="ExactSolution">
+            <E VAR="u" VALUE="(1-exp(-0.963740544195769*x)*cos(2*PI*y))" />
+            <E VAR="v"
+                VALUE="(-0.963740544195769/(2*PI))*exp(-0.963740544195769*x)*sin(2*PI*y)" />
+            <E VAR="p" VALUE="0.5*(1-exp(-2*0.963740544195769*x))" />
+        </FUNCTION>
+
+    </CONDITIONS>
+
+    <GEOMETRY DIM="2" SPACE="2">
+        <VERTEX COMPRESSED="B64Z-LittleEndian" BITSIZE="64">eJxjYEAGD/aj0gwMjFjlf9jDRJhQ5D/Yo8szY5VHmM+C1fwPcP2s2N0Hl2dDk18/MZnnf1aFLUyEHUpHO9a8Cq6u2YtuPweqfgz3caK5H2r+XpgIF5o8uvu4scoj/MeD1X5E+PGiut8WXZ4Pq36E+fxY5RHuE0CTdze//3OXWwLcf4JQ+vG+paKKtTm26OYLQemmsp+NV/7326CbL4zi/py9UPPh8QMANMBNrQAA</VERTEX>
+        <EDGE COMPRESSED="B64Z-LittleEndian" BITSIZE="64">eJx1kskSgjAQBUFckDUIigv6/5/pgelLVyWXV/3yMjWTpCj2dQjtQkvxoBzchlZickflx9CTuAk9i6l3UT6F1uJrRhupz7di5unE1Ou1Dw/iSX1PyiX5twyPYnTOKHUX+dzjXczcD/XPP1jFvOtTTN8v+bz3W8wcH/n8g03MHF/5rJ+Y+f7dGwRp</EDGE>
+        <ELEMENT>
+            <Q COMPRESSED="B64Z-LittleEndian" BITSIZE="64">eJx1zskOgyAABmFsRaq0bl3t4vs/ppfh4CTl8iVDAn8I+1PhAY9/eo0RG92XnvCErd4tvcOMZ/1zwR4H7YrqI056p1Gf8aq9Sf2Gd+0t+x/4xJf2tuoLvrW3U//gV3uz+g9X7d0A9rcC/gAA</Q>
+        </ELEMENT>
+        <COMPOSITE>
+            <C ID="0"> Q[0-11] </C>
+            <C ID="1"> E[23,25,27,29] </C>
+            <C ID="2"> E[3,6,9,12] </C>
+            <C ID="3"> E[0,11,13,21-22,30] </C>
+        </COMPOSITE>
+        <DOMAIN> C[0] </DOMAIN>
+    </GEOMETRY>
+</NEKTAR>

solvers/IncNavierStokesSolver/Tests/KovaFlow_expFromFile_par.tst

+<?xml version="1.0" encoding="utf-8"?>
+<test>
+    <description>Kovasznay Flow, expansion from restart file, par(2)</description>
+    <executable>IncNavierStokesSolver</executable>
+    <parameters>--use-metis KovaFlow_expFromFile.xml</parameters>
+    <processes>2</processes>
+    <files>
+        <file description="Session File">KovaFlow_expFromFile.xml</file>
+	<file description="Session File">KovaFlow_m8.rst</file>
+    </files>
+    <metrics>
+        <metric type="L2" id="1">
+            <value variable="u" tolerance="1e-12">3.22767e-07</value>
+            <value variable="v" tolerance="1e-12">1.68898e-06</value>
+	    <value variable="p" tolerance="1e-12">9.87849e-06</value>
+        </metric>
+        <metric type="Linf" id="2">
+            <value variable="u" tolerance="1e-12">4.31288e-07</value>
+            <value variable="v" tolerance="1e-12">2.23918e-06</value>
+	    <value variable="p" tolerance="1e-12">3.43342e-05</value>
+        </metric>
+    </metrics>
+</test>