Commit 2e5a08dd authored by Michael Turner's avatar Michael Turner
Browse files

update to master

parents f08a7d6c 2db0ae69
......@@ -39,6 +39,7 @@ v4.4.0
- Fix bug in the calculation of the RHS magnitude in CG solver (!721)
- Fix bug in CMake Homebrew and MacPorts detection for OS X (!729)
- Fix bug in FieldUtils when using half mode expansions (!734)
- Fix bug in CMake PETSc detection for Ubuntu 16.04/Debian 9 (!735)
**ADRSolver:**
- Add a projection equation system for C^0 projections (!675)
......@@ -55,6 +56,7 @@ v4.4.0
- Fix linearised advection for full 3D cases (!708)
- Added a weak pressure formulation following Guermond & Shen (!713)
- Added a convective like outflow boundary condition from Dong (!713)
- Added the ability to specifiy Womersley boundary conditions for pulsatile flow (!472)
**FieldConvert:**
- Allow equi-spaced output for 1D and 2DH1D fields (!613)
......@@ -92,6 +94,12 @@ v4.4.0
- Change variable names in mcf file to make more sense (!736)
- Fix issues in varopti module so that in can be compiled without meshgen on
(!736)
- Replace LAPACK Eigenvalue calculation with handwritten function in
varopti (!738)
- Improved node-colouring algorithm for better load-balancing
in varopti (!738)
- Simplified calculation of the energy functional in varopti for improved
performance (!738)
**FieldConvert:**
- Move all modules to a new library, FieldUtils, to support post-processing
......@@ -100,6 +108,9 @@ v4.4.0
- Add module to add composite ID of elements as a field (!674)
- Add reader for Nek5000 field files (!680)
**Tester:**
- Fix output not displayed on segfault or system error (!745)
v4.3.5
------
**Library:**
......@@ -128,6 +139,12 @@ v4.3.4
**IncNavierStokesSolver:**
- Fix 2nd order time-integration for VCSMapping (!687)
v4.3.4
------
**Library:**
- Fix performance issue with `v_ExtractDataToCoeffs` for post-processing of large
simulations (!672)
v4.3.3
------
**Library**:
......
......@@ -139,6 +139,8 @@ OPTION(NEKTAR_BUILD_PACKAGES "Build Nektar++ binary packages" OFF)
MARK_AS_ADVANCED(NEKTAR_BUILD_PACKAGES)
OPTION(NEKTAR_TEST_ALL "Include full set of regression tests to this build." OFF)
OPTION(NEKTAR_TEST_USE_HOSTFILE "Use a hostfile to explicitly specify number of
slots." OFF)
# Meshing utilities and library
IF (NOT WIN32)
......@@ -231,7 +233,6 @@ INCLUDE (ThirdPartyArpack)
INCLUDE (ThirdPartyMPI)
INCLUDE (ThirdPartyPETSc)
INCLUDE (ThirdPartyVTK)
INCLUDE (ThirdPartyQT4)
INCLUDE (ThirdPartySMV)
INCLUDE (ThirdPartyOCE)
INCLUDE (ThirdPartyTetGen)
......
......@@ -225,6 +225,13 @@ show :
else ()
set (PETSC_LIBRARY_VEC "NOTFOUND" CACHE INTERNAL "Cleared" FORCE) # There is no libpetscvec
petsc_find_library (SINGLE petsc)
# Debian 9/Ubuntu 16.04 uses _real and _complex extensions when using libraries in /usr/lib/petsc.
if (NOT PETSC_LIBRARY_SINGLE)
petsc_find_library (SINGLE petsc_real)
endif()
if (NOT PETSC_LIBRARY_SINGLE)
petsc_find_library (SINGLE petsc_complex)
endif()
foreach (pkg SYS VEC MAT DM KSP SNES TS ALL)
set (PETSC_LIBRARIES_${pkg} "${PETSC_LIBRARY_SINGLE}")
endforeach ()
......
......@@ -8,7 +8,7 @@
# Try to find system Loki headers. Hint /opt/local/include for MacPorts
# (although there is no Portfile for Loki currently).
FIND_PATH(LOKI_INCLUDE_DIR loki/Typelist.h PATHS /opt/local/include)
FIND_PATH(LOKI_INCLUDE_DIR loki/Singleton.h PATHS /opt/local/include)
IF (LOKI_INCLUDE_DIR)
SET(BUILD_LOKI OFF)
......
# QT4
OPTION(NEKTAR_USE_QT4
"Use QT4 for graphical tools and utilities." OFF)
IF( NEKTAR_USE_QT4 )
FIND_PACKAGE(Qt4 COMPONENTS QtCore QtGui REQUIRED)
INCLUDE(${QT_USE_FILE})
ADD_DEFINITIONS(${QT_DEFINITIONS})
IF (QT4_FOUND)
MESSAGE(STATUS "Found QT4: ${QT_INCLUDE_DIR}")
ENDIF(QT4_FOUND)
ENDIF( NEKTAR_USE_QT4 )
......@@ -872,6 +872,83 @@ the session file:
\item \inltt{SVVDiffCoeff}: sets the SVV diffusion coefficient (default value = 0.1)
\end{itemize}
\subsection{Womersley Boundary Condition}
It is possible to define the time-dependent Womersley velocity profile
for pulsatile flow in a pipe. The modulation of the velocity profile
is based on the desired peak or centerline velocity which can be
represented by a Fourier expansion $U_{max}=A(\omega_n)e^{i\omega_n
t}$ where $A$ are the Fourier modes and $\omega $ the frequency. The
womersely solution is then defined as:
$$ w(r,t) = A_0(1-(r/R)^2) + \sum_{n=1}^N
\tilde{A_n}[1-\frac{J_0(i^{3/2}\alpha_n r/R)}{J_0(i^{3/2}
\alpha)}]e^{i\omega_n t} $$
where the womersley number $\alpha$ is defined:
$$ \alpha_n = R\sqrt{\frac{2\pi n}{T\nu}}$$
and $\tilde{A_n}$ ($n=1:N$)are the Fourier coefficients scaled in the
following way:
$$ \tilde{A_n} = 2A_n/[1 - \frac{1}{J_0(i^{3/2}\alpha)}] $$
The Womersley velocity profile is implemented in the following way:
\begin{lstlisting}[style=XMLStyle]
<REGION REF="0">
<D VAR="u" USERDEFINEDTYPE="Womersley:WomParams.xml" VALUE="0" />
<D VAR="v" USERDEFINEDTYPE="Womersley:WomParams.xml" VALUE="0" />
<D VAR="w" USERDEFINEDTYPE="Womersley:WomParams.xml" VALUE="0" />
<N VAR="p" USERDEFINEDTYPE="H" VALUE="0" />
</REGION>
\end{lstlisting}
A file containing the Fourier coefficients, $\tilde{A}$, must be in
the directory where the solver is called from. The name of the file is
defined by the string given in the attribute \inltt{USERDEFINEDTYPE}
after the ``:'' and contains the real and imaginary coefficients. This
file has the format
\begin{lstlisting}[style=XMLStyle]
<NEKTAR>
<WOMERSLEYBC>
<WOMPARAMS>
<W PROPERTY="Radius" VALUE="0.5" />
<W PROPERTY="Period" VALUE="1.0" />
<W PROPERTY="axisnormal" VALUE="0.0,0.0,1.0" />
<W PROPERTY="axispoint" VALUE="0.0,0.0,0.0" />
</WOMPARAMS>
<FOURIERCOEFFS>
<F ID="0"> 0.600393641193, 0.0 </F>
<F ID="1"> -0.277707172935, 0.0767582715413 </F>
<F ID="2"> -0.0229953131146, 0.0760936232478 </F>
<F ID="3"> 0.00858135174058, 0.017089888642 </F>
<F ID="4"> 0.0140332527651, 0.0171575122496 </F>
<F ID="5"> 0.0156970122129, -0.00547357750345 </F>
<F ID="6"> 0.00473626554238, -0.00498786519876 </F>
<F ID="7"> 0.00204434981523, -0.00614566561937 </F>
<F ID="8"> -0.000274697215201, 0.000153571881197 </F>
<F ID="9"> -0.000148037910774, 2.68919619581e-05 </F>
</FOURIERCOEFFS>
</WOMERSLEYBC>
</NEKTAR>
\end{lstlisting}
Each value of $\tilde{A}$ is provided in the \inltt{FOURIERCOEFFS}
section and provided as separate entries containing the real and
imaginary components, i.e. the mean component provided above is
$0.600393641193,0.0$.
Similarly in the \inltt{WOMPARAMS} section the key parameters of the boundary condition are also provided as:
\begin{itemize}
\item \inltt{RADIUS} is the radius of the boundary.
\item \inltt{PERIOD} is the cycle time period,
\item \inltt{AXISNORMAL} defines the normal direction to the boundary,
\item \inltt{AXISPOINT} defines a coordinate in the boundary centre,
\end{itemize}
\subsection{Forcing}
\subsubsection{MovingBody}
......
......@@ -495,7 +495,7 @@ where the step time is used as variable. For example, the function
</FUNCTION>
\end{lstlisting}
For .pts files, the time consuming computation of interpolation weights in only
For .pts files, the time consuming computation of interpolation weights is 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.
......@@ -674,4 +674,4 @@ will be the y-axis and the z-axis.
%%% Local Variables:
%%% mode: latex
%%% TeX-master: "../user-guide"
%%% End:
\ No newline at end of file
%%% End:
......@@ -247,6 +247,7 @@ namespace Nektar
m_functionMapNameToInstanceType["atan"] = E_ATAN;
m_functionMapNameToInstanceType["atan2"] = E_ATAN2;
m_functionMapNameToInstanceType["ang"] = E_ANG;
m_functionMapNameToInstanceType["bessel"] = E_BESSEL;
m_functionMapNameToInstanceType["ceil"] = E_CEIL;
m_functionMapNameToInstanceType["cos"] = E_COS;
m_functionMapNameToInstanceType["cosh"] = E_COSH;
......@@ -285,6 +286,8 @@ namespace Nektar
m_function2[E_ATAN2] = atan2;
m_function2[E_ANG] = ang;
m_function2[E_RAD] = rad;
m_function2[E_BESSEL] = boost::math::cyl_bessel_j;
// there is no entry to m_function that correspond to awgn function.
// this is made in purpose. This function need not be pre-evaluated once!
}
......@@ -780,6 +783,9 @@ namespace Nektar
case E_ANG:
stack.push_back ( makeStep<EvalAng>( stateIndex, stateIndex, stateIndex+1 ) );
return std::make_pair(false,0);
case E_BESSEL:
stack.push_back ( makeStep<EvalBessel>( stateIndex, stateIndex, stateIndex+1 ) );
return std::make_pair(false,0);
case E_CEIL:
stack.push_back ( makeStep<EvalCeil>( stateIndex, stateIndex ) );
return std::make_pair(false,0);
......
......@@ -45,6 +45,8 @@
#include <boost/random/variate_generator.hpp> // for variate_generator
#include <boost/random/normal_distribution.hpp>
#include <boost/math/special_functions/bessel.hpp>
#define BOOST_SPIRIT_THREADSAFE
#if( BOOST_VERSION / 100 % 1000 >= 36 )
#include <boost/spirit/include/classic_core.hpp>
......@@ -457,7 +459,6 @@ namespace Nektar
std::map<int, OneArgFunc> m_function;
std::map<int, TwoArgFunc> m_function2;
// ======================================================
// Internal representation of evaluation step
// ======================================================
......@@ -480,7 +481,7 @@ namespace Nektar
E_ABS, E_ASIN, E_ACOS, E_ATAN, E_ATAN2, E_ANG,
E_CEIL, E_COS, E_COSH, E_EXP, E_FABS, E_FLOOR,
E_LOG, E_LOG10, E_POW, E_RAD, E_SIN, E_SINH,
E_SQRT, E_TAN, E_TANH, E_SIGN, E_AWGN
E_SQRT, E_TAN, E_TANH, E_SIGN, E_AWGN, E_BESSEL
};
......@@ -643,6 +644,12 @@ namespace Nektar
virtual void run_many(ci n) { for(int i=0;i<n;i++) state[storeIdx*n+i] = ang( state[argIdx1*n+i], state[argIdx2*n+i] ); }
virtual void run_once() { state[storeIdx] = ang( state[argIdx1], state[argIdx2] ); }
};
struct EvalBessel: public EvaluationStep
{
EvalBessel(rgt rn, vr s, cvr c, cvr p, cvr v, ci i, ci l, ci r): EvaluationStep(rn,i,l,r,s,c,p,v) {}
virtual void run_many(ci n) { for(int i=0;i<n;i++) state[storeIdx*n+i] = boost::math::cyl_bessel_j( state[argIdx1*n+i], state[argIdx2*n+i] ); }
virtual void run_once() { state[storeIdx] = boost::math::cyl_bessel_j( state[argIdx1], state[argIdx2] ); }
};
struct EvalCeil: public EvaluationStep
{
EvalCeil(rgt rn, vr s, cvr c, cvr p, cvr v, ci i, ci l, ci r): EvaluationStep(rn,i,l,r,s,c,p,v) {}
......
......@@ -9,6 +9,11 @@
#include <float.h>
#include <complex>
#include <LibUtilities/BasicConst/NektarUnivTypeDefs.hpp>
/// Maximum number of iterations in polynomial defalation routine Jacobz
#define STOP 30
......@@ -2977,10 +2982,45 @@ namespace Polylib {
}
}
/**
\brief
Calcualte the bessel function of the first kind with complex double input y.
Taken from Numerical Recipies in C
Returns a complex double
*/
std::complex<Nektar::NekDouble> ImagBesselComp(int n,std::complex<Nektar::NekDouble> y)
{
std::complex<Nektar::NekDouble> z (1.0,0.0);
std::complex<Nektar::NekDouble> zbes (1.0,0.0);
std::complex<Nektar::NekDouble> zarg;
Nektar::NekDouble tol = 1e-15;
int maxit = 10000;
int i = 1;
zarg = -0.25*y*y;
while (abs(z) > tol && i <= maxit){
z = z*(1.0/i/(i+n)*zarg);
if (abs(z) <= tol) break;
zbes = zbes + z;
i++;
}
zarg = 0.5*y;
for (i=1;i<=n;i++){
zbes = zbes*zarg;
}
return zbes;
}
} // end of namespace
......
#ifndef H_PLYLIB
#include <LibUtilities/BasicConst/NektarUnivTypeDefs.hpp>
#include <LibUtilities/LibUtilitiesDeclspec.h>
#include <complex>
/*
* LIBRARY ROUTINES FOR POLYNOMIAL CALCULUS AND INTERPOLATION
......@@ -158,7 +160,9 @@ namespace Polylib {
LIB_UTILITIES_EXPORT void jacobd (const int, const double *, double *, const int ,
const double, const double);
LIB_UTILITIES_EXPORT std::complex<Nektar::NekDouble> ImagBesselComp(const int,
std::complex<Nektar::NekDouble> );
} // end of namespace
......
......@@ -106,6 +106,7 @@ IF( NEKTAR_SOLVER_INCNAVIERSTOKES )
ADD_NEKTAR_TEST(ChannelSpongeLNSE)
ADD_NEKTAR_TEST(ChanFlow_Standard_BodyForce)
ADD_NEKTAR_TEST(Cyl_AdaptiveSFD)
ADD_NEKTAR_TEST(Womersley_PipeFlow)
ADD_NEKTAR_TEST(CylFlow_Mov_mapping)
ADD_NEKTAR_TEST(PPF_R10000_ModifiedArnoldi_Shift)
ADD_NEKTAR_TEST(PPF_R15000_ModifiedArnoldi_Shift)
......
......@@ -44,6 +44,18 @@
#include <LocalRegions/Expansion2D.h>
#include <LocalRegions/Expansion3D.h>
#include <LibUtilities/Polylib/Polylib.h>
#include <LibUtilities/BasicUtils/FileSystem.h>
#include <LibUtilities/BasicUtils/PtsIO.h>
#include <algorithm>
#include <complex>
#include <iostream>
#include <fstream>
#include <sstream>
#include <tinyxml.h>
#include <LibUtilities/BasicUtils/ParseUtils.hpp>
using namespace std;
namespace Nektar
......@@ -234,6 +246,48 @@ namespace Nektar
}
}
// Set up maping for womersley BC - and load variables
for (int i = 0; i < m_fields.num_elements(); ++i)
{
for(int n = 0; n < m_fields[i]->GetBndConditions().num_elements(); ++n)
{
if(boost::istarts_with(m_fields[i]->GetBndConditions()[n]->GetUserDefined(),"Womersley"))
{
m_womersleyParams[n] = MemoryManager<WomersleyParams>::AllocateSharedPtr(m_spacedim);
#if 0
m_session->LoadParameter("Period",m_womersleyParams[n]->m_period);
m_session->LoadParameter("Radius",m_womersleyParams[n]->m_radius);
NekDouble n0,n1,n2;
m_session->LoadParameter("n0",n0);
m_session->LoadParameter("n1",n1);
m_session->LoadParameter("n2",n2);
m_womersleyParams[n]->m_axisnormal[0] = n0;
m_womersleyParams[n]->m_axisnormal[1] = n1;
m_womersleyParams[n]->m_axisnormal[2] = n2;
NekDouble x0,x1,x2;
m_session->LoadParameter("x0",x0);
m_session->LoadParameter("x1",x1);
m_session->LoadParameter("x2",x2);
m_womersleyParams[n]->m_axispoint[0] = x0;
m_womersleyParams[n]->m_axispoint[1] = x1;
m_womersleyParams[n]->m_axispoint[2] = x2;
#endif
// Read in fourier coeffs
SetUpWomersley(n,
m_fields[i]->GetBndConditions()[n]->GetUserDefined());
m_fields[i]->GetBoundaryToElmtMap(m_fieldsBCToElmtID[i],m_fieldsBCToTraceID[i]);
}
}
}
// Set up Field Meta Data for output files
m_fieldMetaDataMap["Kinvis"] = boost::lexical_cast<std::string>(m_kinvis);
m_fieldMetaDataMap["TimeStep"] = boost::lexical_cast<std::string>(m_timestep);
......@@ -345,7 +399,10 @@ namespace Nektar
varName = m_session->GetVariable(i);
m_fields[i]->EvaluateBoundaryConditions(time, varName);
}
else if(boost::istarts_with(m_fields[i]->GetBndConditions()[n]->GetUserDefined(),"Womersley"))
{
SetWomersleyBoundary(i,n);
}
}
// Set Radiation conditions if required
......@@ -449,7 +506,7 @@ namespace Nektar
Array<OneD, Array<OneD, NekDouble> > normals;
Array<OneD, NekDouble> Bphys,Bcoeffs;
int fieldid = m_velocity[0];
int fldid = m_velocity[0];
for(cnt = n = 0; n < BndConds[0].num_elements(); ++n)
{
......@@ -457,9 +514,9 @@ namespace Nektar
{
for(i = 0; i < BndExp[0][n]->GetExpSize(); ++i,cnt++)
{
elmtid = m_fieldsBCToElmtID[fieldid][cnt];
elmtid = m_fieldsBCToElmtID[fldid][cnt];
elmt = m_fields[0]->GetExp(elmtid);
boundary = m_fieldsBCToTraceID[fieldid][cnt];
boundary = m_fieldsBCToTraceID[fldid][cnt];
normals = elmt->GetSurfaceNormal(boundary);
......@@ -500,8 +557,253 @@ namespace Nektar
/**
* Add an additional forcing term programmatically.
* Womersley boundary condition defintion
*/
void IncNavierStokes::SetWomersleyBoundary(const int fldid, const int bndid)
{
ASSERTL1(m_womersleyParams.count(bndid) == 1, "Womersley parameters for this boundary have not been set up");
WomersleyParamsSharedPtr WomParam = m_womersleyParams[bndid];
std::complex<NekDouble> za, zar, zJ0, zJ0r, zq, zvel, zJ0rJ0;
int i,j,k;
int M = WomParam->m_wom_vel_r.size();
NekDouble R = WomParam->m_radius;
NekDouble T = WomParam->m_period;
Array<OneD, NekDouble > normals = WomParam->m_axisnormal;
Array<OneD, NekDouble > x0 = WomParam->m_axispoint;
// Womersley Number
NekDouble alpha = R*sqrt(2*M_PI/T/m_kinvis);
NekDouble r,kt;
std::complex<NekDouble> z1 (1.0,0.0);
std::complex<NekDouble> zi (0.0,1.0);
std::complex<NekDouble> comp_conj (-1.0,1.0); //complex conjugate
Array<OneD, MultiRegions::ExpListSharedPtr> BndExp;
BndExp = m_fields[fldid]->GetBndCondExpansions();
StdRegions::StdExpansionSharedPtr elmt;
StdRegions::StdExpansionSharedPtr bc;
int cnt=0;
int elmtid,offset, boundary,nfq;
Array<OneD, NekDouble> Bvals,w;
//Loop over all expansions
for(i = 0; i < BndExp[bndid]->GetExpSize(); ++i,cnt++)
{
// Get element id and offset
elmtid = m_fieldsBCToElmtID[fldid][cnt];
elmt = m_fields[fldid]->GetExp(elmtid);
offset = m_fields[fldid]->GetPhys_Offset(elmtid);
// Get Boundary and trace expansion
bc = BndExp[bndid]->GetExp(i);
boundary = m_fieldsBCToTraceID[fldid][cnt];
nfq=bc->GetTotPoints();
w = m_fields[fldid]->UpdatePhys() + offset;
Array<OneD, NekDouble> x(nfq,0.0);
Array<OneD, NekDouble> y(nfq,0.0);
Array<OneD, NekDouble> z(nfq,0.0);
Array<OneD, NekDouble> wbc(nfq,0.0);
bc->GetCoords(x,y,z);
// Add edge values (trace) into the wbc
elmt->GetTracePhysVals(boundary,bc,w,wbc);
//Compute womersley solution
for (j=0;j<nfq;j++)
{
//NOTE: only need to calculate these two once, could
//be stored or precomputed?
r = sqrt((x[j]-x0[0])*(x[j]-x0[0]) +
(y[j]-x0[1])*(y[j]-x0[1]) +
(z[j]-x0[2])*(z[j]-x0[2]))/R;
wbc[j] = WomParam->m_wom_vel_r[0]*(1. - r*r); // Compute Poiseulle Flow
for (k=1; k<M; k++)
{
kt = 2.0 * M_PI * k * m_time / T;
za = alpha * sqrt((NekDouble)k/2.0) * comp_conj;
zar = r * za;
zJ0 = Polylib::ImagBesselComp(0,za);
zJ0r = Polylib::ImagBesselComp(0,zar);
zJ0rJ0 = zJ0r / zJ0;
zq = std::exp(zi * kt) * std::complex<NekDouble>(
WomParam->m_wom_vel_r[k],
WomParam->m_wom_vel_i[k]);
zvel = zq * (z1 - zJ0rJ0);
wbc[j] = wbc[j] + zvel.real();
}
}
// Multiply w by normal to get u,v,w component of velocity
Vmath::Smul(nfq,normals[fldid],wbc,1,wbc,1);
Bvals = BndExp[bndid]->UpdateCoeffs()+
BndExp[bndid]->GetCoeff_Offset(i);
// Push back to Coeff space
bc->FwdTrans(wbc,Bvals);
}
}
void IncNavierStokes::SetUpWomersley(const int bndid, std::string womStr)
{
std::string::size_type indxBeg = womStr.find_first_of(':') + 1;
string filename = womStr.substr(indxBeg,string::npos);
std::complex<NekDouble> coef;
#if 1
TiXmlDocument doc(filename);
bool loadOkay = doc.LoadFile();
ASSERTL0(loadOkay,(std::string("Failed to load file: ") +
filename).c_str());
TiXmlHandle docHandle(&doc);
int err; /// Error value returned by TinyXML.
TiXmlElement *nektar = doc.FirstChildElement("NEKTAR");
ASSERTL0(nektar, "Unable to find NEKTAR tag in file.");
TiXmlElement *wombc = nektar->FirstChildElement("WOMERSLEYBC");
ASSERTL0(wombc, "Unable to find WOMERSLEYBC tag in file.");
// read womersley parameters
TiXmlElement *womparam = wombc->FirstChildElement("WOMPARAMS");
ASSERTL0(womparam, "Unable to find WOMPARAMS tag in file.");
// Input coefficients