Commit be533508 authored by Giacomo Castiglioni's avatar Giacomo Castiglioni

removed smooth AV!

parent d6e40656
......@@ -444,8 +444,7 @@ Under the two following sections it is possible to define the initial conditions
\section{Examples}
\subsection{Shock capturing}
Compressible flows can be characterised by abrupt changes in density within the flow domain often referred to as shocks. These discontinuities lead to numerical instabilities (Gibbs phenomena). This problem is prevented by locally adding a diffusion term to the equations to damp the numerical oscillations.
An artificial diffusion term is introduced locally to the Euler equations to deal with flow discontinuity and the consequential numerical oscillations. Two models are implemented, a non-smooth and a smooth artificial viscosity model.
Compressible flows can be characterised by abrupt changes in density within the flow domain often referred to as shocks. These discontinuities can lead to numerical instabilities (Gibbs phenomena). This problem is prevented by locally adding a diffusion term to the equations to damp the numerical oscillations.
\subsubsection{Non-smooth artificial viscosity model}
For the non-smooth artificial viscosity model the added artificial viscosity is constant in each element and discontinuous between the elements. The Euler system is augmented by an added Laplacian term on right hand side of equation \ref{eq:euler} \cite{persson2006sub}.
......@@ -480,9 +479,9 @@ where $s_0 = s_\kappa - 4.25\;log_{10}(p)$.
To enable the non-smooth viscosity model, the following line has to be added to the \inltt{SOLVERINFO} section:
\begin{lstlisting}[style=XmlStyle]
<SOLVERINFO>
<SOLVERINFO>
<I PROPERTY="ShockCaptureType" VALUE="NonSmooth" />
<SOLVERINFO>
<SOLVERINFO>
\end{lstlisting}
The diffusivity and the sensor can be controlled by the following parameters:
\begin{lstlisting}[style=XmlStyle]
......@@ -502,45 +501,6 @@ The diffusivity and the sensor can be controlled by the following parameters:
\end{center}
\end{figure}
\subsubsection{Smooth artificial viscosity model}
For the smooth artificial viscosity model an extra PDE for the artificial viscosity is appended to the Euler system
\begin{equation}\label{eq:eulerplusvis}\begin{split}
\frac{\partial \epsilon}{\partial t} &= \nabla\cdot \left(\nabla \epsilon\right) + \frac{1}{\tau}\left(\frac{h}{p}\lambda_{max}S_\kappa - \epsilon\right)\ \ \ \ \ \ \textrm{on}\ \ \ \ \ \ \ \Omega\\
\frac{\partial \epsilon}{\partial n} &= 0\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \textrm{on}\ \ \ \ \ \ \ \Gamma
\end{split}
\end{equation}
where $S_\kappa$ is a normalised sensor value and serves as a forcing term for the artificial viscosity. A smooth artificial viscosity distribution is obtained.\\
\\
To enable the smooth viscosity model, the following line has to be added to the \inltt{SOLVERINFO} section:
\begin{lstlisting}[style=XmlStyle]
<SOLVERINFO>
<I PROPERTY="ShockCaptureType" VALUE="Smooth" />
<SOLVERINFO>
\end{lstlisting}
Furthermore, the extra viscosity variable \inltt{eps} has to be added to the variable list:
\begin{lstlisting}[style=XmlStyle]
<VARIABLES>
<V ID="0"> rho </V>
<V ID="1"> rhou </V>
<V ID="2"> rhov </V>
<V ID="4"> E </V>
<V ID="5"> eps </V>
</VARIABLES>
\end{lstlisting}
A similar addition has to be made for the boundary conditions and initial conditions. The tests that have been run started with a uniform homogeneous boundary condition and initial condition.
The following parameters can be set in the xml session file:
\begin{lstlisting}[style=XmlStyle]
<PARAMETERS>
<P> Skappa = -1.3 </P>
<P> Kappa = 0.2 </P>
<P> mu0 = 1.0 </P>
<P> FH = 3 </P>
<P> FL = 0.01*FH </P>
<P> C1 = 0.03 </P>
<P> C2 = 5/3*C1 </P>
</PARAMETERS>
\end{lstlisting}
where for now \inltt{FH} and \inltt{FL} are used to tune which range of the sensor is used as a forcing term and \inltt{C1} and \inltt{C2} are fixed constants which can be played around with to make the model more diffusive or not. However these constants are generally fixed.
\subsection{Variable polynomial order}
A sensor based $p$-adaptive algorithm is implemented to optimise the computational cost and accuracy.
The DG scheme allows one to use different polynomial orders since the fluxes over the elements are determined using a Riemann solver and there is now further coupling between the elements. Furthermore, the initial $p$-adaptive algorithm uses the same sensor as the shock capturing algorithm to identify the smoothness of the local solution so it rather straightforward to implement both algorithms at the same time.\\
......
......@@ -130,11 +130,6 @@ namespace Nektar
int numConvFields = nConvectiveFields;
if (m_shockCaptureType == "Smooth")
{
numConvFields = nConvectiveFields - 1;
}
for (j = 0; j < nDim; ++j)
{
for (i = 0; i < numConvFields; ++i)
......
///////////////////////////////////////////////////////////////////////////////
//
// File: SmoothShockCapture.cpp
//
// For more information, please see: http://www.nektar.info
//
// The MIT License
//
// Copyright (c) 2006 Division of Applied Mathematics, Brown University (USA),
// Department of Aeronautics, Imperial College London (UK), and Scientific
// Computing and Imaging Institute, University of Utah (USA).
//
// Permission is hereby granted, free of charge, to any person obtaining a
// copy of this software and associated documentation files (the "Software"),
// to deal in the Software without restriction, including without limitation
// the rights to use, copy, modify, merge, publish, distribute, sublicense,
// and/or sell copies of the Software, and to permit persons to whom the
// Software is furnished to do so, subject to the following conditions:
//
// The above copyright notice and this permission notice shall be included
// in all copies or substantial portions of the Software.
//
// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
// OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
// THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
// DEALINGS IN THE SOFTWARE.
//
// Description: Smooth artificial diffusion for shock capture
//
///////////////////////////////////////////////////////////////////////////////
#include "SmoothShockCapture.h"
using namespace std;
namespace Nektar
{
std::string SmoothShockCapture::className = GetArtificialDiffusionFactory().
RegisterCreatorFunction("Smooth",
SmoothShockCapture::create,
"Smooth artificial diffusion for shock capture.");
SmoothShockCapture::SmoothShockCapture(
const LibUtilities::SessionReaderSharedPtr& pSession,
const Array<OneD, MultiRegions::ExpListSharedPtr>& pFields,
const int spacedim)
: ArtificialDiffusion(pSession, pFields, spacedim)
{
ASSERTL0(m_fields.num_elements() == spacedim + 3,
"Not enough variables for smooth shock capturing; "
"make sure you have added eps to variable list.");
m_session->LoadParameter ("FL", m_FacL, 0.0);
m_session->LoadParameter ("FH", m_FacH, 0.0);
m_session->LoadParameter ("hFactor", m_hFactor, 1.0);
m_session->LoadParameter ("C1", m_C1, 3.0);
m_session->LoadParameter ("C2", m_C2, 5.0);
m_session->LoadParameter ("mu0", m_mu0, 1.0);
m_session->LoadParameter ("SensorOffset", m_offset, 1);
ROOTONLY_NEKERROR(Nektar::ErrorUtil::ewarning,
"h/p Lambda scaling not implemented for SmoothShockCapture."
"There seems to be something wrong with the boundary conditions "
"as well.")
}
void SmoothShockCapture::v_DoArtificialDiffusion(
const Array<OneD, const Array<OneD, NekDouble> > &inarray,
Array<OneD, Array<OneD, NekDouble> > &outarray)
{
int i;
int nvariables = inarray.num_elements();
int npoints = m_fields[0]->GetNpoints();
Array<OneD, Array<OneD, NekDouble> > outarrayDiff(nvariables);
for (i = 0; i < nvariables; ++i)
{
outarrayDiff[i] = Array<OneD, NekDouble>(npoints, 0.0);
}
m_diffusion->Diffuse(nvariables, m_fields, inarray, outarrayDiff);
const Array<OneD, int> ExpOrder = m_fields[0]->EvalBasisNumModesMaxPerExp();
NekDouble pOrder = Vmath::Vmax(ExpOrder.num_elements(), ExpOrder, 1);
Array <OneD, NekDouble > a_vel (npoints, 0.0);
Array <OneD, NekDouble > u_abs (npoints, 0.0);
Array <OneD, NekDouble > tmp (npoints, 0.0);
m_varConv->GetSoundSpeed(inarray, a_vel);
m_varConv->GetAbsoluteVelocity(inarray, u_abs);
Vmath::Vadd(npoints, a_vel, 1, u_abs, 1, tmp, 1);
NekDouble max_wave_sp = Vmath::Vmax(npoints, tmp, 1);
NekDouble fac = m_C2*max_wave_sp*pOrder;
Vmath::Smul(npoints,
fac,
outarrayDiff[nvariables-1], 1,
outarrayDiff[nvariables-1], 1);
for (i = 0; i < nvariables; ++i)
{
Vmath::Vadd(npoints,
outarray[i], 1,
outarrayDiff[i], 1,
outarray[i], 1);
}
Array<OneD, Array<OneD, NekDouble> > outarrayForcing(nvariables);
for (i = 0; i < nvariables; ++i)
{
outarrayForcing[i] = Array<OneD, NekDouble>(npoints, 0.0);
}
GetForcingTerm(inarray, outarrayForcing);
for (i = 0; i < nvariables; ++i)
{
// Add Forcing Term
Vmath::Vadd(npoints,
outarray[i], 1,
outarrayForcing[i], 1,
outarray[i], 1);
}
}
void SmoothShockCapture::v_GetArtificialViscosity(
const Array<OneD, Array<OneD, NekDouble> > &physfield,
Array<OneD, NekDouble > &mu)
{
int nvariables = physfield.num_elements();
NekDouble Phi0 = (m_FacH+m_FacL)/2;
NekDouble DeltaPhi = m_FacH-Phi0;
for (int e = 0; e < mu.num_elements(); e++)
{
if (physfield[nvariables-1][e] <= (Phi0 - DeltaPhi))
{
mu[e] = 0;
}
else if(physfield[nvariables-1][e] >= (Phi0 + DeltaPhi))
{
mu[e] = m_mu0;
}
else if(abs(physfield[nvariables-1][e]-Phi0) < DeltaPhi)
{
mu[e] = m_mu0/2*(1+sin(M_PI*
(physfield[nvariables-1][e]-Phi0)/(2*DeltaPhi)));
}
}
}
void SmoothShockCapture::GetForcingTerm(
const Array<OneD, const Array<OneD, NekDouble> > &inarray,
Array<OneD, Array<OneD, NekDouble> > outarrayForcing)
{
const int nPts = m_fields[0]->GetTotPoints();
const int nvariables = m_fields.num_elements();
const int nElements = m_fields[0]->GetExpSize();
Array<OneD, NekDouble> Sensor(nPts, 0.0);
Array<OneD, NekDouble> SensorKappa(nPts, 0.0);
Array <OneD, NekDouble > Lambda(nPts, 0.0);
Array <OneD, NekDouble > soundspeed(nPts, 0.0);
Array <OneD, NekDouble > absVelocity(nPts, 0.0);
NekDouble tau, pOrder;
Array<OneD,int> pOrderElmt = m_fields[0]->EvalBasisNumModesMaxPerExp();
// Thermodynamic related quantities
m_varConv->GetSoundSpeed(inarray, soundspeed);
m_varConv->GetAbsoluteVelocity(inarray, absVelocity);
m_varConv->GetSensor(m_fields[0], inarray, Sensor, SensorKappa, m_offset);
// Determine the maximum wavespeed
Vmath::Vadd(nPts, absVelocity, 1, soundspeed, 1, Lambda, 1);
NekDouble LambdaMax = Vmath::Vmax(nPts, Lambda, 1);
int PointCount = 0;
for (int e = 0; e < nElements; e++)
{
int nQuadPointsElement = m_fields[0]->GetExp(e)->GetTotPoints();
for (int n = 0; n < nQuadPointsElement; n++)
{
pOrder = pOrderElmt[e];
// order 1.0e-06
tau = 1.0 / (m_C1*pOrder*LambdaMax);
outarrayForcing[nvariables-1][n + PointCount] =
1 / tau * (m_hFactor * LambdaMax /
pOrder *
SensorKappa[n + PointCount] -
inarray[nvariables-1][n + PointCount]);
}
PointCount += nQuadPointsElement;
}
}
}
///////////////////////////////////////////////////////////////////////////////
//
// File: SmoothShockCapture.h
//
// For more information, please see: http://www.nektar.info
//
// The MIT License
//
// Copyright (c) 2006 Division of Applied Mathematics, Brown University (USA),
// Department of Aeronautics, Imperial College London (UK), and Scientific
// Computing and Imaging Institute, University of Utah (USA).
//
// Permission is hereby granted, free of charge, to any person obtaining a
// copy of this software and associated documentation files (the "Software"),
// to deal in the Software without restriction, including without limitation
// the rights to use, copy, modify, merge, publish, distribute, sublicense,
// and/or sell copies of the Software, and to permit persons to whom the
// Software is furnished to do so, subject to the following conditions:
//
// The above copyright notice and this permission notice shall be included
// in all copies or substantial portions of the Software.
//
// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
// OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
// THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
// DEALINGS IN THE SOFTWARE.
//
// Description: Smooth artificial diffusion for shock capture
//
///////////////////////////////////////////////////////////////////////////////
#ifndef NEKTAR_SOLVERS_COMPRESSIBLEFLOWSOLVER_ARTIFICIALDIFFUSION_SMOOTH
#define NEKTAR_SOLVERS_COMPRESSIBLEFLOWSOLVER_ARTIFICIALDIFFUSION_SMOOTH
#include "ArtificialDiffusion.h"
namespace Nektar
{
/**
* @brief Smooth artificial diffusion for shock capture for compressible flow
* problems.
*/
class SmoothShockCapture : public ArtificialDiffusion
{
public:
friend class MemoryManager<SmoothShockCapture>;
/// Creates an instance of this class
static ArtificialDiffusionSharedPtr create(
const LibUtilities::SessionReaderSharedPtr& pSession,
const Array<OneD, MultiRegions::ExpListSharedPtr>& pFields,
const int spacedim)
{
ArtificialDiffusionSharedPtr p = MemoryManager<SmoothShockCapture>::
AllocateSharedPtr(pSession, pFields, spacedim);
return p;
}
///Name of the class
static std::string className;
protected:
virtual void v_DoArtificialDiffusion(
const Array<OneD, const Array<OneD, NekDouble> > &inarray,
Array<OneD, Array<OneD, NekDouble> > &outarray);
virtual void v_GetArtificialViscosity(
const Array<OneD, Array<OneD, NekDouble> > &physfield,
Array<OneD, NekDouble > &mu);
private:
SmoothShockCapture(const LibUtilities::SessionReaderSharedPtr& pSession,
const Array<OneD, MultiRegions::ExpListSharedPtr>& pFields,
const int spacedim);
virtual ~SmoothShockCapture(void){};
void GetForcingTerm(
const Array<OneD, const Array<OneD, NekDouble> > &inarray,
Array<OneD, Array<OneD, NekDouble> > outarrayForcing);
/// Parameters
NekDouble m_FacL;
NekDouble m_FacH;
NekDouble m_hFactor;
NekDouble m_C1;
NekDouble m_C2;
NekDouble m_mu0;
int m_offset;
};
}
#endif
......@@ -10,7 +10,6 @@ IF( NEKTAR_SOLVER_COMPRESSIBLE_FLOW )
SET(CompressibleFlowSolverSource
./ArtificialDiffusion/ArtificialDiffusion.cpp
./ArtificialDiffusion/NonSmoothShockCapture.cpp
./ArtificialDiffusion/SmoothShockCapture.cpp
./BoundaryConditions/CFSBndCond.cpp
./BoundaryConditions/ExtrapOrder0BC.cpp
./BoundaryConditions/IsentropicVortexBC.cpp
......
Markdown is supported
0%
or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment