Commit 13aa1f16 authored by Gianmarco Mengaldo's avatar Gianmarco Mengaldo
Browse files

Added FR diffusion for compressible NS.

parent 7ba01654
......@@ -7,6 +7,7 @@ SET(SOLVER_UTILS_SOURCES
Diffusion/DiffusionLDG.cpp
Diffusion/DiffusionLDGNS.cpp
Diffusion/DiffusionLFR.cpp
Diffusion/DiffusionLFRNS.cpp
Driver.cpp
DriverArnoldi.cpp
DriverModifiedArnoldi.cpp
......@@ -32,6 +33,7 @@ SET(SOLVER_UTILS_HEADERS
Diffusion/DiffusionLDG.h
Diffusion/DiffusionLDGNS.h
Diffusion/DiffusionLFR.h
Diffusion/DiffusionLFRNS.h
Driver.h
DriverArnoldi.h
DriverModifiedArnoldi.h
......
......@@ -61,8 +61,7 @@ namespace Nektar
const Array<OneD, Array<OneD, NekDouble> > &inarray,
Array<OneD, Array<OneD, NekDouble> > &outarray)
{
cout<<setprecision(16);
int i, j, k, num;
int i, j, k;
int nDim = fields[0]->GetCoordim(0);
int nPts = fields[0]->GetTotPoints();
int nCoeffs = fields[0]->GetNcoeffs();
......@@ -118,7 +117,7 @@ namespace Nektar
fields[i]->BwdTrans (qcoeffs, qfield[j][i]);
}
}
// Compute u from q_{\eta} and q_{\xi}
// Obtain numerical fluxes
v_NumFluxforVector(fields, inarray, qfield, flux[0]);
......
......@@ -34,6 +34,8 @@
///////////////////////////////////////////////////////////////////////////////
#include <SolverUtils/Diffusion/DiffusionLDGNS.h>
#include <iostream>
#include <iomanip>
namespace Nektar
{
......@@ -51,6 +53,16 @@ namespace Nektar
Array<OneD, MultiRegions::ExpListSharedPtr> pFields)
{
m_session = pSession;
m_session->LoadParameter ("Gamma", m_gamma, 1.4);
m_session->LoadParameter ("GasConstant", m_gasConstant, 287.058);
m_session->LoadParameter ("Twall", m_Twall, 300.15);
m_session->LoadSolverInfo("ViscosityType", m_ViscosityType,
"Constant");
m_session->LoadParameter ("mu", m_mu, 1.78e-05);
m_session->LoadParameter ("thermalConductivity",
m_thermalConductivity, 0.0257);
m_session->LoadParameter ("rhoInf", m_rhoInf, 1.225);
m_session->LoadParameter ("pInf", m_pInf, 101325);
}
/**
......@@ -71,7 +83,8 @@ namespace Nektar
const Array<OneD, MultiRegions::ExpListSharedPtr> &fields,
const Array<OneD, Array<OneD, NekDouble> > &inarray,
Array<OneD, Array<OneD, NekDouble> > &outarray)
{
{
cout<<setprecision(16);
int i, j;
int nDim = fields[0]->GetCoordim(0);
int nScalars = inarray.num_elements();
......@@ -247,12 +260,7 @@ namespace Nektar
const Array<OneD, MultiRegions::ExpListSharedPtr> &fields,
const Array<OneD, Array<OneD, NekDouble> > &inarray,
Array<OneD, Array<OneD, NekDouble> > &penaltyfluxO1)
{
// They should be taken from CFS
const NekDouble gamma = 1.4;
const NekDouble R = 287.05;
const NekDouble Twall = 0.012;
{
int cnt;
int i, j, e;
int id1, id2;
......@@ -266,8 +274,7 @@ namespace Nektar
Array<OneD, NekDouble> tmp1(nTracePts, 0.0);
Array<OneD, NekDouble> tmp2(nTracePts, 0.0);
Array<OneD, NekDouble> Tw(nTracePts, Twall);
Array<OneD, NekDouble> Tw(nTracePts, m_Twall);
Array< OneD, Array<OneD, NekDouble > > scalarVariables(nScalars);
Array< OneD, Array<OneD, NekDouble > > uplus(nScalars);
......@@ -305,9 +312,6 @@ namespace Nektar
GetPhys_Offset(fields[0]->GetTraceMap()->
GetBndCondTraceToGlobalTraceMap(cnt++));
// Should I reinforce velocity bcs at the wall?
// They don't give too different answers!
// =====================================================
// Reinforcing bcs for velocity in case of Wall bcs
if (fields[i]->GetBndConditions()[j]->
GetUserDefined() ==
......@@ -316,9 +320,7 @@ namespace Nektar
Vmath::Zero(nBndEdgePts,
&scalarVariables[i][id2], 1);
}
// =====================================================
// =====================================================
// Imposing velocity bcs if not Wall
else if (fields[i]->GetBndConditions()[j]->
GetBoundaryConditionType() ==
......@@ -331,7 +333,6 @@ namespace Nektar
UpdatePhys())[id1], 1,
&scalarVariables[i][id2], 1);
}
// =====================================================
// For Dirichlet boundary condition: uflux = u_bcs
if (fields[i]->GetBndConditions()[j]->
......@@ -377,10 +378,12 @@ namespace Nektar
GetBndCondExpansions().num_elements();
for (j = 0; j < nBndRegions; ++j)
{
//cout<<"bcRegion = "<< j << endl;
nBndEdges = fields[nScalars]->
GetBndCondExpansions()[j]->GetExpSize();
for (e = 0; e < nBndEdges; ++e)
{
//cout<<"bcEdge = "<< e << endl;
nBndEdgePts = fields[nScalars]->
GetBndCondExpansions()[j]->GetExp(e)->GetNumPoints(0);
......@@ -391,12 +394,11 @@ namespace Nektar
GetPhys_Offset(fields[0]->GetTraceMap()->
GetBndCondTraceToGlobalTraceMap(cnt++));
/*
// Imposing Temperature Twall at the wall
if (fields[i]->GetBndConditions()[j]->
GetUserDefined() ==
SpatialDomains::eWallViscous)
{
{
Vmath::Vcopy(nBndEdgePts,
&Tw[0], 1,
&scalarVariables[nScalars-1][id2], 1);
......@@ -407,30 +409,37 @@ namespace Nektar
GetBoundaryConditionType() ==
SpatialDomains::eDirichlet)
{
*/
// Divide E by rho
Vmath::Vdiv(nBndEdgePts,
&(fields[nScalars]->
GetBndCondExpansions()[j]->
GetPhys())[id1], 1,
GetPhys())[id1], 1,
&(fields[0]->
GetBndCondExpansions()[j]->
GetPhys())[id1], 1,
&scalarVariables[nScalars-1][id2], 1);
/*
// Subtract kinetic energy to E/rho
Vmath::Vsub(nBndEdgePts,
&scalarVariables[nScalars-1][id2], 1,
&tmp2[id2], 1,
&scalarVariables[nScalars-1][id2], 1);
*/
// Multiply by constant factor (gamma-1)/R
// (Note: I should extract R and gamma)
Vmath::Smul(nBndEdgePts, (gamma - 1)/R,
Vmath::Smul(nBndEdgePts, (m_gamma - 1)/m_gasConstant,
&scalarVariables[nScalars-1][id2], 1,
&scalarVariables[nScalars-1][id2], 1);
//}
}
/*
for (int bb = 0; bb < nBndEdgePts; ++bb)
{
cout<<"T-bcs = "<<scalarVariables[nScalars-1][id2+bb]<<endl;
}
int num;
cin>>num;
*/
// For Dirichlet boundary condition: uflux = u_bcs
if (fields[nScalars]->GetBndConditions()[j]->
GetBoundaryConditionType() ==
......@@ -450,6 +459,13 @@ namespace Nektar
&uplus[nScalars-1][id2], 1,
&penaltyfluxO1[nScalars-1][id2], 1);
}
/*
for (int bb = 0; bb < nBndEdgePts; ++bb)
{
cout<<"PenaltyFlux = "<<penaltyfluxO1[nScalars-1][id2+bb]<<endl;
}
cin>>num;
*/
}
}
}
......@@ -567,9 +583,11 @@ namespace Nektar
// Extract the physical values of the solution at the boundaries
fields[var]->ExtractTracePhys(qfield, qtemp);
//cout<<"var = "<< var << endl;
// Loop on the boundary regions to apply appropriate bcs
for (i = 0; i < nBndRegions; ++i)
{
//cout<<"bcRegion = "<< i << endl;
// Number of boundary regions related to region 'i'
nBndEdges = fields[var]->
GetBndCondExpansions()[i]->GetExpSize();
......@@ -577,6 +595,7 @@ namespace Nektar
// Weakly impose bcs by modifying flux values
for (e = 0; e < nBndEdges; ++e)
{
//cout<<"bcEdge = "<< e << endl;
nBndEdgePts = fields[var]->
GetBndCondExpansions()[i]->GetExp(e)->GetNumPoints(0);
......@@ -588,8 +607,7 @@ namespace Nektar
GetBndCondTraceToGlobalTraceMap(cnt++));
// In case of Dirichlet bcs:
// uflux = g_D
// qflux = q+
// uflux = gD
if(fields[var]->GetBndConditions()[i]->
GetBoundaryConditionType() == SpatialDomains::eDirichlet)
{
......@@ -615,6 +633,14 @@ namespace Nektar
&penaltyflux[id2], 1);
*/
}
/*
for (int bb = 0; bb < nBndEdgePts; ++bb)
{
cout<<"2nd - bcs = "<<penaltyflux[id2+bb]<<endl;
}
int num;
cin>>num;
*/
}
}
}
......
......@@ -55,8 +55,17 @@ namespace Nektar
protected:
DiffusionLDGNS();
Array<OneD, Array<OneD, NekDouble> > m_traceNormals;
LibUtilities::SessionReaderSharedPtr m_session;
Array<OneD, Array<OneD, NekDouble> > m_traceNormals;
LibUtilities::SessionReaderSharedPtr m_session;
NekDouble m_gamma;
NekDouble m_gasConstant;
NekDouble m_Twall;
std::string m_ViscosityType;
NekDouble m_mu;
NekDouble m_thermalConductivity;
NekDouble m_rhoInf;
NekDouble m_pInf;
virtual void v_InitObject(
LibUtilities::SessionReaderSharedPtr pSession,
......
This diff is collapsed.
///////////////////////////////////////////////////////////////////////////////
//
// File: DiffusionLFRNS.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: LFRNS diffusion class.
//
///////////////////////////////////////////////////////////////////////////////
#ifndef NEKTAR_SOLVERUTILS_DIFFUSIONLFRNS
#define NEKTAR_SOLVERUTILS_DIFFUSIONLFRNS
#include <SolverUtils/Diffusion/Diffusion.h>
namespace Nektar
{
namespace SolverUtils
{
class DiffusionLFRNS : public Diffusion
{
public:
static DiffusionSharedPtr create(std::string diffType)
{
return DiffusionSharedPtr(new DiffusionLFRNS(diffType));
}
static std::string type[];
Array<OneD, Array<OneD, NekDouble> > m_Q2D_e0;
Array<OneD, Array<OneD, NekDouble> > m_Q2D_e1;
Array<OneD, Array<OneD, NekDouble> > m_Q2D_e2;
Array<OneD, Array<OneD, NekDouble> > m_Q2D_e3;
Array<OneD, Array<OneD, NekDouble> > m_dGL_xi1;
Array<OneD, Array<OneD, NekDouble> > m_dGR_xi1;
Array<OneD, Array<OneD, NekDouble> > m_dGL_xi2;
Array<OneD, Array<OneD, NekDouble> > m_dGR_xi2;
Array<OneD, Array<OneD, NekDouble> > m_dGL_xi3;
Array<OneD, Array<OneD, NekDouble> > m_dGR_xi3;
DNekMatSharedPtr m_Ixm;
DNekMatSharedPtr m_Ixp;
protected:
DiffusionLFRNS(std::string diffType);
Array<OneD, Array<OneD, NekDouble> > m_traceNormals;
Array<OneD, Array<OneD, Array<OneD,NekDouble> > > m_tanbasis;
LibUtilities::SessionReaderSharedPtr m_session;
std::string m_diffType;
virtual void v_InitObject(
LibUtilities::SessionReaderSharedPtr pSession,
Array<OneD, MultiRegions::ExpListSharedPtr> pFields);
virtual void v_SetupMetrics(
LibUtilities::SessionReaderSharedPtr pSession,
Array<OneD, MultiRegions::ExpListSharedPtr> pFields);
virtual void v_SetupCFunctions(
LibUtilities::SessionReaderSharedPtr pSession,
Array<OneD, MultiRegions::ExpListSharedPtr> pFields);
virtual void v_SetupInterpolationMatrices(
LibUtilities::SessionReaderSharedPtr pSession,
Array<OneD, MultiRegions::ExpListSharedPtr> pFields);
virtual void v_Diffuse(
const int nConvective,
const Array<OneD, MultiRegions::ExpListSharedPtr> &fields,
const Array<OneD, Array<OneD, NekDouble> > &inarray,
Array<OneD, Array<OneD, NekDouble> > &outarray);
virtual void v_NumericalFluxO1(
const Array<OneD, MultiRegions::ExpListSharedPtr> &fields,
const Array<OneD, Array<OneD, NekDouble> > &inarray,
Array<OneD, Array<OneD, Array<OneD, NekDouble> > >
&numericalFluxO1);
virtual void v_WeakPenaltyO1(
const Array<OneD, MultiRegions::ExpListSharedPtr> &fields,
const Array<OneD, Array<OneD, NekDouble> > &inarray,
Array<OneD, Array<OneD, NekDouble> > &penaltyfluxO1);
virtual void v_NumericalFluxO2(
const Array<OneD, MultiRegions::ExpListSharedPtr> &fields,
const Array<OneD, Array<OneD, NekDouble> > &ufield,
Array<OneD, Array<OneD, Array<OneD, NekDouble> > >&qfield,
Array<OneD, Array<OneD, NekDouble> > &qflux);
virtual void v_WeakPenaltyO2(
const Array<OneD, MultiRegions::ExpListSharedPtr> &fields,
const int var,
const int dir,
const Array<OneD, const NekDouble> &qfield,
Array<OneD, NekDouble> &penaltyflux);
virtual void v_DerCFlux_1D(
const int nConvectiveFields,
const Array<OneD, MultiRegions::ExpListSharedPtr>&fields,
const Array<OneD, const NekDouble> &flux,
const Array<OneD, const NekDouble> &iFlux,
Array<OneD, NekDouble> &derCFlux);
virtual void v_DerCFlux_2D(
const int nConvectiveFields,
const int direction,
const Array<OneD, MultiRegions::ExpListSharedPtr>&fields,
const Array<OneD, const NekDouble> &flux,
const Array<OneD, NekDouble> &iFlux,
Array<OneD, NekDouble> &derCFlux);
virtual void v_DivCFlux_2D(
const int nConvectiveFields,
const Array<OneD, MultiRegions::ExpListSharedPtr>&fields,
const Array<OneD, const NekDouble> &fluxX1,
const Array<OneD, const NekDouble> &fluxX2,
const Array<OneD, const NekDouble> &numericalFlux,
Array<OneD, NekDouble> &divCFlux);
};
}
}
#endif
......@@ -2205,32 +2205,39 @@ namespace Nektar
case MultiRegions::eDiscontinuous:
{
if (DiffusionType == "LDG")
if (DiffusionType == "LDG" || DiffusionType == "LDGNS")
{
out << "\tDiffusion Type : LDG" <<endl;
}
else if (DiffusionType == "LFRDG")
else if (DiffusionType == "LFRDG" || DiffusionType == "LFRDGNS")
{
out << "\tDiffusion Type : LFRDG" <<endl;
}
else if (DiffusionType == "LFRSD")
else if (DiffusionType == "LFRSD" || DiffusionType == "LFRSDNS")
{
out << "\tDiffusion Type : LFRSD" <<endl;
}
else if (DiffusionType == "LFRHU")
else if (DiffusionType == "LFRHU" || DiffusionType == "LFRHUNS")
{
out << "\tDiffusion Type : LFRHU" <<endl;
}
else if (DiffusionType == "LFRcmin")
else if (DiffusionType == "LFRcmin" || DiffusionType == "LFRcminNS")
{
out << "\tDiffusion Type : LFR c = c-min" <<endl;
}
else if (DiffusionType == "LFRcinf")
else if (DiffusionType == "LFRcinf" || DiffusionType == "LFRcinfNS")
{
out << "\tDiffusion Type : LFR c = c-infinity" <<endl;
}
break;
}
case MultiRegions::eMixed_CG_Discontinuous:
{
break;
}
default:
break;
}
}
}
......
......@@ -117,6 +117,9 @@ namespace Nektar
m_session->LoadParameter ("mu", m_mu, 1.78e-05);
m_session->LoadParameter ("thermalConductivity",
m_thermalConductivity, 0.0257);
m_Cp = m_gamma / (m_gamma - 1.0) * m_gasConstant;
m_Prandtl = m_Cp * m_mu / m_thermalConductivity;
// Type of advection class to be used
......@@ -1090,7 +1093,7 @@ namespace Nektar
&derivativesO1[0][2][0], 1,
&tmp2[0], 1);
// STx = u * Sxx + v * Sxy + (K / mu) * dT/dx
// STx = u * Sxx + v * Sxy + K * dT/dx
Vmath::Vadd(nPts, &STx[0], 1, &tmp1[0], 1, &STx[0], 1);
Vmath::Vadd(nPts, &STx[0], 1, &tmp2[0], 1, &STx[0], 1);
......@@ -1113,7 +1116,7 @@ namespace Nektar
&derivativesO1[1][2][0], 1,
&tmp2[0], 1);
// STy = v * Syy + u * Sxy + (K / mu) * dT/dy
// STy = v * Syy + u * Sxy + K * dT/dy
Vmath::Vadd(nPts, &STy[0], 1, &tmp1[0], 1, &STy[0], 1);
Vmath::Vadd(nPts, &STy[0], 1, &tmp2[0], 1, &STy[0], 1);
}
......@@ -1200,7 +1203,7 @@ namespace Nektar
&derivativesO1[2][2][0], 1,
&tmp3[0], 1);
// STy = w * Szz + u * Sxz + v * Syz + K * dT/dz
// STz = w * Szz + u * Sxz + v * Syz + K * dT/dz
Vmath::Vadd(nPts, &STz[0], 1, &tmp1[0], 1, &STz[0], 1);
Vmath::Vadd(nPts, &STz[0], 1, &tmp2[0], 1, &STz[0], 1);
Vmath::Vadd(nPts, &STz[0], 1, &tmp3[0], 1, &STz[0], 1);
......
......@@ -103,6 +103,9 @@ namespace Nektar
std::string m_ViscosityType;
NekDouble m_mu;
NekDouble m_thermalConductivity;
NekDouble m_Cp;
NekDouble m_Prandtl;
CompressibleFlowSystem(
const LibUtilities::SessionReaderSharedPtr& pSession);
......
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