Commit 19715c63 authored by Dave Moxey's avatar Dave Moxey
Browse files

Merge remote-tracking branch 'upstream/master' into feature/FCHomogeneousStretch

parents 829bdb12 7036708b
Changelog
=========
v4.4.0
------
**IncNavierStokesSolver:**
- Add ability to simulate additional scalar fields (!624)
**NekMesh:**
- Modify curve module to allow for spline input (!628)
v4.3.3
------
**Library**:
- 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)
**CompressibleFlowSolver**:
- Fix issue with residual output (!647)
**Packaging**:
- Fix NekMesh dependencies for DEB package (!650)
v4.3.2
------
**Library**:
......@@ -10,9 +35,17 @@ v4.3.2
- 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)
- Fix `IO_CheckSteps` to avoid missing first checkpoint (!639)
- Fix bug in iterative solver where only root process would ASSERT when
exceeding the maximum number of iterations (!636)
**FieldConvert**:
- Fix appearence of duplicate messages when running in parallel (!626)
- Fix issue with efficiency when using large number of 3DH1D planes (!627)
- Add module for combining average fields (!620)
- Fix wall shear stress processing module for parallel execution (!635)
**Packaging**:
- Fixes for DEB package dependencies (!630)
......
Subproject commit 244b5128da8a9c9baedf71230db81344c790b03b
Subproject commit a20044f232e15e3fe2754a90215f8646db241147
......@@ -34,7 +34,7 @@ $\epsilon \nabla^2 u + \mathbf{V}\nabla u = f$ &
\inltt{SteadyAdvectionDiffusion} & 2D only & Continuous/Discontinuous \\
$\epsilon \nabla^2 u + \lambda u = f$ &
\inltt{SteadyDiffusionReaction} & 2D only & Continuous/Discontinuous \\
$\epsilon \nabla^2 u \mathbf{V}\nabla u + \lambda u = f$ &
$\epsilon \nabla^2 u + \mathbf{V}\nabla u + \lambda u = f$ &
\inltt{SteadyAdvectionDiffusionReaction} & 2D only &
Continuous/Discontinuous \\
$ \dfrac{\partial u}{\partial t} + \mathbf{V}\nabla u = f$ &
......
......@@ -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]
......@@ -1103,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}
......
......@@ -96,6 +96,8 @@ Specifically, FieldConvert has these additional functionalities
\item \inltt{C0Projection}: Computes the C0 projection of a given output file;
\item \inltt{QCriterion}: Computes the Q-Criterion for a given output file;
\item \inltt{addFld}: Sum two .fld files;
\item \inltt{combineAvg}: Combine two \nekpp binary output (.chk or .fld) field file containing averages of fields (and
possibly also Reynolds stresses) into single file;
\item \inltt{concatenate}: Concatenate a \nekpp binary output (.chk or .fld) field file 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;
......@@ -203,6 +205,21 @@ the result either in Tecplot, Paraview or VisIt.
%
%
%
\subsection{Combine two .fld files containing time averages: \textit{combineAvg} module}
To combine two .fld files obtained through the AverageFields or ReynoldsStresses filters,
use the \inltt{combineAvg} module of FieldConvert
%
\begin{lstlisting}[style=BashInputStyle]
FieldConvert -m combineAvg:fromfld=file1.fld file1.xml file2.fld \
file3.fld
\end{lstlisting}
%
\inltt{file3.fld} can be processed in a similar way as described
in section \ref{s:utilities:fieldconvert:sub:convert} to visualise
the result either in Tecplot, Paraview or VisIt.
%
%
%
\subsection{Concatenate two files: \textit{concatenate} module}
To concatenate \inltt{file1.fld} and \inltt{file2.fld} into \inltt{file-conc.fld}
one can run the following command
......
......@@ -32,7 +32,7 @@ mesh-generator \gmsh to the XML file format.
converted into base64. This is identified for each section by the attribute
\inltt{COMPRESSED="B64Z-LittleEndian''}. To output
in ascii format add the module option ``:xml:uncompress'' to the .xml file,
i.e. \\ \inltt{ MeshConvert file.msh newfile.xml:xml:uncompress}
i.e. \\ \inltt{ \mc file.msh newfile.xml:xml:uncompress}
\end{notebox}
\section{Exporting a mesh from \gmsh}
......@@ -290,7 +290,7 @@ converted into base64. If you wish to see an ascii output you need to
specify the output module option \inltt{uncompress} by executing:
\begin{lstlisting}[style=BashInputStyle]
MeshConvert Mesh.msh output.xml:xml:uncompress
NekMesh Mesh.msh output.xml:xml:uncompress
\end{lstlisting}
In the rest of these subsections, we discuss the various processing modules
......@@ -302,13 +302,13 @@ To extract composite surfaces 2 and 3 from a mesh use the module \inltt{extract
module:
%
\begin{lstlisting}[style=BashInputStyle]
MeshConvert -m extract:surf=2,3 Mesh.xml output.xml
NekMesh -m extract:surf=2,3 Mesh.xml output.xml
\end{lstlisting}
%
If you also wish to have the boundaries of the extracted surface detected add the \inltt{detectbnd} option
%
\begin{lstlisting}[style=BashInputStyle]
MeshConvert -m extract:surf=2,3:detectbnd Mesh.xml output.xml
NekMesh -m extract:surf=2,3:detectbnd Mesh.xml output.xml
\end{lstlisting}
\subsection{Negative Jacobian detection}
......
......@@ -48,7 +48,7 @@ The description below explains how the \inltt{GEOMETRY} section is laid out in
uncompressed ascii format. From Jan 2016 the distribution uses the compressed
format for each of the above sections. To convert a compressed xml file into
ascii format use \\
\inltt{ MeshConvert file.msh newfile.xml:xml:uncompress}
\inltt{ NekMesh file.msh newfile.xml:xml:uncompress}
\end{notebox}
......
......@@ -41,6 +41,10 @@
#include <boost/optional.hpp>
#include <LibUtilities/LibUtilitiesDeclspec.h>
#if defined(NEKTAR_USE_MPI)
#include <mpi.h>
#endif
#ifndef _WIN32
#include <execinfo.h>
#endif
......@@ -63,7 +67,7 @@ namespace ErrorUtil
{
efatal,
ewarning
};
};
class NekError : public std::runtime_error
{
......@@ -86,6 +90,18 @@ namespace ErrorUtil
#endif
msg;
// Default rank is zero. If MPI used and initialised, populate with
// the correct rank. Messages are only printed on rank zero.
int rank = 0;
#if defined(NEKTAR_USE_MPI)
int flag;
MPI_Initialized(&flag);
if(flag)
{
MPI_Comm_rank(MPI_COMM_WORLD,&rank);
}
#endif
std::string btMessage("");
#if defined(NEKTAR_FULLDEBUG)
#ifndef _WIN32
......@@ -103,10 +119,13 @@ namespace ErrorUtil
free(btStrings);
#endif
#endif
switch(type)
switch (type)
{
case efatal:
if( outStream )
case efatal:
if (!rank)
{
if (outStream)
{
(*outStream) << btMessage;
(*outStream) << "Fatal : " << baseMsg << std::endl;
......@@ -114,14 +133,22 @@ namespace ErrorUtil
else
{
std::cerr << btMessage;
std::cerr << std::endl << "Fatal : " << baseMsg << std::endl;
std::cerr << std::endl << "Fatal : " << baseMsg
<< std::endl;
}
throw NekError(baseMsg);
break;
case ewarning:
if( outStream )
}
#if defined(NEKTAR_USE_MPI)
if (flag)
{
MPI_Barrier(MPI_COMM_WORLD);
}
#endif
throw NekError(baseMsg);
break;
case ewarning:
if (!rank)
{
if (outStream)
{
(*outStream) << btMessage;
(*outStream) << "Warning: " << baseMsg << std::endl;
......@@ -131,10 +158,10 @@ namespace ErrorUtil
std::cerr << btMessage;
std::cerr << "Warning: " << baseMsg << std::endl;
}
break;
default:
std::cerr << "Unknown warning type: " << baseMsg << std::endl;
}
break;
default:
std::cerr << "Unknown warning type: " << baseMsg << std::endl;
}
}
......
......@@ -402,7 +402,8 @@ void PtsField::CalcW_Linear(const int physPtIdx, const NekDouble coord)
for (i = 0; i < npts - 1; ++i)
{
if ((m_pts[0][i] <= coord) && (coord <= m_pts[0][i + 1]))
if ((m_pts[0][i] <= (coord + NekConstants::kNekZeroTol)) &&
(coord <= (m_pts[0][i + 1]+ NekConstants::kNekZeroTol)) )
{
NekDouble pdiff = m_pts[0][i + 1] - m_pts[0][i];
......@@ -499,7 +500,8 @@ void PtsField::CalcW_Quadratic(const int physPtIdx, const NekDouble coord)
for (i = 0; i < npts - 1; ++i)
{
if ((m_pts[0][i] <= coord) && (coord <= m_pts[0][i + 1]))
if ((m_pts[0][i] <= (coord + NekConstants::kNekZeroTol)) &&
(coord <= (m_pts[0][i + 1]+ NekConstants::kNekZeroTol)) )
{
NekDouble pdiff = m_pts[0][i + 1] - m_pts[0][i];
NekDouble h1, h2, h3;
......
......@@ -302,12 +302,12 @@ namespace Nektar
int retval = MPI_Sendrecv(pSendData.get(),
(int) pSendData.num_elements(),
MPI_DOUBLE,
pRecvProc,
pSendProc,
0,
pRecvData.get(),
(int) pRecvData.num_elements(),
MPI_DOUBLE,
pSendProc,
pRecvProc,
0,
m_comm,
&status);
......
......@@ -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++)
......
......@@ -1217,7 +1217,7 @@ namespace Nektar
<< ntot << "\" NumberOfCells=\""
<< ntotminus << "\">" << endl;
outfile << " <Points>" << endl;
outfile << " <DataArray type=\"Float32\" "
outfile << " <DataArray type=\"Float64\" "
<< "NumberOfComponents=\"3\" format=\"ascii\">" << endl;
outfile << " ";
for (i = 0; i < ntot; ++i)
......
......@@ -367,9 +367,24 @@ namespace Nektar
int expansion,
int istrip)
{
// If there is only one plane (e.g. HalfMode), we write a 2D plane.
if (m_planes.num_elements() == 1)
{
m_planes[0]->WriteVtkPieceHeader(outfile, expansion);
return;
}
// If we are using Fourier points, output extra plane to fill domain
int outputExtraPlane = 0;
if ( m_homogeneousBasis->GetBasisType() == LibUtilities::eFourier
&& m_homogeneousBasis->GetPointsType() ==
LibUtilities::eFourierEvenlySpaced)
{
outputExtraPlane = 1;
}
int i, j;
int nquad0 = (*m_exp)[expansion]->GetNumPoints(0);
int nquad1 = m_planes.num_elements();
int nquad1 = m_planes.num_elements() + outputExtraPlane;
int ntot = nquad0*nquad1;
int ntotminus = (nquad0-1)*(nquad1-1);
......@@ -379,6 +394,20 @@ namespace Nektar
coords[2] = Array<OneD,NekDouble>(ntot);
GetCoords(expansion,coords[0],coords[1],coords[2]);
if (outputExtraPlane)
{
// Copy coords[0] and coords[1] to extra plane
Array<OneD,NekDouble> tmp;
Vmath::Vcopy (nquad0, coords[0], 1,
tmp = coords[0] + (nquad1-1)*nquad0, 1);
Vmath::Vcopy (nquad0, coords[1], 1,
tmp = coords[1] + (nquad1-1)*nquad0, 1);
// Fill coords[2] for extra plane
NekDouble z = coords[2][nquad0*m_planes.num_elements()-1] +
(coords[2][nquad0] - coords[2][0]);
Vmath::Fill(nquad0, z, tmp = coords[2] + (nquad1-1)*nquad0, 1);
}
outfile << " <Piece NumberOfPoints=\""
<< ntot << "\" NumberOfCells=\""
<< ntotminus << "\">" << endl;
......
......@@ -109,7 +109,7 @@ namespace Nektar
(*m_exp).push_back(m_planes[0]->GetExp(j));
}
for(n = 1; n < m_homogeneousBasis->GetNumPoints(); ++n)
for(n = 1; n < m_planes.num_elements(); ++n)
{
m_planes[n] = MemoryManager<ExpList2D>::AllocateSharedPtr(*plane_zero,False);
for(j = 0; j < nel; ++j)
......@@ -371,10 +371,19 @@ namespace Nektar
return;
}
// If we are using Fourier points, output extra plane to fill domain
int outputExtraPlane = 0;
if ( m_homogeneousBasis->GetBasisType() == LibUtilities::eFourier
&& m_homogeneousBasis->GetPointsType() ==
LibUtilities::eFourierEvenlySpaced)
{
outputExtraPlane = 1;
}
int i,j,k;
int nq0 = (*m_exp)[expansion]->GetNumPoints(0);
int nq1 = (*m_exp)[expansion]->GetNumPoints(1);
int nq2 = m_planes.num_elements();
int nq2 = m_planes.num_elements() + outputExtraPlane;
int ntot = nq0*nq1*nq2;
int ntotminus = (nq0-1)*(nq1-1)*(nq2-1);
......@@ -384,6 +393,20 @@ namespace Nektar
coords[2] = Array<OneD,NekDouble>(ntot);
GetCoords(expansion,coords[0],coords[1],coords[2]);
if (outputExtraPlane)
{
// Copy coords[0] and coords[1] to extra plane
Array<OneD,NekDouble> tmp;
Vmath::Vcopy (nq0*nq1, coords[0], 1,
tmp = coords[0] + (nq2-1)*nq0*nq1, 1);
Vmath::Vcopy (nq0*nq1, coords[1], 1,
tmp = coords[1] + (nq2-1)*nq0*nq1, 1);
// Fill coords[2] for extra plane
NekDouble z = coords[2][nq0*nq1*m_planes.num_elements()-1] +
(coords[2][nq0*nq1] - coords[2][0]);
Vmath::Fill(nq0*nq1, z, tmp = coords[2] + (nq2-1)*nq0*nq1, 1);
}
NekDouble DistStrip;
m_session->LoadParameter("DistStrip", DistStrip, 0);
// Reset the z-coords for homostrips
......@@ -396,7 +419,7 @@ namespace Nektar
<< ntot << "\" NumberOfCells=\""
<< ntotminus << "\">" << endl;
outfile << " <Points>" << endl;
outfile << " <DataArray type=\"Float32\" "
outfile << " <DataArray type=\"Float64\" "
<< "NumberOfComponents=\"3\" format=\"ascii\">" << endl;
outfile << " ";
for (i = 0; i < ntot; ++i)
......
......@@ -367,7 +367,7 @@ namespace Nektar
<< ntot << "\" NumberOfCells=\""
<< ntotminus << "\">" << endl;
outfile << " <Points>" << endl;
outfile << " <DataArray type=\"Float32\" "
outfile << " <DataArray type=\"Float64\" "
<< "NumberOfComponents=\"3\" format=\"ascii\">" << endl;
outfile << " ";
for (i = 0; i < ntot; ++i)
......
......@@ -872,8 +872,37 @@ namespace Nektar
int nq = (*m_exp)[expansion]->GetTotPoints();
int npoints_per_plane = m_planes[0]->GetTotPoints();
// If we are using Fourier points, output extra plane to fill domain
int outputExtraPlane = 0;
Array<OneD, NekDouble> extraPlane;
if ( m_homogeneousBasis->GetBasisType() == LibUtilities::eFourier
&& m_homogeneousBasis->GetPointsType() ==
LibUtilities::eFourierEvenlySpaced)
{
outputExtraPlane = 1;
// Get extra plane data
if (m_StripZcomm->GetSize() == 1)
{
extraPlane = m_phys + m_phys_offset[expansion];
}
else
{
// Determine to and from rank for communication
int size = m_StripZcomm->GetSize();
int rank = m_StripZcomm->GetRank();
int fromRank = (rank+1) % size;
int toRank = (rank == 0) ? size-1 : rank-1;
// Communicate using SendRecv
extraPlane = Array<OneD, NekDouble>(nq);
Array<OneD, NekDouble> send (nq,
m_phys + m_phys_offset[expansion]);
m_StripZcomm->SendRecv(toRank, send,
fromRank, extraPlane);
}
}
// printing the fields of that zone
outfile << " <DataArray type=\"Float32\" Name=\""
outfile << " <DataArray type=\"Float64\" Name=\""
<< var << "\">" << endl;
outfile << " ";
for (int n = 0; n < m_planes.num_elements(); ++n)
......@@ -884,6 +913,14 @@ namespace Nektar
outfile << (fabs(phys[i]) < NekConstants::kNekZeroTol ? 0 : phys[i]) << " ";
}
}
if (outputExtraPlane)
{
for(i = 0; i < nq; ++i)
{
outfile << (fabs(extraPlane[i]) < NekConstants::kNekZeroTol ?
0 : extraPlane[i]) << " ";
}
}
outfile << endl;
outfile << " </DataArray>" << endl;
}
......
......@@ -734,7 +734,7 @@ namespace Nektar
int npoints_per_line = m_lines[0]->GetTotPoints();
// printing the fields of that zone
outfile << " <DataArray type=\"Float32\" Name=\""
outfile << " <DataArray type=\"Float64\" Name=\""
<< var << "\">" << endl;
outfile << " ";
for (int n = 0; n < m_lines.num_elements(); ++n)
......
......@@ -192,7 +192,9 @@ namespace Nektar
&pLocToGloMap):
m_linSysKey(pKey),
m_expList(pExpList),
m_robinBCInfo(m_expList.lock()->GetRobinBCInfo())
m_robinBCInfo(m_expList.lock()->GetRobinBCInfo()),
m_verbose(m_expList.lock()->GetSession()->
DefinesCmdLineArgument("verbose"))
{
}
......
......@@ -129,6 +129,8 @@ namespace Nektar
const boost::weak_ptr<ExpList> m_expList;
/// Robin boundary info
const std::map<int, RobinBCInfoSharedPtr> m_robinBCInfo;
// Provide verbose output
bool m_verbose;
virtual int v_GetNumBlocks ();
virtual DNekScalMatSharedPtr v_GetBlock (unsigned int n);
......
......@@ -65,7 +65,7 @@ namespace Nektar
const Array<OneD,const NekDouble> &pInput,
Array<OneD, NekDouble> &pOutput,
const AssemblyMapSharedPtr &locToGloMap,
const int pNumDir = 0);
const int pNumDir);
};
}
}
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment