Commit 7c085bc8 authored by Dave Moxey's avatar Dave Moxey
Browse files

Merge branch 'master' into fix/FC3DH1Defficiency

parents 1711e3d9 fda40e6f
Changelog
=========
v4.3.2
------
**Library**:
- Add small optimisation for DriverAdaptive (!618)
- Updated FFTW build to use the compiler used for building Nektar++ (!629)
- Fix numbering bug in periodic boundary conditions (!631)
- Print error message for invalid equation also in release version (!634)
- HistoryPoints filter now uses closest plane to requested z-coordinate and
output is produced in physical space (!621).
- Fix minor performance issue with time integration schemes (!632)
- Fix FilterCheckpoint filter to be consistent with `IO_CheckSteps` (!633)
**FieldConvert**:
- Fix appearence of duplicate messages when running in parallel (!626)
**Packaging**:
- Fixes for DEB package dependencies (!630)
v4.3.1
------
**Library**:
......
......@@ -38,7 +38,7 @@ IF (NEKTAR_USE_FFTW)
BINARY_DIR ${TPBUILD}/fftw-3.2.2
TMP_DIR ${TPBUILD}/fftw-3.2.2-tmp
INSTALL_DIR ${TPDIST}
CONFIGURE_COMMAND ${TPSRC}/fftw-3.2.2/configure --prefix=${TPDIST} --quiet --enable-shared --disable-dependency-tracking
CONFIGURE_COMMAND CC=${CMAKE_C_COMPILER} ${TPSRC}/fftw-3.2.2/configure --prefix=${TPDIST} --quiet --enable-shared --disable-dependency-tracking
)
SET(FFTW_LIBRARY fftw3 CACHE FILEPATH
......
......@@ -592,11 +592,15 @@ the advection term using the pressure inverse mass matrix. It can be used just i
\item \inltt{SpectralVanishingViscosity}: activates a stabilization technique
which increases the viscosity on the highest Fourier frequencies of a Quasi-3D approach.
which increases the viscosity on the modes with the highest frequencies.
\begin{lstlisting}[style=XMLStyle]
<I PROPERTY="SpectralVanishingViscosity" VALUE="True"/>
\end{lstlisting}
In a Quasi-3D simulation, this will affect both the Fourier and the spectral/hp expansions.
To activate them independently, use \inltt{SpectralVanishingViscositySpectralHP}
and \inltt{SpectralVanishingViscosityHomo1D}.
\item \inltt{DEALIASING}: activates the 3/2 padding rule on the advection term
of a Quasi-3D simulation.
\begin{lstlisting}[style=XMLStyle]
......
......@@ -80,7 +80,7 @@ namespace Nektar
{
m_expr_id = -1;
std::string msg(std::string("Equation::Equation() fails on expression [") + m_expr + std::string("]\n"));
ASSERTL1(false, msg);
ASSERTL0(false, msg);
throw e;
return;
}
......@@ -88,7 +88,7 @@ namespace Nektar
{
m_expr_id = -1;
std::string msg(std::string("Equation::Equation() fails on expression [") + m_expr + std::string("]\n"));
ASSERTL1(false, msg);
ASSERTL0(false, msg);
throw e;
return;
}
......
......@@ -180,7 +180,10 @@ namespace Nektar
// and determine the full pathname to the file to write out.
// Any existing file/directory which is in the way is removed.
std::string filename = SetUpOutput(outFile);
SetUpFieldMetaData(outFile, fielddefs, fieldmetadatamap);
if (m_comm->GetSize() > 1)
{
SetUpFieldMetaData(outFile, fielddefs, fieldmetadatamap);
}
// Create the file (partition)
TiXmlDocument doc;
......
......@@ -258,6 +258,19 @@ namespace Nektar
}
}
static bool PoolCreated(std::string whichPool)
{
bool value = false;
typename ValueContainerPool::iterator x;
x = m_ValueContainerPool.find(whichPool);
if (x != m_ValueContainerPool.end())
{
value = true;
}
return value;
}
static void EnableManagement(std::string whichPool = "")
{
typename FlagContainerPool::iterator x;
......
......@@ -1614,8 +1614,8 @@ namespace Nektar
{
t_new[0] += V(0,j)*t_old[j];
}
i_start = 1;
}
i_start = 1;
}
for(i = i_start; i < m_numsteps; i++)
......
......@@ -498,7 +498,8 @@ namespace Nektar
// we find it, set as Dirichlet with the vertex id gId.
if (pIt->first == meshVertId)
{
graph[0][meshVertId] = gId < 0 ? graphVertId++ : gId;
gId = gId < 0 ? graphVertId++ : gId;
graph[0][meshVertId] = gId;
for (i = 0; i < pIt->second.size(); ++i)
{
......@@ -522,7 +523,8 @@ namespace Nektar
if (found)
{
graph[0][pIt->first] = gId < 0 ? graphVertId++ : gId;
gId = gId < 0 ? graphVertId++ : gId;
graph[0][pIt->first] = gId;
for (i = 0; i < pIt->second.size(); ++i)
{
......
......@@ -761,20 +761,12 @@ namespace Nektar
int planes_offset = 0;
Array<OneD, NekDouble> coeff_tmp;
std::map<int,int>::iterator it;
int IDoffset = 0;
// introduce a 2 plane offset for single mode case so can
// be post-processed or used in MultiMode expansion.
if(m_homogeneousBasis->GetBasisType() == LibUtilities::eFourierSingleMode)
{
IDoffset = 2;
}
// Build map of plane IDs lying on this processor.
std::map<int,int> homoZids;
for (i = 0; i < m_planes.num_elements(); ++i)
{
homoZids[m_transposition->GetPlaneID(i)+IDoffset] = i;
homoZids[m_transposition->GetPlaneID(i)] = i;
}
if(fielddef->m_numHomogeneousDir)
......
......@@ -293,9 +293,14 @@ void DriverAdaptive::v_Execute(ostream &out)
//
// @todo This could be made better by replacing individual matrices
// within the linear system.
LibUtilities::NekManager<MultiRegions::GlobalLinSysKey,
MultiRegions::GlobalLinSys>::
ClearManager(std::string("GlobalLinSys"));
if (LibUtilities::NekManager<MultiRegions::GlobalLinSysKey,
MultiRegions::GlobalLinSys>::
PoolCreated(std::string("GlobalLinSys")))
{
LibUtilities::NekManager<MultiRegions::GlobalLinSysKey,
MultiRegions::GlobalLinSys>::
ClearManager(std::string("GlobalLinSys"));
}
int chkNumber = m_equ[0]->GetCheckpointNumber();
int chkSteps = m_equ[0]->GetCheckpointSteps();
......
......@@ -2267,8 +2267,6 @@ namespace Nektar
else if (m_multipleModes)
{
AddSummaryItem(s, "ModeType", "Multiple Modes");
AddSummaryItem(s, "Selected Mode",
boost::lexical_cast<string>(m_NumMode));
}
}
else if(m_HomogeneousType == eHomogeneous2D)
......
......@@ -549,8 +549,6 @@ namespace Nektar
int m_HomoDirec; ///< number of homogenous directions
int m_NumMode; ///< Mode to use in case of single mode analysis
/// Initialises EquationSystem class members.
SOLVER_UTILS_EXPORT EquationSystem( const LibUtilities::SessionReaderSharedPtr& pSession);
......
......@@ -63,12 +63,11 @@ FilterCheckpoint::FilterCheckpoint(
}
// OutputFrequency
it = pParams.find("OutputFrequency");
ASSERTL0(it != pParams.end(), "Missing parameter 'OutputFrequency'.");
LibUtilities::Equation equ(m_session, it->second);
m_outputFrequency = floor(equ.Evaluate());
m_outputIndex = 0;
m_index = 0;
m_fld = MemoryManager<LibUtilities::FieldIO>
::AllocateSharedPtr(pSession->GetComm());
......@@ -85,14 +84,14 @@ void FilterCheckpoint::v_Initialise(
{
m_index = 0;
m_outputIndex = 0;
v_Update(pFields, time);
}
void FilterCheckpoint::v_Update(
const Array<OneD, const MultiRegions::ExpListSharedPtr> &pFields,
const NekDouble &time)
{
m_index++;
if (m_index % m_outputFrequency > 0)
if (m_index++ % m_outputFrequency > 0)
{
return;
}
......
......@@ -36,6 +36,7 @@
#include <LibUtilities/Memory/NekMemoryManager.hpp>
#include <iomanip>
#include <SolverUtils/Filters/FilterHistoryPoints.h>
#include <MultiRegions/ExpList3DHomogeneous1D.h>
#include <boost/format.hpp>
......@@ -93,13 +94,25 @@ FilterHistoryPoints::FilterHistoryPoints(
it = pParams.find("OutputPlane");
if (it == pParams.end())
{
m_outputPlane = 0;
m_outputPlane = -1;
}
else
{
LibUtilities::Equation equ(m_session, it->second);
m_outputPlane = floor(equ.Evaluate());
}
it = pParams.find("WaveSpace");
if (it == pParams.end())
{
m_waveSpace = false;
}
else
{
std::string sOption = it->second.c_str();
m_waveSpace = ( boost::iequals(sOption,"true")) ||
( boost::iequals(sOption,"yes"));
}
}
// Points
......@@ -132,29 +145,53 @@ void FilterHistoryPoints::v_Initialise(
m_index = 0;
m_historyList.clear();
vector<unsigned int> planeIDs;
// Read history points
Array<OneD, NekDouble> gloCoord(3,0.0);
int dim = pFields[0]->GetGraph()->GetSpaceDimension();
if (m_isHomogeneous1D)
{
dim++;
}
int i = 0;
while (!m_historyPointStream.fail())
{
m_historyPointStream >> gloCoord[0]
>> gloCoord[1]
>> gloCoord[2];
if(m_isHomogeneous1D) // overwrite with plane z
if (!m_historyPointStream.fail())
{
NekDouble Z = (pFields[0]->GetHomogeneousBasis()
->GetZ())[m_outputPlane];
if(fabs(gloCoord[2]-Z) > NekConstants::kVertexTheSameDouble)
// Overwrite gloCoord[2] for 3DH1D using m_outputPlane if it is
// defined, or a nearby plane otherwise
if(m_isHomogeneous1D)
{
cout << "Reseting History point from " << gloCoord[2]
<< " to " << Z << endl;
int nplanes = pFields[0]->GetHomogeneousBasis()
->GetZ().num_elements();
NekDouble lhom = pFields[0]->GetHomoLen();
int plane;
if (m_outputPlane == -1)
{
// Pick plane immediately before the point
plane = floor((gloCoord[2]*nplanes)/lhom);
}
else
{
plane = m_outputPlane;
}
NekDouble Z = (pFields[0]->GetHomogeneousBasis()
->GetZ())[plane];
Z = (Z+1)*lhom/2;
if(fabs(gloCoord[2]-Z) > NekConstants::kVertexTheSameDouble)
{
cout << "Reseting History point from z = " << gloCoord[2]
<< " to z = " << Z << endl;
}
gloCoord[2] = Z;
planeIDs.push_back(plane);
}
gloCoord[2] = Z;
}
if (!m_historyPointStream.fail())
{
SpatialDomains::PointGeomSharedPtr vert
= MemoryManager<SpatialDomains::PointGeom>
::AllocateSharedPtr(dim, i, gloCoord[0],
......@@ -169,7 +206,7 @@ void FilterHistoryPoints::v_Initialise(
// Determine the unique process responsible for each history point
// For points on a partition boundary, must select a single process
LibUtilities::CommSharedPtr vComm = pFields[0]->GetComm();
int vRank = vComm->GetRank();
int vRank = vComm->GetRowComm()->GetRank();
int vHP = m_historyPoints.size();
Array<OneD, int> procList(vHP, -1 );
Array<OneD, int> idList (vHP, -1 );
......@@ -188,8 +225,16 @@ void FilterHistoryPoints::v_Initialise(
gloCoord[2]);
// Determine the expansion and local coordinates
idList[i] = pFields[0]->GetExpIndex(gloCoord,locCoords,
NekConstants::kGeomFactorsTol);
if (m_isHomogeneous1D)
{
idList[i] = pFields[0]->GetPlane(0)->GetExpIndex(gloCoord,locCoords,
NekConstants::kNekZeroTol);
}
else
{
idList[i] = pFields[0]->GetExpIndex(gloCoord,locCoords,
NekConstants::kNekZeroTol);
}
// Save Local coordinates for later
LocCoords.push_back(locCoords);
......@@ -222,6 +267,7 @@ void FilterHistoryPoints::v_Initialise(
// If multiple processes find they are the nearest (e.g. point lies
// on a partition boundary, we will choose the process of highest
// rank.
m_planeIDs = Array<OneD, int> (planeIDs.size(),-1);
for (i = 0; i < vHP; ++i)
{
if (dist_loc[i] == dist[i])
......@@ -239,6 +285,7 @@ void FilterHistoryPoints::v_Initialise(
// process ID
if (idList[i] != -1)
{
procList[i] = vRank;
if(m_isHomogeneous1D)
{
int j;
......@@ -246,7 +293,7 @@ void FilterHistoryPoints::v_Initialise(
= pFields[0]->GetZIDs();
for(j = 0; j < IDs.num_elements(); ++j)
{
if(IDs[j] == m_outputPlane)
if(IDs[j] == planeIDs[i])
{
break;
}
......@@ -254,14 +301,9 @@ void FilterHistoryPoints::v_Initialise(
if(j != IDs.num_elements())
{
m_outputPlane = j;
procList[i] = vRank;
m_planeIDs[i] = j;
}
}
else
{
procList[i] = vRank;
}
}
}
......@@ -336,7 +378,7 @@ void FilterHistoryPoints::v_Initialise(
if(m_isHomogeneous1D)
{
m_outputStream << ") at points:";
m_outputStream << ") at points:" << endl;
}
else
{
......@@ -358,7 +400,10 @@ void FilterHistoryPoints::v_Initialise(
if(m_isHomogeneous1D)
{
m_outputStream << "(in Wavespace)" << endl;
if (m_waveSpace)
{
m_outputStream << "# (in Wavespace)" << endl;
}
}
}
v_Update(pFields, time);
......@@ -382,7 +427,6 @@ void FilterHistoryPoints::v_Update(const Array<OneD, const MultiRegions::ExpList
int numFields = pFields.num_elements();
LibUtilities::CommSharedPtr vComm = pFields[0]->GetComm();
Array<OneD, NekDouble> data(numPoints*numFields, 0.0);
Array<OneD, NekDouble> gloCoord(3, 0.0);
std::list<std::pair<SpatialDomains::PointGeomSharedPtr, Array<OneD, NekDouble> > >::iterator x;
Array<OneD, NekDouble> physvals;
......@@ -399,18 +443,94 @@ void FilterHistoryPoints::v_Update(const Array<OneD, const MultiRegions::ExpList
{
locCoord = (*x).second;
expId = (*x).first->GetVid();
NekDouble value;
int plane = m_planeIDs[m_historyLocalPointMap[k]];
physvals = pFields[j]->GetPlane(m_outputPlane)->UpdatePhys() + pFields[j]->GetPhys_Offset(expId);
// transform elemental data if required.
if(pFields[j]->GetPhysState() == false)
if (m_waveSpace)
{
ASSERTL0( pFields[j]->GetWaveSpace() == true,
"HistoryPoints in wavespace require that solution is in wavespace");
}
if ( pFields[j]->GetWaveSpace() == false || m_waveSpace)
{
pFields[j]->GetPlane(m_outputPlane)->GetExp(expId)->BwdTrans(pFields[j]->GetPlane(m_outputPlane)->GetCoeffs() + pFields[j]->GetCoeff_Offset(expId),physvals);
if (plane != -1)
{
physvals = pFields[j]->GetPlane(plane)->
UpdatePhys() + pFields[j]->GetPhys_Offset(expId);
// transform elemental data if required.
if(pFields[j]->GetPhysState() == false)
{
pFields[j]->GetPlane(plane)->GetExp(expId)->
BwdTrans(pFields[j]->GetPlane(plane)->
GetCoeffs() + pFields[j]->
GetCoeff_Offset(expId),physvals);
}
// Interpolate data
value = pFields[j]->GetPlane(plane)->GetExp(expId)->
StdPhysEvaluate(locCoord,physvals);
}
}
else
{
// Create vector with eIDs across all planes
std::vector<unsigned int> eIDs;
int nPlanes = pFields[j]->GetZIDs().num_elements();
int elmtsPerPlane = pFields[j]->GetExpSize()/nPlanes;
// interpolate point can do with zero plane methods
data[m_historyLocalPointMap[k]*numFields+j] = pFields[j]->GetExp(expId)->StdPhysEvaluate(locCoord,physvals);
for ( int n = 0; n < nPlanes; n++)
{
eIDs.push_back(expId + n*elmtsPerPlane);
}
// Create new 3DH1D expansion with one element per plane
MultiRegions::ExpList3DHomogeneous1DSharedPtr tmp =
boost::dynamic_pointer_cast<MultiRegions::
ExpList3DHomogeneous1D>(pFields[j]);
ASSERTL0(tmp,"Failed to type cast expansion");
MultiRegions::ExpList3DHomogeneous1DSharedPtr exp =
MemoryManager<MultiRegions::
ExpList3DHomogeneous1D>::
AllocateSharedPtr(*tmp, eIDs);
// Fill phys array of new expansion and apply HomoBwdTrans
for ( int n = 0; n < nPlanes; n++)
{
Array<OneD, NekDouble> toPhys =
exp->GetPlane(n)->UpdatePhys();
if(pFields[j]->GetPhysState())
{
int nq = exp->GetPlane(0)->GetTotPoints();
Array<OneD, NekDouble> fromPhys =
pFields[j]->GetPlane(n)->GetPhys() +
pFields[j]->GetPhys_Offset(expId);
Vmath::Vcopy(nq, fromPhys, 1, toPhys, 1);
}
else
{
Array<OneD, NekDouble> fromCoeffs =
pFields[j]->GetPlane(n)->GetCoeffs() +
pFields[j]->GetCoeff_Offset(expId);
exp->GetPlane(n)->GetExp(0)->
BwdTrans(fromCoeffs, toPhys);
}
}
exp->HomogeneousBwdTrans(exp->GetPhys(), exp->UpdatePhys());
// Interpolate data
if (plane != -1)
{
physvals = exp->GetPlane(plane)->UpdatePhys();
value = exp->GetPlane(plane)->GetExp(0)->
StdPhysEvaluate(locCoord,physvals);
}
}
// store data
if (plane != -1)
{
data[m_historyLocalPointMap[k]*numFields+j] = value;
}
}
}
else
......
......@@ -80,9 +80,11 @@ class FilterHistoryPoints : public Filter
SpatialDomains::PointGeomVector m_historyPoints;
unsigned int m_index;
unsigned int m_outputFrequency;
/// plane to take history point from if using a homogeneous1D expansion
/// plane to take history point from if using a homogeneous1D expansion
unsigned int m_outputPlane;
Array<OneD, int> m_planeIDs;
bool m_isHomogeneous1D;
bool m_waveSpace;
std::string m_outputFile;
std::ofstream m_outputStream;
std::stringstream m_historyPointStream;
......
......@@ -54,6 +54,8 @@ if (NEKTAR_BUILD_UTILITIES)
endif()
if (NEKTAR_BUILD_LIBRARY AND DPKG)
set(NEK_DEP "openmpi-bin")
### DEBIAN PACKAGES ############################################
# Note that the formatting of the DESCRIPTION field is VERY specific. It
# must start with a new line, each line must start with a space, and there
......@@ -65,8 +67,9 @@ if (NEKTAR_BUILD_LIBRARY AND DPKG)
DESCRIPTION "
This library provides core routines including linear algebra and integration
with ThirdParty libraries."
DEPENDS "${NEK_DEP}"
INSTALL_LIBS "${libnektar++-utilities_LIBS}")
set(NEK_DEP "libnektar++-utilities (= ${NEKTAR_VERSION})")
set(NEK_DEP "${NEK_DEP}, libnektar++-utilities (= ${NEKTAR_VERSION})")
add_nektar_package(
NAME libnektar++-stdregions
......@@ -473,6 +476,7 @@ set(CPACK_SOURCE_IGNORE_FILES
"/\\\\.gitmodules"
"/build/"
"/builds/"
"/docs/tutorial"
"/ThirdParty/"
"/Testing/"
"/library/Demos/MultiRegions/ExtraDemos/"
......
......@@ -495,77 +495,59 @@ void AdjointAdvection::v_SetBaseFlow(
* coefficient storage.
* @param infile Filename to read.
*/
void AdjointAdvection::ImportFldBase(std::string pInfile,
Array<OneD, MultiRegions::ExpListSharedPtr>& pFields, int slice)
void AdjointAdvection::ImportFldBase(
std::string pInfile,
Array<OneD, MultiRegions::ExpListSharedPtr> &pFields,
int pSlice)
{
std::vector<LibUtilities::FieldDefinitionsSharedPtr> FieldDef;
std::vector<std::vector<NekDouble> > FieldData;
std::vector<std::vector<NekDouble> > FieldData;
int nqtot = pFields[0]->GetTotPoints();
int nqtot = m_baseflow[0].num_elements();
Array<OneD, NekDouble> tmp_coeff(pFields[0]->GetNcoeffs(), 0.0);
//Get Homogeneous
LibUtilities::FieldIOSharedPtr fld =
MemoryManager<LibUtilities::FieldIO>::AllocateSharedPtr(
m_session->GetComm());
fld->Import(pInfile, FieldDef, FieldData);
int numexp = pFields[0]->GetExpSize();
Array<OneD,int> ElementGIDs(numexp);
int nvar = m_session->GetVariables().size();
int s;