...
 
Commits (38)
......@@ -14,3 +14,6 @@ ThirdParty
*.pdf
# Kdevelop project files
*.kdev4
# Sublime project files
*.sublime-project
*.sublime-workspace
......@@ -132,6 +132,7 @@ v5.0.0
- Fix surface extraction, added regression test (!994)
- Fix 2D meshing running out of memory due to missing else (!1012)
- Add support for .msh v4.1 file input (!1054)
- Added penalty term to LDG and LDGNS, slight generalization of LDG (!1080)
**FieldConvert**:
- Add input module for Semtex field files (!777)
......@@ -157,6 +158,7 @@ v5.0.0
- Fixed scaling for compressed xml, fixed error printout for mesh only (!1040)
- Add field conversion from Halfmode to SingleMode (!1032)
- Fix double precision output in .dat format (!1059)
- Add phase sampling feature in FilterFieldConvert (!1068)
**IncNavierStokesSolver**
- Replace steady-state check based on difference of norms by check based on
......
......@@ -10,6 +10,12 @@ IF (NEKTAR_BUILD_PYTHON)
CMAKE_DEPENDENT_OPTION(NEKTAR_USE_PYTHON3
"If true, prefer to use Python 3." OFF "NEKTAR_BUILD_PYTHON" OFF)
# Set the Python3 status flag as the opposite of the Python3 flag by
# default on first run to ensure Python is searched for.
IF (NOT DEFINED NEKTAR_PYTHON3_STATUS)
SET(NEKTAR_PYTHON3_STATUS NOT ${NEKTAR_USE_PYTHON3} CACHE INTERNAL "")
ENDIF()
IF (NOT NEKTAR_PYTHON3_STATUS STREQUAL NEKTAR_USE_PYTHON3)
# Unset any existing python executable/library settings so that we can
# rediscover correct python version if v2/3 settings changed.
......@@ -30,9 +36,10 @@ IF (NEKTAR_BUILD_PYTHON)
# Find Python
FIND_PACKAGE(PythonInterp ${PYTHONVER} REQUIRED)
FIND_PACKAGE(PythonLibsNew ${PYTHONVER} REQUIRED)
INCLUDE_DIRECTORIES(${PYTHON_INCLUDE_DIRS})
ADD_DEFINITIONS(-DWITH_PYTHON)
# Save include dir in Cache for subsequent configures.
SET(PYTHON_INCLUDE_DIRS ${PYTHON_INCLUDE_DIRS} CACHE INTERNAL "" FORCE)
# Include headers from root directory for config file.
# Now try to find Boost::Python. For now we are relying entirely on
......@@ -81,6 +88,8 @@ IF (NEKTAR_BUILD_PYTHON)
SET(NEKTAR_PYTHON3_STATUS ${NEKTAR_USE_PYTHON3} CACHE INTERNAL "")
ENDIF()
INCLUDE_DIRECTORIES(${PYTHON_INCLUDE_DIRS})
ADD_DEFINITIONS(-DWITH_PYTHON)
MESSAGE(STATUS "Found Python: ${PYTHON_EXECUTABLE}")
MESSAGE(STATUS "Found Boost.Python: ${BOOST_PYTHON_LIB}")
......
......@@ -11,14 +11,14 @@ The ADRSolver is designed to solve partial differential equations of the form:
\begin{equation}
\alpha \dfrac{\partial u}{\partial t} + \lambda u + \nu \nabla u + \epsilon \nabla \cdot (D \nabla u) = f
\end{equation}
in either discontinuous or continuous projections of the solution field.
For a full list of the equations which are supported, and the capabilities of each equation,
in either discontinuous or continuous projections of the solution field.
For a full list of the equations which are supported, and the capabilities of each equation,
see the table below.
\begin{table}[h!]
\begin{center}
\tiny
\renewcommand\arraystretch{2.2}
\renewcommand\arraystretch{2.2}
\begin{tabular}{llll}
\toprule
\textbf{Equation to solve} & \textbf{EquationType} & \textbf{Dimensions} &
......@@ -26,22 +26,22 @@ see the table below.
\midrule
$u = f$ &
\inltt{Projection} & All & Continuous/Discontinuous \\
$\nabla^2 u = 0$ &
$\nabla^2 u = 0$ &
\inltt{Laplace} & All & Continuous/Discontinuous \\
$\nabla^2 u = f$ &
$\nabla^2 u = f$ &
\inltt{Poisson} & All & Continuous/Discontinuous \\
$\nabla^2 u + \lambda u = f$ &
$\nabla^2 u + \lambda u = f$ &
\inltt{Helmholtz} & All & Continuous/Discontinuous \\
$\epsilon \nabla^2 u + \mathbf{V}\nabla u = f$ &
$\epsilon \nabla^2 u + \mathbf{V}\nabla u = f$ &
\inltt{SteadyAdvectionDiffusion} & 2D only & Continuous/Discontinuous \\
$\epsilon \nabla^2 u + \lambda u = f$ &
$\epsilon \nabla^2 u + \lambda u = f$ &
\inltt{SteadyDiffusionReaction} & 2D only & Continuous/Discontinuous \\
$\epsilon \nabla^2 u + \mathbf{V}\nabla u + \lambda u = f$ &
\inltt{SteadyAdvectionDiffusionReaction} & 2D only &
$\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$ &
\inltt{UnsteadyAdvection} & All & Continuous/Discontinuous \\
$\dfrac{\partial u}{\partial t} = \epsilon \nabla^2 u$ &
$\dfrac{\partial u}{\partial t} = \epsilon \nabla^2 u$ &
\inltt{UnsteadyDiffusion} & All & Continuous/Discontinuous \\
$\dfrac{\partial u}{\partial t} = \epsilon \nabla^2 u + R(u)$ &
\inltt{UnsteadyReactionDiffusion} & All & Continuous \\
......@@ -64,7 +64,7 @@ ADRSolver session.xml
\section{Session file configuration}
The type of equation which is to be solved is specified through the EquationType
The type of equation which is to be solved is specified through the EquationType
SOLVERINFO option in the session file. This can be set as in table \ref{t:ADR1}.
At present, the Steady non-symmetric solvers cannot be used in parallel. \\
......@@ -76,7 +76,7 @@ The solver info are listed below:
\item \textbf{TimeIntegrationMethod}: The following types of time integration methods have been tested with each solver:
\begin{center}
\footnotesize
\renewcommand\arraystretch{1.2}
\renewcommand\arraystretch{1.2}
\begin{tabular}{lcccc}
\toprule
\textbf{EqType} & \textbf{Explicit} & \textbf{Diagonally Implicit} &
......@@ -116,7 +116,7 @@ The solver info are listed below:
\end{itemize}
\item \textbf{DiffusionType}:
\begin{itemize}
\item \inltt{LDG}.
\item \inltt{LDG} (The penalty term is proportional to an optional parameter \inltt{LDGc11} which is by default set to one; proportionality to polynomial order can be manually imposed by setting the parameter \inltt{LDGc11} equal to $p^2$).
\end{itemize}
\item \textbf{UpwindType}:
\begin{itemize}
......@@ -129,13 +129,13 @@ The solver info are listed below:
The following parameters can be specified in the \inltt{PARAMETERS} section of
the session file:
\begin{itemize}
\item \inltt{epsilon}: sets the diffusion coefficient $\epsilon$.\\
\item \inltt{epsilon}: sets the diffusion coefficient $\epsilon$.\\
\textit{Can be used} in: SteadyDiffusionReaction, SteadyAdvectionDiffusionReaction, UnsteadyDiffusion, UnsteadyAdvectionDiffusion. \\
\textit{Default value}: 0.
\item \inltt{d00}, \inltt{d11}, \inltt{d22}: sets the diagonal entries of the
diffusion tensor $D$. \\
\textit{Can be used in}: UnsteadyDiffusion \\
\textit{Default value}: All set to 1 (i.e. identity matrix).
\textit{Default value}: All set to 1 (i.e. identity matrix).
\item \inltt{lambda}: sets the reaction coefficient $\lambda$. \\
\textit{Can be used in}: SteadyDiffusionReaction, Helmholtz, SteadyAdvectionDiffusionReaction\\
\textit{Default value}: 0.
......@@ -215,7 +215,7 @@ and then define the actual advection function as
Two boundary regions are defined, one at each end of the domain, and periodicity
is enforced
\begin{lstlisting}[style=XMLStyle]
\begin{lstlisting}[style=XMLStyle]
<BOUNDARYREGIONS>
<B ID="0"> C[1] </B>
<B ID="1"> C[2] </B>
......@@ -268,7 +268,7 @@ We consider the elliptic partial differential equation:
where $\nabla^2$ is the Laplacian and $\lambda$ is a real positive constant.
\subsubsection{Input file}
\subsubsection{Input file}
The input for this example is given in the example file
\inlsh{Helmholtz2D\_modal.xml}
......@@ -300,7 +300,7 @@ basis with $7$ modes (maximum polynomial order is $6$).
\end{lstlisting}
Only one parameter is needed for this problem. In this example $\lambda = 1$ and
the Continuous Galerkin Method is used as projection scheme to solve the
the Continuous Galerkin Method is used as projection scheme to solve the
Helmholtz equation, so we need to specify the following parameters and solver
information.
\begin{lstlisting}[style=XMLStyle]
......@@ -331,7 +331,7 @@ region is then assigned an appropriate boundary condition.
<D VAR="u" VALUE="sin(PI*x)*sin(PI*y)" />
</REGION>
<REGION REF="1">
<R VAR="u" VALUE="sin(PI*x)*sin(PI*y)-PI*sin(PI*x)*cos(PI*y)"
<R VAR="u" VALUE="sin(PI*x)*sin(PI*y)-PI*sin(PI*x)*cos(PI*y)"
PRIMCOEFF="1" />
</REGION>
<REGION REF="2">
......@@ -343,7 +343,7 @@ region is then assigned an appropriate boundary condition.
</BOUNDARYCONDITIONS>
\end{lstlisting}
We know that for $f = -(\lambda + 2 \pi^2)sin(\pi x)cos(\pi y)$, the exact
We know that for $f = -(\lambda + 2 \pi^2)sin(\pi x)cos(\pi y)$, the exact
solution of the two-dimensional Helmholtz equation is $u = sin(\pi x)cos(\pi
y)$. These functions are defined specified to initialise the problem and verify
the correct solution is obtained by evaluating the $L_2$ and $L_{inf}$ errors.
......@@ -363,7 +363,7 @@ the correct solution is obtained by evaluating the $L_2$ and $L_{inf}$ errors.
ADRSolver Test_Helmholtz2D_modal.xml
\end{lstlisting}
This execution should print out a summary of input file, the $L_2$ and
This execution should print out a summary of input file, the $L_2$ and
$L_{inf}$ errors and the time spent on the calculation.
\subsubsection{Post-processing}
......@@ -417,9 +417,9 @@ Schmidt number (the relative thickness of the momentum to mass transfer boundary
layer) is sufficiently large. The analytical solution for the non-dimensional
mass transfer at the wall is given by:
\begin{align*}
S h(z) = \dfrac{2^{4/3}(Pe R/z)^{1/3}}{g^{1/3}\Gamma(4/3)} ,
S h(z) = \dfrac{2^{4/3}(Pe R/z)^{1/3}}{g^{1/3}\Gamma(4/3)} ,
\end{align*}
where $z$ is the streamwise coordinate, $R$ the pipe radius, $\Gamma(4/3)$ an incomplete
where $z$ is the streamwise coordinate, $R$ the pipe radius, $\Gamma(4/3)$ an incomplete
Gamma function and $Pe$ the Peclet number given by:
\begin{align*}
Pe = \dfrac{2 U R}{\epsilon}
......@@ -538,7 +538,7 @@ flow), hence we can define this through an analytical function as follows:
\end{lstlisting}
We assume that the initial domain concentration is uniform everywhere and the
same as the inlet. This is defined by,
same as the inlet. This is defined by,
\begin{lstlisting}[style=XMLStyle]
<FUNCTION NAME="InitialConditions">
<E VAR="u" VALUE="1" />
......
......@@ -272,7 +272,7 @@ Note that only \inltt{WeakDG} is fully supported, the other operators work only
\item \inltt{DiffusionType} is the diffusion operator we want to use
for the Navier-Stokes equations:
\begin{itemize}
\item \inltt{LDGNS} (LDG);
\item \inltt{LDGNS} (LDG with primitive variables. The penalty term is inversely proportional to the element size, proportional to the local viscosity for the momentum equations and to the thermal conductivity for the energy equation, and proportional to an optional parameter \inltt{LDGNSc11} which is by default set to one; proportionality to polynomial order can be manually imposed by setting the parameter \inltt{LDGNSc11} equal to $p^2$);
\item \inltt{LFRDGNS} (Flux-Reconstruction recovering nodal DG scheme);
\item \inltt{LFRSDNS} (Flux-Reconstruction recovering a spectral difference (SD) scheme);
\item \inltt{LFRHUNS} (Flux-Reconstruction recovering Huynh (G2) scheme);
......
......@@ -35,6 +35,62 @@ In the following we document the filters implemented. Note that some filters are
solver-specific and will therefore only work for a given subset of the available
solvers.
\subsection{Phase sampling}
\begin{notebox}
This feature is currently only supported for filters derived from the
FieldConvert filter: AverageFields, MovingAverage, ReynoldsStresses.
\end{notebox}
When analysing certain time-dependent problems, it might be of interest to
activate a filter in a specific physical phase and with a certain period
(for instance, to carry out phase averaging).
The simulation time can be written as
$t = m \mathcal{T} + n_{\mathcal{T}} \mathcal{T}$,
where $m$ is an integer representing the number of periods $\mathcal{T}$
elapsed, and $0 \leq n_{\mathcal{T}} \leq 1$ is the phase.
This feature is not a filter in itself and it is activated by adding the
parameters below to the filter of interest:
\begin{center}
\begin{tabularx}{0.99\textwidth}{lllX}
\toprule
\textbf{Option name} & \textbf{Required} & \textbf{Default} &
\textbf{Description} \\
\midrule
\inltt{PhaseAverage} & \cmark & &
Feature activation\\
\inltt{PhaseAveragePeriod} & \cmark & &
Period $\mathcal{T}$\\
\inltt{PhaseAveragePhase} & \cmark & &
Phase $n_{\mathcal{T}}$.\\
\bottomrule
\end{tabularx}
\end{center}
For instance, to activate phase averaging with a period of $\mathcal{T}=10$
at phase $n_{\mathcal{T}}=0.5$:
\begin{lstlisting}[style=XMLStyle,gobble=2]
<FILTER TYPE="FilterName">
<PARAM NAME="Param1"> Value1 </PARAM>
<PARAM NAME="Param2"> Value2 </PARAM>
<PARAM NAME="PhaseAverage"> True </PARAM>
<PARAM NAME="PhaseAveragePeriod"> 10 </PARAM>
<PARAM NAME="PhaseAveragePhase"> 0.5 </PARAM>
</FILTER>
\end{lstlisting}
Since this feature monitors $n_{\mathcal{T}}$ every \inltt{SampleFrequency},
for best results it is recommended to set \inltt{SampleFrequency}$=1$. \\
The maximum error in sampling phase is $n_{\mathcal{T}, \text{tol}} =
\frac{\Delta t}{ 2\mathcal{T}}\cdot$ \inltt{SampleFrequency}, which is displayed
at the beginning of the simulation if the solver is run with the verbose
\inltt{-v} option.\\
The number of periods elapsed is calculated based on simulation time. Caution
is therefore recommended when modifying time information in the restart field,
because if the new time does not correspond to the same phase, the feature
will produce erroneous results.
\subsection{Aerodynamic forces}\label{filters:AeroForces}
\begin{notebox}
......
......@@ -1214,13 +1214,13 @@ void Mapping::v_UpdateBCs( const NekDouble time)
if( BndConds[n]->GetUserDefined() =="" ||
BndConds[n]->GetUserDefined() =="MovingBody")
{
BndExp[n]->FwdTrans_BndConstrained(BndExp[n]->GetPhys(),
BndExp[n]->UpdateCoeffs());
if (m_fields[i]->GetExpType() == MultiRegions::e3DH1D)
{
BndExp[n]->HomogeneousFwdTrans(BndExp[n]->GetCoeffs(),
BndExp[n]->UpdateCoeffs());
BndExp[n]->SetWaveSpace(false);
}
BndExp[n]->FwdTrans_BndConstrained(
BndExp[n]->GetPhys(), BndExp[n]->UpdateCoeffs());
}
}
}
......
......@@ -705,6 +705,12 @@ namespace Nektar
// Now fill in all other Dirichlet coefficients.
for(int i = 0; i < m_bndCondExpansions.num_elements(); ++i)
{
if (m_bndConditions[i]->GetBoundaryConditionType() ==
SpatialDomains::ePeriodic)
{
continue;
}
Array<OneD, NekDouble>& coeffs = m_bndCondExpansions[i]->UpdateCoeffs();
for(int j = 0; j < (m_bndCondExpansions[i])->GetNcoeffs(); ++j)
......@@ -734,7 +740,13 @@ namespace Nektar
Array<OneD, NekDouble>& coeffs = m_bndCondExpansions[nreg]->UpdateCoeffs();
for(int j = 0; j < nreg; ++j)
{
{
if (m_bndConditions[j]->GetBoundaryConditionType() ==
SpatialDomains::ePeriodic)
{
continue;
}
bndcnt += m_bndCondExpansions[j]->GetNcoeffs();
}
......
......@@ -217,42 +217,7 @@ namespace Nektar
m_useFFT, false,
PlanesBndCondExp, comm);
}
EvaluateBoundaryConditions(0.0, variable);
}
void DisContField3DHomogeneous1D::EvaluateBoundaryConditions(
const NekDouble time,
const std::string varName)
{
int n;
const Array<OneD, const NekDouble> z = m_homogeneousBasis->GetZ();
Array<OneD, NekDouble> local_z(m_planes.num_elements());
for (n = 0; n < m_planes.num_elements(); n++)
{
local_z[n] = z[m_transposition->GetPlaneID(n)];
}
for (n = 0; n < m_planes.num_elements(); ++n)
{
m_planes[n]->EvaluateBoundaryConditions(
time, varName, 0.5*m_lhom*(1.0+local_z[n]));
}
// Fourier transform coefficient space boundary values
// This will only be undertaken for time dependent
// boundary conditions unless time == 0.0 which is the
// case when the method is called from the constructor.
for (n = 0; n < m_bndCondExpansions.num_elements(); ++n)
{
if (time == 0.0 ||
m_bndConditions[n]->IsTimeDependent() )
{
m_bndCondExpansions[n]->HomogeneousFwdTrans(
m_bndCondExpansions[n]->GetCoeffs(),
m_bndCondExpansions[n]->UpdateCoeffs());
}
}
v_EvaluateBoundaryConditions(0.0, variable);
}
void DisContField3DHomogeneous1D::v_HelmSolve(
......@@ -315,9 +280,127 @@ namespace Nektar
const NekDouble x2_in,
const NekDouble x3_in)
{
boost::ignore_unused(x2_in, x3_in);
boost::ignore_unused(x2_in,x3_in);
int i;
int npoints;
int nbnd = m_bndCondExpansions.num_elements();
MultiRegions::ExpListSharedPtr locExpList;
for (i = 0; i < nbnd; ++i)
{
if (time == 0.0 || m_bndConditions[i]->IsTimeDependent())
{
locExpList = m_bndCondExpansions[i];
npoints = locExpList->GetNpoints();
Array<OneD, NekDouble> x0(npoints, 0.0);
Array<OneD, NekDouble> x1(npoints, 0.0);
Array<OneD, NekDouble> x2(npoints, 0.0);
Array<OneD, NekDouble> valuesFile(npoints, 1.0), valuesExp(npoints, 1.0);
locExpList->GetCoords(x0, x1, x2);
if (m_bndConditions[i]->GetBoundaryConditionType()
== SpatialDomains::eDirichlet)
{
SpatialDomains::DirichletBCShPtr bcPtr = std::static_pointer_cast<
SpatialDomains::DirichletBoundaryCondition>(
m_bndConditions[i]);
std::string filebcs = bcPtr->m_filename;
std::string exprbcs = bcPtr->m_expr;
if (filebcs != "")
{
ExtractFileBCs(filebcs, bcPtr->GetComm(), varName, locExpList);
valuesFile = locExpList->GetPhys();
}
if (exprbcs != "")
{
LibUtilities::Equation condition =
std::static_pointer_cast<SpatialDomains::
DirichletBoundaryCondition >(
m_bndConditions[i])->m_dirichletCondition;
condition.Evaluate(x0, x1, x2, time, valuesExp);
}
EvaluateBoundaryConditions(time, varName);
Vmath::Vmul(npoints, valuesExp, 1, valuesFile, 1,
locExpList->UpdatePhys(), 1);
// set wave space to false since have set up phys values
locExpList->SetWaveSpace(false);
locExpList->FwdTrans_BndConstrained(
locExpList->GetPhys(),
locExpList->UpdateCoeffs());
}
else if (m_bndConditions[i]->GetBoundaryConditionType()
== SpatialDomains::eNeumann)
{
SpatialDomains::NeumannBCShPtr bcPtr = std::static_pointer_cast<
SpatialDomains::NeumannBoundaryCondition>(
m_bndConditions[i]);
std::string filebcs = bcPtr->m_filename;
if (filebcs != "")
{
ExtractFileBCs(filebcs, bcPtr->GetComm(), varName, locExpList);
}
else
{
LibUtilities::Equation condition = std::
static_pointer_cast<SpatialDomains::
NeumannBoundaryCondition>(
m_bndConditions[i])->m_neumannCondition;
condition.Evaluate(x0, x1, x2, time,
locExpList->UpdatePhys());
locExpList->IProductWRTBase(locExpList->GetPhys(),
locExpList->UpdateCoeffs());
}
}
else if (m_bndConditions[i]->GetBoundaryConditionType()
== SpatialDomains::eRobin)
{
SpatialDomains::RobinBCShPtr bcPtr = std::static_pointer_cast<
SpatialDomains::RobinBoundaryCondition>(
m_bndConditions[i]);
std::string filebcs = bcPtr->m_filename;
if (filebcs != "")
{
ExtractFileBCs(filebcs, bcPtr->GetComm(), varName, locExpList);
}
else
{
LibUtilities::Equation condition = std::
static_pointer_cast<SpatialDomains::
RobinBoundaryCondition>(
m_bndConditions[i])->m_robinFunction;
condition.Evaluate(x0, x1, x2, time,
locExpList->UpdatePhys());
}
locExpList->IProductWRTBase(locExpList->GetPhys(),
locExpList->UpdateCoeffs());
}
else if (m_bndConditions[i]->GetBoundaryConditionType()
== SpatialDomains::ePeriodic)
{
continue;
}
else
{
ASSERTL0(false, "This type of BC not implemented yet");
}
}
}
}
std::shared_ptr<ExpList> &DisContField3DHomogeneous1D::
......
......@@ -801,7 +801,7 @@ namespace Nektar
}
void ExpList::FwdTrans_BndConstrained(
void ExpList::v_FwdTrans_BndConstrained(
const Array<OneD, const NekDouble>& inarray,
Array<OneD, NekDouble> &outarray)
{
......
......@@ -1297,6 +1297,10 @@ namespace Nektar
const Array<OneD,const NekDouble> &inarray,
Array<OneD,NekDouble> &outarray);
virtual void v_FwdTrans_BndConstrained(
const Array<OneD,const NekDouble> &inarray,
Array<OneD,NekDouble> &outarray);
virtual void v_SmoothField(Array<OneD,NekDouble> &field);
virtual void v_IProductWRTBase(
......@@ -1746,6 +1750,17 @@ namespace Nektar
v_FwdTrans_IterPerExp(inarray,outarray);
}
/**
*
*/
inline void ExpList::FwdTrans_BndConstrained (
const Array<OneD, const NekDouble> &inarray,
Array<OneD,NekDouble> &outarray)
{
v_FwdTrans_BndConstrained(inarray,outarray);
}
/**
*
*/
......
......@@ -483,7 +483,9 @@ namespace Nektar
/**
* Forward transform element by element
*/
void ExpListHomogeneous1D::v_FwdTrans_IterPerExp(const Array<OneD, const NekDouble> &inarray, Array<OneD, NekDouble> &outarray)
void ExpListHomogeneous1D::v_FwdTrans_IterPerExp(const Array<OneD,
const NekDouble> &inarray,
Array<OneD, NekDouble> &outarray)
{
int cnt = 0, cnt1 = 0;
Array<OneD, NekDouble> tmparray;
......@@ -500,6 +502,29 @@ namespace Nektar
HomogeneousFwdTrans(outarray,outarray);
}
}
/**
* Forward transform element by element with boundaries constrained
*/
void ExpListHomogeneous1D::v_FwdTrans_BndConstrained(const Array<OneD,
const NekDouble> &inarray,
Array<OneD, NekDouble> &outarray)
{
int cnt = 0, cnt1 = 0;
Array<OneD, NekDouble> tmparray;
//spectral element FwdTrans plane by plane
for(int n = 0; n < m_planes.num_elements(); ++n)
{
m_planes[n]->FwdTrans_BndConstrained(inarray+cnt, tmparray = outarray + cnt1);
cnt += m_planes[n]->GetTotPoints();
cnt1 += m_planes[n]->GetNcoeffs();
}
if(!m_WaveSpace)
{
HomogeneousFwdTrans(outarray,outarray);
}
}
/**
* Backward transform
......
......@@ -190,7 +190,10 @@ namespace Nektar
virtual void v_FwdTrans_IterPerExp(const Array<OneD,const NekDouble> &inarray,
Array<OneD, NekDouble> &outarray);
virtual void v_FwdTrans_BndConstrained(const Array<OneD,const NekDouble> &inarray,
Array<OneD, NekDouble> &outarray);
virtual void v_BwdTrans(const Array<OneD,const NekDouble> &inarray,
Array<OneD, NekDouble> &outarray,
CoeffState coeffstate);
......
......@@ -43,16 +43,16 @@ namespace Nektar
static DiffusionFactory instance;
return instance;
}
void Diffusion::InitObject(
const LibUtilities::SessionReaderSharedPtr pSession,
Array<OneD, MultiRegions::ExpListSharedPtr> pFields)
{
v_InitObject(pSession, pFields);
}
void Diffusion::Diffuse(
const int nConvectiveFields,
const std::size_t nConvectiveFields,
const Array<OneD, MultiRegions::ExpListSharedPtr> &fields,
const Array<OneD, Array<OneD, NekDouble> > &inarray,
Array<OneD, Array<OneD, NekDouble> > &outarray,
......
......@@ -52,23 +52,28 @@ namespace Nektar
namespace SolverUtils
{
typedef std::function<void (
const int,
const int,
const Array<OneD, Array<OneD, NekDouble> >&,
Array<OneD, Array<OneD, NekDouble> >&,
Array<OneD, Array<OneD, NekDouble> >&)> DiffusionFluxVecCB;
const Array<OneD, Array<OneD, NekDouble> > &,
const Array<OneD, Array<OneD, Array<OneD, NekDouble> > >&,
Array<OneD, Array<OneD, Array<OneD, NekDouble> > >&)>
DiffusionFluxVecCB;
typedef std::function<void (
const Array<OneD, Array<OneD, NekDouble> >&,
Array<OneD, Array<OneD, Array<OneD, NekDouble> > >&,
Array<OneD, Array<OneD, Array<OneD, NekDouble> > >&)>
Array<OneD, Array<OneD, Array<OneD, NekDouble> > >&)>
DiffusionFluxVecCBNS;
typedef std::function<void (
const Array<OneD, Array<OneD, NekDouble> >&,
const Array<OneD, Array<OneD, NekDouble> >&,
Array<OneD, Array<OneD, NekDouble> >&)>
DiffusionFluxPenaltyNS;
typedef std::function<void (
const Array<OneD, Array<OneD, NekDouble> >&,
Array<OneD, NekDouble >&)>
DiffusionArtificialDiffusion;
class Diffusion
{
public:
......@@ -79,34 +84,33 @@ namespace Nektar
SOLVER_UTILS_EXPORT void InitObject(
LibUtilities::SessionReaderSharedPtr pSession,
Array<OneD, MultiRegions::ExpListSharedPtr> pFields);
SOLVER_UTILS_EXPORT void Diffuse(
const int nConvectiveFields,
const std::size_t nConvectiveFields,
const Array<OneD, MultiRegions::ExpListSharedPtr> &fields,
const Array<OneD, Array<OneD, NekDouble> > &inarray,
Array<OneD, Array<OneD, NekDouble> > &outarray,
const Array<OneD, Array<OneD, NekDouble> > &pFwd = NullNekDoubleArrayofArray,
const Array<OneD, Array<OneD, NekDouble> > &pBwd = NullNekDoubleArrayofArray);
SOLVER_UTILS_EXPORT void FluxVec(
Array<OneD, Array<OneD, Array<OneD, NekDouble> > >
&fluxvector);
template<typename FuncPointerT, typename ObjectPointerT>
template<typename FuncPointerT, typename ObjectPointerT>
void SetFluxVector(FuncPointerT func, ObjectPointerT obj)
{
m_fluxVector = std::bind(
func, obj, std::placeholders::_1, std::placeholders::_2,
std::placeholders::_3, std::placeholders::_4,
std::placeholders::_5);
std::placeholders::_3);
}
void SetFluxVectorVec(DiffusionFluxVecCB fluxVector)
void SetFluxVector(DiffusionFluxVecCB fluxVector)
{
m_fluxVector = fluxVector;
}
template<typename FuncPointerT, typename ObjectPointerT>
template<typename FuncPointerT, typename ObjectPointerT>
void SetFluxVectorNS(FuncPointerT func, ObjectPointerT obj)
{
m_fluxVectorNS = std::bind(
......@@ -114,16 +118,22 @@ namespace Nektar
std::placeholders::_3);
}
void SetFluxVectorNS(DiffusionFluxVecCBNS fluxVector)
{
m_fluxVectorNS = fluxVector;
}
template<typename FuncPointerT, typename ObjectPointerT>
void SetArtificialDiffusionVector(FuncPointerT func, ObjectPointerT obj)
void SetFluxPenaltyNS(FuncPointerT func, ObjectPointerT obj)
{
m_ArtificialDiffusionVector = std::bind(
func, obj, std::placeholders::_1, std::placeholders::_2);
m_fluxPenaltyNS = std::bind(
func, obj, std::placeholders::_1, std::placeholders::_2,
std::placeholders::_3);
}
void SetFluxVectorNS(DiffusionFluxVecCBNS fluxVector)
void SetFluxPenaltyNS(DiffusionFluxPenaltyNS flux)
{
m_fluxVectorNS = fluxVector;
m_fluxPenaltyNS = flux;
}
inline void SetHomoDerivs(Array<OneD, Array<OneD, NekDouble> > &deriv)
......@@ -135,11 +145,11 @@ namespace Nektar
{
return v_GetFluxTensor();
}
protected:
DiffusionFluxVecCB m_fluxVector;
DiffusionFluxVecCBNS m_fluxVectorNS;
DiffusionArtificialDiffusion m_ArtificialDiffusionVector;
DiffusionFluxPenaltyNS m_fluxPenaltyNS;
virtual void v_InitObject(
LibUtilities::SessionReaderSharedPtr pSession,
......@@ -147,9 +157,9 @@ namespace Nektar
{
boost::ignore_unused(pSession, pFields);
};
virtual void v_Diffuse(
const int nConvectiveFields,
const std::size_t nConvectiveFields,
const Array<OneD, MultiRegions::ExpListSharedPtr> &fields,
const Array<OneD, Array<OneD, NekDouble> > &inarray,
Array<OneD, Array<OneD, NekDouble> > &outarray,
......@@ -161,17 +171,17 @@ namespace Nektar
{
boost::ignore_unused(deriv);
}
virtual Array<OneD, Array<OneD, Array<OneD, NekDouble> > > &v_GetFluxTensor()
{
static Array<OneD, Array<OneD, Array<OneD, NekDouble> > > tmp;
return tmp;
}
};
};
/// A shared pointer to an EquationSystem object
typedef std::shared_ptr<Diffusion> DiffusionSharedPtr;
/// Datatype of the NekFactory used to instantiate classes derived
/// from the Diffusion class.
typedef LibUtilities::NekFactory<std::string, Diffusion, std::string> DiffusionFactory;
......
......@@ -72,7 +72,7 @@ namespace Nektar
};
/**
* @brief Diffusion3DHomogeneous1D uses the 2D WeakDG approach
* @brief Diffusion3DHomogeneous1D uses the 2D WeakDG approach
* to compute the diffusion term looping on the planes in the z
* direction and adding the flux in z direction at the end.
*/
......@@ -80,12 +80,12 @@ namespace Nektar
{
// Strip trailing string "3DHomogeneous1D" to determine 2D diffusion
// type, and create a diffusion object for the plane.
string name = diffType.substr(0, diffType.length()-15);
m_planeDiff = GetDiffusionFactory().CreateInstance(name, name);
m_diffType = diffType.substr(0, diffType.length()-15);
m_planeDiff = GetDiffusionFactory().CreateInstance(m_diffType, m_diffType);
}
/**
* @brief Initiliase Diffusion3DHomogeneous1D objects and store
* @brief Initiliase Diffusion3DHomogeneous1D objects and store
* them before starting the time-stepping.
*
* @param pSession Pointer to session reader.
......@@ -114,7 +114,26 @@ namespace Nektar
m_homoLen = pFields[0]->GetHomoLen();
m_trans = pFields[0]->GetTransposition();
m_planeCounter = 0;
m_planeDiff->SetFluxVectorNS(m_fluxVectorNS);
if (m_diffType == "LDG")
{
// Set viscous flux for LDG
m_planeDiff->SetFluxVector(m_fluxVector);
}
else if (m_diffType == "LDGNS")
{
// Set viscous flux for LDGNS
m_planeDiff->SetFluxVectorNS(m_fluxVectorNS);
// Set penalty flux
m_planeDiff->SetFluxPenaltyNS(m_fluxPenaltyNS);
}
else if (m_diffType == "LFRDGNS" ||
m_diffType == "LFRHUNS" ||
m_diffType == "LFRSDNS" )
{
// Set viscous flux for FR cases
m_planeDiff->SetFluxVectorNS(m_fluxVectorNS);
}
m_fieldsPlane = Array<OneD, MultiRegions::ExpListSharedPtr>
(nConvectiveFields);
......@@ -169,7 +188,7 @@ namespace Nektar
* using an LDG interface flux and the the flux in the third direction.
*/
void Diffusion3DHomogeneous1D::v_Diffuse(
const int nConvectiveFields,
const std::size_t nConvectiveFields,
const Array<OneD, MultiRegions::ExpListSharedPtr> &fields,
const Array<OneD, Array<OneD, NekDouble> > &inarray,
Array<OneD, Array<OneD, NekDouble> > &outarray,
......@@ -181,22 +200,21 @@ namespace Nektar
Array<OneD, NekDouble> tmp(m_numPoints), tmp2;
Array<OneD, Array<OneD, NekDouble> > viscHComp;
const int nPointsTot = fields[0]->GetNpoints();
int i, j;
NekDouble beta;
if (m_fluxVectorNS)
{
viscHComp = Array<OneD, Array<OneD, NekDouble> >(nConvectiveFields);
for (i = 0; i < nConvectiveFields - 1; ++i)
for (int i = 0; i < nConvectiveFields - 1; ++i)
{
fields[0]->PhysDeriv(2, inarray[i], m_homoDerivStore[i]);
viscHComp[i] = Array<OneD, NekDouble>(m_numPoints);
}
}
for (i = 0; i < m_numPlanes; ++i)
for (int i = 0; i < m_numPlanes; ++i)
{
// Set up memory references for fields, inarray and outarray for
// this plane.
......@@ -212,7 +230,7 @@ namespace Nektar
m_outarrayPlane[j] = Array<OneD, NekDouble>(
m_numPointsPlane, tmp2 = outarray[j] + m_planePos[i]);
}
if (m_fluxVectorNS)
{
......@@ -220,15 +238,43 @@ namespace Nektar
}
m_planeDiff->Diffuse(nConvectiveFields,
m_fieldsPlane,
m_inarrayPlane,
m_outarrayPlane);
if (m_diffType == "LDGNS")
{
// Store plane Fwd/Bwd traces
std::size_t nTracePts = m_fieldsPlane[0]->GetTrace()
->GetTotPoints();
std::size_t nScalar = m_inarrayPlane.num_elements();
Array<OneD, Array<OneD, NekDouble> > Fwd(nScalar);
Array<OneD, Array<OneD, NekDouble> > Bwd(nScalar);
{
for(std::size_t k = 0; k < nScalar; ++k)
{
Fwd[k] = Array<OneD, NekDouble>(nTracePts, 0.0);
Bwd[k] = Array<OneD, NekDouble>(nTracePts, 0.0);
m_fieldsPlane[k]->GetFwdBwdTracePhys(
m_inarrayPlane[k], Fwd[k], Bwd[k]);
}
}
m_planeDiff->Diffuse(nConvectiveFields,
m_fieldsPlane,
m_inarrayPlane,
m_outarrayPlane,
Fwd,
Bwd);
}
else
{
m_planeDiff->Diffuse(nConvectiveFields,
m_fieldsPlane,
m_inarrayPlane,
m_outarrayPlane);
}
if (m_fluxVectorNS)
{
Array<OneD, Array<OneD, Array<OneD, NekDouble> > > &viscTensor = m_planeDiff->GetFluxTensor();
Array<OneD, Array<OneD, Array<OneD, NekDouble> > >
&viscTensor = m_planeDiff->GetFluxTensor();
// Extract H (viscTensor[2])
for (int j = 0; j < nConvectiveFields - 1; ++j)
......@@ -239,24 +285,25 @@ namespace Nektar
}
}
}
if (m_fluxVectorNS)
{
for (j = 0; j < nConvectiveFields - 1; ++j)
for (int j = 0; j < nConvectiveFields - 1; ++j)
{
fields[j+1]->PhysDeriv(2, viscHComp[j], tmp);
Vmath::Vadd(nPointsTot, outarray[j+1], 1, tmp, 1, outarray[j+1], 1);
Vmath::Vadd(nPointsTot, outarray[j+1], 1, tmp, 1,
outarray[j+1], 1);
}
}
else
{
for (j = 0; j < nConvectiveFields; ++j)
for (int j = 0; j < nConvectiveFields; ++j)
{
fields[j]->HomogeneousFwdTrans(inarray[j], tmp);
for (i = 0; i < m_numPlanes; ++i)
for (int i = 0; i < m_numPlanes; ++i)
{
beta = 2*M_PI*m_trans->GetK(i)/m_homoLen;
beta *= beta;
......
......@@ -51,7 +51,7 @@ namespace Nektar
new Diffusion3DHomogeneous1D(diffType));
}
static std::string type[];
protected:
Diffusion3DHomogeneous1D(std::string diffType);
......@@ -59,10 +59,10 @@ namespace Nektar
std::string m_diffType;
SolverUtils::DiffusionSharedPtr m_planeDiff;
NekDouble m_homoLen;
int m_numPoints;
int m_numPointsPlane;
int m_numPlanes;
int m_planeCounter;
std::size_t m_numPoints;
std::size_t m_numPointsPlane;
std::size_t m_numPlanes;
std::size_t m_planeCounter;
Array<OneD, unsigned int> m_planes;
Array<OneD, unsigned int> m_planePos;
Array<OneD, Array<OneD, NekDouble> > m_homoDerivStore;
......@@ -76,14 +76,14 @@ namespace Nektar
LibUtilities::SessionReaderSharedPtr pSession,
Array<OneD, MultiRegions::ExpListSharedPtr> pFields);
virtual void v_Diffuse(
const int nConvective,
const std::size_t nConvective,
const Array<OneD, MultiRegions::ExpListSharedPtr> &fields,
const Array<OneD, Array<OneD, NekDouble> > &inarray,
Array<OneD, Array<OneD, NekDouble> > &outarray,
const Array<OneD, Array<OneD, NekDouble> > &pFwd = NullNekDoubleArrayofArray,
const Array<OneD, Array<OneD, NekDouble> > &pBwd = NullNekDoubleArrayofArray);
};
};
}
}
#endif
......@@ -51,59 +51,63 @@ namespace Nektar
boost::ignore_unused(diffType);
return DiffusionSharedPtr(new DiffusionLDG());
}
static std::string type;
protected:
DiffusionLDG();
std::string m_shockCaptureType;
std::string m_shockCaptureType;
/// Coefficient of penalty term
NekDouble m_C11;
Array<OneD, Array<OneD, NekDouble> > m_traceNormals;
LibUtilities::SessionReaderSharedPtr m_session;
virtual void v_InitObject(
LibUtilities::SessionReaderSharedPtr pSession,
Array<OneD, MultiRegions::ExpListSharedPtr> pFields);
virtual void v_Diffuse(
const int nConvective,
const std::size_t nConvective,
const Array<OneD, MultiRegions::ExpListSharedPtr> &fields,
const Array<OneD, Array<OneD, NekDouble> > &inarray,
Array<OneD, Array<OneD, NekDouble> > &outarray,
const Array<OneD, Array<OneD, NekDouble> > &pFwd = NullNekDoubleArrayofArray,
const Array<OneD, Array<OneD, NekDouble> > &pBwd = NullNekDoubleArrayofArray);
virtual void v_NumFluxforScalar(
void NumFluxforScalar(
const Array<OneD, MultiRegions::ExpListSharedPtr> &fields,
const Array<OneD, Array<OneD, NekDouble> > &ufield,
Array<OneD, Array<OneD, Array<OneD, NekDouble> > >&uflux,
const Array<OneD, Array<OneD, NekDouble> > &pFwd = NullNekDoubleArrayofArray,
const Array<OneD, Array<OneD, NekDouble> > &pBwd = NullNekDoubleArrayofArray);
virtual void v_WeakPenaltyforScalar(
void ApplyScalarBCs(
const Array<OneD, MultiRegions::ExpListSharedPtr> &fields,
const int var,
const std::size_t var,
const Array<OneD, const NekDouble> &ufield,
const Array<OneD, const NekDouble> &uplus,
const Array<OneD, const NekDouble> &Fwd,
const Array<OneD, const NekDouble> &Bwd,
Array<OneD, NekDouble> &penaltyflux);
virtual void v_NumFluxforVector(
void NumFluxforVector(
const Array<OneD, MultiRegions::ExpListSharedPtr> &fields,
const Array<OneD, Array<OneD, NekDouble> > &ufield,
Array<OneD, Array<OneD, Array<OneD, NekDouble> > >&qfield,
Array<OneD, Array<OneD, NekDouble> > &qflux);
virtual void v_WeakPenaltyforVector(
void ApplyVectorBCs(
const Array<OneD, MultiRegions::ExpListSharedPtr> &fields,
const int var,
const int dir,
const std::size_t var,
const std::size_t dir,
const Array<OneD, const NekDouble> &qfield,
const Array<OneD, const NekDouble> &qtemp,
Array<OneD, NekDouble> &penaltyflux,
NekDouble C11);
};
const Array<OneD, const NekDouble> &qFwd,
const Array<OneD, const NekDouble> &qBwd,
Array<OneD, NekDouble> &penaltyflux);
};
}
}
#endif
This source diff could not be displayed because it is too large. You can view the blob instead.
......@@ -48,32 +48,32 @@ namespace Nektar
{
return DiffusionSharedPtr(new DiffusionLFR(diffType));
}
static std::string type[];
Array<OneD, NekDouble> m_jac;
Array<OneD, Array<OneD, NekDouble> > m_gmat;
Array<OneD, Array<OneD, NekDouble> > m_Q2D_e0;
Array<OneD, Array<OneD, NekDouble> > m_Q2D_e1;
Array<OneD, Array<OneD, NekDouble> > m_Q2D_e2;
Array<OneD, Array<OneD, NekDouble> > m_Q2D_e3;
Array<OneD, Array<OneD, NekDouble> > m_dGL_xi1;
Array<OneD, Array<OneD, NekDouble> > m_Q2D_e0;
Array<OneD, Array<OneD, NekDouble> > m_Q2D_e1;
Array<OneD, Array<OneD, NekDouble> > m_Q2D_e2;
Array<OneD, Array<OneD, NekDouble> > m_Q2D_e3;
Array<OneD, Array<OneD, NekDouble> > m_dGL_xi1;
Array<OneD, Array<OneD, NekDouble> > m_dGR_xi1;
Array<OneD, Array<OneD, NekDouble> > m_dGL_xi2;
Array<OneD, Array<OneD, NekDouble> > m_dGL_xi2;
Array<OneD, Array<OneD, NekDouble> > m_dGR_xi2;
Array<OneD, Array<OneD, NekDouble> > m_dGL_xi3;
Array<OneD, Array<OneD, NekDouble> > m_dGL_xi3;
Array<OneD, Array<OneD, NekDouble> > m_dGR_xi3;
DNekMatSharedPtr m_Ixm;
DNekMatSharedPtr m_Ixp;
protected:
DiffusionLFR(std::string diffType);
Array<OneD, Array<OneD, NekDouble> > m_traceNormals;
LibUtilities::SessionReaderSharedPtr m_session;
Array<OneD, Array<OneD, Array<OneD, NekDouble> > > m_IF1;
Array<OneD, Array<OneD, Array<OneD, NekDouble> > > m_DU1;
Array<OneD, Array<OneD, Array<OneD, NekDouble> > > m_DFC1;
......@@ -84,49 +84,49 @@ namespace Nektar
Array<OneD, Array<OneD, Array<OneD, NekDouble> > > m_DFC2;
Array<OneD, Array<OneD, NekDouble> > m_divFD;
Array<OneD, Array<OneD, NekDouble> > m_divFC;
Array<OneD, Array<OneD, Array<OneD, NekDouble> > > m_tmp1;
Array<OneD, Array<OneD, Array<OneD, NekDouble> > > m_tmp2;
std::string m_diffType;
virtual void v_InitObject(
LibUtilities::SessionReaderSharedPtr pSession,
Array<OneD, MultiRegions::ExpListSharedPtr> pFields);
virtual void v_SetupMetrics(
LibUtilities::SessionReaderSharedPtr pSession,
Array<OneD, MultiRegions::ExpListSharedPtr> pFields);
virtual void v_SetupCFunctions(
LibUtilities::SessionReaderSharedPtr pSession,
Array<OneD, MultiRegions::ExpListSharedPtr> pFields);
virtual void v_Diffuse(
const int nConvectiveFields,
const std::size_t nConvectiveFields,
const Array<OneD, MultiRegions::ExpListSharedPtr> &fields,
const Array<OneD, Array<OneD, NekDouble> > &inarray,
Array<OneD, Array<OneD, NekDouble> > &outarray,
const Array<OneD, Array<OneD, NekDouble> > &pFwd = NullNekDoubleArrayofArray,
const Array<OneD, Array<OneD, NekDouble> > &pBwd = NullNekDoubleArrayofArray);
virtual void v_NumFluxforScalar(
const Array<OneD, MultiRegions::ExpListSharedPtr> &fields,
const Array<OneD, Array<OneD, NekDouble> > &ufield,
Array<OneD, Array<OneD, Array<OneD, NekDouble> > >&uflux);
virtual void v_WeakPenaltyforScalar(
const Array<OneD, MultiRegions::ExpListSharedPtr> &fields,
const int var,
const Array<OneD, const NekDouble> &ufield,
Array<OneD, NekDouble> &penaltyflux);
virtual void v_NumFluxforVector(
const Array<OneD, MultiRegions::ExpListSharedPtr> &fields,
const Array<OneD, Array<OneD, NekDouble> > &ufield,
Array<OneD, Array<OneD, Array<OneD, NekDouble> > >&qfield,
Array<OneD, Array<OneD, NekDouble> > &qflux);
virtual void v_WeakPenaltyforVector(
const Array<OneD, MultiRegions::ExpListSharedPtr> &fields,
const int var,
......@@ -134,30 +134,30 @@ namespace Nektar
const Array<OneD, const NekDouble> &qfield,
Array<OneD, NekDouble> &penaltyflux,
NekDouble C11);
virtual void v_DerCFlux_1D(
const int nConvectiveFields,
const Array<OneD, MultiRegions::ExpListSharedPtr>&fields,
const Array<OneD, const NekDouble> &flux,
const Array<OneD, const NekDouble> &flux,
const Array<OneD, const NekDouble> &iFlux,
Array<OneD,