Commit c37bfdcf authored by Douglas Serson's avatar Douglas Serson

Move shock capture to separate classes

parent 378a0eb0
......@@ -115,11 +115,6 @@ namespace Nektar
m_fluxVectorNS = fluxVector;
}
inline void SetRiemannSolver(RiemannSolverSharedPtr riemann)
{
m_riemann = riemann;
}
inline void SetHomoDerivs(Array<OneD, Array<OneD, NekDouble> > &deriv)
{
v_SetHomoDerivs(deriv);
......@@ -133,7 +128,6 @@ namespace Nektar
protected:
DiffusionFluxVecCB m_fluxVector;
DiffusionFluxVecCBNS m_fluxVectorNS;
RiemannSolverSharedPtr m_riemann;
DiffusionArtificialDiffusion m_ArtificialDiffusionVector;
virtual void v_InitObject(
......
......@@ -115,28 +115,9 @@ namespace Nektar
m_planeCounter = 0;
m_planeDiff->SetFluxVectorNS(m_fluxVectorNS);
if (m_riemann)
{
// Set Riemann solver and flux vector callback for this plane.
m_planeDiff->SetRiemannSolver(m_riemann);
// Override Riemann solver scalar and vector callbacks.
map<string, RSScalarFuncType>::iterator it1;
map<string, RSScalarFuncType> scalars = m_riemann->GetScalars();
for (it1 = scalars.begin(); it1 != scalars.end(); ++it1)
{
boost::shared_ptr<HomoRSScalar> tmp =
MemoryManager<HomoRSScalar>
::AllocateSharedPtr(it1->second, m_numPlanes);
m_riemann->SetScalar(it1->first, &HomoRSScalar::Exec, tmp);
}
}
m_fieldsPlane = Array<OneD, MultiRegions::ExpListSharedPtr>
(nConvectiveFields);
if (m_fluxVectorNS)
{
m_inarrayPlane = Array<OneD, Array<OneD, NekDouble> >
......
......@@ -228,7 +228,7 @@ namespace Nektar
int i, j;
int nTracePts = fields[0]->GetTrace()->GetTotPoints();
int nvariables = fields.num_elements();
int nDim = uflux.num_elements();
int nDim = fields[0]->GetCoordim(0);;
Array<OneD, NekDouble > Fwd (nTracePts);
Array<OneD, NekDouble > Bwd (nTracePts);
......
///////////////////////////////////////////////////////////////////////////////
//
// File: ArtificialDiffusion.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).
//
// License for the specific language governing rights and limitations under
// 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: Abstract base class for compressible solver artificial diffusion
// used for shock capturing artificial diffusion.
//
///////////////////////////////////////////////////////////////////////////////
#include "ArtificialDiffusion.h"
using namespace std;
namespace Nektar
{
ArtificialDiffusionFactory& GetArtificialDiffusionFactory()
{
typedef Loki::SingletonHolder<ArtificialDiffusionFactory,
Loki::CreateUsingNew,
Loki::NoDestroy,
Loki::SingleThreaded> Type;
return Type::Instance();
}
ArtificialDiffusion::ArtificialDiffusion(
const LibUtilities::SessionReaderSharedPtr& pSession,
const Array<OneD, MultiRegions::ExpListSharedPtr>& pFields,
const int spacedim)
: m_session(pSession),
m_fields(pFields)
{
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 ("Skappa", m_Skappa, -2.048);
m_session->LoadParameter ("Kappa", m_Kappa, 0.0);
// Create auxiliary object to convert variables
m_varConv = MemoryManager<VariableConverter>::AllocateSharedPtr(
m_session, spacedim);
m_diffusion = SolverUtils::GetDiffusionFactory()
.CreateInstance("LDG", "LDG");
m_diffusion->SetArtificialDiffusionVector(
&ArtificialDiffusion::GetArtificialViscosity, this);
m_diffusion->InitObject (m_session, m_fields);
}
/**
*
*/
void ArtificialDiffusion::DoArtificialDiffusion(
const Array<OneD, const Array<OneD, NekDouble> > &inarray,
Array<OneD, Array<OneD, NekDouble> > &outarray)
{
v_DoArtificialDiffusion(inarray, outarray);
}
/**
*
*/
void ArtificialDiffusion::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);
for (i = 0; i < nvariables; ++i)
{
Vmath::Vadd(npoints,
outarray[i], 1,
outarrayDiff[i], 1,
outarray[i], 1);
}
}
/**
*
*/
void ArtificialDiffusion::GetArtificialViscosity(
const Array<OneD, Array<OneD, NekDouble> > &physfield,
Array<OneD, NekDouble > &mu)
{
v_GetArtificialViscosity(physfield, mu);
}
}
///////////////////////////////////////////////////////////////////////////////
//
// File: ArtificialDiffusion.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).
//
// License for the specific language governing rights and limitations under
// 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: Abstract base class for compressible solver artificial diffusion
// used for shock capturing artificial diffusion.
//
///////////////////////////////////////////////////////////////////////////////
#ifndef NEKTAR_SOLVERS_COMPRESSIBLEFLOWSOLVER_ARTIFICIALDIFFUSION_BASE
#define NEKTAR_SOLVERS_COMPRESSIBLEFLOWSOLVER_ARTIFICIALDIFFUSION_BASE
#include <string>
#include <CompressibleFlowSolver/Misc/VariableConverter.h>
#include <LibUtilities/BasicUtils/NekFactory.hpp>
#include <LibUtilities/BasicUtils/SharedArray.hpp>
#include <MultiRegions/ExpList.h>
#include <SolverUtils/Diffusion/Diffusion.h>
namespace Nektar
{
// Forward declaration
class ArtificialDiffusion;
/// A shared pointer to a boundary condition object
typedef boost::shared_ptr<ArtificialDiffusion> ArtificialDiffusionSharedPtr;
/// Declaration of the boundary condition factory
typedef LibUtilities::NekFactory<std::string, ArtificialDiffusion,
const LibUtilities::SessionReaderSharedPtr&,
const Array<OneD, MultiRegions::ExpListSharedPtr>&,
const int > ArtificialDiffusionFactory;
/// Declaration of the boundary condition factory singleton
ArtificialDiffusionFactory& GetArtificialDiffusionFactory();
/**
* @class ArtificialDiffusion
* @brief Encapsulates the artificial diffusion used in shock capture
*/
class ArtificialDiffusion
{
public:
virtual ~ArtificialDiffusion() {}
/// Apply the boundary condition
void DoArtificialDiffusion(
const Array<OneD, const Array<OneD, NekDouble> > &inarray,
Array<OneD, Array<OneD, NekDouble> > &outarray);
void GetArtificialViscosity(
const Array<OneD, Array<OneD, NekDouble> > &physfield,
Array<OneD, NekDouble > &mu);
protected:
/// Session reader
LibUtilities::SessionReaderSharedPtr m_session;
/// Array of fields
Array<OneD, MultiRegions::ExpListSharedPtr> m_fields;
/// Auxiliary object to convert variables
VariableConverterSharedPtr m_varConv;
/// LDG Diffusion operator
SolverUtils::DiffusionSharedPtr m_diffusion;
/// Parameters
NekDouble m_FacL;
NekDouble m_FacH;
NekDouble m_hFactor;
NekDouble m_C1;
NekDouble m_C2;
NekDouble m_mu0;
NekDouble m_Skappa;
NekDouble m_Kappa;
/// Constructor
ArtificialDiffusion(const LibUtilities::SessionReaderSharedPtr& pSession,
const Array<OneD, MultiRegions::ExpListSharedPtr>& pFields,
const int spacedim);
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)=0;
};
}
#endif
///////////////////////////////////////////////////////////////////////////////
//
// File: NonSmoothShockCapture.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).
//
// License for the specific language governing rights and limitations under
// 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: NonSmooth artificial diffusion for shock capture
//
///////////////////////////////////////////////////////////////////////////////
#include "NonSmoothShockCapture.h"
using namespace std;
namespace Nektar
{
std::string NonSmoothShockCapture::className = GetArtificialDiffusionFactory().
RegisterCreatorFunction("NonSmooth",
NonSmoothShockCapture::create,
"NonSmooth artificial diffusion for shock capture.");
NonSmoothShockCapture::NonSmoothShockCapture(
const LibUtilities::SessionReaderSharedPtr& pSession,
const Array<OneD, MultiRegions::ExpListSharedPtr>& pFields,
const int spacedim)
: ArtificialDiffusion(pSession, pFields, spacedim)
{
}
void NonSmoothShockCapture::v_GetArtificialViscosity(
const Array<OneD, Array<OneD, NekDouble> > &physfield,
Array<OneD, NekDouble > &mu)
{
const int nElements = m_fields[0]->GetExpSize();
int PointCount = 0;
int nTotQuadPoints = m_fields[0]->GetTotPoints();
Array<OneD, NekDouble> Sensor (nTotQuadPoints, 0.0);
Array<OneD, NekDouble> SensorKappa(nTotQuadPoints, 0.0);
m_varConv->GetSensor(m_fields[0], physfield, Sensor, SensorKappa);
for (int e = 0; e < nElements; e++)
{
// Threshold value specified in C. Biottos thesis. Based on a 1D
// shock tube problem S_k = log10(1/p^4). See G.E. Barter and
// D.L. Darmofal. Shock Capturing with PDE-based artificial
// diffusion for DGFEM: Part 1 Formulation, Journal of Computational
// Physics 229 (2010) 1810-1827 for further reference
// Adjustable depending on the coarsness of the mesh. Might want to
// move this variable into the session file
int nQuadPointsElement = m_fields[0]->GetExp(e)->GetTotPoints();
for (int n = 0; n < nQuadPointsElement; n++)
{
NekDouble mu_0 = m_mu0;
if (Sensor[n+PointCount] < (m_Skappa-m_Kappa))
{
mu[n+PointCount] = 0;
}
else if (Sensor[n+PointCount] >= (m_Skappa-m_Kappa)
&& Sensor[n+PointCount] <= (m_Skappa+m_Kappa))
{
mu[n+PointCount] = mu_0 * (0.5 * (1 + sin(
M_PI * (Sensor[n+PointCount] -
m_Skappa - m_Kappa) /
(2*m_Kappa))));
}
else if (Sensor[n+PointCount] > (m_Skappa+m_Kappa))
{
mu[n+PointCount] = mu_0;
}
}
PointCount += nQuadPointsElement;
}
}
}
///////////////////////////////////////////////////////////////////////////////
//
// File: NonSmoothShockCapture.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).
//
// License for the specific language governing rights and limitations under
// 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: NonSmooth artificial diffusion for shock capture
//
///////////////////////////////////////////////////////////////////////////////
#ifndef NEKTAR_SOLVERS_COMPRESSIBLEFLOWSOLVER_ARTIFICIALDIFFUSION_NONSMOOTH
#define NEKTAR_SOLVERS_COMPRESSIBLEFLOWSOLVER_ARTIFICIALDIFFUSION_NONSMOOTH
#include "ArtificialDiffusion.h"
namespace Nektar
{
/**
* @brief Wall boundary conditions for compressible flow problems.
*/
class NonSmoothShockCapture : public ArtificialDiffusion
{
public:
friend class MemoryManager<NonSmoothShockCapture>;
/// 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<NonSmoothShockCapture>::
AllocateSharedPtr(pSession, pFields, spacedim);
return p;
}
///Name of the class
static std::string className;
protected:
virtual void v_GetArtificialViscosity(
const Array<OneD, Array<OneD, NekDouble> > &physfield,
Array<OneD, NekDouble > &mu);
private:
NonSmoothShockCapture(const LibUtilities::SessionReaderSharedPtr& pSession,
const Array<OneD, MultiRegions::ExpListSharedPtr>& pFields,
const int spacedim);
virtual ~NonSmoothShockCapture(void){};
};
}
#endif
///////////////////////////////////////////////////////////////////////////////
//
// 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).
//
// License for the specific language governing rights and limitations under
// 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.");
}
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->GetPressure(inarray, tmp);
m_varConv->GetSoundSpeed(inarray, tmp, 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();
int nPts = m_fields[0]->GetTotPoints();
Array<OneD, NekDouble > sensor (nPts, 0.0);
Array<OneD, NekDouble > SensorKappa(nPts, 0.0);
// Calculate sensor
m_varConv->GetSensor(m_fields[0], physfield, sensor, SensorKappa);
NekDouble ThetaH = m_FacH;
NekDouble ThetaL = m_FacL;
NekDouble Phi0 = (ThetaH+ThetaL)/2;
NekDouble DeltaPhi = ThetaH-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,