Commit 07a0458e authored by Dave Moxey's avatar Dave Moxey

Merge branch 'master' into fix/mesh-partition

parents fa1f7400 952804d8
......@@ -334,8 +334,8 @@ year={2011}
publisher={Birkh{\"a}user}
}
@article{DoKa14,
title={A robust and accurate outflow boundary condition for incompressible flow
simulations on severely-truncated unbounded domains},
title={A robust and accurate outflow boundary condition for incompressible
flow simulations on severely-truncated unbounded domains},
author={S. Dong and G. E. Karniadakis and C. Chryssostomidis},
journal={Journal of Computational Physics},
volume={261},
......@@ -352,4 +352,4 @@ year={2011}
pages = {157--172},
year = 2007,
publisher = {IOS Press}
}
\ No newline at end of file
}
......@@ -1573,6 +1573,71 @@ It is possible to visualise the wall shear stress distribution by running the Fl
\includegraphics[width=12cm]{Figures/WSS.png}
\caption{Non-dimensional wall shear stress distribution.}
\end{center}
\end{figure}
\subsection{finite-strip modeling of flow past flexible cables}
As a computationally efficient model, strip theory-based modeling technique has been proposed previously to predict vortex-induced vibration (VIV) for higher Reynolds number flows. In the strip theory-based model, the fluid flow solution is obtained on a series of 2D computational planes (also called as “strips”) along the riser’s axis direction. These strips then are coupled with each other through structural dynamic model of the riser, and then VIV response prediction is achieved by the strip-structure interactions.In the 2D strip theory, it is assumed that the flow is purely two-dimensional without spanwise correlation, which allows the problems to be split into various 2D planes. A consequence of 2D strip solution under this assumption is that it is unable to reflect the influence of spanwise wake turbulence on the structural dynamics. In order to overcome this shortcoming, we proposed a new module in the framework of Nektar++, in which a spanwise scale is locally allocated to each one of the strips, so that the spanwise velocity correlation is reconstructed in the flow field within each strips. In particular, this model lets the fluid domain to be divided in $N$ strips with thickness ratio of $L_{z}/D$ and evenly distributed along the spanwise ($z$) direction. The gap between the neighboring strips, represented by $L_{g}$, satisfies relation $L_{c}=N(L_{z}+L_{g})$. Since the strip in this model has finite scale in the $z$-direction, we named it as finite strip to distinguish from traditional 2D strip plane. Next, the flow dynamics within each individual strips are modeled by viscous incompressible Navier-Stokes equations, while a tensioned beam model is employed to govern the dynamics of the flexible structures. In this example, we will show how to perform a finite-strip model to predict the vortex-induced vibration responses of flexible cables. Let us consider a vortex-induced vibration of a slender cable with an aspect ratio of $L_z/D$=4$\pi$, which is immersed in uniform flows at Re=100.
\subsubsection{Input File}
The cable with a mass ratio (defined as the ratio of the total oscillating mass to the mass of displaced fluid) of 1 has diameter $D=1$, the 2D mesh is composed of 284 quadrilateral elements. The spanwise direction is split in 16 strips with thickness ratio of $L_c/D$=$\pi$/8 and one pair of complex Fourier modes for each one of the strips. We will use a sixth order polynomial expansion for the spectral element and the input file for this example is \inlsh{CylFlow\_HomoStrip.xml}.
\begin{lstlisting}[style=XMLStyle]
<E COMPOSITE="C[73]" NUMMODES="6" TYPE="MODIFIED" FIELDS="u,v,w,p" />
\end{lstlisting}
To use the finite strip routines we need just to insert a flag of \inlsh{"HomoStrip"} in the solver information as below, in addition, we need to specify the types of vibration and support ends for the cables. In this case, the vibration type is specified as \inlsh{VALUE="CONSTRAINED"}, which means that the cable's vibration is constrained only in the crossflow direction. Other options include \inlsh{VALUE="FREE"} and \inlsh{"FORCED"}, respectively corresponding to the free vibrations in both streamwise and crossflow directions and forced vibration by specified functions given in input file. For the support ends of the cable, another option of \inlsh{VALUE="PINNED-PINNED"} is available for the simulations, which satisfies the condition of zero values of displacements on the support ends.
\subsubsection{Solver information:~}
\begin{lstlisting}[style=XMLStyle]
<SOLVERINFO>
<I PROPERTY="HomoStrip" VALUE="True"/>
<I PROPERTY="VibrationType" VALUE="CONSTRAINED"/>
<I PROPERTY="SupportType" VALUE="FREE-FREE"/>
</SOLVERINFO>
\end{lstlisting}
\subsubsection{Parameters}
All the simulation parameters are specified in the section as follows.
\begin{lstlisting}[style=XMLStyle]
<PARAMETERS>
<P> LZ = PI/8 </P> <!--thickness ratio-->
<P> LC = 4*PI </P> <!--aspect ratio-->
<P> A = 0.025 </P>
<P> omega = 1.0 </P>
<P> PROC_Z = 16 </P>
<P> Strip_Z = 16 </P> <!--number of the strips-->
<P> DistStrip = PI/4 </P> <!--distance of the strips-->
<P> StructStiff = 0.02 </P>
<P> StructRho = 2.0 </P>
<P> CableTension = 8.82 </P>
<P> BendingStiff = 0.0 </P>
<P> FictDamp = 0.0 </P>
<P> FictMass = 3.0 </P>
</PARAMETERS>
\end{lstlisting}
\subsubsection{Running the solver}
In this example we will run the solver in parallel. We can specify the number of the strips by providing an additional flag to the solver, –nsz. In this example, we will run 16 strips, therefore it would be specified as –nsz 16. The solver can now be run as follows
\begin{lstlisting}[style=BashInputStyle]
mpirun -np 16 IncNavierStokesSolver CylFlow_HomoStrip.xml --npz 16 --nsz 16
\end{lstlisting}
The simulation results are illustrated in spanwise vorticity contours in Figure \ref{f:incns:finite-strip-modeling}. The wake response of the cable appears as standing wave patter in the earlier stage and then it transitions into traveling wave response, as shown in this figure.
\begin{figure}
\begin{center}
\includegraphics[width=7cm]{Figures/strip-16-time-100.png}
\includegraphics[width=7cm]{Figures/strip-16-time-600.png}
\caption{Spanwise vorticity contours in standing wave and traveling wave patterns predicted in finite strip modeling.}
\label{f:incns:finite-strip-modeling}
\end{center}
\end{figure}
\subsection{2D direct stability analysis of the channel flow}
......
This diff is collapsed.
......@@ -2,12 +2,4 @@
\input{meshconvert}
\input{fieldconvert}
\input{fldtovtk}
\input{fldtotecplot}
\input{xmltovtk}
\input{probefld}
\input{fieldconvert}
\ No newline at end of file
......@@ -298,7 +298,7 @@ the Incompressible Navier-Stokes solver,
Finally, multi-variable functions such as initial conditions and analytic
solutions may be specified for use in, or comparison with, simulations. These
may be specified using expressions (\inltt{<E>}) or imported from a file
(\inltt{<F>}) using the Nektar++ .fld or .pts file formats
(\inltt{<F>}) using the Nektar++ FLD file format
\begin{lstlisting}[style=XMLStyle]
<FUNCTION NAME="ExactSolution">
......@@ -309,31 +309,11 @@ may be specified using expressions (\inltt{<E>}) or imported from a file
</FUNCTION>
\end{lstlisting}
An .fld file is a solution file (in other words an .fld renamed as .rst) where
A restart file is a solution file (in other words an .fld renamed as .rst) where
the field data is specified. The expansion order used to generate the .rst file
must be the same as that for the simulation.
.pts files contain scattered point data which needs to be interpolated to the field.
For further information on the file format and the different interpolation schemes, see
section~\ref{s:utilities:fieldconvert:sub:interppointdatatofld}.
All filenames must be specified relative to the location of the .xml file.
must be the same as that for the simulation. The filename must be specified
relative to the location of the .xml file.
With the additional argument \inltt{TIMEDEPENDENT="1"}, different files can be
loaded for each timestep. The filenames are defined using
\href{http://www.boost.org/doc/libs/1_56_0/libs/format/doc/format.html#syntax}{boost::format syntax}
where the step time is used as variable. For example, the function
\inltt{Baseflow} would load the files \inltt{U0V0\_1.00000000E-05.fld},
\inltt{U0V0\_2.00000000E-05.fld} and so on.
\begin{lstlisting}[style=XMLStyle]
<FUNCTION NAME="Baseflow">
<F VAR="U0,V0" TIMEDEPENDENT="1"FILE="U0V0_%14.8E.fld"/>
</FUNCTION>
\end{lstlisting}
For .pts files, the time consuming computation of interpolation weights in only
performed for the first timestep. The weights are stored and reused in all subsequent steps,
which is why all consecutive .pts files must use the same ordering, number and location of
data points.
Other examples of this input feature can be the insertion of a forcing term,
\begin{lstlisting}[style=XMLStyle]
......
......@@ -875,6 +875,10 @@ namespace Nektar
}
}
/** \brief Get the normals along specficied face
* Get the face normals interplated to a points0 x points 0
* type distribution
**/
void PrismExp::v_ComputeFaceNormal(const int face)
{
const SpatialDomains::GeomFactorsSharedPtr &geomFactors =
......
......@@ -197,7 +197,7 @@ namespace Nektar
}
}
/**
/**
*
*/
ExpansionType ExpList::GetExpType(void)
......@@ -560,7 +560,7 @@ namespace Nektar
* array of size \f$N_{\mathrm{eof}}\f$.
*/
void ExpList::v_FwdTrans_IterPerExp(const Array<OneD, const NekDouble> &inarray,
Array<OneD, NekDouble> &outarray)
Array<OneD, NekDouble> &outarray)
{
Array<OneD,NekDouble> f(m_ncoeffs);
......@@ -1187,7 +1187,7 @@ namespace Nektar
* \f$Q_{\mathrm{tot}}\f$.
*/
void ExpList::v_BwdTrans_IterPerExp(const Array<OneD, const NekDouble> &inarray,
Array<OneD, NekDouble> &outarray)
Array<OneD, NekDouble> &outarray)
{
Array<OneD, NekDouble> tmp;
for (int i = 0; i < m_collections.size(); ++i)
......@@ -1235,7 +1235,7 @@ namespace Nektar
NekDouble tol,
bool returnNearestElmt)
{
NekDouble resid;
NekDouble nearpt = 1e6;
if (GetNumElmts() == 0)
{
......@@ -1255,7 +1255,7 @@ namespace Nektar
{
if ((*m_exp)[i]->GetGeom()->ContainsPoint(gloCoords,
locCoords,
tol, resid))
tol, nearpt))
{
w.SetX(gloCoords[0]);
w.SetY(gloCoords[1]);
......@@ -1291,7 +1291,7 @@ namespace Nektar
// retrieve local coordinate of point
(*m_exp)[min_id]->GetGeom()->GetLocCoords(gloCoords,
locCoords);
locCoords);
return min_id;
}
else
......@@ -1304,56 +1304,64 @@ namespace Nektar
{
static int start = 0;
int min_id = 0;
NekDouble resid_min = 1e6;
NekDouble nearpt_min = 1e6;
Array<OneD, NekDouble> savLocCoords(locCoords.num_elements());
// restart search from last found value
for (int i = start; i < (*m_exp).size(); ++i)
{
if ((*m_exp)[i]->GetGeom()->ContainsPoint(gloCoords, locCoords,
tol, resid))
if ((*m_exp)[i]->GetGeom()->ContainsPoint(gloCoords,
locCoords,
tol, nearpt))
{
start = i;
return i;
}
else
{
if(resid < resid_min)
if(nearpt < nearpt_min)
{
min_id = i;
resid_min = resid;
Vmath::Vcopy(locCoords.num_elements(),savLocCoords,1,locCoords,1);
nearpt_min = nearpt;
Vmath::Vcopy(locCoords.num_elements(),locCoords,1,savLocCoords,1);
}
}
}
for (int i = 0; i < start; ++i)
{
if ((*m_exp)[i]->GetGeom()->ContainsPoint(gloCoords, locCoords,
tol, resid))
if ((*m_exp)[i]->GetGeom()->ContainsPoint(gloCoords,
locCoords,
tol, nearpt))
{
start = i;
return i;
}
else
{
if(resid < resid_min)
if(nearpt < nearpt_min)
{
min_id = i;
resid_min = resid;
Vmath::Vcopy(locCoords.num_elements(),savLocCoords,1,locCoords,1);
nearpt_min = nearpt;
Vmath::Vcopy(locCoords.num_elements(),
locCoords,1,savLocCoords,1);
}
}
}
std::string msg = "Failed to find point in element to tolerance of "
+ boost::lexical_cast<std::string>(resid)
+ " using nearest point found";
WARNINGL0(true,msg.c_str());
std::string msg = "Failed to find point within element to tolerance of "
+ boost::lexical_cast<std::string>(tol)
+ " using local point ("
+ boost::lexical_cast<std::string>(locCoords[0]) +","
+ boost::lexical_cast<std::string>(locCoords[1]) +","
+ boost::lexical_cast<std::string>(locCoords[1])
+ ") in element: "
+ boost::lexical_cast<std::string>(min_id);
WARNINGL1(false,msg.c_str());
if(returnNearestElmt)
{
Vmath::Vcopy(locCoords.num_elements(),locCoords,1,savLocCoords,1);
Vmath::Vcopy(locCoords.num_elements(),savLocCoords,1,locCoords,1);
return min_id;
}
else
......@@ -1817,7 +1825,6 @@ namespace Nektar
ASSERTL0(false,
"This method is not defined or valid for this class type");
LibUtilities::TranspositionSharedPtr trans;
return trans;
}
......@@ -1834,7 +1841,6 @@ namespace Nektar
ASSERTL0(false,
"This method is not defined or valid for this class type");
Array<OneD, unsigned int> NoModes(1);
return NoModes;
}
......@@ -1843,7 +1849,6 @@ namespace Nektar
ASSERTL0(false,
"This method is not defined or valid for this class type");
Array<OneD, unsigned int> NoModes(1);
return NoModes;
}
......@@ -1936,6 +1941,7 @@ namespace Nektar
void ExpList::GeneralGetFieldDefinitions(std::vector<LibUtilities::FieldDefinitionsSharedPtr> &fielddef,
int NumHomoDir,
int NumHomoStrip,
Array<OneD, LibUtilities::BasisSharedPtr> &HomoBasis,
std::vector<NekDouble> &HomoLen,
std::vector<unsigned int> &HomoZIDs,
......@@ -2025,8 +2031,16 @@ namespace Nektar
if(elementIDs.size() > 0)
{
LibUtilities::FieldDefinitionsSharedPtr fdef = MemoryManager<LibUtilities::FieldDefinitions>::AllocateSharedPtr(shape, elementIDs, basis, UniOrder, numModes,fields, NumHomoDir, HomoLen, HomoZIDs, HomoYIDs);
fielddef.push_back(fdef);
for(int i = 0; i < NumHomoStrip; ++i)
{
LibUtilities::FieldDefinitionsSharedPtr fdef =
MemoryManager<LibUtilities::FieldDefinitions>::
AllocateSharedPtr(shape, elementIDs, basis,
UniOrder, numModes,fields,
NumHomoDir, HomoLen, HomoZIDs,
HomoYIDs);
fielddef.push_back(fdef);
}
}
}
}
......@@ -2395,7 +2409,6 @@ namespace Nektar
"This method is not defined or valid for this class type");
}
void ExpList::v_GetBCValues(Array<OneD, NekDouble> &BndVals,
const Array<OneD, NekDouble> &TotField,
int BndID)
......
......@@ -64,13 +64,13 @@ namespace Nektar
class GlobalMatrix;
enum Direction
{
eX,
eY,
eZ,
eS,
eN
};
{
eX,
eY,
eZ,
eS,
eN
};
enum ExpansionType
{
......@@ -81,21 +81,21 @@ namespace Nektar
e3DH2D,
e3D,
eNoType
};
};
MultiRegions::Direction const DirCartesianMap[] =
{
eX,
eY,
eZ
};
{
eX,
eY,
eZ
};
/// A map between global matrix keys and their associated block
/// matrices.
typedef map<GlobalMatrixKey,DNekScalBlkMatSharedPtr> BlockMatrixMap;
/// A shared pointer to a BlockMatrixMap.
typedef boost::shared_ptr<BlockMatrixMap> BlockMatrixMapShPtr;
/// Base class for all multi-elemental spectral/hp expansions.
class ExpList: public boost::enable_shared_from_this<ExpList>
......@@ -134,7 +134,7 @@ namespace Nektar
/// Returns the type of the expansion
MULTI_REGIONS_EXPORT void SetExpType(ExpansionType Type);
/// Evaulates the maximum number of modes in the elemental basis
/// order over all elements
inline int EvalBasisNumModesMax(void) const;
......@@ -160,12 +160,12 @@ namespace Nektar
/// Returns the total number of qudature points scaled by
/// the factor scale on each 1D direction
inline int Get1DScaledTotPoints(const NekDouble scale) const;
/// Sets the wave space to the one of the possible configuration
/// true or false
inline void SetWaveSpace(const bool wavespace);
///Set Modified Basis for the stability analysis
inline void SetModifiedBasis(const bool modbasis);
......@@ -321,8 +321,8 @@ namespace Nektar
Array<OneD, NekDouble> &coord_0,
Array<OneD, NekDouble> &coord_1 = NullNekDouble1DArray,
Array<OneD, NekDouble> &coord_2 = NullNekDouble1DArray);
// Homogeneous transforms
// Homogeneous transforms
inline void HomogeneousFwdTrans(
const Array<OneD, const NekDouble> &inarray,
Array<OneD, NekDouble> &outarray,
......@@ -342,7 +342,7 @@ namespace Nektar
const Array<OneD, NekDouble> &inarray2,
Array<OneD, NekDouble> &outarray,
CoeffState coeffstate = eLocal);
inline void GetBCValues(
Array<OneD, NekDouble> &BndVals,
const Array<OneD, NekDouble> &TotField,
......@@ -353,7 +353,7 @@ namespace Nektar
Array<OneD, const NekDouble> &V2,
Array<OneD, NekDouble> &outarray,
int BndID);
/// Apply geometry information to each expansion.
MULTI_REGIONS_EXPORT void ApplyGeomInfo();
......@@ -397,7 +397,7 @@ namespace Nektar
}
void WriteVtkPieceHeader(std::ofstream &outfile, int expansion,
int istrip)
int istrip = 0)
{
v_WriteVtkPieceHeader(outfile, expansion, istrip);
}
......@@ -627,8 +627,8 @@ namespace Nektar
inline void PhysDeriv(
Direction edir,
const Array<OneD, const NekDouble> &inarray,
Array<OneD, NekDouble> &out_d);
Array<OneD, NekDouble> &out_d);
/// This function discretely evaluates the derivative of a function
/// \f$f(\boldsymbol{x})\f$ on the domain consisting of all
/// elements of the expansion.
......@@ -738,6 +738,7 @@ namespace Nektar
MULTI_REGIONS_EXPORT void GeneralGetFieldDefinitions(
std::vector<LibUtilities::FieldDefinitionsSharedPtr> &fielddef,
int NumHomoDir = 0,
int NumHomoStrip = 1,
Array<OneD, LibUtilities::BasisSharedPtr> &HomoBasis =
LibUtilities::NullBasisSharedPtr1DArray,
std::vector<NekDouble> &HomoLen =
......@@ -821,15 +822,15 @@ namespace Nektar
const boost::shared_ptr<ExpList> &fromExpList,
const Array<OneD, const NekDouble> &fromCoeffs,
Array<OneD, NekDouble> &toCoeffs);
//Extract data in fielddata into the m_coeffs_list for the 3D stability analysis (base flow is 2D)
MULTI_REGIONS_EXPORT void ExtractDataToCoeffs(
LibUtilities::FieldDefinitionsSharedPtr &fielddef,
std::vector<NekDouble> &fielddata,
std::string &field,
Array<OneD, NekDouble> &coeffs);
/// Returns a shared pointer to the current object.
boost::shared_ptr<ExpList> GetSharedThisPtr()
......@@ -975,7 +976,7 @@ namespace Nektar
NekOptimize::GlobalOptParamSharedPtr m_globalOptParam;
BlockMatrixMapShPtr m_blockMat;
//@todo should this be in ExpList or ExpListHomogeneous1D.cpp
// it's a bool which determine if the expansion is in the wave space (coefficient space)
// or not
......@@ -1050,7 +1051,7 @@ namespace Nektar
Array<OneD, NekDouble> &Upwind);
virtual boost::shared_ptr<ExpList> &v_GetTrace();
virtual boost::shared_ptr<AssemblyMapDG> &v_GetTraceMap();
virtual const Array<OneD, const int> &v_GetTraceBndMap();
......@@ -1134,11 +1135,11 @@ namespace Nektar
const Array<OneD,const NekDouble> &inarray,
Array<OneD, NekDouble> &outarray,
CoeffState coeffstate);
virtual void v_BwdTrans_IterPerExp(
const Array<OneD,const NekDouble> &inarray,
Array<OneD,NekDouble> &outarray);
virtual void v_FwdTrans(
const Array<OneD,const NekDouble> &inarray,
Array<OneD, NekDouble> &outarray,
......@@ -1154,11 +1155,11 @@ namespace Nektar
const Array<OneD, const NekDouble> &inarray,
Array<OneD, NekDouble> &outarray,
CoeffState coeffstate);
virtual void v_IProductWRTBase_IterPerExp(
const Array<OneD,const NekDouble> &inarray,
Array<OneD, NekDouble> &outarray);
virtual void v_GeneralMatrixOp(
const GlobalMatrixKey &gkey,
const Array<OneD,const NekDouble> &inarray,
......@@ -1169,7 +1170,7 @@ namespace Nektar
Array<OneD, NekDouble> &coord_0,
Array<OneD, NekDouble> &coord_1,
Array<OneD, NekDouble> &coord_2 = NullNekDouble1DArray);
virtual void v_PhysDeriv(
const Array<OneD, const NekDouble> &inarray,
Array<OneD, NekDouble> &out_d0,
......@@ -1451,7 +1452,7 @@ namespace Nektar
{
m_WaveSpace = wavespace;
}
/**
*
*/
......@@ -1520,7 +1521,7 @@ namespace Nektar
v_IProductWRTBase(inarray,outarray, coeffstate);
}
/**
/**
*
*/
inline void ExpList::IProductWRTBase_IterPerExp(
......@@ -1540,8 +1541,8 @@ namespace Nektar
{
v_FwdTrans(inarray,outarray,coeffstate);
}
/**
/**
*
*/
inline void ExpList::FwdTrans_IterPerExp (
......@@ -1569,8 +1570,8 @@ namespace Nektar
{
v_BwdTrans(inarray,outarray,coeffstate);
}
/**
/**
*
*/
inline void ExpList::BwdTrans_IterPerExp (
......@@ -1644,7 +1645,7 @@ namespace Nektar
{
v_GetCoords(coord_0,coord_1,coord_2);
}
/**
*
*/
......@@ -1655,7 +1656,7 @@ namespace Nektar
{
v_PhysDeriv(inarray,out_d0,out_d1,out_d2);
}
/**
*
*/
......@@ -1673,8 +1674,8 @@ namespace Nektar
Array<OneD, NekDouble> &out_d)
{
v_PhysDeriv(edir, inarray,out_d);
}