Commit 1aae6b06 by Michael Turner

Merge remote-tracking branch 'upstream/master' into feature/projection-meshing

parents 92868911 e2bd965e
......@@ -17,6 +17,8 @@ v5.0.0
- Added native support for csv files in addititon to pts (!760 !835)
- Utilize LAPACK_DIR env variable to find the native blas/lapack install (!827)
- Remove StdExpansion use from MultiRegion (use Expansions instead). (!831)
- Move steady state check and CFL output from solvers to SolverUtils (!832)
- Remove DG advection implementation from EquationSystem (!832)
- Simplify RawType typedefs (!840)
- Remove unused files from BasicUtils (!841)
- Remove checks for old boost versions which are no longer supported (!841)
......@@ -24,6 +26,7 @@ v5.0.0
- Added support for using the distance to a specific region (e.g. outlet) in the
function definitions for the Absorption Forcing (!769)
- Improve performance of DisContField2D::v_ExtractTracePhys (!824)
- Fix small bug in Jacobian Energy (!857)
**NekMesh**:
- Add feature to read basic 2D geo files as CAD (!731)
......@@ -48,8 +51,16 @@ v5.0.0
- Enable output to multiple files (!844)
- Allow using xml file without expansion tag in FieldConvert (!849)
**IncNavierStokesSolver**
- Replace steady-state check based on difference of norms by check based on
norm of the difference, to be consistent with the compressible solver (!832)
**CompressibleFlowSolver**
- Add 3D regression tests (!567)
- Introduce forcing for quasi-1D Euler simulations (!771)
- Allow performing axi-symmetric Euler simulations (!771)
- Add ability to use an exponential filtering for stabilization with
seg, quad and hex elements (!771, !862)
**Documentation**:
- Added the developer-guide repository as a submodule (!751)
......@@ -59,6 +70,15 @@ v4.4.2
**Library**
- Fix evaluation of points (e.g. HistoryPoints, Interpolation to pts) close to
the interface of two elements (!836)
- Fix deadlock in Hdf5 with homogeneous expansions (!858)
**NekMesh**
- Fix missing periodic boundary meshing and boundary layer mesh adjustment
configurations in 2D (!859)
**Documentation**:
- Fix sign of the viscous term in the velocity correction scheme equations in
the user guide (!856)
v4.4.1
------
......
......@@ -66,7 +66,7 @@ zero. If we now integrate the 1st, 2nd and last term in equation
(\ref{eqn.weakp}) by parts we can obtain the weak pressure equation
\begin{align}
\int_{\Omega} \nabla q \cdot \nabla p^{n+1} &= \int_{\Omega} q\, \nabla \cdot \left ( \frac{\partial \mathbf{u}}{\partial t}^{n+1} + \mathbf{N}(\mathbf{u})^{n+1} \right ) \nonumber \\
&- \int_{\partial \Omega} q \left ( \frac{\partial \mathbf{u}}{\partial t}^{n+1} + \mathbf{N}(\mathbf{u})^{n+1} - \nu \nabla \times \nabla \times \mathbf{u}^{n+1} \right ) \cdot \mathbf{n}
&- \int_{\partial \Omega} q \left ( \frac{\partial \mathbf{u}}{\partial t}^{n+1} + \mathbf{N}(\mathbf{u})^{n+1} + \nu \nabla \times \nabla \times \mathbf{u}^{n+1} \right ) \cdot \mathbf{n}
\label{eqn.weakp1}
\end{align}
where $\partial \Omega$ is the boundary of the domain and we have used
......@@ -98,7 +98,7 @@ which to decouple the system we impose that $\nabla \cdot
\int_{\Omega} \nabla q \cdot \nabla p^{n+1} &=
\int_{\Omega} q \, \nabla \cdot \left (-\frac{\hat{\mathbf{u}}}{\Delta t}
+ \mathbf{N}(\mathbf{u})^{*,n+1} \right ) \nonumber \\
&- \int_{\partial \Omega} q \left ( \frac{\partial \mathbf{u}}{\partial t}^{n+1} + \mathbf{N}(\mathbf{u})^{*.n+1} - \nu (\nabla \times \nabla \times \mathbf{u})^{*,n+1} \right ) \cdot \mathbf{n}
&- \int_{\partial \Omega} q \left ( \frac{\partial \mathbf{u}}{\partial t}^{n+1} + \mathbf{N}(\mathbf{u})^{*.n+1} + \nu (\nabla \times \nabla \times \mathbf{u})^{*,n+1} \right ) \cdot \mathbf{n}
\label{eqn.weakpfinal}
\end{align}
We note this can be recast into an equivalent strong form of the
......@@ -109,7 +109,7 @@ which to decouple the system we impose that $\nabla \cdot
\end{equation}
with consistent Neumann boundary conditions prescribed as
\begin{equation}
\frac{\partial p}{\partial n}^{n+1}= - \Bigl[ \frac{\partial \mathbf{u}}{\partial t}^{n+1} - \nu (\nabla \times \nabla \times \mathbf{u})^{*,n+1} + \mathbf{N}^{*,n+1}\Bigr]\cdot \mathbf{n}
\frac{\partial p}{\partial n}^{n+1}= - \Bigl[ \frac{\partial \mathbf{u}}{\partial t}^{n+1} + \nu (\nabla \times \nabla \times \mathbf{u})^{*,n+1} + \mathbf{N}^{*,n+1}\Bigr]\cdot \mathbf{n}
\label{eqn.pressurebcs}
\end{equation}
......
......@@ -98,7 +98,7 @@ void ProcessJacobianEnergy::Process(po::variables_map &vm)
StdRegions::StdExpansionSharedPtr Elmt = exp->GetExp(i);
const StdRegions::StdExpansion * sep = &( *Elmt );
const LocalRegions::Expansion * lep = dynamic_cast<const LocalRegions::Expansion*>( lep );
const LocalRegions::Expansion * lep = dynamic_cast<const LocalRegions::Expansion*>( sep );
int nquad = Elmt->GetTotPoints();
int coeffoffset = exp->GetCoeff_Offset(i);
......
......@@ -522,7 +522,7 @@ void FieldIOHdf5::v_Write(const std::string &outFile,
H5::DataSpaceSharedPtr homz_space = H5::DataSpace::OneD(nTotHomZ);
H5::DataSetSharedPtr homz_dset =
root->CreateDataSet("HOMOGENEOUSZIDS", homz_type, homz_space);
ASSERTL1(data_dset,
ASSERTL1(homz_dset,
prfx.str() + "cannot create HOMOGENEOUSZIDS dataset.");
}
......@@ -830,6 +830,27 @@ void FieldIOHdf5::v_Write(const std::string &outFile,
ids_dset->Write(
fielddefs[nFields - 1]->m_elementIDs, ids_fspace, writePL);
data_dset->Write(fielddata[nFields - 1], data_fspace, writePL);
if (order_dset)
{
order_dset->Write(numModesPerDirVar[nFields - 1],
order_fspace, writePL);
}
if (homy_dset)
{
homy_dset->Write(homoYIDs[nFields - 1], homy_fspace, writePL);
}
if (homz_dset)
{
homz_dset->Write(homoZIDs[nFields - 1], homz_fspace, writePL);
}
if (homs_dset)
{
homs_dset->Write(homoSIDs[nFields - 1], homs_fspace, writePL);
}
}
m_comm->Block();
......
......@@ -63,6 +63,92 @@ AdvectionSystem::~AdvectionSystem()
void AdvectionSystem::v_InitObject()
{
UnsteadySystem::v_InitObject();
m_session->LoadParameter("IO_CFLSteps", m_cflsteps, 0);
}
/**
*
*/
bool AdvectionSystem::v_PostIntegrate(int step)
{
bool result = UnsteadySystem::v_PostIntegrate(step);
if(m_cflsteps && !((step+1)%m_cflsteps))
{
int elmtid;
NekDouble cfl = GetCFLEstimate(elmtid);
if(m_comm->GetRank() == 0)
{
if( m_HomogeneousType == eNotHomogeneous)
{
cout << "CFL: ";
}
else
{
cout << "CFL (zero plane): ";
}
cout << cfl << " (in elmt " << elmtid << ")" << endl;
}
}
return result;
}
/**
*
*/
Array<OneD, NekDouble> AdvectionSystem::GetElmtCFLVals(void)
{
int nelmt = m_fields[0]->GetExpSize();
const Array<OneD, int> expOrder = GetNumExpModesPerExp();
const NekDouble cLambda = 0.2; // Spencer book pag. 317
Array<OneD, NekDouble> stdVelocity(nelmt, 0.0);
stdVelocity = v_GetMaxStdVelocity();
Array<OneD, NekDouble> cfl(nelmt, 0.0);
for(int el = 0; el < nelmt; ++el)
{
cfl[el] = m_timestep*(stdVelocity[el] * cLambda *
(expOrder[el]-1) * (expOrder[el]-1));
}
return cfl;
}
/**
*
*/
NekDouble AdvectionSystem::GetCFLEstimate(int &elmtid)
{
int n_element = m_fields[0]->GetExpSize();
Array<OneD, NekDouble> cfl = GetElmtCFLVals();
elmtid = Vmath::Imax(n_element,cfl,1);
NekDouble CFL,CFL_loc;
CFL = CFL_loc = cfl[elmtid];
m_comm->AllReduce(CFL,LibUtilities::ReduceMax);
// unshuffle elmt id if data is not stored in consecutive order.
elmtid = m_fields[0]->GetExp(elmtid)->GetGeom()->GetGlobalID();
if(CFL != CFL_loc)
{
elmtid = -1;
}
m_comm->AllReduce(elmtid,LibUtilities::ReduceMax);
// express element id with respect to plane
if(m_HomogeneousType == eHomogeneous1D)
{
elmtid = elmtid%m_fields[0]->GetPlane(0)->GetExpSize();
}
return CFL;
}
}
......
......@@ -54,14 +54,32 @@ public:
SOLVER_UTILS_EXPORT virtual void v_InitObject();
/// Returns the advection object held by this instance.
AdvectionSharedPtr GetAdvObject()
SOLVER_UTILS_EXPORT AdvectionSharedPtr GetAdvObject()
{
return m_advObject;
}
SOLVER_UTILS_EXPORT Array<OneD, NekDouble> GetElmtCFLVals(void);
SOLVER_UTILS_EXPORT NekDouble GetCFLEstimate(int &elmtid);
protected:
/// Advection term
SolverUtils::AdvectionSharedPtr m_advObject;
SOLVER_UTILS_EXPORT virtual bool v_PostIntegrate(int step);
SOLVER_UTILS_EXPORT virtual Array<OneD, NekDouble> v_GetMaxStdVelocity()
{
ASSERTL0(false,
"v_GetMaxStdVelocity is not implemented by the base class.");
Array<OneD, NekDouble> dummy(1);
return dummy;
}
private:
/// dump cfl estimate
int m_cflsteps;
};
/// Shared pointer to an AdvectionSystem class
......
......@@ -54,10 +54,10 @@ namespace Nektar
/// Calculate the larger time-step mantaining the problem stable.
SOLVER_UTILS_EXPORT NekDouble GetTimeStep(
const Array<OneD, const Array<OneD, NekDouble> > &inarray);
/// CFL safety factor (comprise between 0 to 1).
NekDouble m_cflSafetyFactor;
protected:
/// Number of time steps between outputting status information.
int m_infosteps;
......@@ -81,6 +81,15 @@ namespace Nektar
/// forward transformed state.
bool m_homoInitialFwd;
/// Tolerance to which steady state should be evaluated at
NekDouble m_steadyStateTol;
/// Check for steady state at step interval
int m_steadyStateSteps;
/// Storage for previous solution for steady-state check
Array<OneD, Array<OneD, NekDouble> > m_previousSolution;
// Steady-state residual file
std::ofstream m_errFile;
std::vector<int> m_intVariables;
std::vector<FilterSharedPtr> m_filters;
......@@ -108,34 +117,11 @@ namespace Nektar
SOLVER_UTILS_EXPORT virtual void v_AppendOutput1D(
Array<OneD, Array<OneD, NekDouble> > &solution1D);
///
SOLVER_UTILS_EXPORT virtual void v_NumericalFlux(
Array<OneD, Array<OneD, NekDouble> > &physfield,
Array<OneD, Array<OneD, NekDouble> > &numflux);
///
SOLVER_UTILS_EXPORT virtual void v_NumericalFlux(
Array<OneD, Array<OneD, NekDouble> > &physfield,
Array<OneD, Array<OneD, NekDouble> > &numfluxX,
Array<OneD, Array<OneD, NekDouble> > &numfluxY );
///
SOLVER_UTILS_EXPORT virtual void v_NumFluxforScalar(
const Array<OneD, Array<OneD, NekDouble> > &ufield,
Array<OneD, Array<OneD, Array<OneD, NekDouble> > > &uflux);
///
SOLVER_UTILS_EXPORT virtual void v_NumFluxforVector(
const Array<OneD, Array<OneD, NekDouble> > &ufield,
Array<OneD, Array<OneD, Array<OneD, NekDouble> > > &qfield,
Array<OneD, Array<OneD, NekDouble> > &qflux);
SOLVER_UTILS_EXPORT virtual NekDouble v_GetTimeStep(
const Array<OneD, const Array<OneD, NekDouble> > &inarray);
SOLVER_UTILS_EXPORT virtual bool v_PreIntegrate(int step);
SOLVER_UTILS_EXPORT virtual bool v_PostIntegrate(int step);
SOLVER_UTILS_EXPORT virtual bool v_SteadyStateCheck(int step);
SOLVER_UTILS_EXPORT virtual bool v_RequireFwdTrans()
{
......@@ -153,21 +139,10 @@ namespace Nektar
private:
///
void WeakPenaltyforScalar(
const int var,
const Array<OneD, const NekDouble> &physfield,
Array<OneD, NekDouble> &penaltyflux,
NekDouble time=0.0);
///
void WeakPenaltyforVector(
const int var,
const int dir,
const Array<OneD, const NekDouble> &physfield,
Array<OneD, NekDouble> &penaltyflux,
NekDouble C11,
NekDouble time=0.0);
void InitializeSteadyState();
bool CheckSteadyState(int step);
};
}
......
......@@ -2436,5 +2436,69 @@ namespace Nektar
// backward transform to physical space
OrthoExp.BwdTrans(orthocoeffs,array);
}
void StdHexExp::v_ExponentialFilter(
Array<OneD, NekDouble> &array,
const NekDouble alpha,
const NekDouble exponent,
const NekDouble cutoff)
{
// Generate an orthogonal expansion
int qa = m_base[0]->GetNumPoints();
int qb = m_base[1]->GetNumPoints();
int qc = m_base[2]->GetNumPoints();
int nmodesA = m_base[0]->GetNumModes();
int nmodesB = m_base[1]->GetNumModes();
int nmodesC = m_base[2]->GetNumModes();
int P = nmodesA - 1;
int Q = nmodesB - 1;
int R = nmodesC - 1;
// Declare orthogonal basis.
LibUtilities::PointsKey pa(qa,m_base[0]->GetPointsType());
LibUtilities::PointsKey pb(qb,m_base[1]->GetPointsType());
LibUtilities::PointsKey pc(qc,m_base[2]->GetPointsType());
LibUtilities::BasisKey Ba(LibUtilities::eOrtho_A, nmodesA, pa);
LibUtilities::BasisKey Bb(LibUtilities::eOrtho_A, nmodesB, pb);
LibUtilities::BasisKey Bc(LibUtilities::eOrtho_A, nmodesC, pc);
StdHexExp OrthoExp(Ba,Bb,Bc);
// Cutoff
int Pcut = cutoff*P;
int Qcut = cutoff*Q;
int Rcut = cutoff*R;
// Project onto orthogonal space.
Array<OneD, NekDouble> orthocoeffs(OrthoExp.GetNcoeffs());
OrthoExp.FwdTrans(array,orthocoeffs);
//
NekDouble fac, fac1, fac2, fac3;
int index = 0;
for(int i = 0; i < nmodesA; ++i)
{
for(int j = 0; j < nmodesB; ++j)
{
for(int k = 0; k < nmodesC; ++k, ++index)
{
//to filter out only the "high-modes"
if(i > Pcut || j > Qcut || k > Rcut)
{
fac1 = (NekDouble) (i - Pcut)/( (NekDouble)(P - Pcut) );
fac2 = (NekDouble) (j - Qcut)/( (NekDouble)(Q - Qcut) );
fac3 = (NekDouble) (k - Rcut)/( (NekDouble)(R - Rcut) );
fac = max( max(fac1, fac2), fac3);
fac = pow(fac, exponent);
orthocoeffs[index] *= exp(-alpha*fac);
}
}
}
}
// backward transform to physical space
OrthoExp.BwdTrans(orthocoeffs,array);
}
}
}
......@@ -285,6 +285,11 @@ namespace Nektar
STD_REGIONS_EXPORT virtual void v_SVVLaplacianFilter(Array<OneD, NekDouble> &array,const StdMatrixKey &mkey);
STD_REGIONS_EXPORT virtual void v_ExponentialFilter(
Array<OneD, NekDouble> &array,
const NekDouble alpha,
const NekDouble exponent,
const NekDouble cutoff);
};
typedef std::shared_ptr<StdHexExp> StdHexExpSharedPtr;
......
......@@ -47,13 +47,14 @@ namespace Nektar
CFLtester::CFLtester(
const LibUtilities::SessionReaderSharedPtr& pSession)
: UnsteadySystem(pSession)
: UnsteadySystem(pSession),
AdvectionSystem(pSession)
{
}
void CFLtester::v_InitObject()
{
UnsteadySystem::v_InitObject();
AdvectionSystem::v_InitObject();
m_velocity = Array<OneD, Array<OneD, NekDouble> >(m_spacedim);
std::vector<std::string> vel;
......@@ -64,6 +65,60 @@ namespace Nektar
GetFunction( "AdvectionVelocity")->Evaluate(vel, m_velocity);
// Type of advection class to be used
switch(m_projectionType)
{
// Discontinuous field
case MultiRegions::eDiscontinuous:
{
// Do not forwards transform initial condition
m_homoInitialFwd = false;
// Define the normal velocity fields
if (m_fields[0]->GetTrace())
{
m_traceVn = Array<OneD, NekDouble>(GetTraceNpoints());
}
string advName;
string riemName;
m_session->LoadSolverInfo(
"AdvectionType", advName, "WeakDG");
m_advObject = SolverUtils::
GetAdvectionFactory().CreateInstance(advName, advName);
m_advObject->SetFluxVector(
&CFLtester::GetFluxVector, this);
m_session->LoadSolverInfo(
"UpwindType", riemName, "Upwind");
m_riemannSolver = SolverUtils::
GetRiemannSolverFactory().CreateInstance(riemName);
m_riemannSolver->SetScalar(
"Vn", &CFLtester::GetNormalVelocity, this);
m_advObject->SetRiemannSolver(m_riemannSolver);
m_advObject->InitObject(m_session, m_fields);
break;
}
// Continuous field
case MultiRegions::eGalerkin:
case MultiRegions::eMixed_CG_Discontinuous:
{
string advName;
m_session->LoadSolverInfo(
"AdvectionType", advName, "NonConservative");
m_advObject = SolverUtils::
GetAdvectionFactory().CreateInstance(advName, advName);
m_advObject->SetFluxVector(
&CFLtester::GetFluxVector, this);
break;
}
default:
{
ASSERTL0(false, "Unsupported projection type.");
break;
}
}
if (m_explicitAdvection)
{
m_ode.DefineOdeRhs (&CFLtester::DoOdeRhs, this);
......@@ -84,47 +139,14 @@ namespace Nektar
Array<OneD, Array<OneD, NekDouble> >&outarray,
const NekDouble time)
{
int j;
int nvariables = inarray.num_elements();
int npoints = GetNpoints();
switch (m_projectionType)
m_advObject->Advect(nvariables, m_fields, m_velocity, inarray,
outarray, time);
for(int j = 0; j < nvariables; ++j)
{
case MultiRegions::eDiscontinuous:
{
int ncoeffs = inarray[0].num_elements();
Array<OneD, Array<OneD, NekDouble> > WeakAdv(nvariables);
WeakAdv[0] = Array<OneD, NekDouble>(ncoeffs*nvariables);
for(j = 1; j < nvariables; ++j)
{
WeakAdv[j] = WeakAdv[j-1] + ncoeffs;
}
WeakDGAdvection(inarray, WeakAdv, true, true);
for(j = 0; j < nvariables; ++j)
{
m_fields[j]->MultiplyByElmtInvMass(WeakAdv[j], WeakAdv[j]);
m_fields[j]->BwdTrans(WeakAdv[j],outarray[j]);
Vmath::Neg(npoints,outarray[j],1);
}
break;
}
case MultiRegions::eGalerkin:
case MultiRegions::eMixed_CG_Discontinuous:
{
// Calculate -V\cdot Grad(u);
for(j = 0; j < nvariables; ++j)
{
AdvectionNonConservativeForm(m_velocity,
inarray[j],
outarray[j]);
Vmath::Neg(npoints, outarray[j], 1);
}
break;
}
Vmath::Neg(npoints,outarray[j],1);
}
}
......@@ -133,7 +155,6 @@ namespace Nektar
Array<OneD, Array<OneD, NekDouble> >&outarray,
const NekDouble time)
{
int j;
int nvariables = inarray.num_elements();
SetBoundaryConditions(time);
......@@ -144,7 +165,7 @@ namespace Nektar
// Just copy over array
int npoints = GetNpoints();
for(j = 0; j < nvariables; ++j)
for(int j = 0; j < nvariables; ++j)
{
Vmath::Vcopy(npoints,inarray[j],1,outarray[j],1);
}
......@@ -155,7 +176,7 @@ namespace Nektar
{
Array<OneD, NekDouble> coeffs(m_fields[0]->GetNcoeffs());
for(j = 0; j < nvariables; ++j)
for(int j = 0; j < nvariables; ++j)
{
m_fields[j]->FwdTrans(inarray[j],coeffs);
m_fields[j]->BwdTrans_IterPerExp(coeffs,outarray[j]);
......@@ -168,60 +189,56 @@ namespace Nektar
}
}
void CFLtester::v_GetFluxVector(
const int i,
Array<OneD, Array<OneD, NekDouble> > &physfield,
Array<OneD, Array<OneD, NekDouble> > &flux)
/**
* @brief Get the normal velocity
*/
Array<OneD, NekDouble> &CFLtester::GetNormalVelocity()
{
ASSERTL1(flux.num_elements() == m_velocity.num_elements(),
"Dimension of flux array and velocity array do not match");
// Number of trace (interface) points
int nTracePts = GetTraceNpoints();
for (int j = 0; j < flux.num_elements(); ++j)
// Auxiliary variable to compute the normal velocity
Array<OneD, NekDouble> tmp(nTracePts);
// Reset the normal velocity
Vmath::Zero(nTracePts, m_traceVn, 1);
for (int i = 0; i < m_velocity.num_elements(); ++i)
{
Vmath::Vmul(GetNpoints(),
physfield[i], 1,
m_velocity[j], 1,
flux[j], 1);
m_fields[0]->ExtractTracePhys(m_velocity[i], tmp);
Vmath::Vvtvp(nTracePts,
m_traceNormals[i], 1,
tmp, 1,
m_traceVn, 1,
m_traceVn, 1);
}
return m_traceVn;
}
void CFLtester::v_NumericalFlux(
Array<OneD, Array<OneD, NekDouble> > &physfield,
Array<OneD, Array<OneD, NekDouble> > &numflux)
void CFLtester::GetFluxVector(
const Array<OneD, Array<OneD, NekDouble> > &physfield,
Array<OneD, Array<OneD, Array<OneD, NekDouble> > > &flux)
{
int i;
int nTraceNumPoints = GetTraceNpoints();
int nvel = m_spacedim;
Array<OneD, NekDouble > Fwd(nTraceNumPoints);
Array<OneD, NekDouble > Bwd(nTraceNumPoints);
Array<OneD, NekDouble > Vn (nTraceNumPoints,0.0);
ASSERTL1(flux[0].num_elements() == m_velocity.num_elements(),
"Dimension of flux array and velocity array do not match");
// Get Edge Velocity - Could be stored if time independent
for (i = 0; i < nvel; ++i)
{
m_fields[0]->ExtractTracePhys(m_velocity[i], Fwd);
Vmath::Vvtvp(nTraceNumPoints,
m_traceNormals[i], 1,
Fwd, 1, Vn, 1, Vn, 1);
}
int nq = physfield[0].num_elements();
for (i = 0; i < numflux.num_elements(); ++i)
for (int i = 0; i < flux.num_elements(); ++i)
{
m_fields[i]->GetFwdBwdTracePhys(physfield[i], Fwd, Bwd);
//evaulate upwinded m_fields[i]
m_fields[i]->GetTrace()->Upwind(Vn, Fwd, Bwd, numflux[i]);
// calculate m_fields[i]*Vn
Vmath::Vmul(nTraceNumPoints, numflux[i], 1, Vn, 1, numflux[i], 1);
for (int j = 0; j < flux[0].num_elements(); ++j)
{
Vmath::Vmul(nq, physfield[i], 1, m_velocity[j], 1,
flux[i][j], 1);
}
}
}
void CFLtester::v_GenerateSummary(SummaryList& s)
{
UnsteadySystem::v_GenerateSummary(s);
AdvectionSystem::v_GenerateSummary(s);
}
......
......@@ -36,7 +36,8 @@
#ifndef NEKTAR_SOLVERS_ADRSOLVER_EQUATIONSYSTEMS_CFLTESTER_H
#define NEKTAR_SOLVERS_ADRSOLVER_EQUATIONSYSTEMS_CFLTESTER_H
#include <SolverUtils/UnsteadySystem.h>
#include <SolverUtils/AdvectionSystem.h>
#include <SolverUtils/RiemannSolvers/RiemannSolver.h>
using namespace Nektar::SolverUtils;
......@@ -78,7 +79,7 @@ namespace Nektar
{724.129308, 1297.665239, 1890.995890, 2477.394447, 3025.267650, 3552.997300, 4093.052410, 4609.237151, 5136.747279, 5661.792457, 6397.687000, 7383.760000, 8431.678000, 9540.598000},
{732.175780, 1305.179484, 1897.615616, 2483.199183, 3030.256925, 3557.347822, 4096.823380, 4612.552660, 5139.604517, 5664.300842, 6771.812000, 7388.760000, 8436.678000, 9546.598000}};
class CFLtester : public UnsteadySystem
class CFLtester : public SolverUtils::AdvectionSystem
{
public:
friend class MemoryManager<CFLtester>;
......@@ -93,8 +94,9 @@ namespace Nektar
virtual ~CFLtester();
protected:
Array<OneD, Array<OneD, NekDouble> > m_velocity;
SolverUtils::RiemannSolverSharedPtr m_riemannSolver;
Array<OneD, Array<OneD, NekDouble> > m_velocity;
Array<OneD, NekDouble> m_traceVn;
CFLtester(const LibUtilities::SessionReaderSharedPtr& pSession);
......@@ -106,11 +108,14 @@ namespace Nektar
Array<OneD, Array<OneD, NekDouble> > &outarray,
const NekDouble time);
virtual void v_InitObject();
/// Get the normal velocity