Commit 6bc3c14d authored by Douglas Serson's avatar Douglas Serson
Browse files

Tidy artificial diffusion

parent 6d337ab7
......@@ -139,8 +139,8 @@ namespace Nektar
{
Array<OneD, NekDouble> muvar(nPts, 0.0);
m_ArtificialDiffusionVector(inarray, muvar);
int numConvFields = nConvectiveFields;
int numConvFields = nConvectiveFields;
if (m_shockCaptureType == "Smooth")
{
......
......@@ -54,18 +54,8 @@ ArtificialDiffusion::ArtificialDiffusion(
const Array<OneD, MultiRegions::ExpListSharedPtr>& pFields,
const int spacedim)
: m_session(pSession),
m_fields(pFields)
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);
m_session->LoadParameter ("SensorOffset", m_offset, 1);
// Create auxiliary object to convert variables
m_varConv = MemoryManager<VariableConverter>::AllocateSharedPtr(
m_session, spacedim);
......
......@@ -91,19 +91,9 @@ class ArtificialDiffusion
/// 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;
int m_offset;
/// Constructor
ArtificialDiffusion(const LibUtilities::SessionReaderSharedPtr& pSession,
ArtificialDiffusion(
const LibUtilities::SessionReaderSharedPtr& pSession,
const Array<OneD, MultiRegions::ExpListSharedPtr>& pFields,
const int spacedim);
......
......@@ -42,8 +42,8 @@ namespace Nektar
std::string NonSmoothShockCapture::className = GetArtificialDiffusionFactory().
RegisterCreatorFunction("NonSmooth",
NonSmoothShockCapture::create,
"NonSmooth artificial diffusion for shock capture.");
NonSmoothShockCapture::create,
"NonSmooth artificial diffusion for shock capture.");
NonSmoothShockCapture::NonSmoothShockCapture(
const LibUtilities::SessionReaderSharedPtr& pSession,
......@@ -51,60 +51,17 @@ NonSmoothShockCapture::NonSmoothShockCapture(
const int spacedim)
: ArtificialDiffusion(pSession, pFields, spacedim)
{
m_session->LoadParameter ("SensorOffset", m_offset, 1);
}
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();
int nTotPoints = 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, m_offset);
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 (m_fields[0]->GetExp(e)->GetNcoeffs() <= m_offset)
{
mu[n+PointCount] = 0;
}
else if (Sensor[n+PointCount] < (m_Skappa-m_Kappa))
{
mu[n+PointCount] = 0;
}
else if (Sensor[n+PointCount] > (m_Skappa+m_Kappa))
{
mu[n+PointCount] = mu_0;
}
else
{
mu[n+PointCount] = mu_0 * (0.5 * (1 + sin(
M_PI * (Sensor[n+PointCount] -
m_Skappa) / (2*m_Kappa))));
}
}
PointCount += nQuadPointsElement;
}
Array<OneD, NekDouble> Sensor (nTotPoints, 0.0);
m_varConv->GetSensor(m_fields[0], physfield, Sensor, mu, m_offset);
}
}
......@@ -74,11 +74,15 @@ class NonSmoothShockCapture : public ArtificialDiffusion
Array<OneD, NekDouble > &mu);
private:
NonSmoothShockCapture(const LibUtilities::SessionReaderSharedPtr& pSession,
NonSmoothShockCapture(
const LibUtilities::SessionReaderSharedPtr& pSession,
const Array<OneD, MultiRegions::ExpListSharedPtr>& pFields,
const int spacedim);
virtual ~NonSmoothShockCapture(void){};
/// Parameters
int m_offset;
};
}
......
......@@ -45,7 +45,8 @@ std::string SmoothShockCapture::className = GetArtificialDiffusionFactory().
SmoothShockCapture::create,
"Smooth artificial diffusion for shock capture.");
SmoothShockCapture::SmoothShockCapture(const LibUtilities::SessionReaderSharedPtr& pSession,
SmoothShockCapture::SmoothShockCapture(
const LibUtilities::SessionReaderSharedPtr& pSession,
const Array<OneD, MultiRegions::ExpListSharedPtr>& pFields,
const int spacedim)
: ArtificialDiffusion(pSession, pFields, spacedim)
......@@ -53,6 +54,14 @@ SmoothShockCapture::SmoothShockCapture(const LibUtilities::SessionReaderSharedPt
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);
}
void SmoothShockCapture::v_DoArtificialDiffusion(
......@@ -127,19 +136,9 @@ void SmoothShockCapture::v_GetArtificialViscosity(
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, m_offset);
NekDouble ThetaH = m_FacH;
NekDouble ThetaL = m_FacL;
NekDouble Phi0 = (ThetaH+ThetaL)/2;
NekDouble DeltaPhi = ThetaH-Phi0;
NekDouble Phi0 = (m_FacH+m_FacL)/2;
NekDouble DeltaPhi = m_FacH-Phi0;
for (int e = 0; e < mu.num_elements(); e++)
{
......@@ -170,13 +169,13 @@ void SmoothShockCapture::GetForcingTerm(
Array<OneD, NekDouble> Sensor(nPts, 0.0);
Array<OneD, NekDouble> SensorKappa(nPts, 0.0);
Array <OneD, NekDouble > Lambda(nPts, 0.0);
Array <OneD, NekDouble > Tau(nPts, 1.0);
Array <OneD, NekDouble > soundspeed(nPts, 0.0);
Array <OneD, NekDouble > pressure(nPts, 0.0);
Array <OneD, NekDouble > absVelocity(nPts, 0.0);
NekDouble tau, pOrder;
Array<OneD,int> pOrderElmt = m_fields[0]->EvalBasisNumModesMaxPerExp();
Array<OneD, NekDouble> pOrder (nPts, 0.0);
// Thermodynamic related quantities
m_varConv->GetPressure(inarray, pressure);
......@@ -197,15 +196,14 @@ void SmoothShockCapture::GetForcingTerm(
for (int n = 0; n < nQuadPointsElement; n++)
{
pOrder[n + PointCount] = pOrderElmt[e];
pOrder = pOrderElmt[e];
// order 1.0e-06
Tau[n + PointCount] =
1.0 / (m_C1*pOrder[n + PointCount]*LambdaMax);
tau = 1.0 / (m_C1*pOrder*LambdaMax);
outarrayForcing[nvariables-1][n + PointCount] =
1 / Tau[n + PointCount] * (m_hFactor * LambdaMax /
pOrder[n + PointCount] *
1 / tau * (m_hFactor * LambdaMax /
pOrder *
SensorKappa[n + PointCount] -
inarray[nvariables-1][n + PointCount]);
}
......
......@@ -86,6 +86,15 @@ class SmoothShockCapture : public ArtificialDiffusion
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;
};
}
......
......@@ -338,130 +338,78 @@ namespace Nektar
Array<OneD, NekDouble> &SensorKappa,
int offset)
{
int e, NumModesElement, nQuadPointsElement;
int nTotQuadPoints = field->GetTotPoints();
int nElements = field->GetExpSize();
Array<OneD, NekDouble> tmp;
// Find solution (SolP) at p = P;
// The input array (physarray) is the solution at p = P;
Array<OneD,int> expOrderElement = field->EvalBasisNumModesMaxPerExp();
Array<OneD,int> ExpOrderElement = field->EvalBasisNumModesMaxPerExp();
Array<OneD, NekDouble> solution = physarray[0];
Array<OneD, NekDouble> SolP (nTotQuadPoints, 0.0);
Array<OneD, NekDouble> SolPmOne(nTotQuadPoints, 0.0);
Array<OneD, NekDouble> SolNorm (nTotQuadPoints, 0.0);
Vmath::Vcopy(nTotQuadPoints, physarray[0], 1, SolP, 1);
int CoeffsCount = 0;
for (e = 0; e < nElements; e++)
for (int e = 0; e < field->GetExpSize(); e++)
{
NumModesElement = ExpOrderElement[e];
int nQuadPointsElement = field->GetExp(e)->GetTotPoints();
int nCoeffsElement = field->GetExp(e)->GetNcoeffs();
int numCutOff = NumModesElement - offset;
Array<OneD, NekDouble> SolPElementPhys (nQuadPointsElement, 0.0);
Array<OneD, NekDouble> SolPElementCoeffs(nCoeffsElement, 0.0);
Array<OneD, NekDouble> SolPmOneElementPhys(nQuadPointsElement, 0.0);
Array<OneD, NekDouble> SolPmOneElementCoeffs(nCoeffsElement, 0.0);
// create vector the save the solution points per element at P = p;
for (int i = 0; i < nQuadPointsElement; i++)
{
SolPElementPhys[i] = SolP[CoeffsCount+i];
}
field->GetExp(e)->FwdTrans(SolPElementPhys,
SolPElementCoeffs);
int numModesElement = expOrderElement[e];
int nElmtPoints = field->GetExp(e)->GetTotPoints();
int physOffset = field->GetPhys_Offset(e);
int nElmtCoeffs = field->GetExp(e)->GetNcoeffs();
int numCutOff = numModesElement - offset;
// create vector to save the solution points per element at P = p;
Array<OneD, NekDouble> elmtPhys (nElmtPoints, solution+physOffset);
// Compute coefficients
Array<OneD, NekDouble> elmtCoeffs(nElmtCoeffs, 0.0);
field->GetExp(e)->FwdTrans(elmtPhys, elmtCoeffs);
// ReduceOrderCoeffs reduces the polynomial order of the solution
// that is represented by the coeffs given as an inarray. This is
// done by projecting the higher order solution onto the orthogonal
// basis and padding the higher order coefficients with zeros.
Array<OneD, NekDouble> reducedElmtCoeffs(nElmtCoeffs, 0.0);
field->GetExp(e)->ReduceOrderCoeffs(numCutOff,
SolPElementCoeffs,
SolPmOneElementCoeffs);
elmtCoeffs,
reducedElmtCoeffs);
field->GetExp(e)->BwdTrans(SolPmOneElementCoeffs,
SolPmOneElementPhys);
Array<OneD, NekDouble> reducedElmtPhys (nElmtPoints, 0.0);
field->GetExp(e)->BwdTrans(reducedElmtCoeffs,
reducedElmtPhys);
for (int i = 0; i < nQuadPointsElement; i++)
{
SolPmOne[CoeffsCount+i] = SolPmOneElementPhys[i];
}
NekDouble SolPmeanNumerator = 0.0;
NekDouble SolPmeanDenumerator = 0.0;
NekDouble numerator = 0.0;
NekDouble denominator = 0.0;
// Determining the norm of the numerator of the Sensor
Array<OneD, NekDouble> difference (nElmtPoints, 0.0);
Vmath::Vsub(nElmtPoints,
elmtPhys, 1,
reducedElmtPhys, 1,
difference, 1);
numerator = Vmath::Dot(nElmtPoints, difference, difference);
Vmath::Vsub(nQuadPointsElement,
SolPElementPhys, 1,
SolPmOneElementPhys, 1,
SolNorm, 1);
denominator = Vmath::Dot(nElmtPoints, elmtPhys, elmtPhys);
Vmath::Vmul(nQuadPointsElement,
SolNorm, 1,
SolNorm, 1,
SolNorm, 1);
NekDouble elmtSensor = sqrt(numerator / denominator);
elmtSensor = log10(elmtSensor);
Vmath::Fill(nElmtPoints, elmtSensor, tmp = Sensor + physOffset, 1);
for (int i = 0; i < nQuadPointsElement; i++)
NekDouble elmtSensorKappa;
if (numModesElement <= offset)
{
SolPmeanNumerator += SolNorm[i];
elmtSensorKappa = 0;
}
Vmath::Vmul(nQuadPointsElement,
SolPElementPhys, 1,
SolPElementPhys, 1,
SolNorm, 1);
for (int i = 0; i < nQuadPointsElement; i++)
else if (elmtSensor < (m_Skappa-m_Kappa))
{
SolPmeanDenumerator += SolNorm[i];
elmtSensorKappa = 0;
}
NekDouble elmtSensor =
sqrt(SolPmeanNumerator) / sqrt(SolPmeanDenumerator);
elmtSensor = log10(elmtSensor);
Vmath::Fill(nQuadPointsElement, elmtSensor,
(&Sensor[CoeffsCount]), 1);
CoeffsCount += nQuadPointsElement;
}
CoeffsCount = 0.0;
for (e = 0; e < nElements; e++)
{
NumModesElement = ExpOrderElement[e];
NekDouble ThetaS = m_mu0;
NekDouble Phi0 = m_Skappa;
NekDouble DeltaPhi = m_Kappa;
nQuadPointsElement = field->GetExp(e)->GetTotPoints();
for (int i = 0; i < nQuadPointsElement; i++)
else if (elmtSensor > (m_Skappa+m_Kappa))
{
if (Sensor[CoeffsCount+i] <= (Phi0 - DeltaPhi))
{
SensorKappa[CoeffsCount+i] = 0;
}
else if(Sensor[CoeffsCount+i] >= (Phi0 + DeltaPhi))
{
SensorKappa[CoeffsCount+i] = ThetaS;
}
else if(abs(Sensor[CoeffsCount+i]-Phi0) < DeltaPhi)
{
SensorKappa[CoeffsCount+i] =
ThetaS / 2 * (1 + sin(M_PI * (Sensor[CoeffsCount+i] -
Phi0) / (2 * DeltaPhi)));
}
elmtSensorKappa = m_mu0;
}
CoeffsCount += nQuadPointsElement;
else
{
elmtSensorKappa =
0.5 * m_mu0 *
(1 + sin(M_PI * (elmtSensor - m_Skappa) / (2*m_Kappa)));
}
Vmath::Fill(nElmtPoints, elmtSensorKappa,
tmp = SensorKappa + physOffset, 1);
}
}
......
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