Commit 4aff7ccd authored by Dave Moxey's avatar Dave Moxey
Browse files

Merge branch 'feature/DriverAdaptive' into 'master'

Feature/driver adaptive

This MR introduces a driver which uses an adaptive polynomial order procedure for 2D and 3DH1D expansions. 

I only tested it for the incompressible solver, but I suppose it might work with other solvers as well. Extending it to 3D should be simple, but I did not attempt that since the incompressible solver does not support variable order in 3D yet. 

A bug in ExtractCoeffsToCoeffs is also fixed by this branch.

See merge request !538
parents 6f3d6c13 371ec699
......@@ -775,7 +775,63 @@ Once these three files (the \inltt{Geometry} in \inlsh{Session.xml.gz}, the non
IncNavierStokesSolver Session.xml.gz Session.xml
\end{lstlisting}
\section{Adaptive polynomial order Session file configuration}
\label{SectionDriverAdaptive}
An adaptive polynomial order procedure is available for 2D and Quasi-3D simulations.
This procedure consists of the following steps:
\begin{itemize}
\item Advance the equations for a determined number of time steps
\item Use the sensor defined in equation \ref{eq:sensor} to
estimate an error measure (the variable used in the sensor can be specified).
The error is defined here as the square of the sensor.
\item Use the error to determine if the order in each element should be
increased by one, decreased by one, or left unaltered.
\item Project the solution in each element to the new polynomial order and use it as an
initial condition to restart the equation, repeating all steps a given number of times.
\end{itemize}
It is important to note that restarting the simulation after the refinement can be an
expensive operation (in a typical case 200 times the cost of a single time step). Therefore, the number
of steps between successive refinements needs to be carefully chosen, since if this value is too
low the procedure becomes inefficient, while if it is too high the refinement might not capture
accurately structures that are convected by the flow.
\subsection{Solver Info}
The use of the adaptive polynomial order procedure
is enforced through the definition of the \inltt{Driver} which has to be
\inltt{Adaptive}.
\subsection{Parameters}
The following parameters can be specified in the \texttt{PARAMETERS} section of the session file:
\begin{itemize}
\item \inltt{NumSteps}: when using the adaptive order procedure, this parameter determines
the number of time steps between successive refinements.
\item \inltt{NumRuns}: this parameter defines the number of times the sequence of solving
the equation and refining is performed. Therefore,
the total number of time steps in the simulation is $NumSteps \times NumRuns$.
\item \inltt{AdaptiveMaxModes}: sets the maximum number of modes (in each direction)
that can be used in an element during the adaptive procedure. The solution will not be
refined beyond this point, even if the error is higher than the tolerance. Default value: 12.
\item \inltt{AdaptiveMinModes}: sets the minimal number of modes (in each direction)
that can be used in an element during the adaptive procedure. Default value: 4.
\item \inltt{AdaptiveUpperTolerance}: defines a constant tolerance. The polynomial order
in an element is increased whenever the error is higher than this value. This can be replaced
by a spatially-varying function, as described below. Default value: $10^{-6}$.
\item \inltt{AdaptiveLowerinModes}: defines a constant tolerance. The polynomial order in an
element is decreased whenever the error is lower than this value. This can also be replaced
by a spatially-varying function. Default value: $10^{-8}$.
\item \inltt{AdaptiveSensorVariable}: integer defining which variable will be considered when calculating
the error. For example, if this parameter is set to $1$ in the Incompressible Navier-Stokes Solver,
the error will be estimated using the $v$ velocity. Default value: $0$.
\end{itemize}
\subsection{Functions}
Spatially varying tolerances can be specified by defining the functions \inltt{AdaptiveLowerinModes} and/or \inltt{AdaptiveUpperTolerance}. In this case, the tolerance in an 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.
\section{Examples}
......
......@@ -1315,6 +1315,8 @@ namespace Nektar
std::map<int, MeshVertex>::iterator vVertIt;
std::map<std::string, std::string>::iterator vAttrIt;
std::vector<unsigned int> idxList;
// Populate lists of elements, edges and vertices required.
for (boost::tie(vertit, vertit_end) = boost::vertices(pGraph);
vertit != vertit_end;
......@@ -1925,10 +1927,8 @@ namespace Nektar
// which belong to this partition.
for (vIt = m_meshComposites.begin(); vIt != m_meshComposites.end(); ++vIt)
{
bool comma = false; // Set to true after first entity output
bool range = false; // True when entity IDs form a range
int last_idx = -2; // Last entity ID output
std::string vCompositeStr = "";
idxList.clear();
for (unsigned int j = 0; j < vIt->second.list.size(); ++j)
{
// Based on entity type, check if in this partition
......@@ -1960,34 +1960,11 @@ namespace Nektar
break;
}
// Condense consecutive entity IDs into ranges
// last_idx initially -2 to avoid error for ID=0
if (last_idx + 1 == vIt->second.list[j])
{
last_idx++;
range = true;
continue;
}
// This entity is not in range, so close previous range with
// last_idx
if (range)
{
vCompositeStr += "-" + boost::lexical_cast<std::string>(last_idx);
range = false;
}
// Output ID, which is either standalone, or will start a
// range.
vCompositeStr += comma ? "," : "";
vCompositeStr += boost::lexical_cast<std::string>(vIt->second.list[j]);
last_idx = vIt->second.list[j];
comma = true;
}
// If last entity is part of a range, it must be output now
if (range)
{
vCompositeStr += "-" + boost::lexical_cast<std::string>(last_idx);
idxList.push_back(vIt->second.list[j]);
}
std::string vCompositeStr = ParseUtils::GenerateSeqString(idxList);
if (vCompositeStr.length() > 0)
{
vComposites[vIt->first] = vIt->second;
......@@ -2001,18 +1978,16 @@ namespace Nektar
}
}
idxList.clear();
std::string vDomainListStr;
bool comma = false;
for (unsigned int i = 0; i < m_domain.size(); ++i)
{
if (vComposites.find(m_domain[i]) != vComposites.end())
{
vDomainListStr += comma ? "," : "";
comma = true;
vDomainListStr += boost::lexical_cast<std::string>(m_domain[i]);
idxList.push_back(m_domain[i]);
}
}
vDomainListStr = "C[" + vDomainListStr + "]";
vDomainListStr = "C[" + ParseUtils::GenerateSeqString(idxList) + "]";
TiXmlText* vDomainList = new TiXmlText(vDomainListStr);
vDomain->LinkEndChild(vDomainList);
......@@ -2055,19 +2030,20 @@ namespace Nektar
vSeqStr = vSeqStr.substr(indxBeg, indxEnd - indxBeg + 1);
std::vector<unsigned int> vSeq;
ParseUtils::GenerateSeqVector(vSeqStr.c_str(), vSeq);
std::string vListStr;
bool comma = false;
idxList.clear();
for (unsigned int i = 0; i < vSeq.size(); ++i)
{
if (vComposites.find(vSeq[i]) != vComposites.end())
{
vListStr += comma ? "," : "";
comma = true;
vListStr += boost::lexical_cast<std::string>(vSeq[i]);
idxList.push_back(vSeq[i]);
}
}
int p = atoi(vItem->Attribute("ID"));
std::string vListStr = ParseUtils::GenerateSeqString(idxList);
if (vListStr.length() == 0)
{
TiXmlElement* tmp = vItem;
......
......@@ -39,6 +39,7 @@
#include <LibUtilities/BasicConst/NektarUnivTypeDefs.hpp>
#include <vector>
#include <sstream>
#include <boost/version.hpp>
#include <LibUtilities/LibUtilitiesDeclspec.h>
......@@ -155,6 +156,52 @@ namespace Nektar
space_p).full;
}
static std::string GenerateSeqString(const std::vector<unsigned int> &elmtids)
{
std::stringstream idStringStream;
bool setdash = true;
unsigned int endval;
if (elmtids.size() == 0)
{
return std::string("");
}
idStringStream << elmtids[0];
for (int i = 1; i < elmtids.size(); ++i)
{
if (elmtids[i] == elmtids[i - 1] + 1)
{
if (setdash)
{
idStringStream << "-";
setdash = false;
}
if (i == elmtids.size() - 1) // last element
{
idStringStream << elmtids[i];
}
else
{
endval = elmtids[i];
}
}
else
{
if (setdash == false) // finish off previous dash sequence
{
idStringStream << endval;
setdash = true;
}
idStringStream << "," << elmtids[i];
}
}
return idStringStream.str();
}
private:
struct SymbolFunctor
......
......@@ -161,7 +161,10 @@ namespace Gs
MPI_Comm_dup(vCommMpi->GetComm(), &vComm.c);
vComm.id = vCommMpi->GetRank();
vComm.np = vCommMpi->GetSize();
return nektar_gs_setup(pId.get(),pId.num_elements(), &vComm, 0, gs_auto, 1);
gs_data* result = nektar_gs_setup(pId.get(),pId.num_elements(),
&vComm, 0, gs_auto, 1);
MPI_Comm_free(&vComm.c);
return result;
#else
return 0;
#endif
......@@ -202,7 +205,9 @@ namespace Gs
static inline void Finalise (gs_data *pGsh)
{
#ifdef NEKTAR_USE_MPI
if (pGsh)
int finalized;
MPI_Finalized(&finalized);
if (pGsh && !finalized)
{
nektar_gs_free(pGsh);
}
......
......@@ -169,7 +169,10 @@ namespace Xxt
MPI_Comm_dup(vCommMpi->GetComm(), &vComm.c);
vComm.id = vCommMpi->GetRank();
vComm.np = vCommMpi->GetSize();
return nektar_crs_setup(pRank, &pId[0], nz, &pAi[0], &pAj[0], &pAr[0], 0, &vComm);
crs_data* result = nektar_crs_setup(pRank, &pId[0], nz, &pAi[0],
&pAj[0], &pAr[0], 0, &vComm);
MPI_Comm_free(&vComm.c);
return result;
#else
return 0;
#endif
......@@ -198,9 +201,11 @@ namespace Xxt
static inline void Finalise (crs_data* pCrs)
{
#ifdef NEKTAR_USE_MPI
if (pCrs)
int finalized;
MPI_Finalized(&finalized);
if (pCrs && !finalized)
{
//nektar_crs_free(pCrs);
nektar_crs_free(pCrs);
}
#endif
}
......
......@@ -1113,12 +1113,14 @@ namespace Nektar
Array<OneD, NekDouble> tmp5(unique_edges, 1.0);
Array<OneD, NekDouble> tmp6(unique_faces, 1.0);
Gs::Gather(tmp4, Gs::gs_add, tmp1);
Gs::Finalise(tmp1);
if (unique_edges > 0)
{
Array<OneD, long> edgeArray(unique_edges, &procEdges[0]);
Gs::gs_data *tmp2 = Gs::Init(edgeArray, vComm);
Gs::Gather(tmp5, Gs::gs_add, tmp2);
Gs::Finalise(tmp2);
}
if (unique_faces > 0)
......@@ -1126,6 +1128,7 @@ namespace Nektar
Array<OneD, long> faceArray(unique_faces, &procFaces[0]);
Gs::gs_data *tmp3 = Gs::Init(faceArray, vComm);
Gs::Gather(tmp6, Gs::gs_add, tmp3);
Gs::Finalise(tmp3);
}
// Finally, fill the partVerts set with all non-Dirichlet
......@@ -1363,6 +1366,7 @@ namespace Nektar
}
Gs::gs_data *tmp = Gs::Init(edgeId, vComm);
Gs::Gather(edgeDof, Gs::gs_min, tmp);
Gs::Finalise(tmp);
for (i=0; i < dofs[1].size(); i++)
{
dofs[1][edgeId[i]] = (int) (edgeDof[i]+0.5);
......@@ -1377,6 +1381,7 @@ namespace Nektar
}
Gs::gs_data *tmp2 = Gs::Init(faceId, vComm);
Gs::Gather(faceDof, Gs::gs_min, tmp2);
Gs::Finalise(tmp2);
for (i=0; i < dofs[2].size(); i++)
{
dofs[2][faceId[i]] = (int) (faceDof[i]+0.5);
......@@ -1842,6 +1847,8 @@ namespace Nektar
*/
AssemblyMapCG::~AssemblyMapCG()
{
Gs::Finalise(m_gsh);
Gs::Finalise(m_bndGsh);
}
/**
......
......@@ -2229,27 +2229,19 @@ namespace Nektar
int i;
int offset = 0;
// check if the same and if so just copy over coeffs
if(fromExpList->GetNcoeffs() == m_ncoeffs)
{
Vmath::Vcopy(m_ncoeffs,fromCoeffs,1,toCoeffs,1);
}
else
for(i = 0; i < (*m_exp).size(); ++i)
{
std::vector<unsigned int> nummodes;
for(i = 0; i < (*m_exp).size(); ++i)
int eid = m_offset_elmt_id[i];
for(int j= 0; j < fromExpList->GetExp(eid)->GetNumBases(); ++j)
{
int eid = m_offset_elmt_id[i];
for(int j= 0; j < fromExpList->GetExp(eid)->GetNumBases(); ++j)
{
nummodes.push_back(fromExpList->GetExp(eid)->GetBasisNumModes(j));
}
nummodes.push_back(fromExpList->GetExp(eid)->GetBasisNumModes(j));
}
(*m_exp)[eid]->ExtractDataToCoeffs(&fromCoeffs[offset], nummodes,0,
&toCoeffs[m_coeff_offset[eid]]);
(*m_exp)[eid]->ExtractDataToCoeffs(&fromCoeffs[offset], nummodes,0,
&toCoeffs[m_coeff_offset[eid]]);
offset += fromExpList->GetExp(eid)->GetNcoeffs();
}
offset += fromExpList->GetExp(eid)->GetNcoeffs();
}
}
......
......@@ -829,11 +829,12 @@ namespace Nektar
{
int i;
int fromNcoeffs_per_plane = fromExpList->GetPlane(0)->GetNcoeffs();
int toNcoeffs_per_plane = m_planes[0]->GetNcoeffs();
Array<OneD, NekDouble> tocoeffs_tmp, fromcoeffs_tmp;
for(i = 0; i < m_planes.num_elements(); ++i)
{
m_planes[i]->ExtractCoeffsToCoeffs(fromExpList->GetPlane(i),fromcoeffs_tmp = fromCoeffs + fromNcoeffs_per_plane*i, tocoeffs_tmp = toCoeffs + m_ncoeffs*i);
m_planes[i]->ExtractCoeffsToCoeffs(fromExpList->GetPlane(i),fromcoeffs_tmp = fromCoeffs + fromNcoeffs_per_plane*i, tocoeffs_tmp = toCoeffs + toNcoeffs_per_plane*i);
}
}
......
......@@ -54,6 +54,7 @@ namespace Nektar
&pLocToGloMap)
: GlobalLinSys(pKey, pExp, pLocToGloMap)
{
m_crsData = 0;
}
GlobalLinSysXxt::~GlobalLinSysXxt()
......
......@@ -14,6 +14,7 @@ SET(SOLVER_UTILS_SOURCES
Diffusion/DiffusionLFR.cpp
Diffusion/DiffusionLFRNS.cpp
Driver.cpp
DriverAdaptive.cpp
DriverArnoldi.cpp
DriverModifiedArnoldi.cpp
DriverStandard.cpp
......@@ -57,6 +58,7 @@ SET(SOLVER_UTILS_HEADERS
Diffusion/DiffusionLFR.h
Diffusion/DiffusionLFRNS.h
Driver.h
DriverAdaptive.h
DriverArnoldi.h
DriverModifiedArnoldi.h
DriverStandard.h
......
This diff is collapsed.
///////////////////////////////////////////////////////////////////////////////
//
// File DriverAdaptive.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: Driver class for adaptive solver
//
///////////////////////////////////////////////////////////////////////////////
#ifndef NEKTAR_SOLVERUTILS_DRIVERADAPTIVE_H
#define NEKTAR_SOLVERUTILS_DRIVERADAPTIVE_H
#include <SolverUtils/Driver.h>
#include <MultiRegions/GlobalLinSys.h>
namespace Nektar
{
namespace SolverUtils
{
/// Base class for the adaptive polynomial order driver.
class DriverAdaptive : public Driver
{
public:
friend class MemoryManager<DriverAdaptive>;
/// Creates an instance of this class
static DriverSharedPtr create(
const LibUtilities::SessionReaderSharedPtr &pSession)
{
DriverSharedPtr p =
MemoryManager<DriverAdaptive>::AllocateSharedPtr(pSession);
p->InitObject();
return p;
}
/// Name of the class
static std::string className;
protected:
/// Constructor
SOLVER_UTILS_EXPORT DriverAdaptive(
const LibUtilities::SessionReaderSharedPtr pSession);
/// Destructor
SOLVER_UTILS_EXPORT virtual ~DriverAdaptive();
/// Second-stage initialisation
SOLVER_UTILS_EXPORT virtual void v_InitObject(ostream &out = cout);
/// Virtual function for solve implementation.
SOLVER_UTILS_EXPORT virtual void v_Execute(ostream &out = cout);
SOLVER_UTILS_EXPORT void ReplaceExpansion(
Array<OneD, MultiRegions::ExpListSharedPtr> &fields,
map<int, int> deltaP);
static std::string driverLookupId;
};
}
} // end of namespace
#endif // NEKTAR_SOLVERUTILS_DRIVERADAPTIVE_H
......@@ -153,7 +153,7 @@ void DriverSteadyState::v_Execute(ostream &out)
// m_steps is set to 1. Then "m_equ[m_nequ - 1]->DoSolve()" will run
// for only one time step
m_equ[m_nequ - 1]->SetStepsToOne();
m_equ[m_nequ - 1]->SetSteps(1);
ofstream m_file("ConvergenceHistory.txt", ios::out | ios::trunc);
Array<OneD, Array<OneD, NekDouble> > q0(NumVar_SFD);
......
......@@ -660,6 +660,8 @@ namespace Nektar
m_session->LoadParameter("NumQuadPointsError",
m_NumQuadPointsError, 0);
m_nchk = 1;
// Zero all physical fields initially
ZeroPhysFields();
}
......
......@@ -359,7 +359,7 @@ namespace Nektar
SOLVER_UTILS_EXPORT inline void CopyToPhysField(const int i,
Array<OneD, NekDouble> &output);
SOLVER_UTILS_EXPORT inline void SetStepsToOne();
SOLVER_UTILS_EXPORT inline void SetSteps(const int steps);
SOLVER_UTILS_EXPORT void ZeroPhysFields();
......@@ -406,6 +406,41 @@ namespace Nektar
/// Perform a case-insensitive string comparison.
SOLVER_UTILS_EXPORT int NoCaseStringCompare(
const string & s1, const string& s2) ;
SOLVER_UTILS_EXPORT int GetCheckpointNumber()
{
return m_nchk;
}
SOLVER_UTILS_EXPORT void SetCheckpointNumber(int num)
{
m_nchk = num;
}
SOLVER_UTILS_EXPORT int GetCheckpointSteps()
{
return m_checksteps;
}
SOLVER_UTILS_EXPORT void SetCheckpointSteps(int num)
{
m_checksteps = num;
}
SOLVER_UTILS_EXPORT void SetTime(
const NekDouble time)
{
m_time = time;
}
SOLVER_UTILS_EXPORT void SetInitialStep(
const int step)
{
m_initialStep = step;
}
/// Evaluates the boundary conditions at the given time.
SOLVER_UTILS_EXPORT void SetBoundaryConditions(NekDouble time);
/// Virtual function to identify if operator is negated in DoSolve
SOLVER_UTILS_EXPORT virtual bool v_NegatedOp();
......@@ -435,6 +470,8 @@ namespace Nektar
std::string m_sessionName;
/// Current time of simulation.
NekDouble m_time;
/// Number of the step where the simulation should begin
int m_initialStep;
/// Finish time of the simulation.
NekDouble m_fintime;
/// Time step size
......@@ -445,6 +482,8 @@ namespace Nektar
std::set<std::string> m_loadedFields;
/// Time between checkpoints.
NekDouble m_checktime;
/// Number of checkpoints written so far
int m_nchk;
/// Number of steps to take.
int m_steps;
/// Number of steps between checkpoints.
......@@ -524,9 +563,6 @@ namespace Nektar
SOLVER_UTILS_EXPORT virtual void v_InitObject();
/// Evaluates the boundary conditions at the given time.
SOLVER_UTILS_EXPORT void SetBoundaryConditions(NekDouble time);
/// Virtual function for initialisation implementation.
SOLVER_UTILS_EXPORT virtual void v_DoInitialise();
......@@ -859,9 +895,9 @@ namespace Nektar
return m_timestep;
}
inline void EquationSystem::SetStepsToOne(void)
inline void EquationSystem::SetSteps(const int steps)