Skip to content

GitLab

Commit ec3ebd52 authored by Chris Cantwell's avatar Chris Cantwell
Browse files

Merge branch 'master' into fix/PetSc_MKL

parents 43aa3700 db898d38
Hide whitespace changes
Inline Side-by-side
CHANGELOG.md
View file @ ec3ebd52
......@@ -34,9 +34,12 @@ v4.4.0
- Enabled MUMPS support in PETSc if a Fortran compiler was found and added 3D
support to the Helmholtz smoother used e.g. in FieldConverts C0Projection
module (!714)
- Fix bug in `Vmath::FillWhiteNoise` which caused `ForcingNoise` to have
a repeated pattern (!718)
- 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)
......@@ -53,12 +56,17 @@ 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)
**CardiacEPSolver:**
- Added a Python translator utility to generate cell models from CellML (!723)
**FieldConvert:**
- Allow equi-spaced output for 1D and 2DH1D fields (!613)
- Update quality metric to include scaled Jacobian output (!695)
- Allow multiple XML files to be specified in InterpField module (!705)
- Fix issues with isocontour module (!719)
- Fix issue with interpolator routine (!746)
**NekMesh:**
- Modify curve module to allow for spline input (!628)
......@@ -85,6 +93,17 @@ v4.4.0
(!712)
- 2D to 3D mesh extrusion module (!715)
- Add new two-dimensional mesher from NACA code or step file (!720)
- Fix inverted boundary layer in 2D (!736)
- More sensible element sizing with boundary layers in 2D (!736)
- 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
......@@ -93,6 +112,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:**
......@@ -121,6 +143,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**:
......
CMakeLists.txt
View file @ ec3ebd52
......@@ -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)
......@@ -230,7 +232,6 @@ INCLUDE (ThirdPartyFFTW)
INCLUDE (ThirdPartyArpack)
INCLUDE (ThirdPartyMPI)
INCLUDE (ThirdPartyVTK)
INCLUDE (ThirdPartyQT4)
INCLUDE (ThirdPartySMV)
INCLUDE (ThirdPartyOCE)
INCLUDE (ThirdPartyTetGen)
......
cmake/FindPETSc.cmake
View file @ ec3ebd52
......@@ -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 ()
......
cmake/ThirdPartyLoki.cmake
View file @ ec3ebd52
......@@ -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)
......
cmake/ThirdPartyQT4.cmake deleted 100644 → 0
View file @ 43aa3700
# 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 )
docs/user-guide/solvers/incompressible-ns.tex
View file @ ec3ebd52
......@@ -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}
......
docs/user-guide/xml/xml-conditions.tex
View file @ ec3ebd52
......@@ -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:
library/FieldUtils/Interpolator.cpp
View file @ ec3ebd52
......@@ -339,8 +339,7 @@ void Interpolator::Interpolate(
if (elmtid >= 0)
{
int offset = m_expInField[0]->GetPhys_Offset(
m_expInField[0]->GetOffset_Elmt_Id(elmtid));
int offset = m_expInField[0]->GetPhys_Offset(elmtid);
for (int f = 0; f < m_expInField.size(); ++f)
{
......@@ -424,8 +423,7 @@ void Interpolator::Interpolate(
if (elmtid >= 0)
{
int offset = m_expInField[0]->GetPhys_Offset(
m_expInField[0]->GetOffset_Elmt_Id(elmtid));
int offset = m_expInField[0]->GetPhys_Offset(elmtid);
for (int f = 0; f < m_expInField.size(); ++f)
{
......
library/LibUtilities/BasicUtils/Thread.h
View file @ ec3ebd52
......@@ -100,19 +100,19 @@ class ThreadJob
{
public:
/// Base constructor
ThreadJob();
LIB_UTILITIES_EXPORT ThreadJob();
/// Base destructor.
virtual ~ThreadJob();
LIB_UTILITIES_EXPORT virtual ~ThreadJob();
/**
* This method will be called when the task is loaded
* onto a worker thread and is ready to run. When Run
* has finished this instance will be destructed.
*/
virtual void Run() = 0;
LIB_UTILITIES_EXPORT virtual void Run() = 0;
/// Set number of worker threads.
void SetWorkerNum(unsigned int num);
LIB_UTILITIES_EXPORT void SetWorkerNum(unsigned int num);
protected:
/**
......@@ -168,7 +168,7 @@ class ThreadManager : public boost::enable_shared_from_this<ThreadManager>
{
public:
/// Destructor.
virtual ~ThreadManager();
LIB_UTILITIES_EXPORT virtual ~ThreadManager();
/**
* @brief Pass a list of tasklets to the master queue.
* @param joblist Vector of ThreadJob pointers.
......@@ -181,7 +181,7 @@ class ThreadManager : public boost::enable_shared_from_this<ThreadManager>
*
* @see SchedType
*/
virtual void QueueJobs(std::vector<ThreadJob*>& joblist) = 0;
LIB_UTILITIES_EXPORT virtual void QueueJobs(std::vector<ThreadJob*>& joblist) = 0;
/**
* @brief Pass a single job to the master queue.
* @param job A pointer to a ThreadJob subclass.
......@@ -190,14 +190,14 @@ class ThreadManager : public boost::enable_shared_from_this<ThreadManager>
* issue then suspend the workers with SetNumWorkers(0) until the jobs
* are queued.
*/
virtual void QueueJob(ThreadJob* job) = 0;
LIB_UTILITIES_EXPORT virtual void QueueJob(ThreadJob* job) = 0;
/**
* @brief Return the number of active workers.
*
* Active workers are threads that are either running jobs
* or are waiting for jobs to be queued.
*/
virtual unsigned int GetNumWorkers() = 0;
LIB_UTILITIES_EXPORT virtual unsigned int GetNumWorkers() = 0;
/**
* @brief Returns the worker number of the executing thread.
*
......@@ -216,7 +216,7 @@ class ThreadManager : public boost::enable_shared_from_this<ThreadManager>
*
* Returns 0 if called by non-thread.
*/
virtual unsigned int GetWorkerNum() = 0;
LIB_UTILITIES_EXPORT virtual unsigned int GetWorkerNum() = 0;
/**
* @brief Sets the number of active workers.
* @param num The number of active workers.
......@@ -227,18 +227,18 @@ class ThreadManager : public boost::enable_shared_from_this<ThreadManager>
* If num is greater than the maximum allowed number of active workers,
* then the maximum value will be used instead.
*/
virtual void SetNumWorkers(const unsigned int num) = 0;
LIB_UTILITIES_EXPORT virtual void SetNumWorkers(const unsigned int num) = 0;
/**
* @brief Sets the number of active workers to the maximum.
*
* Sets the number of active workers to the maximum available.
*/
virtual void SetNumWorkers() = 0;
LIB_UTILITIES_EXPORT virtual void SetNumWorkers() = 0;
/**
* @brief Gets the maximum available number of threads.
* @return The maximum number of workers.
*/
virtual unsigned int GetMaxNumWorkers() = 0;
LIB_UTILITIES_EXPORT virtual unsigned int GetMaxNumWorkers() = 0;
/**
* @brief Waits until all queued jobs are finished.
*
......@@ -261,7 +261,7 @@ class ThreadManager : public boost::enable_shared_from_this<ThreadManager>
* number of worker threads, implementations should increase the number
* of active workers by 1 on entering Wait().
*/
virtual void Wait() = 0;
LIB_UTILITIES_EXPORT virtual void Wait() = 0;
/**
* @brief Controls how many jobs are sent to each worker at a time.
*
......@@ -271,17 +271,17 @@ class ThreadManager : public boost::enable_shared_from_this<ThreadManager>
* @see SchedType
* @see SetSchedType()
*/
virtual void SetChunkSize(unsigned int chnk) = 0;
LIB_UTILITIES_EXPORT virtual void SetChunkSize(unsigned int chnk) = 0;
/**
* @brief Sets the current scheduling algorithm.
* @see SetChunkSize()
*/
virtual void SetSchedType(SchedType s) = 0;
LIB_UTILITIES_EXPORT virtual void SetSchedType(SchedType s) = 0;
/**
* @brief Indicates whether the code is in a worker thread or not.
* @return True if the caller is in a worker thread.
*/
virtual bool InThread() = 0;
LIB_UTILITIES_EXPORT virtual bool InThread() = 0;
/**
* @brief A calling threads holds until all active threads call this
* method.
......@@ -294,16 +294,16 @@ class ThreadManager : public boost::enable_shared_from_this<ThreadManager>
* is altered after a thread has called this method. It is only safe to
* call SetNumWorkers() when no threads are holding.
*/
virtual void Hold() = 0;
LIB_UTILITIES_EXPORT virtual void Hold() = 0;
/**
* @brief Returns a description of the type of threading.
*
* E.g. "Threading with Boost"
*/
virtual const std::string& GetType() const = 0;
LIB_UTILITIES_EXPORT virtual const std::string& GetType() const = 0;
/// ThreadManager implementation.
virtual bool IsInitialised();
LIB_UTILITIES_EXPORT virtual bool IsInitialised();
inline int GetThrFromPartition(int pPartition)
{
......@@ -354,17 +354,17 @@ class ThreadMaster
THREADMANAGER_MAX
};
/// Constructor
ThreadMaster();
LIB_UTILITIES_EXPORT ThreadMaster();
/// Destructor
~ThreadMaster();
LIB_UTILITIES_EXPORT ~ThreadMaster();
/// Sets what ThreadManagers will be created in CreateInstance.
void SetThreadingType(const std::string &p_type);
LIB_UTILITIES_EXPORT void SetThreadingType(const std::string &p_type);
/// Gets the ThreadManager associated with string s.
ThreadManagerSharedPtr& GetInstance(const ThreadManagerName t);
LIB_UTILITIES_EXPORT ThreadManagerSharedPtr& GetInstance(const ThreadManagerName t);
/// Creates an instance of a ThreadManager (which one is determined by
/// a previous call to SetThreadingType) and associates it with
/// the string s.
ThreadManagerSharedPtr CreateInstance(const ThreadManagerName t,
LIB_UTILITIES_EXPORT ThreadManagerSharedPtr CreateInstance(const ThreadManagerName t,
unsigned int nThr);
};
......
library/LibUtilities/BasicUtils/Vmath.cpp
View file @ ec3ebd52
......@@ -135,30 +135,44 @@ namespace Vmath
template LIB_UTILITIES_EXPORT Nektar::NekDouble ran2 (long* idum);
/// \brief Fills a vector with white noise.
template<class T> void FillWhiteNoise( int n, const T eps, T *x, const int incx, int outseed)
template<class T> void FillWhiteNoise( int n, const T eps, T *x,
const int incx, int outseed)
{
// Protect the static vars here and in ran2
boost::mutex::scoped_lock l(mutex);
// Define static variables for generating random numbers
static int iset = 0;
static T gset;
static long seed = 0;
// Bypass seed if outseed was specified
if( outseed != 9999)
{
seed = long(outseed);
}
while( n-- )
{
static int iset = 0;
static T gset;
long seed = long(outseed);
T fac, rsq, v1, v2;
if (iset == 0) {
do {
v1 = 2.0 * ran2<T> (&seed) - 1.0;
v2 = 2.0 * ran2<T> (&seed) - 1.0;
rsq = v1*v1 + v2*v2;
} while (rsq >= 1.0 || rsq == 0.0);
fac = sqrt(-2.0 * log (rsq) / rsq);
gset = v1 * fac;
iset = 1;
*x = eps * v2 * fac;
} else {
iset = 0;
*x = eps * gset;
if (iset == 0)
{
do
{
v1 = 2.0 * ran2<T> (&seed) - 1.0;
v2 = 2.0 * ran2<T> (&seed) - 1.0;
rsq = v1*v1 + v2*v2;
} while (rsq >= 1.0 || rsq == 0.0);
fac = sqrt(-2.0 * log (rsq) / rsq);
gset = v1 * fac;
iset = 1;
*x = eps * v2 * fac;
}
else
{
iset = 0;
*x = eps * gset;
}
x += incx;
}
......
library/LibUtilities/BasicUtils/Vmath.hpp
View file @ ec3ebd52
......@@ -54,7 +54,8 @@ namespace Vmath
/// \brief Fills a vector with white noise.
template<class T> LIB_UTILITIES_EXPORT void FillWhiteNoise( int n, const T eps, T *x, const int incx, int seed = 0);
template<class T> LIB_UTILITIES_EXPORT void FillWhiteNoise(
int n, const T eps, T *x, const int incx, int seed = 9999);
/// \brief Multiply vector z = x*y
template<class T> LIB_UTILITIES_EXPORT void Vmul( int n, const T *x, const int incx, const T *y,
......
library/LibUtilities/BasicUtils/VmathArray.hpp
View file @ ec3ebd52
......@@ -53,7 +53,9 @@
Fill(n,alpha,&x[0],incx);
}
template<class T> void FillWhiteNoise( int n, const T eps, Array<OneD, T> &x, const int incx, int outseed)
template<class T> void FillWhiteNoise(
int n, const T eps, Array<OneD, T> &x,
const int incx, int outseed = 9999)
{
ASSERTL1(n*incx <= x.num_elements()+x.GetOffset(),"Out of bounds");
......
library/LibUtilities/Interpreter/AnalyticExpressionEvaluator.cpp
View file @ ec3ebd52
......@@ -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);
......
library/LibUtilities/Interpreter/AnalyticExpressionEvaluator.hpp
View file @ ec3ebd52
......@@ -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) {}
......
library/LibUtilities/Polylib/Polylib.cpp
View file @ ec3ebd52
......@@ -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;