Commit 8cdcfa17 authored by Simon Clifford's avatar Simon Clifford

Merge remote-tracking branch 'origin/master' into feature/threading

parents 68496b2c fa4ca62a
FIND_PATH(MKL_INCLUDE_DIR mkl_cblas.h /usr/include
/usr/local/include
/opt/intel/mkl/8.1.1/include
$ENV{MKL_HOME}/include)
FIND_PATH(MKL_INCLUDE_DIR mkl_cblas.h
HINTS $ENV{MKL_HOME}/include $ENV{MKLROOT}/include
)
IF( CMAKE_SYSTEM_PROCESSOR STREQUAL "x86_64" )
GET_FILENAME_COMPONENT(MKL_LIB_PATH ${MKL_INCLUDE_DIR}/../lib/64 ABSOLUTE)
SET(INTEL_LIBDIR "intel64")
ELSE ( CMAKE_SYSTEM_PROCESSOR STREQUAL "x86_64" )
GET_FILENAME_COMPONENT(MKL_LIB_PATH ${MKL_INCLUDE_DIR}/../lib/32 ABSOLUTE)
SET(INTEL_LIBDIR "ia32")
ENDIF( CMAKE_SYSTEM_PROCESSOR STREQUAL "x86_64" )
FIND_LIBRARY( MKL_LAPACK NAMES mkl_lapack PATHS ${MKL_LIB_PATH} )
FIND_LIBRARY(MKL_LAPACK
NAMES mkl_lapack
HINTS $ENV{MKL_HOME}/lib/${INTEL_LIBDIR} $ENV{MKLROOT}/lib/${INTEL_LIBDIR}
)
IF ( MKL_LAPACK_FOUND )
# old MKL versions
IF ( CMAKE_SYSTEM_PROCESSOR STREQUAL "x86_64" )
FIND_LIBRARY( MKL NAMES mkl_ia64 PATHS ${MKL_LIB_PATH} )
ELSE ( CMAKE_SYSTEM_PROCESSOR STREQUAL "x86_64" )
FIND_LIBRARY( MKL NAMES mkl_ia32 PATHS ${MKL_LIB_PATH} )
ENDIF( CMAKE_SYSTEM_PROCESSOR STREQUAL "x86_64" )
FIND_LIBRARY( MKL_GUIDE NAMES guide PATHS ${MKL_LIB_PATH} )
SET (MKL ${MKL} ${MKL_LAPACK} ${MKL_GUIDE} )
# old MKL versions
IF ( CMAKE_SYSTEM_PROCESSOR STREQUAL "x86_64" )
FIND_LIBRARY(MKL
NAMES mkl_ia64
HINTS $ENV{MKL_HOME}/lib/${INTEL_LIBDIR}
$ENV{MKLROOT}/lib/${INTEL_LIBDIR}
)
ELSE ( CMAKE_SYSTEM_PROCESSOR STREQUAL "x86_64" )
FIND_LIBRARY(MKL
NAMES mkl_ia32
HINTS $ENV{MKL_HOME}/lib/${INTEL_LIBDIR}
$ENV{MKLROOT}/lib/${INTEL_LIBDIR}
)
ENDIF( CMAKE_SYSTEM_PROCESSOR STREQUAL "x86_64" )
FIND_LIBRARY(MKL_GUIDE
NAMES guide
HINTS $ENV{MKL_HOME}/lib/${INTEL_LIBDIR}
$ENV{MKLROOT}/lib/${INTEL_LIBDIR}
)
SET (MKL ${MKL} ${MKL_LAPACK} ${MKL_GUIDE} )
ELSE (MKL_LAPACK_FOUND )
# newer MKL version
SET (MKL_LAPACK "")
FIND_LIBRARY( MKL_INTEL NAMES mkl_intel_lp64 PATHS ${MKL_LIB_PATH} )
FIND_LIBRARY( MKL_SEQUENTIAL NAMES mkl_sequential PATHS ${MKL_LIB_PATH} )
FIND_LIBRARY( MKL_CORE NAMES mkl_core PATHS ${MKL_LIB_PATH} )
SET (MKL ${MKL_INTEL} ${MKL_SEQUENTIAL} ${MKL_CORE} )
# newer MKL version
SET (MKL_LAPACK "")
FIND_LIBRARY(MKL_INTEL
NAMES mkl_intel_lp64
HINTS $ENV{MKL_HOME}/lib/${INTEL_LIBDIR}
$ENV{MKLROOT}/lib/${INTEL_LIBDIR}
)
FIND_LIBRARY(MKL_SEQUENTIAL
NAMES mkl_sequential
HINTS $ENV{MKL_HOME}/lib/${INTEL_LIBDIR}
$ENV{MKLROOT}/lib/${INTEL_LIBDIR}
)
FIND_LIBRARY(MKL_CORE
NAMES mkl_core
HINTS $ENV{MKL_HOME}/lib/${INTEL_LIBDIR}
$ENV{MKLROOT}/lib/${INTEL_LIBDIR}
)
SET (MKL ${MKL_INTEL} ${MKL_SEQUENTIAL} ${MKL_CORE} )
ENDIF ( MKL_LAPACK_FOUND )
......@@ -35,15 +64,23 @@ SET( MKL_BLAS_INCLUDE_FILE ${MKL_INCLUDE_DIR}/mkl_blas.h )
SET( MKL_LAPACK_INCLUDE_FILE ${MKL_INCLUDE_DIR}/mkl_lapack.h )
IF (MKL_INCLUDE_DIR)
SET(MKL_FOUND ON)
SET(MKL_FOUND ON)
ENDIF (MKL_INCLUDE_DIR)
IF (MKL_FOUND OR MKL_INTEL_FOUND )
IF (NOT MKL_FIND_QUIETLY)
MESSAGE(STATUS "Found MKL: ${MKL_INCLUDE_DIR}")
ENDIF (NOT MKL_FIND_QUIETLY)
IF (NOT MKL_FIND_QUIETLY)
MESSAGE(STATUS "Found MKL: ${MKL_INCLUDE_DIR}")
ENDIF (NOT MKL_FIND_QUIETLY)
ELSE(MKL_FOUND OR MKL_INTEL_FOUND )
IF (MKL_FIND_REQUIRED)
MESSAGE(FATAL_ERROR "Could not find MKL")
ENDIF (MKL_FIND_REQUIRED)
IF (MKL_FIND_REQUIRED)
MESSAGE(FATAL_ERROR "Could not find MKL")
ENDIF (MKL_FIND_REQUIRED)
ENDIF (MKL_FOUND OR MKL_INTEL_FOUND )
MARK_AS_ADVANCED(MKL_INCLUDE_DIR)
MARK_AS_ADVANCED(MKL_LAPACK)
MARK_AS_ADVANCED(MKL)
MARK_AS_ADVANCED(MKL_GUIDE)
MARK_AS_ADVANCED(MKL_INTEL)
MARK_AS_ADVANCED(MKL_SEQUENTIAL)
MARK_AS_ADVANCED(MKL_CORE)
......@@ -41,7 +41,7 @@ IF (NEKTAR_USE_FFTW)
CONFIGURE_COMMAND ${TPSRC}/fftw-3.2.2/configure --prefix=${TPDIST} --quiet --enable-shared --disable-dependency-tracking
)
SET(FFTW_LIBRARY fftw CACHE FILEPATH
SET(FFTW_LIBRARY fftw3 CACHE FILEPATH
"FFTW library" FORCE)
SET(FFTW_INCLUDE_DIR ${TPDIST}/include CACHE FILEPATH
"FFTW include" FORCE)
......
......@@ -4,7 +4,7 @@
The aim of APESolver is to predict aerodynamic sound generation. Through
the application of a splitting technique, the flow-induced acoustic field is
totally decoupled from the underlying incompressible hydrodynamic field. The
acoustic perturbation equations proposed by Ewert and Shroeder are employed as
acoustic perturbation equations (APE-1/APE-4) proposed by Ewert and Shroeder are employed as
the governing equations of the acoustic field and they assure stable
aeroacoustic simulation due to the suppression of the term related to the
production of perturbed vorticity. These equations are similar to the linearised
......@@ -12,16 +12,21 @@ perturbed compressible equations, but while in the original formulation the flow
decomposition is based on solenoidal vortical perturbations as well as
irrotational acoustic perturbations, in this case perturbations are assumed to
be exclusively of acoustic nature.
\begin{align*}
\frac{\partial \mathbf{u}'}{\partial t}
+\nabla(\mathbf{u}' \cdot \mathbf{U})+\frac{1}{\rho_0}\nabla p' &= 0 \\
\frac{\partial p'}{\partial t}
+\nabla \cdot (\gamma P \mathbf{u}' + p'\mathbf{U})&=-\frac{DP'}{Dt}
\end{align*}
where $(\mathbf{U},P)$ represents the base flow, $(\mathbf{u}',p')$
the perturbations and $\mathrm{D}/\mathrm{D}t$ the material derivative.
$P'=P-p_{\infty}$ is the acoustic source term, with $p_\infty$
the pressure at a reference value.
\begin{subequations}
\begin{align*}
\frac{\partial p'}{\partial t}
+ \overline{c}^2 \frac{\partial \overline{\rho} u'_i}{\partial x_i}
+ \overline{c}^2 \frac{\partial \overline{u}_i p' / \overline{c}^2}{\partial x_i}
&= \overline{c}^2 q_c
\\
\frac{\partial u'_i}{\partial t}
+ \frac{\partial \overline{u}_j u'_j}{\partial x_i}
+ \frac{\partial p' / \overline{\rho}}{\partial x_i}
&= 0
\end{align*}
\end{subequations}
where $(\overline{u}_i,\overline{p}, \overline{\rho}, \overline{c}^2 = \gamma \overline{p} / \overline{\rho} )$ represents the base flow and $(u'_i,p')$ the perturbations.
$\overline{c}^2 q_c$ is the acoustic source term.
\section{Usage}
\begin{lstlisting}[style=BashInputStyle]
......@@ -40,15 +45,13 @@ Currently, only \inltt{APEUpwind} supported.
\subsection{Parameters}
\begin{itemize}
\item \inltt{Rho0}: Density
\item \inltt{Gamma}: Ratio of specific heats
\item \inltt{Pinfinity}: Ambient pressure
\end{itemize}
\subsection{Functions}
\begin{itemize}
\item \inltt{BaseFlow} Baseflow $(\mathbf{U},P)$ defined by the variables \inltt{U0,V0,W0,P0}
\item \inltt{Source} Source term $P'=P-p_{\infty}$
\item \inltt{BaseFlow} Baseflow $(\overline{u}_i, \overline{p}, \overline{\rho})$ defined by the variables \inltt{u0,v0,w0,p0,rho0}
\item \inltt{Source} Source term $\overline{c}^2 q_c$
\item \inltt{InitialConditions}
\end{itemize}
......@@ -57,14 +60,14 @@ Currently, only \inltt{APEUpwind} supported.
\subsection{Aeroacoustic Wave Propagation}
In this section we explain how to set up a simple simulation of aeroacoustics in
Nektar++. We will study the propagation of an acoustic wave in the simple case
where the base flow is $\mathbf{U}=0, P=p_{\infty}=10^6$. The geometry consists
where the base flow is $\overline{u}_i = 0, \, \overline{p}=p_{\infty}=10^6, \, \overline{\rho} = \rho_0 = 1.204$. The geometry consists
of $64$ quadrilateral elements, as shown in Fig.~\ref{f:apesolver:geometry}.
\begin{figure}
\centering
\includegraphics[width=0.5\linewidth]{Figures/APE_Geometry.png}
\caption{Geometry used for the example case of modelling propagation of
acoustic waves where $\mathbf{U}=0, P=p_{\infty}=10^6$}
acoustic waves where $\overline{u}_i = 0, \, \overline{p}=p_{\infty}=10^6, \, \overline{\rho} = \rho_0 = 1.204$}
\label{f:apesolver:geometry}
\end{figure}
......@@ -86,8 +89,8 @@ ambient pressure.
<P> TimeStep = 0.00001 </P>
<P> NumSteps = 150 </P>
<P> FinTime = TimeStep*NumSteps </P>
<P> Rho0 = 1.204 </P> <!-- Incompressible density -->
<P> Gamma = 1.4 </P> <!-- Ratio of specific heats -->
<P> Rho0 = 1.204 </P> <!-- Incompressible density -->
<P> Gamma = 1.4 </P> <!-- Ratio of specific heats -->
<P> Pinfinity = 100000 </P> <!-- Ambient pressure -->
\end{lstlisting}
......@@ -99,9 +102,10 @@ eventual source terms using the following functions:
\begin{lstlisting}[style=XmlStyle]
<FUNCTION NAME="Baseflow">
<E VAR="U0" VALUE="0" />
<E VAR="V0" VALUE="0" />
<E VAR="P0" VALUE="Pinfinity" />
<E VAR="u0" VALUE="0" />
<E VAR="v0" VALUE="0" />
<E VAR="p0" VALUE="Pinfinity" />
<E VAR="rho0" VALUE="Rho0" />
</FUNCTION>
<FUNCTION NAME="Source">
......
......@@ -8,7 +8,7 @@ A useful tool implemented in Nektar++ is the incompressible Navier Stokes solver
\begin{subequations}
\begin{equation}
\frac{\partial \mathbf{V}}{\partial t} + \mathbf{V} \cdot \nabla \mathbf{V} = -\nabla p + \nu \nabla^2 \mathbf{V} + f
\frac{\partial \mathbf{V}}{\partial t} + \mathbf{V} \cdot \nabla \mathbf{V} = -\nabla p + \nu \nabla^2 \mathbf{V} + \mathbf{f}
\end{equation}
\begin{equation}
......@@ -22,19 +22,33 @@ where $V$ is the velocity, $p$ the pressure and $\nu$ the kinematic viscosity. T
\begin{enumerate}
\item Calculate a first intermediate velocity field evaluating the advection term explicitly and combining it with the solution at previous time-steps:
\begin{equation}
\frac{\tilde{\mathbf{V}}-\sum_{q=0}^{J-1}\frac{\alpha_q}{\gamma_0}\mathbf{V}^{n-q}}{\Delta t}=\sum_{q=0}^{J-1}\frac{\beta_q}{\gamma_0}[-(\mathbf{V}\cdot\nabla)\mathbf{V}]^{n-q} + f^{n+1}\end{equation}
\frac{\tilde{\mathbf{V}}-\sum_{q=0}^{J-1}\frac{\alpha_q}{\gamma_0}\mathbf{V}^{n-q}}{\Delta t}=\sum_{q=0}^{J-1}\frac{\beta_q}{\gamma_0}[-(\mathbf{V}\cdot\nabla)\mathbf{V}]^{n-q} + \mathbf{f}^{n+1}\end{equation}
\item Solve a Poisson equation to obtain the pressure solution at the new time level:
\begin{equation}
\Delta p^{n+1}=\Bigl(\frac{\gamma_0}{\Delta t}\Bigr)\nabla\cdot\tilde{\mathbf{V}}
\end{equation}
with consistent boundary conditions:
with consistent Neumann boundary conditions prescribed as
\begin{equation}
\frac{\partial p}{\partial n}^{n+1}= - \Bigl[ \frac{\partial \mathbf{V}}{\partial t}^{n+1} + \nu \sum_{q=0}^{J-1}\beta_q(\nabla \times \nabla \times \mathbf{V})^{n-q} + \sum_{q=0}^{J-1}\beta_q [(\mathbf{V} \cdot \nabla)\mathbf{V}]^{n-q}\Bigr]\cdot n
\frac{\partial p}{\partial n}^{n+1}= - \Bigl[ \frac{\partial \mathbf{V}}{\partial t}^{n+1} + \nu \sum_{q=0}^{J-1}\beta_q(\nabla \times \nabla \times \mathbf{V})^{n-q} + \sum_{q=0}^{J-1}\beta_q [(\mathbf{V} \cdot \nabla)\mathbf{V}]^{n-q}\Bigr]\cdot \mathbf{n}
\end{equation}
and on the outflow boundary, either setting $ p^{n+1}=0$ or using high order Dirichlet boundary condition, which is given as\cite{DoKa14}
\begin{equation}
p^{n+1}= \nu \sum_{q=0}^{J-1}\nabla\mathbf{V}^{n-q}\cdot \mathbf{n}-\frac{1}{2}
\mid \sum_{q=0}^{J-1}\mathbf{V}^{n-q} \mid^2 S_o(\mathbf{n}\cdot \sum_{q=0}^{J-1}\mathbf{V}^{n-q})-\mathbf{f}_b^{n+1}\cdot \mathbf{n}
\end{equation}
with a step function defined by $S_o(n\cdot
\mathbf{V})=\frac{1}{2}(1-\tanh\frac{n\cdot\mathbf{V}}{\mathbf{V_0}\delta})$,
where $\mathbf{V_0}$ is the characteristic velocity scale and $\delta$
is a non-dimensional positive constant chosen to be sufficiently
small. $\mathbf{f}_b$ is the forcing term in this case the analytical
conditions can be given but if these are not known explicitly, it is
set to zero, i.e. $\mathbf{f}_b=0$.
\item Calculate a second intermediate velocity field:
......@@ -48,9 +62,17 @@ with consistent boundary conditions:
\Bigl(\Delta-\frac{\gamma_0}{\Delta t \nu}\Bigr)\mathbf{V}^{n+1}=-\Bigl(\frac{\gamma_0}{\Delta t \nu}\Bigr)\tilde{\tilde{\mathbf{V}}}
\end{equation}
with the outflow boundary conditions either set as $n\cdot\nabla\mathbf{V}^{n+1}=0$ or as high order Neumann boundary conditions prescribed by\cite{DoKa14}:
\begin{equation}
\mathbf{n}\cdot\nabla\mathbf{V}^{n+1}=\frac{1}{\nu}\Bigl[p^{n+1}\mathbf{n}+\frac{1}{2}
\mid\sum_{q=0}^{J-1} \mathbf{V}^{n-q} \mid^2 S_o(\mathbf{n}\cdot
\sum_{q=0}^{J-1}\mathbf{V}^{n-q})\mathbf{n}-\nu\sum_{q=0}^{J-1}(\nabla\cdot\mathbf{V}^{n-q})\mathbf{n}+\mathbf{f}_b^{n+1}\Bigr]
\end{equation}
\end{enumerate}
Here, $J$ is the integration order and $\gamma_0$, $alpha_q$, $\beta_q$
Here, $J$ is the integration order and $\gamma_0$, $\alpha_q$, $\beta_q$
are the stiffly stable time integration coefficients given in the table below.
\begin{equation}
......@@ -66,7 +88,7 @@ are the stiffly stable time integration coefficients given in the table below.
\end{array}
\end{equation}
This multi-step implicit-explicit splitting scheme decouples the velocity field $mathbf{V}$ from the pressure $p$, leading to an
This multi-step implicit-explicit splitting scheme decouples the velocity field $\mathbf{V}$ from the pressure $p$, leading to an
explicit treatment of the advection term and an implicit treatment of the pressure and the diffusion term.
\subsection{Direct solver (coupled approach)}
......@@ -848,6 +870,83 @@ The result should look similar to that shown in Figure~\ref{f:incns:kovaflow}.
\end{center}
\end{figure}
\subsection{Kovasznay Flow 2D using high-order outflow boundary conditions}
\label{s:incns:kovasznay2D_HOBC}
In this example, we solve the same case of 2D Kovasznay flow on
severely-truncated computational domain but using a high-order outflow
boundary condition, which is much more accurate and robust for
unbounded flows~\cite{DoKa14}. The solver information and parameters
used here are similar to the previous one. What only we need to modify
in the input file is just the boundary condition type upon the outlet
region shown as below
\begin{lstlisting}[style=XMLStyle]
<BOUNDARYCONDITIONS>
<REGION REF="0">
<D VAR="u" VALUE="1-exp(KovLam*x)*cos(2*PI*y)" />
<D VAR="v" VALUE="(KovLam/2/PI)*exp(KovLam*x)*sin(2*PI*y)" />
<N VAR="p" USERDEFINEDTYPE="H" VALUE="0" />
</REGION>
<REGION REF="1">
<N VAR="u" USERDEFINEDTYPE="HOutflow"
VALUE="-Kinvis*KovLam*exp(KovLam*x)*cos(2*PI*y)
- 0.5*(1-exp(2*KovLam*x))-0.5*(((1-exp(KovLam*x)*cos(2*PI*y))
*(1-exp(KovLam*x)*cos(2*PI*y))+(KovLam/(2*PI)*exp(KovLam*x)
*sin(2*PI*y))*(KovLam/(2*PI)*exp(KovLam*x)*sin(2*PI*y))))
*(0.5*(1.0-tanh((1-exp(KovLam*x)*cos(2*PI*y))*20)))" />
<N VAR="v" USERDEFINEDTYPE="HOutflow"
VALUE="Kinvis*KovLam*KovLam/(2*PI)*exp(KovLam*x)*sin(2*PI*y)" />
<D VAR="p" USERDEFINEDTYPE="HOutflow"
VALUE="-Kinvis*KovLam*exp(KovLam*x)*cos(2*PI*y)
- 0.5*(1-exp(2*KovLam*x)) -0.5*(((1-exp(KovLam*x)*cos(2*PI*y))
*(1-exp(KovLam*x)*cos(2*PI*y))+(KovLam/(2*PI)*exp(KovLam*x)
*sin(2*PI*y))*(KovLam/(2*PI)*exp(KovLam*x)*sin(2*PI*y))))
*(0.5*(1.0-tanh((1-exp(KovLam*x)*cos(2*PI*y))*20)))" />
</REGION>
<REGION REF="2">
<N VAR="u" VALUE="0" />
<D VAR="v" VALUE="0" />
<N VAR="p" VALUE="0" />
</REGION>
</BOUNDARYCONDITIONS>
\end{lstlisting}
We note that in this example the ``VALUE'' property is set based on
the analytic solution but this is not typically known and so often a
VALUE of zero will be specified.
Instead of loading an initial condition from a specified file, we initialized the flow fields in this example by using following expressions
\begin{lstlisting}[style=XMLStyle]
<FUNCTION NAME="InitialConditions">
<E VAR="u" VALUE="(1-exp(KovLam*x)*cos(2*PI*y))" />
<E VAR="v" VALUE="(KovLam/(2*PI))*exp(KovLam*x)*sin(2*PI*y)" />
<E VAR="p" VALUE="0.5*(1-exp(2*KovLam*x))" />
</FUNCTION>
\end{lstlisting}
\subsubsection{Running the simulation}
We then launch the simulation by the same solver as that in the previous example
\begin{lstlisting}[style=BashInputStyle]
IncNavierStokesSolver KovaFlow_m8_short_HOBC.xml
\end{lstlisting}
The solution with errors displayed as below
\begin{lstlisting}[style=BashInputStyle]
L 2 error (variable u) : 2.51953e-08
L inf error (variable u) : 9.56014e-09
L 2 error (variable v) : 1.10694e-08
L inf error (variable v) : 9.47464e-08
L 2 error (variable p) : 5.59175e-08
L inf error (variable p) : 2.93085e-07
\end{lstlisting}
The physical solution visualized in velocity profiles is also illustrated in
Figure~\ref{f:incns:kovaflow2d_hobc}.
\begin{figure}
\begin{center}
\includegraphics[width=7cm]{Figures/KF2DCVP8HOBC_U.png}
\includegraphics[width=7cm]{Figures/KF2DCVP8HOBC_V.png}
\caption{Velocity profiles for the Kovasznay Flow in truncated domain (2D).}
\label{f:incns:kovaflow2d_hobc}
\end{center}
\end{figure}
\subsection{Steady Kovasznay Oseen Flow using the direct solver}
\label{s:incns:kovasznay2Ddirect}
......
......@@ -267,3 +267,13 @@ in the human arterial system},
year={1997},
publisher={Cambridge Univ Press}
}
@article{DoKa14,
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},
pages={95--136},
year={2014}
}
......@@ -404,6 +404,9 @@ namespace Nektar
ImportMultiFldFileIDs(infile,filenames, elementIDs_OnPartitions,
fieldmetadatamap);
// Load metadata
ImportFieldMetaData(infile,fieldmetadatamap);
if(ElementIDs == NullInt1DArray) //load all fields
{
for(int i = 0; i < filenames.size(); ++i)
......
......@@ -949,8 +949,8 @@ namespace Nektar
x->SetAttribute("ID", vVertIt->first);
std::stringstream vCoords;
vCoords.precision(12);
vCoords << std::setw(15) << vVertIt->second.x
<< std::setw(15) << vVertIt->second.y
vCoords << std::setw(15) << vVertIt->second.x << " "
<< std::setw(15) << vVertIt->second.y << " "
<< std::setw(15) << vVertIt->second.z << " ";
y = new TiXmlText(vCoords.str());
x->LinkEndChild(y);
......
......@@ -451,5 +451,10 @@ IF( NEKTAR_USE_BLAS_LAPACK )
ENDIF( NEKTAR_USE_BLAS_LAPACK )
IF( NEKTAR_USE_PETSC )
TARGET_LINK_LIBRARIES(LibUtilities LINK_PRIVATE ${PETSC_LIBRARIES})
ADD_DEPENDENCIES(LibUtilities petsc-3.5.2)
ENDIF( NEKTAR_USE_PETSC )
INSTALL(FILES ${ExpressionTemplates} DESTINATION ${NEKTAR_INCLUDE_DIR}/ExpressionTemplates COMPONENT dev)
INSTALL(DIRECTORY ./ DESTINATION ${NEKTAR_INCLUDE_DIR}/LibUtilities COMPONENT dev FILES_MATCHING PATTERN "*.h" PATTERN "*.hpp")
......@@ -33,6 +33,10 @@
//
///////////////////////////////////////////////////////////////////////////////
#ifdef NEKTAR_USING_PETSC
#include "petscsys.h"
#endif
#include <LibUtilities/Communication/CommMpi.h>
#include <LibUtilities/BasicUtils/SharedArray.hpp>
......@@ -66,6 +70,10 @@ namespace Nektar
MPI_Comm_size( m_comm, &m_size );
MPI_Comm_rank( m_comm, &m_rank );
#ifdef NEKTAR_USING_PETSC
PetscInitializeNoArguments();
#endif
m_type = "Parallel MPI";
}
......@@ -107,6 +115,9 @@ namespace Nektar
*/
void CommMpi::v_Finalise()
{
#ifdef NEKTAR_USING_PETSC
PetscFinalize();
#endif
MPI_Finalize();
}
......
......@@ -68,8 +68,8 @@ namespace Nektar
for(int i=2;i<m_N;i++)
{
m_FFTW_w[i] = m_FFTW_w[0];
m_FFTW_w_inv[i] = m_FFTW_w_inv[0];
m_FFTW_w[i] = m_FFTW_w[0]*2;
m_FFTW_w_inv[i] = m_FFTW_w_inv[0]/2;
}
}
......
......@@ -563,11 +563,11 @@ namespace Nektar
{
for(i = 0; i < numPoints; ++i)
{
m_bdata[ 2*p *numPoints+i] = cos(p*M_PI*z[i]);
m_bdata[(2*p+1)*numPoints+i] = -sin(p*M_PI*z[i]);
m_bdata[ 2*p *numPoints+i] = cos(p*M_PI* (z[i]+1) );
m_bdata[(2*p+1)*numPoints+i] = -sin(p*M_PI* (z[i]+1) );
m_dbdata[ 2*p *numPoints+i] = -p*M_PI*sin(p*M_PI*z[i]);
m_dbdata[(2*p+1)*numPoints+i] = -p*M_PI*cos(p*M_PI*z[i]);
m_dbdata[ 2*p *numPoints+i] = -p*M_PI*sin(p*M_PI* (z[i]+1) );
m_dbdata[(2*p+1)*numPoints+i] = -p*M_PI*cos(p*M_PI* (z[i]+1) );
}
}
......@@ -579,11 +579,11 @@ namespace Nektar
for(i = 0; i < numPoints; ++i)
{
m_bdata[i] = cos(M_PI*z[i]);
m_bdata[numPoints+i] = -sin(M_PI*z[i]);
m_bdata[i] = cos(M_PI* (z[i]+1) );
m_bdata[numPoints+i] = -sin(M_PI* (z[i]+1) );
m_dbdata[i] = -M_PI*sin(M_PI*z[i]);
m_dbdata[numPoints+i] = -M_PI*cos(M_PI*z[i]);
m_dbdata[i] = -M_PI*sin(M_PI* (z[i]+1) );
m_dbdata[numPoints+i] = -M_PI*cos(M_PI* (z[i]+1) );
}
for (p=1; p < numModes/2; ++p)
......
This diff is collapsed.
......@@ -972,6 +972,7 @@ namespace Nektar
LibUtilities::PointsKey points0;
LibUtilities::PointsKey points1;
Array<OneD, NekDouble> faceJac(nqtot);
Array<OneD, NekDouble> normals(vCoordDim*nqtot,0.0);
// Extract Jacobian along face and recover local derivatives
......@@ -986,6 +987,7 @@ namespace Nektar
normals[j] = -df[2][j]*jac[j];
normals[nqtot+j] = -df[5][j]*jac[j];
normals[2*nqtot+j] = -df[8][j]*jac[j];
faceJac[j] = jac[j];
}
points0 = ptsKeys[0];
......@@ -1006,6 +1008,7 @@ namespace Nektar
-df[4][tmp]*jac[tmp];
normals[2*nqtot+j+k*nq0] =
-df[7][tmp]*jac[tmp];
faceJac[j+k*nq0] = jac[tmp];
}
}
......@@ -1027,6 +1030,7 @@ namespace Nektar
(df[3][tmp]+df[5][tmp])*jac[tmp];
normals[2*nqtot+j+k*nq1] =
(df[6][tmp]+df[8][tmp])*jac[tmp];
faceJac[j+k*nq1] = jac[tmp];
}
}
......@@ -1048,6 +1052,7 @@ namespace Nektar
df[4][tmp]*jac[tmp];
normals[2*nqtot+j+k*nq0] =
df[7][tmp]*jac[tmp];
faceJac[j+k*nq0] = jac[tmp];
}
}
......@@ -1069,6 +1074,7 @@ namespace Nektar
-df[3][tmp]*jac[tmp];
normals[2*nqtot+j+k*nq1] =
-df[6][tmp]*jac[tmp];
faceJac[j+k*nq1] = jac[tmp];
}
}
......@@ -1084,7 +1090,7 @@ namespace Nektar
Array<OneD, NekDouble> work (nq_face, 0.0);
// Interpolate Jacobian and invert
LibUtilities::Interp2D(points0, points1, jac,
LibUtilities::Interp2D(points0, points1, faceJac,
tobasis0.GetPointsKey(),
tobasis1.GetPointsKey(),
work);
......
......@@ -618,6 +618,7 @@ namespace Nektar
LibUtilities::PointsKey points0;
LibUtilities::PointsKey points1;
Array<OneD, NekDouble> faceJac(nqtot);
Array<OneD, NekDouble> normals(vCoordDim*nqtot,0.0);
// Extract Jacobian along face and recover local derivatives
......@@ -632,6 +633,7 @@ namespace Nektar
normals[j] = -df[2][j]*jac[j];
normals[nqtot+j] = -df[5][j]*jac[j];
normals[2*nqtot+j] = -df[8][j]*jac[j];
faceJac[j] = jac[j];
}
points0 = ptsKeys[0];
......@@ -652,6 +654,7 @@ namespace Nektar
-df[4][tmp]*jac[tmp];
normals[2*nqtot+j+k*nq0] =
-df[7][tmp]*jac[tmp];
faceJac[j+k*nq0] = jac[tmp];
}
}
......@@ -673,6 +676,7 @@ namespace Nektar
(df[3][tmp]+df[5][tmp])*jac[tmp];
normals[2*nqtot+j+k*nq1] =
(df[6][tmp]+df[8][tmp])*jac[tmp];
faceJac[j+k*nq1] = jac[tmp];
}
}
......