Skip to content

GitLab

Commit 8a729181 authored by Chris Cantwell's avatar Chris Cantwell
Browse files

Merge branch 'master' into fix/FCVerbose

parents c11845b3 ffa1489f
Hide whitespace changes
Inline Side-by-side
CHANGELOG.md
View file @ 8a729181
Changelog
=========
v4.4.0
------
**Library:**
- Add support for variable polynomial order for 3D simulations with continuous
Galerkin discretisation (!604)
**IncNavierStokesSolver:**
- Add ability to simulate additional scalar fields (!624)
**NekMesh:**
- Modify curve module to allow for spline input (!628)
**FieldConvert:**
- Add module to stretch homogeneous direction (!609)
v4.3.3
------
**FieldConvert**
**Library**:
- Fix filters when using adaptive driver to avoid output being overwritten after
each adaptive update (!588)
- Minor fix to suppress Xxt output unless `--verbose` is specified (!642)
- Fix of DirectFull solver in case where only Neumann boundary conditions
are imposed. (!655)
**FieldConvert**:
- Fix to avoid repeated import of field file (!649)
- Fix issue with C^0 projection (!644)
- Fix verbose output when using --procid (!648)
**CompressibleFlowSolver**:
- Fix issue with residual output (!647)
- Issues with 1D Euler solver fixed (!565)
- Fix deadlocking issue with boundary conditions (!657)
**Packaging**:
- Fix NekMesh dependencies for DEB package (!650)
>>>>>>> master
v4.3.2
------
**Library**:
......
cmake/ThirdPartyBoost.cmake
View file @ 8a729181
......@@ -8,9 +8,9 @@
#If the user has not set BOOST_ROOT, look in a couple common places first.
MESSAGE(STATUS "Searching for Boost:")
SET(MIN_VER "1.52.0")
SET(MIN_VER "1.56.0")
SET(NEEDED_BOOST_LIBS thread iostreams date_time filesystem system
program_options regex timer)
program_options regex timer chrono)
SET(Boost_DEBUG 0)
SET(Boost_NO_BOOST_CMAKE ON)
IF( BOOST_ROOT )
......@@ -67,7 +67,7 @@ IF (THIRDPARTY_BUILD_BOOST)
# Only build the libraries we need
SET(BOOST_LIB_LIST --with-system --with-iostreams --with-filesystem
--with-program_options --with-date_time --with-thread
--with-regex --with-timer)
--with-regex --with-timer --with-chrono)
IF (NOT WIN32)
# We need -fPIC for 64-bit builds
......@@ -160,6 +160,9 @@ IF (THIRDPARTY_BUILD_BOOST)
ENDIF(THIRDPARTY_BUILD_ZLIB)
# Set up CMake variables
SET(Boost_CHRONO_LIBRARY boost_chrono)
SET(Boost_CHRONO_LIBRARY_DEBUG boost_chrono)
SET(Boost_CHRONO_LIBRARY_RELEASE boost_chrono)
SET(Boost_DATE_TIME_LIBRARY boost_date_time)
SET(Boost_DATE_TIME_LIBRARY_DEBUG boost_date_time)
SET(Boost_DATE_TIME_LIBRARY_RELEASE boost_date_time)
......@@ -189,7 +192,7 @@ IF (THIRDPARTY_BUILD_BOOST)
SET(Boost_CONFIG_INCLUDE_DIR ${TPINC})
SET(Boost_LIBRARY_DIRS ${TPSRC}/dist/lib)
SET(Boost_CONFIG_LIBRARY_DIR ${TPLIB})
SET(Boost_LIBRARIES boost_date_time boost_filesystem boost_iostreams boost_program_options boost_regex boost_system boost_thread boost_timer)
SET(Boost_LIBRARIES boost_chrono boost_date_time boost_filesystem boost_iostreams boost_program_options boost_regex boost_system boost_thread boost_timer)
LINK_DIRECTORIES(${Boost_LIBRARY_DIRS})
STRING(REPLACE ";" ", " NEEDED_BOOST_LIBS_STRING "${NEEDED_BOOST_LIBS}")
......
docs/user-guide/solvers/incompressible-ns.tex
View file @ 8a729181
......@@ -1107,6 +1107,61 @@ 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{Advecting extra passive scalar fields}
In some cases, it might be useful to advect passive scalar fields with the velocity obtained from the
solution of the Navier-Stokes equation. For example, in study of mass transfer or heat transfer problems
where getting analytical expression for advection velocity is not possible, the transport (advection-diffusion)
equation needs to be solved along with the Navier-Stokes equation to get the scalar concentration or
temperature distribution in the flow field.
In the input file, the extra field variables that are being advected need to be defined after the variables
representing the velocity components. The pressure needs to be at the end of the list. For example, for
a 2D simulation the expansions and variables would be defined as
\begin{lstlisting}[style=XMLStyle]
<EXPANSIONS>
<E COMPOSITE="C[0]" NUMMODES="5" FIELDS="u,v,c1,c2,p" TYPE="MODIFIED" />
</EXPANSIONS>
<VARIABLES>
<V ID="0"> u </V>
<V ID="1"> v </V>
<V ID="2"> c1 </V>
<V ID="3"> c2 </V>
<V ID="4"> p </V>
</VARIABLES>
\end{lstlisting}
where u, v are the velocity components, c1 and c2 are extra fields that are being advected and p is the pressure.
In addition, diffusion coefficients for each extra variable can be specified by adding a function
\inltt{DiffusionCoefficient}
\begin{lstlisting}[style=XMLStyle]
<FUNCTION NAME="DiffusionCoefficient">
<E VAR="c1" VALUE="0.1" />
<E VAR="c2" VALUE="0.01" />
</FUNCTION>
\end{lstlisting}
Boundary conditions for the extra fields are set up in the same way as the velocity and pressure
\begin{lstlisting}[style=XMLStyle]
<BOUNDARYCONDITIONS>
<REGION REF="0">
<D VAR="u" VALUE="0" />
<D VAR="v" VALUE="0" />
<D VAR="c1" VALUE="1" />
<D VAR="c2" VALUE="0" />
<N VAR="p" USERDEFINEDTYPE="H" VALUE="0" />
</REGION>
</BOUNDARYCONDITIONS>
\end{lstlisting}
It should be noted that if the diffusion coefficient is too small, the transport equation becomes advection dominated.
This leads to small grid spacing required to resolve all physical scales associated with the transport equation
(the ratio of resolution required for transport to Navier Stokes equation scales with $Sc^{3/4}$, where $Sc$ is the
Schmidt number = kinematic viscosity/diffusion coefficient). Hence, small diffusion coefficient might lead to spurious
oscillations if the mesh spacing is not small enough.
\section{Examples}
\subsection{Kovasznay Flow 2D}
......
docs/user-guide/utilities/fieldconvert.tex
View file @ 8a729181
......@@ -102,6 +102,7 @@ possibly also Reynolds stresses) into single file;
\item \inltt{equispacedoutput}: Write data as equi-spaced output using simplices to represent the data for connecting points;
\item \inltt{extract}: Extract a boundary field;
\item \inltt{homplane}: Extract a plane from 3DH1D expansions;
\item \inltt{homstretch}: Stretch a 3DH1D expansion by an integer factor;
\item \inltt{innerproduct}: take the inner product between one or a series of fields with another field (or series of fields).
\item \inltt{interpfield}: Interpolates one field to another, requires fromxml, fromfld to be defined;
\item \inltt{interppointdatatofld}: Interpolates given discrete data using a finite difference approximation to a fld file given an xml file;
......@@ -327,6 +328,21 @@ The output file \inltt{file-plane.fld} can be processed in a similar
way as described in section \ref{s:utilities:fieldconvert:sub:convert}
to visualise it either in Tecplot or in Paraview.
\subsection{Stretch a 3DH1D expansion: \textit{homstretch} module}
To stretch a 3DH1D expansion in the z-direction, use the command:
\begin{lstlisting}[style=BashInputStyle]
FieldConvert -m homstretch:factor=value file.xml file.fld file-stretch.fld
\end{lstlisting}
The number of modes in the resulting field can be chosen using the command-line
parameter \inltt{output-points-hom-z}. Note that the output for
this module should always be a \inltt{.fld} file and this should not
be used in combination with other modules using a single command.
The output file \inltt{file-stretch.fld} can be processed in a similar
way as described in section \ref{s:utilities:fieldconvert:sub:convert}
to visualise it either in Tecplot or in Paraview.
\subsection{Inner Product of a single or series of fields with respect to a single or series of fields: \textit{innerproduct} module}
You can take the inner product of one field with another field using
......
library/LibUtilities/BasicUtils/PtsField.cpp
View file @ 8a729181
......@@ -36,169 +36,16 @@
#include <LibUtilities/BasicUtils/PtsField.h>
using namespace std;
namespace Nektar
{
namespace LibUtilities
{
/**
* @brief Compute the weights for an interpolation of the field values to physical points
*
* @param physCoords coordinates of the physical points
* @param coord_id id of the coordinate to use for interpolation.
*
* Set coord_id to -1 to use n-D interpolation for an n-dimensional field.
* The most suitable algorithm is chosen automatically.
*/
void PtsField::CalcWeights(
const Array<OneD, Array<OneD, NekDouble> > &physCoords,
short coordId)
{
ASSERTL1(physCoords.num_elements() >= m_dim,
"physCoords is smaller than number of dimesnions");
int nPhysPts = physCoords[0].num_elements();
int lastProg = 0;
m_weights = Array<OneD, Array<OneD, float> >(nPhysPts);
m_neighInds = Array<OneD, Array<OneD, unsigned int> >(nPhysPts);
// interpolate points and transform
for (int i = 0; i < nPhysPts; ++i)
{
Array<OneD, NekDouble> physPt(m_dim);
for (int j = 0; j < m_dim; ++j)
{
physPt[j] = physCoords[j][i];
}
if (m_dim == 1 || coordId >= 0)
{
if (m_dim == 1)
{
coordId = 0;
}
if (m_pts[0].num_elements() <= 2)
{
CalcW_Linear(i, physPt[coordId]);
}
else
{
CalcW_Quadratic(i, physPt[coordId]);
}
}
else
{
CalcW_Shepard(i, physPt);
}
int progress = int(100 * i / nPhysPts);
if (m_progressCallback && progress > lastProg)
{
m_progressCallback(i, nPhysPts);
lastProg = progress;
}
}
}
/**
* @brief Compute weights and perform the interpolate of field values to physical points
*
* @param physCoords coordinates of the physical points
* @param intFields interpolated field at the physical points
* @param coord_id id of the coordinate to use for interpolation.
*
* Set coord_id to -1 to use n-D interpolation for an n-dimensional field.
* The most suitable algorithm is chosen automatically.
*/
void PtsField::Interpolate(
const Array< OneD, Array< OneD, NekDouble > > &physCoords,
Array<OneD, Array<OneD, NekDouble> > &intFields,
short int coordId)
{
CalcWeights(physCoords, coordId);
Interpolate(intFields);
}
/**
* @brief Perform the interpolate of field values to physical points
*
* @param intFields interpolated field at the physical points
*
* The weights must have already been computed by @CalcWeights or set by
* @SetWeights.
*/
void PtsField::Interpolate(Array<OneD, Array<OneD, NekDouble> > &intFields)
{
ASSERTL1(m_weights[0].num_elements() == m_neighInds[0].num_elements(),
"weights / neighInds mismatch")
int nFields = m_fieldNames.size();
int nPhysPts = m_weights.num_elements();
// interpolate points and transform
intFields = Array<OneD, Array<OneD, NekDouble> >(nFields);
for (int i = 0; i < nFields; ++i)
{
intFields[i] = Array<OneD, NekDouble>(nPhysPts);
for (int j = 0; j < nPhysPts; ++j)
{
intFields[i][j] = 0.0;
int nPts = m_weights[j].num_elements();
for (int k = 0; k < nPts; ++k)
{
unsigned int nIdx = m_neighInds[j][k];
intFields[i][j] += m_weights[j][k] * m_pts[m_dim + i][nIdx];
}
}
}
}
/**
* @brief Set the interpolation weights for an interpolation
*
* @param weights Interpolation weights for each neighbour.
* Structure: m_weights[physPtIdx][neighbourIdx]
* @param neighbourInds Indices of the relevant neighbours for each physical point.
* Structure: m_neighInds[ptIdx][neighbourIdx]
*/
void PtsField::SetWeights(
const Array< OneD, Array< OneD, float > > &weights,
const Array< OneD, Array< OneD, unsigned int > > &neighbourInds)
{
ASSERTL0(weights.num_elements() == neighbourInds.num_elements(),
"weights and neighbourInds not of same number of physical points")
m_weights = weights;
m_neighInds = neighbourInds;
}
/**
* @brief Get the interpolation weights and corresponding neighbour indices
*
* @param weights Interpolation weights for each neighbour.
* Structure: m_weights[physPtIdx][neighbourIdx]
* @param neighbourInds Indices of the relevant neighbours for each physical point.
* Structure: m_neighInds[ptIdx][neighbourIdx]
*/
void PtsField::GetWeights(
Array< OneD, Array< OneD, float > > &weights,
Array< OneD, Array< OneD, unsigned int > > &neighbourInds) const
{
weights = m_weights;
neighbourInds = m_neighInds;
}
PtsField::PtsField(const int dim, const Array< OneD, Array< OneD, NekDouble > > &pts):
m_dim(dim),
m_pts(pts),
m_ptsType(ePtsFile)
PtsField::PtsField(const int dim,
const Array<OneD, Array<OneD, NekDouble> > &pts)
: m_dim(dim), m_pts(pts), m_ptsType(ePtsFile)
{
for (int i = 0; i < GetNFields(); ++i)
{
......@@ -206,8 +53,6 @@ PtsField::PtsField(const int dim, const Array< OneD, Array< OneD, NekDouble > >
}
}
/**
* @brief Set the connectivity data for ePtsTetBlock and ePtsTriBlock
*
......@@ -237,46 +82,40 @@ void PtsField::SetConnectivity(const vector< Array< OneD, int > > &conn)
m_ptsConn = conn;
}
void PtsField::SetDim(const int ptsDim)
{
m_dim = ptsDim;
}
int PtsField::GetDim() const
{
return m_dim;
}
int PtsField::GetNFields() const
{
return m_fieldNames.size();
return m_pts.num_elements() - m_dim;
}
vector<std::string> PtsField::GetFieldNames() const
{
return m_fieldNames;
}
std::string PtsField::GetFieldName(const int i) const
{
return m_fieldNames[i];
}
void PtsField::SetFieldNames(const vector<std::string> fieldNames)
{
ASSERTL0(fieldNames.size() == m_pts.num_elements() - m_dim,
"Number of given fieldNames does not match the number of stored fields");
"Number of given fieldNames does not match the number of stored "
"fields");
m_fieldNames = fieldNames;
}
void PtsField::AddField(const Array< OneD, NekDouble > &pts,
const string fieldName)
{
......@@ -298,31 +137,54 @@ void PtsField::AddField(const Array< OneD, NekDouble > &pts,
m_fieldNames.push_back(fieldName);
}
void PtsField::AddPoints(const Array<OneD, const Array<OneD, NekDouble> > &pts)
{
ASSERTL1(pts.num_elements() == m_pts.num_elements(),
"number of variables mismatch");
// TODO: dont copy, dont iterate
for (int i = 0; i < m_pts.num_elements(); ++i)
{
Array<OneD, NekDouble> tmp(m_pts[i].num_elements() + pts[i].num_elements());
for (int j = 0; j < m_pts[i].num_elements(); ++j)
{
tmp[j] = m_pts[i][j];
}
for (int j = 0; j < pts[i].num_elements(); ++j)
{
tmp[m_pts[i].num_elements() + j] = pts[i][j];
}
m_pts[i] = tmp;
}
}
int PtsField::GetNpoints() const
{
return m_pts[0].num_elements();
}
NekDouble PtsField::GetPointVal(const int fieldInd, const int ptInd) const
{
return m_pts[fieldInd][ptInd];
}
void PtsField::SetPointVal(const int fieldInd,
const int ptInd,
const NekDouble val)
{
m_pts[fieldInd][ptInd] = val;
}
void PtsField::GetPts(Array< OneD, Array< OneD, NekDouble > > &pts) const
{
pts = m_pts;
}
Array< OneD, NekDouble > PtsField::GetPts(const int fieldInd) const
{
return m_pts[fieldInd];
}
void PtsField::SetPts(Array< OneD, Array< OneD, NekDouble > > &pts)
{
ASSERTL1(pts.num_elements() == m_pts.num_elements(),
......@@ -331,13 +193,11 @@ void PtsField::SetPts(Array< OneD, Array< OneD, NekDouble > > &pts)
m_pts = pts;
}
vector<int> PtsField::GetPointsPerEdge() const
{
return m_nPtsPerEdge;
}
int PtsField::GetPointsPerEdge(const int i) const
{
return m_nPtsPerEdge[i];
......@@ -352,20 +212,18 @@ int PtsField::GetPointsPerEdge(const int i) const
*/
void PtsField::SetPointsPerEdge(const vector< int > nPtsPerEdge)
{
ASSERTL0(m_ptsType == ePtsLine || m_ptsType == ePtsPlane ||
m_ptsType == ePtsBox,
ASSERTL0(
m_ptsType == ePtsLine || m_ptsType == ePtsPlane || m_ptsType == ePtsBox,
"SetPointsPerEdge only supported for ePtsLine, ePtsPlane and ePtsBox.");
m_nPtsPerEdge = nPtsPerEdge;
}
PtsType PtsField::GetPtsType() const
{
return m_ptsType;
}
void PtsField::SetPtsType(const PtsType type)
{
m_ptsType = type;
......@@ -380,253 +238,5 @@ void PtsField::SetBoxSize(const vector< NekDouble> boxSize)
{
m_boxSize = boxSize;
}
/**
* @brief Compute interpolation weights for a 1D physical point using linear
* interpolation.
*
* @param physPtIdx The index of the physical point in its storage array
* @param coord The coordinate of the physical point
*/
void PtsField::CalcW_Linear(const int physPtIdx, const NekDouble coord)
{
int npts = m_pts[0].num_elements();
int i;
int numPts = 2;
m_neighInds[physPtIdx] = Array<OneD, unsigned int> (numPts);
m_weights[physPtIdx] = Array<OneD, float> (numPts, 0.0);
for (i = 0; i < npts - 1; ++i)
{
if ((m_pts[0][i] <= coord) && (coord <= m_pts[0][i + 1]))
{
NekDouble pdiff = m_pts[0][i + 1] - m_pts[0][i];
m_neighInds[physPtIdx][0] = i;
m_neighInds[physPtIdx][1] = i + 1;
m_weights[physPtIdx][0] = (m_pts[0][i + 1] - coord) / pdiff;
m_weights[physPtIdx][1] = (coord - m_pts[0][i]) / pdiff;
break;
}
}
ASSERTL0(i != npts - 1, "Failed to find coordinate " +
boost::lexical_cast<string>(coord) +
" within provided input points");
};
/**
* @brief Compute interpolation weights for a physical point using a modified
* Shepard algorithm.
*
* @param physPtIdx The index of the physical point in its storage array
* @param physPt The coordinates of the physical point
*
* The algorithm is based on Shepard, D. (1968). A two-dimensional interpolation
* function for irregularly-spaced data. Proceedings of the 1968 23rd ACM
* National Conference. pp. 517–524.
*
* In order to save memory, for n dimesnions, only 2^n points are considered.
* Contrary to Shepard, we use a fixed number of points with fixed weighting
* factors 1/d^n.
*/
void PtsField::CalcW_Shepard(const int physPtIdx,
const Array<OneD, NekDouble> &physPt)
{
// find nearest neighbours
vector<PtsPoint> neighbourPts;
int numPts = pow(float(2), m_dim);
numPts = min(numPts, int(m_pts[0].num_elements() / 2));
FindNeighbours(physPt, neighbourPts, numPts);
m_neighInds[physPtIdx] = Array<OneD, unsigned int> (numPts);
for (int i = 0; i < numPts; i++)
{
m_neighInds[physPtIdx][i] = neighbourPts.at(i).m_idx;
}
m_weights[physPtIdx] = Array<OneD, float> (numPts, 0.0);
// In case d < kVertexTheSameDouble ( d^2 < kNekSqrtTol), use the exact
// point and return
for (int i = 0; i < numPts; ++i)
{
if (neighbourPts[i].m_distSq <= NekConstants::kNekSqrtTol)
{
m_weights[physPtIdx][i] = 1.0;
return;
}
}
NekDouble wSum = 0.0;
for (int i = 0; i < numPts; ++i)
{
m_weights[physPtIdx][i] = 1 / pow(double(neighbourPts[i].m_distSq),
double(m_dim / float(2)));
wSum += m_weights[physPtIdx][i];
}
for (int i = 0; i < numPts; ++i)
{
m_weights[physPtIdx][i] = m_weights[physPtIdx][i] / wSum;
}
ASSERTL0(Vmath::Nnan(numPts, m_weights[physPtIdx], 1) == 0, "NaN found in weights");