Commit 99f0de91 authored by Dirk Ekelschot's avatar Dirk Ekelschot
Browse files

Added the smooth artificial viscosity model

parent 19ba212a
......@@ -195,7 +195,7 @@ namespace Nektar
{
Array<OneD, NekDouble> tmp(m_numPoints), tmp2;
int nVel = advVel.num_elements();
// Call solver's flux vector function to compute the flux vector on
// the entire domain.
m_fluxVector(inarray, m_fluxVecStore);
......
......@@ -847,7 +847,7 @@ namespace Nektar
}
// Computing the interface flux at each trace point
m_riemann->Solve(Fwd, Bwd, numflux);
m_riemann->Solve(m_spaceDim, Fwd, Bwd, numflux);
// Divergence of the flux (computing the RHS)
switch(nDimensions)
......
......@@ -85,7 +85,7 @@ namespace Nektar
int nCoeffs = fields[0]->GetNcoeffs();
int nTracePointsTot = fields[0]->GetTrace()->GetTotPoints();
int i, j;
Array<OneD, Array<OneD, NekDouble> > tmp(nConvectiveFields);
Array<OneD, Array<OneD, Array<OneD, NekDouble> > > fluxvector(
nConvectiveFields);
......@@ -132,7 +132,7 @@ namespace Nektar
fields[i]->GetFwdBwdTracePhys(inarray[i], Fwd[i], Bwd[i]);
}
m_riemann->Solve(Fwd, Bwd, numflux);
m_riemann->Solve(m_spaceDim, Fwd, Bwd, numflux);
// Evaulate <\phi, \hat{F}\cdot n> - OutField[i]
for(i = 0; i < nConvectiveFields; ++i)
......
......@@ -61,12 +61,12 @@ namespace Nektar
const Array<OneD, Array<OneD, NekDouble> >&,
Array<OneD, Array<OneD, Array<OneD, NekDouble> > >&,
Array<OneD, Array<OneD, Array<OneD, NekDouble> > >&)>
DiffusionFluxVecCBNS;
DiffusionFluxVecCBNS;
typedef boost::function<void (
const Array<OneD, Array<OneD, NekDouble> >&,
Array<OneD, NekDouble >&)>
DiffusionArtificialDiffusion;
DiffusionArtificialDiffusion;
class Diffusion
{
......
......@@ -95,6 +95,7 @@ namespace Nektar
* @param flux Resultant flux along trace space.
*/
void RiemannSolver::Solve(
int nDim,
const Array<OneD, const Array<OneD, NekDouble> > &Fwd,
const Array<OneD, const Array<OneD, NekDouble> > &Bwd,
Array<OneD, Array<OneD, NekDouble> > &flux)
......@@ -123,13 +124,13 @@ namespace Nektar
rotateToNormal (Fwd, normals, m_rotStorage[0]);
rotateToNormal (Bwd, normals, m_rotStorage[1]);
v_Solve (m_rotStorage[0], m_rotStorage[1],
v_Solve (nDim, m_rotStorage[0], m_rotStorage[1],
m_rotStorage[2]);
rotateFromNormal(m_rotStorage[2], normals, flux);
}
else
{
v_Solve(Fwd, Bwd, flux);
v_Solve(nDim, Fwd, Bwd, flux);
}
}
......
......@@ -63,6 +63,7 @@ namespace Nektar
{
public:
SOLVER_UTILS_EXPORT void Solve(
int nDim,
const Array<OneD, const Array<OneD, NekDouble> > &Fwd,
const Array<OneD, const Array<OneD, NekDouble> > &Bwd,
Array<OneD, Array<OneD, NekDouble> > &flux);
......@@ -128,6 +129,8 @@ namespace Nektar
{
return m_params;
}
int m_spacedim;
protected:
/// Indicates whether the Riemann solver requires a rotation to be
......@@ -149,6 +152,7 @@ namespace Nektar
SOLVER_UTILS_EXPORT RiemannSolver();
virtual void v_Solve(
int nDim,
const Array<OneD, const Array<OneD, NekDouble> > &Fwd,
const Array<OneD, const Array<OneD, NekDouble> > &Bwd,
Array<OneD, Array<OneD, NekDouble> > &flux) = 0;
......
......@@ -82,6 +82,7 @@ namespace Nektar
* @param flux Resulting flux.
*/
void UpwindLDGSolver::v_Solve(
int nDim,
const Array<OneD, const Array<OneD, NekDouble> > &Fwd,
const Array<OneD, const Array<OneD, NekDouble> > &Bwd,
Array<OneD, Array<OneD, NekDouble> > &flux)
......
......@@ -58,6 +58,7 @@ namespace Nektar
UpwindLDGSolver();
virtual void v_Solve(
int nDim,
const Array<OneD, const Array<OneD, NekDouble> > &Fwd,
const Array<OneD, const Array<OneD, NekDouble> > &Bwd,
Array<OneD, Array<OneD, NekDouble> > &flux);
......
......@@ -82,6 +82,7 @@ namespace Nektar
* @param flux Resulting flux.
*/
void UpwindSolver::v_Solve(
int nDim,
const Array<OneD, const Array<OneD, NekDouble> > &Fwd,
const Array<OneD, const Array<OneD, NekDouble> > &Bwd,
Array<OneD, Array<OneD, NekDouble> > &flux)
......
......@@ -58,6 +58,7 @@ namespace Nektar
UpwindSolver();
virtual void v_Solve(
int nDim,
const Array<OneD, const Array<OneD, NekDouble> > &Fwd,
const Array<OneD, const Array<OneD, NekDouble> > &Bwd,
Array<OneD, Array<OneD, NekDouble> > &flux);
......
......@@ -69,6 +69,7 @@ namespace Nektar
eWall_Forces,
eWall,
eWallViscous,
eArtificialViscosity,
eSymmetry,
eRinglebFlow,
eTimeDependent,
......@@ -99,6 +100,7 @@ namespace Nektar
known_type["MG"] = eMG;
known_type["Wall"] = eWall;
known_type["WallViscous"] = eWallViscous;
known_type["ArtificialVisc"] = eArtificialViscosity;
known_type["Q-inflow"] = eQinflow;
known_type["Terminal"] = eTerminal;
known_type["R-terminal"] = eRterminal;
......
......@@ -103,10 +103,14 @@ namespace Nektar
NekDouble m_gasConstant;
NekDouble m_Twall;
std::string m_ViscosityType;
std::string m_shockCaptureType;
NekDouble m_mu;
NekDouble m_Skappa;
NekDouble m_Kappa;
NekDouble m_mu0;
NekDouble m_FacL;
NekDouble m_FacH;
NekDouble m_eps_max;
NekDouble m_thermalConductivity;
NekDouble m_Cp;
NekDouble m_Prandtl;
......@@ -124,6 +128,9 @@ namespace Nektar
void GetFluxVector(
const Array<OneD, Array<OneD, NekDouble> > &physfield,
Array<OneD, Array<OneD, Array<OneD, NekDouble> > > &flux);
void GetFluxVectorPDESC(
const Array<OneD, Array<OneD, NekDouble> > &physfield,
Array<OneD, Array<OneD, Array<OneD, NekDouble> > > &flux);
void GetFluxVectorDeAlias(
const Array<OneD, Array<OneD, NekDouble> > &physfield,
Array<OneD, Array<OneD, Array<OneD, NekDouble> > > &flux);
......@@ -131,6 +138,10 @@ namespace Nektar
const Array<OneD, Array<OneD, NekDouble> > &physfield,
Array<OneD, Array<OneD, Array<OneD, NekDouble> > > &derivatives,
Array<OneD, Array<OneD, Array<OneD, NekDouble> > > &viscousTensor);
void GetViscousFluxVectorPDESC(
const Array<OneD, Array<OneD, NekDouble> > &physfield,
Array<OneD, Array<OneD, Array<OneD, NekDouble> > > &derivativesO1,
Array<OneD, Array<OneD, Array<OneD, NekDouble> > > &viscousTensor);
void GetViscousFluxVectorDeAlias(
const Array<OneD, Array<OneD, NekDouble> > &physfield,
Array<OneD, Array<OneD, Array<OneD, NekDouble> > > &derivatives,
......@@ -143,6 +154,10 @@ namespace Nektar
int bcRegion,
int cnt,
Array<OneD, Array<OneD, NekDouble> > &physarray);
void ArtificialViscosityBC(
int bcRegion,
int cnt,
Array<OneD, Array<OneD, NekDouble> > &physarray);
void SymmetryBC(
int bcRegion,
int cnt,
......@@ -177,11 +192,18 @@ namespace Nektar
const Array<OneD, const Array<OneD, NekDouble> >&physfield,
const Array<OneD, const Array<OneD, NekDouble> >&velocity,
Array<OneD, NekDouble> &pressure);
void GetEnthalpy(
const Array<OneD, const Array<OneD, NekDouble> > &physfield,
Array<OneD, NekDouble> &pressure,
Array<OneD, NekDouble> &enthalpy);
void GetEntropy(
const Array<OneD, const Array<OneD, NekDouble> > &physfield,
const Array<OneD, const NekDouble> &pressure,
const Array<OneD, const NekDouble> &temperature,
Array<OneD, NekDouble> &entropy);
void GetSmoothArtificialViscosity(
const Array<OneD, Array<OneD, NekDouble> > &physfield,
Array<OneD, NekDouble > &eps_bar);
void GetDynamicViscosity(
const Array<OneD, const NekDouble> &temperature,
Array<OneD, NekDouble >&mu);
......@@ -190,7 +212,11 @@ namespace Nektar
Array<OneD, NekDouble> &stdV);
void GetSensor(
const Array<OneD, const Array<OneD, NekDouble> > &physarray,
Array<OneD, NekDouble> &Sensor);
Array<OneD, NekDouble> &Sensor,
Array<OneD, NekDouble> &SensorKappa);
void GetElementDimensions(
Array<OneD, Array<OneD, NekDouble> > &outarray,
Array<OneD, NekDouble > &hmin);
void GetAbsoluteVelocity(
const Array<OneD, const Array<OneD, NekDouble> > &inarray,
Array<OneD, NekDouble> &Vtot);
......@@ -200,7 +226,9 @@ namespace Nektar
void SetVarPOrderElmt(
const Array<OneD, const Array<OneD, NekDouble> > &physfield,
Array<OneD, NekDouble > &PolyOrder);
void GetForcingTerm(
const Array<OneD, const Array<OneD, NekDouble> > &inarray,
Array<OneD, Array<OneD, NekDouble> > outarrayForcing);
virtual NekDouble v_GetTimeStep(
const Array<OneD, const Array<OneD, NekDouble> > &inarray);
......
......@@ -261,7 +261,8 @@ namespace Nektar
}
Array<OneD, NekDouble> sensor(npts,0.0);
GetSensor(physfield, sensor);
Array<OneD, NekDouble> SensorKappa(npts,0.0);
GetSensor(physfield, sensor, SensorKappa);
m_fields[0]->FwdTrans(sensor, outarray);
}
......
......@@ -342,7 +342,8 @@ namespace Nektar
}
Array<OneD, NekDouble> sensor(npts,0.0);
GetSensor(physfield, sensor);
Array<OneD, NekDouble> SensorKappa(npts,0.0);
GetSensor(physfield, sensor, SensorKappa);
m_fields[0]->FwdTrans(sensor, outarray);
}
......
......@@ -80,6 +80,10 @@ namespace Nektar
{
ASSERTL0(false, "Implicit CFE not set up.");
}
m_checkpointFuncs["Sensor"] = boost::bind(&NavierStokesCFE::CPSensor, this, _1, _2);
m_checkpointFuncs["SensorKappa"] = boost::bind(&NavierStokesCFE::CPSensorKappa, this, _1, _2);
m_checkpointFuncs["EpsSmooth"] = boost::bind(&NavierStokesCFE::CPSmoothArtificialViscosity, this, _1, _2);
}
......@@ -135,69 +139,205 @@ namespace Nektar
int nvariables = inarray.num_elements();
int npoints = GetNpoints();
Array<OneD, Array<OneD, NekDouble> > advVel(m_spacedim);
Array<OneD, Array<OneD, NekDouble> > outarrayAdv(nvariables);
Array<OneD, Array<OneD, NekDouble> > outarrayDiff(nvariables);
Array<OneD, Array<OneD, NekDouble> > inarrayTemp(nvariables-1);
Array<OneD, Array<OneD, NekDouble> > inarrayDiff(nvariables-1);
for (i = 0; i < nvariables; ++i)
{
outarrayAdv[i] = Array<OneD, NekDouble>(npoints, 0.0);
outarrayDiff[i] = Array<OneD, NekDouble>(npoints, 0.0);
}
for (i = 0; i < nvariables-1; ++i)
{
inarrayTemp[i] = Array<OneD, NekDouble>(npoints, 0.0);
inarrayDiff[i] = Array<OneD, NekDouble>(npoints, 0.0);
}
// Advection term in physical rhs form
m_advection->Advect(nvariables, m_fields, advVel, inarray, outarrayAdv);
for (i = 0; i < nvariables; ++i)
if(m_shockCaptureType == "Off")
{
Vmath::Neg(npoints, outarrayAdv[i], 1);
}
// Extract pressure and temperature
Array<OneD, NekDouble > pressure (npoints, 0.0);
Array<OneD, NekDouble > temperature(npoints, 0.0);
GetPressure(inarray, pressure);
GetTemperature(inarray, pressure, temperature);
int i;
int nvariables = inarray.num_elements();
int npoints = GetNpoints();
Array<OneD, Array<OneD, NekDouble> > advVel(m_spacedim);
Array<OneD, Array<OneD, NekDouble> > outarrayAdv(nvariables);
Array<OneD, Array<OneD, NekDouble> > outarrayDiff(nvariables);
Array<OneD, Array<OneD, NekDouble> > inarrayTemp(nvariables-1);
Array<OneD, Array<OneD, NekDouble> > inarrayDiff(nvariables-1);
for (i = 0; i < nvariables; ++i)
{
outarrayAdv[i] = Array<OneD, NekDouble>(npoints, 0.0);
outarrayDiff[i] = Array<OneD, NekDouble>(npoints, 0.0);
}
for (i = 0; i < nvariables-1; ++i)
{
inarrayTemp[i] = Array<OneD, NekDouble>(npoints, 0.0);
inarrayDiff[i] = Array<OneD, NekDouble>(npoints, 0.0);
}
// Advection term in physical rhs form
m_advection->Advect(nvariables, m_fields, advVel, inarray, outarrayAdv);
/*
cout << "ADVECTION = " << Vmath::Vmax(outarrayAdv[0].num_elements(), &outarrayAdv[0][0], 1) << endl;
cout << "ADVECTION = " << Vmath::Vmax(outarrayAdv[0].num_elements(), &outarrayAdv[1][0], 1) << endl;
cout << "ADVECTION = " << Vmath::Vmax(outarrayAdv[0].num_elements(), &outarrayAdv[2][0], 1) << endl;
cout << "ADVECTION = " << Vmath::Vmax(outarrayAdv[0].num_elements(), &outarrayAdv[3][0], 1) << endl;
cout << endl;*/
// Extract velocities
for (i = 1; i < nvariables-1; ++i)
{
Vmath::Vdiv(npoints,
inarray[i], 1,
inarray[0], 1,
inarrayTemp[i-1], 1);
for (i = 0; i < nvariables; ++i)
{
Vmath::Neg(npoints, outarrayAdv[i], 1);
}
// Extract pressure and temperature
Array<OneD, NekDouble > pressure (npoints, 0.0);
Array<OneD, NekDouble > temperature(npoints, 0.0);
GetPressure(inarray, pressure);
GetTemperature(inarray, pressure, temperature);
// Extract velocities
for (i = 1; i < nvariables-1; ++i)
{
Vmath::Vdiv(npoints,
inarray[i], 1,
inarray[0], 1,
inarrayTemp[i-1], 1);
}
// Copy velocities into new inarrayDiff
for (i = 0; i < nvariables-2; ++i)
{
Vmath::Vcopy(npoints, inarrayTemp[i], 1, inarrayDiff[i], 1);
}
// Copy temperature into new inarrayDiffusion
Vmath::Vcopy(npoints,
temperature, 1,
inarrayDiff[nvariables-2], 1);
// Diffusion term in physical rhs form
m_diffusion->Diffuse(nvariables, m_fields, inarrayDiff, outarrayDiff);
for (i = 0; i < nvariables; ++i)
{
Vmath::Vadd(npoints,
outarrayAdv[i], 1,
outarrayDiff[i], 1,
outarray[i], 1);
}
}
// Copy velocities into new inarrayDiff
for (i = 0; i < nvariables-2; ++i)
if(m_shockCaptureType == "Smooth")
{
Vmath::Vcopy(npoints, inarrayTemp[i], 1, inarrayDiff[i], 1);
}
// Copy temperature into new inarrayDiffusion
Vmath::Vcopy(npoints,
temperature, 1,
inarrayDiff[nvariables-2], 1);
// Diffusion term in physical rhs form
m_diffusion->Diffuse(nvariables, m_fields, inarrayDiff, outarrayDiff);
Array<OneD, Array<OneD, NekDouble> > advVel;
Array<OneD, Array<OneD, NekDouble> > outarrayAdv(nvariables);
Array<OneD, Array<OneD, NekDouble> > outarrayDiff(nvariables);
Array<OneD, Array<OneD, NekDouble> > outarrayForcing(nvariables);
for (i = 0; i < nvariables; ++i)
// contains u,v,w
Array<OneD, Array<OneD, NekDouble> > inarrayTemp(nvariables-1);
Array<OneD, Array<OneD, NekDouble> > inarrayDiff(nvariables+1);
for (i = 0; i < nvariables; ++i)
{
outarrayAdv[i] = Array<OneD, NekDouble>(npoints, 0.0);
outarrayDiff[i] = Array<OneD, NekDouble>(npoints, 0.0);
outarrayForcing[i] = Array<OneD, NekDouble>(npoints, 0.0);
}
for (i = 0; i < nvariables-1; ++i)
{
inarrayTemp[i] = Array<OneD, NekDouble>(npoints, 0.0);
}
for (i = 0; i < nvariables+1; ++i)
{
inarrayDiff[i] = Array<OneD, NekDouble>(npoints, 0.0);
}
// Advection term in physical rhs form
m_advection->Advect(nvariables, m_fields, advVel, inarray, outarrayAdv);
/*cout << "ADVECTION = " << Vmath::Vmax(outarrayAdv[0].num_elements(), &outarrayAdv[0][0], 1) << endl;
cout << "ADVECTION = " << Vmath::Vmax(outarrayAdv[0].num_elements(), &outarrayAdv[1][0], 1) << endl;
cout << "ADVECTION = " << Vmath::Vmax(outarrayAdv[0].num_elements(), &outarrayAdv[2][0], 1) << endl;
cout << "ADVECTION = " << Vmath::Vmax(outarrayAdv[0].num_elements(), &outarrayAdv[3][0], 1) << endl;
cout << "ADVECTION = " << Vmath::Vmax(outarrayAdv[0].num_elements(), &outarrayAdv[4][0], 1) << endl;*/
/*cout << Vmath::Vmin(outarray[1].num_elements(), outarray[1], 1)<<
" " << Vmath::Vmax(outarray[1].num_elements(), outarray[1], 1) << endl;
cout << Vmath::Vmin(outarray[2].num_elements(), outarray[2], 1) <<
" " << Vmath::Vmax(outarray[2].num_elements(), outarray[2], 1) << endl;*/
for (i = 0; i < nvariables; ++i)
{
Vmath::Neg(npoints, outarrayAdv[i], 1);
}
// Extract pressure and temperature
Array<OneD, NekDouble > pressure (npoints, 0.0);
Array<OneD, NekDouble > temperature(npoints, 0.0);
Array<OneD, NekDouble > enthalpy (npoints, 0.0);
GetPressure(inarray, pressure);
GetTemperature(inarray, pressure, temperature);
GetEnthalpy(inarray, pressure, enthalpy);
// Extract velocities
for (i = 1; i < nvariables-2; ++i)
{
Vmath::Vdiv(npoints,
inarray[i], 1,
inarray[0], 1,
inarrayTemp[i-1], 1);
}
// Copy velocities into new inarrayDiff
for (i = 0; i < nvariables-3; ++i)
{
Vmath::Vcopy(npoints, inarrayTemp[i], 1, inarrayDiff[i], 1);
}
// Copy temperature into new inarrayDiffusion
Vmath::Vcopy(npoints,
temperature, 1,
inarrayDiff[nvariables-3], 1);
// Copy density into new inarrayDiffusion
Vmath::Vcopy(npoints,
inarray[0], 1,
inarrayDiff[nvariables-2], 1);
// Copy artificial viscosity coefficient into new inarrayDiffusion
Vmath::Vcopy(npoints,
inarray[nvariables-1], 1,
inarrayDiff[nvariables-1], 1);
// Copy enthalpy into new inarrayDiffusion
Vmath::Vcopy(npoints,
enthalpy, 1,
inarrayDiff[nvariables], 1);
// Diffusion term in physical rhs form
m_diffusion->Diffuse(nvariables, m_fields, inarrayDiff, outarrayDiff);
/*cout << "DIFFUSION = " << Vmath::Vmax(outarrayDiff[0].num_elements(), &outarrayDiff[0][0], 1) << endl;
cout << "DIFFUSION = " << Vmath::Vmax(outarrayDiff[0].num_elements(), &outarrayDiff[1][0], 1) << endl;
cout << "DIFFUSION = " << Vmath::Vmax(outarrayDiff[0].num_elements(), &outarrayDiff[2][0], 1) << endl;
cout << "DIFFUSION = " << Vmath::Vmax(outarrayDiff[0].num_elements(), &outarrayDiff[3][0], 1) << endl;
cout << "DIFFUSION = " << Vmath::Vmax(outarrayDiff[0].num_elements(), &outarrayDiff[4][0], 1) << endl;*/
GetForcingTerm(inarray, outarrayForcing);
//cout << "FORCING TERM?! " << Vmath::Vmin(outarrayForcing[0].num_elements(),
// outarrayForcing[nvariables-1], 1) <<
// " " << Vmath::Vmax(outarrayForcing[0].num_elements(),
// outarrayForcing[nvariables-1], 1) << endl;
for (i = 0; i < nvariables; ++i)
{
Vmath::Vadd(npoints,
outarrayAdv[i], 1,
outarrayDiff[i], 1,
outarray[i], 1);
// Add Forcing Term
Vmath::Vadd(npoints,
outarray[i], 1,
outarrayForcing[i], 1,
outarray[i], 1);
}
}
if (m_shockCaptureType == "NonSmooth")
{
Vmath::Vadd(npoints,
outarrayAdv[i], 1,
outarrayDiff[i], 1,
outarray[i], 1);
ASSERTL0(false, "NS with non-smooth shock capturing not yet implemented");
}
//cout << endl;
}
void NavierStokesCFE::DoOdeProjection(
......@@ -260,6 +400,13 @@ namespace Nektar
WallViscousBC(n, cnt, inarray);
}
// artificial Condition
if (m_fields[nvariables-1]->GetBndConditions()[n]->GetUserDefined() ==
SpatialDomains::eArtificialViscosity)
{
ArtificialViscosityBC(n, cnt, inarray);
}
// Symmetric Boundary Condition
if (m_fields[0]->GetBndConditions()[n]->GetUserDefined() ==
SpatialDomains::eSymmetry)
......@@ -294,4 +441,75 @@ namespace Nektar
cnt += m_fields[0]->GetBndCondExpansions()[n]->GetExpSize();
}
}
void NavierStokesCFE::CPSensorKappa(
const Array<OneD, const Array<OneD, NekDouble> > &inarray,
Array<OneD, NekDouble> &outarray)
{