Commit d78f9f39 authored by Dave Moxey's avatar Dave Moxey
Browse files

Merge branch 'feature/APE-BCs' into 'master'

new BCs for the APESolver

See merge request !782
parents e1c397fa f2045ad3
......@@ -83,6 +83,10 @@ v5.0.0
- Add ability to use an exponential filtering for stabilization with
seg, quad and hex elements (!771, !862)
**APESolver:**
- Added two new boundary conditions to the APE system: RiemannInvariantBC
and WhiteNoise (!782)
**Documentation**:
- Added the developer-guide repository as a submodule (!751)
......
......@@ -51,11 +51,57 @@ Currently, only \inltt{APEUpwind} supported.
\subsection{Functions}
\begin{itemize}
\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}
\subsection{Boundary Conditions}
In addition to plain Dirichlet and Neumann boundary conditions, the APESolver features a wall boundray condition, a non-reflecting boundary and a white noise boundary condition.
\begin{itemize}
\item Rigid (Slip-) Wall Boundary Condition
\begin{lstlisting}[style=XmlStyle]
<BOUNDARYCONDITIONS>
<REGION REF="0">
<D VAR="p" USERDEFINEDTYPE="Wall" VALUE="0" />
<D VAR="u" USERDEFINEDTYPE="Wall" VALUE="0" />
<D VAR="v" USERDEFINEDTYPE="Wall" VALUE="0" />
<D VAR="w" USERDEFINEDTYPE="Wall" VALUE="0" />
</REGION>
</BOUNDARYCONDITIONS>
\end{lstlisting}
This BC imposes zero wall-normal perturbation velocity in a way that is more robust than using a Dirichlet boundary condition directly.
\item Non-Reflecting Boundary Condition
\begin{lstlisting}[style=XmlStyle]
<BOUNDARYCONDITIONS>
<REGION REF="0">
<D VAR="p" USERDEFINEDTYPE="RiemannInvariantBC"/>
<D VAR="u" USERDEFINEDTYPE="RiemannInvariantBC"/>
<D VAR="v" USERDEFINEDTYPE="RiemannInvariantBC"/>
<D VAR="w" USERDEFINEDTYPE="RiemannInvariantBC"/>
</REGION>
</BOUNDARYCONDITIONS>
\end{lstlisting}
The Riemann-Invariant BC approximates a non-reflecting (r.g. Farfield) boundary condition by setting incoming invariants to zero.
\item White Noise Boundary Condition
\begin{lstlisting}[style=XmlStyle]
<BOUNDARYCONDITIONS>
<REGION REF="0">
<D VAR="p" USERDEFINEDTYPE="Wall" VALUE="10" />
<D VAR="u" USERDEFINEDTYPE="Wall" VALUE="10" />
<D VAR="v" USERDEFINEDTYPE="Wall" VALUE="10" />
<D VAR="w" USERDEFINEDTYPE="Wall" VALUE="10" />
</REGION>
</BOUNDARYCONDITIONS>
\end{lstlisting}
The white noise BC imposes a stochastic, uniform pressure at the boundary. The implementation uses a Mersenne-Twister pseudo random number generator to generate white Gaussian noise.
The standard deviation $\sigma$ of the pressure is specified by the \inltt{VALUE} attribute.
\end{itemize}
\section{Examples}
\subsection{Aeroacoustic Wave Propagation}
In this section we explain how to set up a simple simulation of aeroacoustics in
......@@ -108,10 +154,6 @@ eventual source terms using the following functions:
<E VAR="rho0" VALUE="Rho0" />
</FUNCTION>
<FUNCTION NAME="Source">
<E VAR="S" VALUE="0" />
</FUNCTION>
<!-- Gaussian pulse located at the origin -->
<FUNCTION NAME="InitialConditions">
<E VAR="p" VALUE="100*exp(-32*(((x)*(x))+((y)*(y))))" />
......
......@@ -36,6 +36,9 @@
#include <iostream>
#include <boost/random/variate_generator.hpp>
#include <boost/random/normal_distribution.hpp>
#include <MultiRegions/AssemblyMap/AssemblyMapDG.h>
#include <APESolver/EquationSystems/APE.h>
#include <LocalRegions/TriExp.h>
......@@ -146,10 +149,10 @@ void APE::v_InitObject()
// Define the normal velocity fields
if (m_fields[0]->GetTrace())
{
m_traceBasefield = Array<OneD, Array<OneD, NekDouble> > (m_spacedim + 2);
for (int i = 0; i < m_spacedim + 2; i++)
m_bfTrace = Array<OneD, Array<OneD, NekDouble> > (m_spacedim + 2);
for (int i = 0; i < m_spacedim + 2; ++i)
{
m_traceBasefield[i] = Array<OneD, NekDouble>(GetTraceNpoints());
m_bfTrace[i] = Array<OneD, NekDouble>(GetTraceNpoints(), 0.0);
}
}
......@@ -167,9 +170,9 @@ void APE::v_InitObject()
m_riemannSolver = SolverUtils::GetRiemannSolverFactory().CreateInstance(
riemName);
m_riemannSolver->SetVector("N", &APE::GetNormals, this);
m_riemannSolver->SetVector("basefield", &APE::GetBasefield, this);
m_riemannSolver->SetVector("basefield", &APE::GetBfTrace, this);
m_riemannSolver->SetAuxVec("vecLocs", &APE::GetVecLocs, this);
m_riemannSolver->SetParam("Gamma", &APE::GetGamma, this);
m_riemannSolver->SetParam("Gamma", &APE::GetGamma, this);
// Set up advection operator
string advName;
......@@ -189,6 +192,9 @@ void APE::v_InitObject()
{
ASSERTL0(false, "Implicit APE not set up.");
}
m_whiteNoiseBC_lastUpdate = -1.0;
m_whiteNoiseBC_p = 0.0;
}
......@@ -363,18 +369,45 @@ void APE::SetBoundaryConditions(Array<OneD, Array<OneD, NekDouble> > &inarray,
Fwd[i] = Array<OneD, NekDouble>(nTracePts);
m_fields[i]->ExtractTracePhys(inarray[i], Fwd[i]);
}
Array<OneD, Array<OneD, NekDouble> > bfFwd = GetBfTrace();
// loop over Boundary Regions
for(int n = 0; n < m_fields[0]->GetBndConditions().num_elements(); ++n)
for (int n = 0; n < m_fields[0]->GetBndConditions().num_elements(); ++n)
{
// Wall Boundary Condition
if (boost::iequals(m_fields[0]->GetBndConditions()[n]->GetUserDefined(),"Wall"))
std::string userDefStr =
m_fields[0]->GetBndConditions()[n]->GetUserDefined();
if (!userDefStr.empty())
{
WallBC(n, cnt, Fwd, inarray);
// Wall Boundary Condition
if (boost::iequals(userDefStr, "Wall"))
{
WallBC(n, cnt, Fwd, inarray);
}
else if (boost::iequals(userDefStr, "WhiteNoise"))
{
WhiteNoiseBC(n, cnt, Fwd, bfFwd, inarray);
}
else if (boost::iequals(userDefStr, "RiemannInvariantBC"))
{
RiemannInvariantBC(n, cnt, Fwd, bfFwd, inarray);
}
else if (boost::iequals(userDefStr, "TimeDependent"))
{
for (int i = 0; i < nvariables; ++i)
{
varName = m_session->GetVariable(i);
m_fields[i]->EvaluateBoundaryConditions(time, varName);
}
}
else
{
string errmsg = "Unrecognised boundary condition: ";
errmsg += userDefStr;
ASSERTL0(false, errmsg.c_str());
}
}
// Time Dependent Boundary Condition (specified in meshfile)
if (m_fields[0]->GetBndConditions()[n]->IsTimeDependent())
else
{
for (int i = 0; i < nvariables; ++i)
{
......@@ -382,7 +415,8 @@ void APE::SetBoundaryConditions(Array<OneD, Array<OneD, NekDouble> > &inarray,
m_fields[i]->EvaluateBoundaryConditions(time, varName);
}
}
cnt +=m_fields[0]->GetBndCondExpansions()[n]->GetExpSize();
cnt += m_fields[0]->GetBndCondExpansions()[n]->GetExpSize();
}
}
......@@ -447,6 +481,200 @@ void APE::WallBC(int bcRegion, int cnt,
}
/**
* @brief Outflow characteristic boundary conditions for compressible
* flow problems.
*/
void APE::RiemannInvariantBC(int bcRegion,
int cnt,
Array<OneD, Array<OneD, NekDouble> > &Fwd,
Array<OneD, Array<OneD, NekDouble> > &BfFwd,
Array<OneD, Array<OneD, NekDouble> > &physarray)
{
int id1, id2, nBCEdgePts;
int nVariables = physarray.num_elements();
const Array<OneD, const int> &traceBndMap = m_fields[0]->GetTraceBndMap();
int eMax = m_fields[0]->GetBndCondExpansions()[bcRegion]->GetExpSize();
for (int e = 0; e < eMax; ++e)
{
nBCEdgePts = m_fields[0]
->GetBndCondExpansions()[bcRegion]
->GetExp(e)
->GetTotPoints();
id1 = m_fields[0]->GetBndCondExpansions()[bcRegion]->GetPhys_Offset(e);
id2 = m_fields[0]->GetTrace()->GetPhys_Offset(traceBndMap[cnt + e]);
// Calculate (v.n)
Array<OneD, NekDouble> Vn(nBCEdgePts, 0.0);
for (int i = 0; i < m_spacedim; ++i)
{
Vmath::Vvtvp(nBCEdgePts,
&Fwd[1 + i][id2], 1,
&m_traceNormals[i][id2], 1,
&Vn[0], 1,
&Vn[0], 1);
}
// Calculate (v0.n)
Array<OneD, NekDouble> Vn0(nBCEdgePts, 0.0);
for (int i = 0; i < m_spacedim; ++i)
{
Vmath::Vvtvp(nBCEdgePts,
&BfFwd[2 + i][id2], 1,
&m_traceNormals[i][id2], 1,
&Vn0[0], 1,
&Vn0[0], 1);
}
for (int i = 0; i < nBCEdgePts; ++i)
{
NekDouble c = sqrt(m_gamma * BfFwd[0][id2 + i] / BfFwd[1][id2 + i]);
NekDouble l0 = Vn0[i] + c;
NekDouble l1 = Vn0[i] - c;
NekDouble h0, h1;
// outgoing
if (l0 > 0)
{
// p/2 + u*c*rho0/2
h0 = Fwd[0][id2 + i] / 2 + Vn[i] * c * BfFwd[1][id2 + i] / 2;
}
// incoming
else
{
h0 = 0.0;
}
// outgoing
if (l1 > 0)
{
// p/2 - u*c*rho0/2
h1 = Fwd[0][id2 + i] / 2 - Vn[i] * c * BfFwd[1][id2 + i] / 2;
}
// incoming
else
{
h1 = 0.0;
}
// compute primitive variables
// p = h0 + h1
// u = ( h0 - h1) / (c*rho0)
Fwd[0][id2 + i] = h0 + h1;
NekDouble VnNew = (h0 - h1) / (c * BfFwd[1][id2 + i]);
// adjust velocity pert. according to new value
for (int j = 0; j < m_spacedim; ++j)
{
Fwd[1 + j][id2 + i] =
Fwd[1 + j][id2 + i] +
(VnNew - Vn[i]) * m_traceNormals[j][id2 + i];
}
}
// Copy boundary adjusted values into the boundary expansion
for (int i = 0; i < nVariables; ++i)
{
Vmath::Vcopy(nBCEdgePts,
&Fwd[i][id2], 1,
&(m_fields[i]
->GetBndCondExpansions()[bcRegion]
->UpdatePhys())[id1], 1);
}
}
}
/**
* @brief Wall boundary conditions for the APE equations.
*/
void APE::WhiteNoiseBC(int bcRegion,
int cnt,
Array<OneD, Array<OneD, NekDouble> > &Fwd,
Array<OneD, Array<OneD, NekDouble> > &BfFwd,
Array<OneD, Array<OneD, NekDouble> > &physarray)
{
int id1, id2, nBCEdgePts;
int nVariables = physarray.num_elements();
const Array<OneD, const int> &traceBndMap = m_fields[0]->GetTraceBndMap();
if (m_rng.count(bcRegion) == 0)
{
m_rng[bcRegion] = boost::mt19937(bcRegion);
}
ASSERTL0(
m_fields[0]->GetBndConditions()[bcRegion]->GetBoundaryConditionType() ==
SpatialDomains::eDirichlet,
"WhiteNoise BCs must be Dirichlet type BCs");
LibUtilities::Equation cond =
std::static_pointer_cast<SpatialDomains::DirichletBoundaryCondition>(
m_fields[0]->GetBndConditions()[bcRegion])
->m_dirichletCondition;
NekDouble sigma = cond.Evaluate();
ASSERTL0(sigma > NekConstants::kNekZeroTol,
"sigma must be greater than zero");
// random velocity perturbation
if (m_whiteNoiseBC_lastUpdate < m_time)
{
m_whiteNoiseBC_lastUpdate = m_time;
boost::normal_distribution<> dist(0, sigma);
m_whiteNoiseBC_p = dist(m_rng[bcRegion]);
}
int eMax = m_fields[0]->GetBndCondExpansions()[bcRegion]->GetExpSize();
for (int e = 0; e < eMax; ++e)
{
nBCEdgePts = m_fields[0]
->GetBndCondExpansions()[bcRegion]
->GetExp(e)
->GetTotPoints();
id1 = m_fields[0]->GetBndCondExpansions()[bcRegion]->GetPhys_Offset(e);
id2 = m_fields[0]->GetTrace()->GetPhys_Offset(traceBndMap[cnt + e]);
Array<OneD, Array<OneD, NekDouble> > tmp(nVariables);
for (int i = 0; i < nVariables; ++i)
{
tmp[i] = Array<OneD, NekDouble>(nBCEdgePts, 0.0);
}
// pressure perturbation
Vmath::Fill(nBCEdgePts, m_whiteNoiseBC_p, &tmp[0][0], 1);
// velocity perturbation
for (int i = 0; i < nBCEdgePts; ++i)
{
NekDouble u = m_whiteNoiseBC_p /
sqrt(m_gamma * BfFwd[0][id2 + i] * BfFwd[1][id2 + i]);
for (int j = 0; j < m_spacedim; ++j)
{
tmp[1 + j][i] = -1.0 * u * m_traceNormals[j][id2 + i];
}
}
// Copy boundary adjusted values into the boundary expansion
for (int i = 0; i < nVariables; ++i)
{
Vmath::Vcopy(nBCEdgePts,
&tmp[i][0], 1,
&(m_fields[i]
->GetBndCondExpansions()[bcRegion]
->UpdatePhys())[id1], 1);
}
}
}
/**
* @brief Compute the advection velocity in the standard space
* for each element of the expansion.
......@@ -602,13 +830,13 @@ const Array<OneD, const Array<OneD, NekDouble> > &APE::GetVecLocs()
/**
* @brief Get the baseflow field.
*/
const Array<OneD, const Array<OneD, NekDouble> > &APE::GetBasefield()
const Array<OneD, const Array<OneD, NekDouble> > &APE::GetBfTrace()
{
for (int i = 0; i < m_spacedim + 2; i++)
{
m_fields[0]->ExtractTracePhys(m_bf[i], m_traceBasefield[i]);
m_fields[0]->ExtractTracePhys(m_bf[i], m_bfTrace[i]);
}
return m_traceBasefield;
return m_bfTrace;
}
......
......@@ -37,6 +37,8 @@
#ifndef NEKTAR_SOLVERS_APESOLVER_EQUATIONSYSTEMS_APE_H
#define NEKTAR_SOLVERS_APESOLVER_EQUATIONSYSTEMS_APE_H
#include <boost/random/mersenne_twister.hpp>
#include <SolverUtils/UnsteadySystem.h>
#include <SolverUtils/AdvectionSystem.h>
#include <SolverUtils/Advection/Advection.h>
......@@ -73,7 +75,7 @@ class APE : public AdvectionSystem
SolverUtils::AdvectionSharedPtr m_advection;
std::vector<SolverUtils::ForcingSharedPtr> m_forcing;
SolverUtils::RiemannSolverSharedPtr m_riemannSolver;
Array<OneD, Array<OneD, NekDouble> > m_traceBasefield;
Array<OneD, Array<OneD, NekDouble> > m_bfTrace;
Array<OneD, Array<OneD, NekDouble> > m_vecLocs;
/// Isentropic coefficient, Ratio of specific heats (APE)
NekDouble m_gamma;
......@@ -111,15 +113,37 @@ class APE : public AdvectionSystem
const Array<OneD, const Array<OneD, NekDouble> > &GetVecLocs();
const Array<OneD, const Array<OneD, NekDouble> > &GetBasefield();
const Array<OneD, const Array<OneD, NekDouble> > &GetBfTrace();
NekDouble GetGamma();
private:
void SetBoundaryConditions(Array<OneD, Array<OneD, NekDouble> > &physarray, NekDouble time);
std::map<int, boost::mt19937> m_rng;
NekDouble m_whiteNoiseBC_lastUpdate;
NekDouble m_whiteNoiseBC_p;
void SetBoundaryConditions(Array<OneD, Array<OneD, NekDouble> > &physarray,
NekDouble time);
void WallBC(int bcRegion,
int cnt,
Array<OneD, Array<OneD, NekDouble> > &Fwd,
Array<OneD, Array<OneD, NekDouble> > &physarray);
void RiemannInvariantBC(int bcRegion,
int cnt,
Array<OneD, Array<OneD, NekDouble> > &Fwd,
Array<OneD, Array<OneD, NekDouble> > &BfFwd,
Array<OneD, Array<OneD, NekDouble> > &physarray);
void WallBC(int bcRegion, int cnt, Array<OneD, Array<OneD, NekDouble> > &Fwd, Array<OneD, Array<OneD, NekDouble> > &physarray);
void WhiteNoiseBC(int bcRegion,
int cnt,
Array<OneD, Array<OneD, NekDouble> > &Fwd,
Array<OneD, Array<OneD, NekDouble> > &BfFwd,
Array<OneD, Array<OneD, NekDouble> > &physarray);
};
}
......
......@@ -8,14 +8,14 @@
</files>
<metrics>
<metric type="L2" id="1">
<value variable="p" tolerance="1e-12">0.102091</value>
<value variable="u" tolerance="1e-12">0.000248662</value>
<value variable="v" tolerance="1e-12">0</value>
<value variable="p" tolerance="1e-12">12.8017</value>
<value variable="u" tolerance="1e-12">0.0314268</value>
<value variable="v" tolerance="1e-12">0.00291624</value>
</metric>
<metric type="Linf" id="2">
<value variable="p" tolerance="1e-12">0.911223</value>
<value variable="u" tolerance="1e-12">0.00221946</value>
<value variable="v" tolerance="1e-12">3.54604e-17</value>
<value variable="p" tolerance="1e-12">57.9942</value>
<value variable="u" tolerance="1e-12">0.180743</value>
<value variable="v" tolerance="1e-12">0.174942</value>
</metric>
</metrics>
</test>
......@@ -780,14 +780,14 @@
</BOUNDARYREGIONS>
<BOUNDARYCONDITIONS>
<REGION REF="0"> <!-- Excitation at the inflow boundary -->
<D VAR="p" USERDEFINEDTYPE="TimeDependent" VALUE="100*(sin(Omega*t))^2" />
<D VAR="u" USERDEFINEDTYPE="TimeDependent" VALUE="100*((sin(Omega*t))^2)/Rho0/((Gamma*Pinfinity/Rho0)^0.5)" />
<D VAR="v" VALUE="0" />
<D VAR="p" USERDEFINEDTYPE="WhiteNoise" VALUE="100" />
<D VAR="u" USERDEFINEDTYPE="WhiteNoise" VALUE="100" />
<D VAR="v" USERDEFINEDTYPE="WhiteNoise" VALUE="100" />
</REGION>
<REGION REF="1"> <!-- Problem homogeneous in y-direction -->
<N VAR="p" VALUE="0" />
<N VAR="u" VALUE="0" />
<N VAR="v" VALUE="0" />
<D VAR="p" USERDEFINEDTYPE="RiemannInvariantBC"/>
<D VAR="u" USERDEFINEDTYPE="RiemannInvariantBC"/>
<D VAR="v" USERDEFINEDTYPE="RiemannInvariantBC"/>
</REGION>
</BOUNDARYCONDITIONS>
<FUNCTION NAME="Baseflow"> <!-- Incompressible base flow -->
......
......@@ -336,9 +336,9 @@
</BOUNDARYREGIONS>
<BOUNDARYCONDITIONS>
<REGION REF="0">
<D VAR="p" VALUE="0"/>
<D VAR="u" VALUE="0"/>
<D VAR="v" VALUE="0"/>
<D VAR="p" USERDEFINEDTYPE="RiemannInvariantBC"/>
<D VAR="u" USERDEFINEDTYPE="RiemannInvariantBC"/>
<D VAR="v" USERDEFINEDTYPE="RiemannInvariantBC"/>
</REGION>
</BOUNDARYCONDITIONS>
<FUNCTION NAME="Baseflow"> <!-- Incompressible base flow -->
......
......@@ -336,9 +336,9 @@
</BOUNDARYREGIONS>
<BOUNDARYCONDITIONS>
<REGION REF="0">
<D VAR="p" VALUE="0"/>
<D VAR="u" VALUE="0"/>
<D VAR="v" VALUE="0"/>
<D VAR="p" USERDEFINEDTYPE="RiemannInvariantBC"/>
<D VAR="u" USERDEFINEDTYPE="RiemannInvariantBC"/>
<D VAR="v" USERDEFINEDTYPE="RiemannInvariantBC"/>
</REGION>
</BOUNDARYCONDITIONS>
<FUNCTION NAME="Baseflow"> <!-- Incompressible base flow -->
......
......@@ -8,14 +8,14 @@
</files>
<metrics>
<metric type="L2" id="1">
<value variable="p" tolerance="1e-4">145.343</value>
<value variable="u" tolerance="1e-7">0.329227</value>
<value variable="v" tolerance="1e-7">0</value>
<value variable="p" tolerance="1e-4">52.0346</value>
<value variable="u" tolerance="1e-7">0.228765</value>
<value variable="v" tolerance="1e-7">0.0598651</value>
</metric>
<metric type="Linf" id="2">
<value variable="p" tolerance="1e-4">119.917</value>
<value variable="u" tolerance="1e-7">0.260865</value>
<value variable="v" tolerance="1e-7">0</value>
<value variable="p" tolerance="1e-4">60.1152</value>
<value variable="u" tolerance="1e-7">0.246703</value>
<value variable="v" tolerance="1e-7">0.109342</value>
</metric>
</metrics>
</test>
......@@ -119,9 +119,9 @@
<D VAR="v" VALUE="0" />
</REGION>
<REGION REF="1"> <!-- Problem homogeneous in y-direction -->
<N VAR="p" VALUE="0" />
<N VAR="u" VALUE="0" />
<N VAR="v" VALUE="0" />
<D VAR="p" USERDEFINEDTYPE="RiemannInvariantBC"/>
<D VAR="u" USERDEFINEDTYPE="RiemannInvariantBC"/>
<D VAR="v" USERDEFINEDTYPE="RiemannInvariantBC"/>
</REGION>
</BOUNDARYCONDITIONS>
<FUNCTION NAME="Baseflow"> <!-- Incompressible base flow -->
......
......@@ -1386,10 +1386,10 @@
</BOUNDARYREGIONS>
<BOUNDARYCONDITIONS>
<REGION REF="0">
<D VAR="p" VALUE="0"/>
<D VAR="u" VALUE="0"/>
<D VAR="v" VALUE="0"/>
<D VAR="w" VALUE="0"/>
<D VAR="p" USERDEFINEDTYPE="RiemannInvariantBC"/>
<D VAR="u" USERDEFINEDTYPE="RiemannInvariantBC"/>
<D VAR="v" USERDEFINEDTYPE="RiemannInvariantBC"/>
<D VAR="w" USERDEFINEDTYPE="RiemannInvariantBC"/>
</REGION>
<REGION REF="1">
<D VAR="p" USERDEFINEDTYPE="Wall" VALUE="0" />
......
......@@ -1385,10 +1385,10 @@
</BOUNDARYREGIONS>
<BOUNDARYCONDITIONS>