Commit f57dfa9d authored by Chris Cantwell's avatar Chris Cantwell
Browse files

Merge branch 'fix/Stability3D' into 'master'

Tidy linearised advection and fix full 3D

This MR simplifies the logic of `LinearisedAdvection::v_Advect`, reducing the dependency on the dimension and handling all components in a single loop. In doing so, I think this fixes the full 3D case, which was not working. Also, the gradient of the base flow is now computed in a separate function and stored in `m_gradBase`, avoiding the need to recompute this every time-step for steady base flows.

The adjoint advection is now derived from LinearisedAdvection, eliminating a lot of duplicate code. It only needs to redefine v_Advect, which now follows a similar structure to the new implementation from LinearisedAdvection.

See merge request !708
parents 41212019 8d010c7f
......@@ -47,6 +47,7 @@ v4.4.0
**IncNavierStokesSolver:**
- Add ability to simulate additional scalar fields (!624)
- Improve performance when using homogeneous dealiasing (!622)
- Fix linearised advection for full 3D cases (!708)
**FieldConvert:**
- Allow equi-spaced output for 1D and 2DH1D fields (!613)
......
......@@ -125,7 +125,8 @@ void Advection::v_InitObject(
*
*/
void Advection::v_SetBaseFlow(
const Array<OneD, Array<OneD, NekDouble> > &inarray)
const Array<OneD, Array<OneD, NekDouble> > &inarray,
const Array<OneD, MultiRegions::ExpListSharedPtr> &fields)
{
ASSERTL0(false,
"A baseflow is not appropriate for this advection type.");
......
......@@ -123,9 +123,11 @@ public:
*
* @param inarray Vector to use as baseflow
*/
inline void SetBaseFlow(const Array<OneD, Array<OneD, NekDouble> >& inarray)
inline void SetBaseFlow(
const Array<OneD, Array<OneD, NekDouble> >& inarray,
const Array<OneD, MultiRegions::ExpListSharedPtr> &fields)
{
v_SetBaseFlow(inarray);
v_SetBaseFlow(inarray, fields);
}
protected:
......@@ -155,7 +157,8 @@ protected:
/// Overrides the base flow used during linearised advection
SOLVER_UTILS_EXPORT virtual void v_SetBaseFlow(
const Array<OneD, Array<OneD, NekDouble> > &inarray);
const Array<OneD, Array<OneD, NekDouble> > &inarray,
const Array<OneD, MultiRegions::ExpListSharedPtr> &fields);
};
/// A shared pointer to an Advection object.
......
......@@ -226,7 +226,7 @@ void DriverSteadyState::v_Execute(ostream &out)
m_equ[m_nequ - 1]->Checkpoint_BaseFlow(m_Check_BaseFlow);
m_Check_BaseFlow++;
A->GetAdvObject()->SetBaseFlow(q0);
A->GetAdvObject()->SetBaseFlow(q0,m_equ[0]->UpdateFields());
DriverModifiedArnoldi::v_Execute(out);
if (m_comm->GetRank() == 0)
......@@ -262,7 +262,7 @@ void DriverSteadyState::v_Execute(ostream &out)
m_equ[m_nequ - 1]->Checkpoint_BaseFlow(m_Check_BaseFlow);
m_Check_BaseFlow++;
A->GetAdvObject()->SetBaseFlow(q0);
A->GetAdvObject()->SetBaseFlow(q0,m_equ[0]->UpdateFields());
DriverModifiedArnoldi::v_Execute(out);
if (m_comm->GetRank() == 0)
......
......@@ -36,148 +36,43 @@
#ifndef NEKTAR_SOLVERS_ADJOINTADVECTION_H
#define NEKTAR_SOLVERS_ADJOINTADVECTION_H
#include <SolverUtils/Advection/Advection.h>
#include <LibUtilities/FFT/NektarFFT.h>
#include <IncNavierStokesSolver/AdvectionTerms/LinearisedAdvection.h>
namespace Nektar
{
/// Advection for the adjoint form of the linearised Navier-Stokes equations
class AdjointAdvection: public SolverUtils::Advection
class AdjointAdvection: public LinearisedAdvection
{
enum FloquetMatType
{
eForwardsCoeff,
eForwardsPhys
};
/// A map between matrix keys and their associated block matrices.
typedef std::map< FloquetMatType, DNekBlkMatSharedPtr> FloquetBlockMatrixMap;
/// A shared pointer to a BlockMatrixMap.
typedef boost::shared_ptr<FloquetBlockMatrixMap> FloquetBlockMatrixMapShPtr;
public:
friend class MemoryManager<AdjointAdvection>;
/// Creates an instance of this class
static SolverUtils::AdvectionSharedPtr create(std::string) {
return MemoryManager<AdjointAdvection>::AllocateSharedPtr();
}
/// Name of class
static std::string className;
protected:
LibUtilities::SessionReaderSharedPtr m_session;
MultiRegions::ProjectionType m_projectionType;
int m_spacedim;
int m_expdim;
/// Storage for base flow
Array<OneD, Array<OneD, NekDouble> > m_baseflow;
//number of slices
int m_slices;
//period length
NekDouble m_period;
//interpolation vector
Array<OneD, Array<OneD, NekDouble> > m_interp;
//auxiliary variables
LibUtilities::NektarFFTSharedPtr m_FFT;
Array<OneD,NekDouble> m_tmpIN;
Array<OneD,NekDouble> m_tmpOUT;
bool m_useFFTW;
/// flag to determine if use single mode or not
bool m_SingleMode;
/// flag to determine if use half mode or not
bool m_HalfMode;
/// flag to determine if use multiple mode or not
bool m_MultipleModes;
bool m_homogen_dealiasing;
MultiRegions::CoeffState m_CoeffState;
DNekBlkMatSharedPtr GetFloquetBlockMatrix(
FloquetMatType mattype,
bool UseContCoeffs = false) const;
DNekBlkMatSharedPtr GenFloquetBlockMatrix(
FloquetMatType mattype,
bool UseContCoeffs = false) const;
FloquetBlockMatrixMapShPtr m_FloquetBlockMat;
AdjointAdvection();
virtual ~AdjointAdvection();
virtual void v_InitObject(
LibUtilities::SessionReaderSharedPtr pSession,
Array<OneD, MultiRegions::ExpListSharedPtr> pFields);
virtual void v_Advect(
const int nConvectiveFields,
const Array<OneD, MultiRegions::ExpListSharedPtr> &fields,
const Array<OneD, Array<OneD, NekDouble> > &advVel,
const Array<OneD, Array<OneD, NekDouble> > &inarray,
Array<OneD, Array<OneD, NekDouble> > &outarray,
const NekDouble &time,
const Array<OneD, Array<OneD, NekDouble> > &pFwd = NullNekDoubleArrayofArray,
const Array<OneD, Array<OneD, NekDouble> > &pBwd = NullNekDoubleArrayofArray);
virtual void v_SetBaseFlow(
const Array<OneD, Array<OneD, NekDouble> > &inarray);
void UpdateBase(
const NekDouble m_slices,
const Array<OneD, const NekDouble> &inarray,
Array<OneD, NekDouble> &outarray,
const NekDouble m_time,
const NekDouble m_period);
void DFT(
const std::string file,
Array<OneD, MultiRegions::ExpListSharedPtr> &pFields,
const NekDouble m_slices);
/// Import Base flow
void ImportFldBase(
std::string pInfile,
Array<OneD, MultiRegions::ExpListSharedPtr> &pFields,
int slice);
private:
///Parameter for homogeneous expansions
enum HomogeneousType
{
eHomogeneous1D,
eHomogeneous2D,
eHomogeneous3D,
eNotHomogeneous
};
/// flag to determine if use or not the FFT for transformations
bool m_useFFT;
enum HomogeneousType m_HomogeneousType;
NekDouble m_LhomX; ///< physical length in X direction (if homogeneous)
NekDouble m_LhomY; ///< physical length in Y direction (if homogeneous)
NekDouble m_LhomZ; ///< physical length in Z direction (if homogeneous)
int m_npointsX; ///< number of points in X direction (if homogeneous)
int m_npointsY; ///< number of points in Y direction (if homogeneous)
int m_npointsZ; ///< number of points in Z direction (if homogeneous)
int m_HomoDirec; ///< number of homogenous directions
int m_NumMode; ///< Mode to use in case of single mode analysis
SpatialDomains::BoundaryConditionsSharedPtr m_boundaryConditions;
public:
friend class MemoryManager<AdjointAdvection>;
/// Creates an instance of this class
static SolverUtils::AdvectionSharedPtr create(std::string)
{
return MemoryManager<AdjointAdvection>::AllocateSharedPtr();
}
/// Name of class
static std::string className;
protected:
AdjointAdvection();
virtual ~AdjointAdvection();
virtual void v_Advect(
const int nConvectiveFields,
const Array<OneD, MultiRegions::ExpListSharedPtr> &fields,
const Array<OneD, Array<OneD, NekDouble> > &advVel,
const Array<OneD, Array<OneD, NekDouble> > &inarray,
Array<OneD, Array<OneD, NekDouble> > &outarray,
const NekDouble &time,
const Array<OneD, Array<OneD, NekDouble> > &pFwd =
NullNekDoubleArrayofArray,
const Array<OneD, Array<OneD, NekDouble> > &pBwd =
NullNekDoubleArrayofArray);
};
}
......
......@@ -36,16 +36,6 @@
#include <IncNavierStokesSolver/AdvectionTerms/LinearisedAdvection.h>
#include <StdRegions/StdSegExp.h>
#include <MultiRegions/ContField3DHomogeneous1D.h>
#include <MultiRegions/ContField3DHomogeneous2D.h>
#include <MultiRegions/ContField1D.h>
#include <MultiRegions/ContField2D.h>
#include <MultiRegions/ContField3D.h>
#include <MultiRegions/DisContField1D.h>
#include <MultiRegions/DisContField2D.h>
#include <MultiRegions/DisContField3DHomogeneous1D.h>
#include <MultiRegions/DisContField3DHomogeneous2D.h>
using namespace std;
namespace Nektar
......@@ -166,6 +156,17 @@ void LinearisedAdvection::v_InitObject(
m_baseflow[i] = Array<OneD, NekDouble>(pFields[i]->GetTotPoints(), 0.0);
}
int nBaseDerivs = (m_halfMode || m_singleMode) ? 2 : m_spacedim;
m_gradBase = Array<OneD, Array<OneD, NekDouble> >(nvar*nBaseDerivs);
for (int i = 0; i < nvar; ++i)
{
for (int j = 0; j < nBaseDerivs; ++j)
{
m_gradBase[i*nBaseDerivs + j ] = Array<OneD, NekDouble>
(pFields[i]->GetTotPoints(), 0.0);
}
}
ASSERTL0(m_session->DefinesFunction("BaseFlow"),
"Base flow must be defined for linearised forms.");
string file = m_session->GetFunctionFilename("BaseFlow", 0);
......@@ -177,6 +178,9 @@ void LinearisedAdvection::v_InitObject(
m_session->LoadParameter("N_slices",m_slices);
if(m_slices>1)
{
ASSERTL0(m_session->GetFunctionType("BaseFlow", 0)
== LibUtilities::eFunctionTypeFile,
"Base flow should be a sequence of files.");
DFT(file,pFields,m_slices);
}
else
......@@ -219,6 +223,11 @@ void LinearisedAdvection::v_InitObject(
}
}
for (int i = 0; i < nvar; ++i)
{
UpdateGradBase(i, pFields[i]);
}
if(m_session->DefinesParameter("period"))
{
m_period=m_session->GetParameter("period");
......@@ -250,8 +259,10 @@ void LinearisedAdvection::v_Advect(
ASSERTL1(nConvectiveFields == inarray.num_elements(),
"Number of convective fields and Inarray are not compatible");
int nPointsTot = fields[0]->GetNpoints();
int ndim = advVel.num_elements();
int nPointsTot = fields[0]->GetNpoints();
int ndim = advVel.num_elements();
int nBaseDerivs = (m_halfMode || m_singleMode) ? 2 : m_spacedim;
int nDerivs = (m_halfMode) ? 2 : m_spacedim;
Array<OneD, Array<OneD, NekDouble> > velocity(ndim);
for(int i = 0; i < ndim; ++i)
......@@ -267,294 +278,88 @@ void LinearisedAdvection::v_Advect(
}
}
Array<OneD, NekDouble> grad0,grad1,grad2;
// Evaluation of the gradient of each component of the base flow
// \nabla U
Array<OneD, NekDouble> grad_base_u0,grad_base_u1,grad_base_u2;
// \nabla V
Array<OneD, NekDouble> grad_base_v0,grad_base_v1,grad_base_v2;
// \nabla W
Array<OneD, NekDouble> grad_base_w0,grad_base_w1,grad_base_w2;
grad0 = Array<OneD, NekDouble> (nPointsTot);
grad_base_u0 = Array<OneD, NekDouble> (nPointsTot);
grad_base_v0 = Array<OneD, NekDouble> (nPointsTot);
grad_base_w0 = Array<OneD, NekDouble> (nPointsTot);
Array<OneD, Array<OneD, NekDouble> > grad (nDerivs);
for( int i = 0; i < nDerivs; ++i)
{
grad[i] = Array<OneD, NekDouble> (nPointsTot);
}
// Evaluation of the base flow for periodic cases
if (m_slices > 1)
{
ASSERTL0(m_session->GetFunctionType("BaseFlow", 0)
== LibUtilities::eFunctionTypeFile,
"Base flow should be a sequence of files.");
for (int i = 0; i < ndim; ++i)
{
UpdateBase(m_slices, m_interp[i], m_baseflow[i],
time, m_period);
UpdateGradBase(i, fields[i]);
}
}
//Evaluate the linearised advection term
switch(ndim)
for( int i = 0; i < ndim; ++i)
{
case 1: // 1D
ASSERTL0(false,"Not set up for 1D");
break;
case 2: //2D
// Calculate gradient
switch(nDerivs)
{
grad1 = Array<OneD, NekDouble> (nPointsTot);
grad_base_u1 = Array<OneD, NekDouble> (nPointsTot);
grad_base_v1 = Array<OneD, NekDouble> (nPointsTot);
//Derivates of the base flow
fields[0]-> PhysDeriv(m_baseflow[0], grad_base_u0, grad_base_u1);
fields[0]-> PhysDeriv(m_baseflow[1], grad_base_v0, grad_base_v1);
//x-equation
fields[0]->PhysDeriv(inarray[0],grad0,grad1);
// Evaluate U du'/dx
Vmath::Vmul (nPointsTot,grad0,1,m_baseflow[0],1,outarray[0],1);
//Evaluate U du'/dx+ V du'/dy
Vmath::Vvtvp(nPointsTot,grad1,1,m_baseflow[1],1,outarray[0],1,
outarray[0],1);
//Evaluate (U du'/dx+ V du'/dy)+u' dU/dx
Vmath::Vvtvp(nPointsTot,grad_base_u0,1,velocity[0],1,outarray[0],1,
outarray[0],1);
//Evaluate (U du'/dx+ V du'/dy +u' dU/dx)+v' dU/dy
Vmath::Vvtvp(nPointsTot,grad_base_u1,1,velocity[1],1,outarray[0],1,
outarray[0],1);
Vmath::Neg(nPointsTot,outarray[0],1);
//y-equation
fields[0]->PhysDeriv(inarray[1],grad0,grad1);
// Evaluate U dv'/dx
Vmath::Vmul (nPointsTot,grad0,1,m_baseflow[0],1,outarray[1],1);
//Evaluate U dv'/dx+ V dv'/dy
Vmath::Vvtvp(nPointsTot,grad1,1,m_baseflow[1],1,outarray[1],1,
outarray[1],1);
//Evaluate (U dv'/dx+ V dv'/dy)+u' dV/dx
Vmath::Vvtvp(nPointsTot,grad_base_v0,1,velocity[0],1,outarray[1],1,
outarray[1],1);
//Evaluate (U dv'/dx+ V dv'/dy +u' dv/dx)+v' dV/dy
Vmath::Vvtvp(nPointsTot,grad_base_v1,1,velocity[1],1,outarray[1],1,
outarray[1],1);
Vmath::Neg(nPointsTot,outarray[1],1);
}
break;
case 3: //3D
{
grad1 = Array<OneD, NekDouble> (nPointsTot);
grad2 = Array<OneD, NekDouble> (nPointsTot);
grad_base_u1 = Array<OneD, NekDouble> (nPointsTot);
grad_base_v1 = Array<OneD, NekDouble> (nPointsTot);
grad_base_w1 = Array<OneD, NekDouble> (nPointsTot,0.0);
grad_base_u2 = Array<OneD, NekDouble> (nPointsTot,0.0);
grad_base_v2 = Array<OneD, NekDouble> (nPointsTot,0.0);
grad_base_w2 = Array<OneD, NekDouble> (nPointsTot,0.0);
// note base flow should be moved to GetBaseFlow method
if(m_halfMode) // can assume W = 0 and d/dz == 0
case 1:
{
// note base flow should be moved to GetBaseFlow method
fields[0]->PhysDeriv(m_baseflow[0],
grad_base_u0, grad_base_u1);
fields[0]->PhysDeriv(m_baseflow[1],
grad_base_v0, grad_base_v1);
fields[i]->PhysDeriv(inarray[i], grad[0]);
}
else if(m_singleMode) // single mode where d/dz = 0
{
fields[0]->PhysDeriv(m_baseflow[0], grad_base_u0,
grad_base_u1);
fields[0]->PhysDeriv(m_baseflow[1], grad_base_v0,
grad_base_v1);
fields[0]->PhysDeriv(m_baseflow[2], grad_base_w0,
grad_base_w1);
}
else if(m_multipleModes)
{
// Differentiate base flow in physical space
bool oldwavespace = fields[0]->GetWaveSpace();
fields[0]->SetWaveSpace(false);
fields[0]->PhysDeriv(m_baseflow[0], grad_base_u0,
grad_base_u1, grad_base_u2);
fields[0]->PhysDeriv(m_baseflow[1], grad_base_v0,
grad_base_v1, grad_base_v2);
fields[0]->PhysDeriv(m_baseflow[2], grad_base_w0,
grad_base_w1, grad_base_w2);
fields[0]->SetWaveSpace(oldwavespace);
}
else
{
ASSERTL0(false, "ERROR: Must be one of half, single or multiple modes");
}
//x-equation
//
// Could probably clean up independent looping if clean up
// base flow derivative evaluation
if(m_multipleModes)
{
fields[0]->PhysDeriv(inarray[0], grad0, grad1, grad2);
// transform gradients into physical Fourier space
fields[0]->HomogeneousBwdTrans(grad0, grad0);
fields[0]->HomogeneousBwdTrans(grad1, grad1);
fields[0]->HomogeneousBwdTrans(grad2, grad2);
}
else
{
if(m_halfMode) //W = 0 so no need for d/dz
{
fields[0]->PhysDeriv(inarray[0], grad0, grad1);
}
else // Single mode
{
fields[0]->PhysDeriv(inarray[0], grad0, grad1, grad2);
}
}
//Evaluate: U du'/dx
Vmath::Vmul (nPointsTot, grad0, 1, m_baseflow[0], 1,
outarray[0], 1);
//Evaluate and add: V du'/dy
Vmath::Vvtvp(nPointsTot, grad1, 1, m_baseflow[1], 1,
outarray[0], 1, outarray[0], 1);
//Evaluate and add: u' dU/dx
Vmath::Vvtvp(nPointsTot,grad_base_u0,1,velocity[0],1,
outarray[0],1,outarray[0],1);
//Evaluate and add: v' dU/dy
Vmath::Vvtvp(nPointsTot,grad_base_u1,1,velocity[1],1,
outarray[0],1,outarray[0],1);
if(!m_halfMode)
{
//Evaluate an add W du'/dz
Vmath::Vvtvp(nPointsTot, grad2, 1, m_baseflow[2],
1,outarray[0], 1, outarray[0], 1);
}
if(m_multipleModes)
{
//Evaluate and add w' dU/dz
Vmath::Vvtvp(nPointsTot,grad_base_u2,1,
velocity[2],1,outarray[0],1,outarray[0],1);
fields[0]->HomogeneousFwdTrans(outarray[0],outarray[0]);
}
Vmath::Neg(nPointsTot,outarray[0],1);
//y-equation------------
if(m_multipleModes)
{
fields[0]->PhysDeriv(inarray[1], grad0, grad1, grad2);
// transform gradients into physical fouier space
fields[0]->HomogeneousBwdTrans(grad0, grad0);
fields[0]->HomogeneousBwdTrans(grad1, grad1);
fields[0]->HomogeneousBwdTrans(grad2, grad2);
}
else
{
if(m_halfMode) //W = 0 so no need for d/dz
{
fields[0]->PhysDeriv(inarray[1], grad0, grad1);
}
else // Single mode
{
fields[0]->PhysDeriv(inarray[1], grad0, grad1, grad2);
}
}
//Evaluate U dv'/dx
Vmath::Vmul (nPointsTot, grad0, 1, m_baseflow[0], 1,
outarray[1], 1);
//Evaluate V dv'/dy
Vmath::Vvtvp(nPointsTot, grad1, 1, m_baseflow[1], 1,
outarray[1], 1, outarray[1], 1);
//Evaluate u' dV/dx
Vmath::Vvtvp(nPointsTot,grad_base_v0,1,velocity[0],1
,outarray[1],1,outarray[1],1);
//Evaluate v' dV/dy
Vmath::Vvtvp(nPointsTot,grad_base_v1,1,velocity[1],1,
outarray[1],1,outarray[1],1);
if(!m_halfMode)
{
//Evaluate W du'/dz
Vmath::Vvtvp(nPointsTot, grad2, 1, m_baseflow[2], 1,
outarray[1], 1, outarray[1], 1);
}
if(m_multipleModes)
{
//Evaluate w' dV/dz
Vmath::Vvtvp(nPointsTot,grad_base_v2,1,velocity[2],1,
outarray[1],1,outarray[1],1);
fields[0]->HomogeneousFwdTrans(outarray[1],outarray[1]);
}
Vmath::Neg(nPointsTot,outarray[1],1);
//z-equation ------------
if(m_multipleModes)
break;
case 2:
{
fields[0]->PhysDeriv(inarray[2], grad0, grad1, grad2);
// transform gradients into physical fouier space
fields[0]->HomogeneousBwdTrans(grad0, grad0);
fields[0]->HomogeneousBwdTrans(grad1, grad1);
fields[0]->HomogeneousBwdTrans(grad2, grad2);
fields[i]->PhysDeriv(inarray[i], grad[0], grad[1]);
}
else
break;
case 3:
{
if(m_halfMode) //W = 0 so no need for d/dz