Skip to content

GitLab

Commit 9556cb2a authored by Dave Moxey's avatar Dave Moxey
Browse files

Merge branch 'feature/cfs_cleaning' of localhost:nektar

parents a7c8cd9d 29f1a54d
Expand all Hide whitespace changes
Inline Side-by-side
library/MultiRegions/ExpList0D.cpp
View file @ 9556cb2a
......@@ -56,26 +56,20 @@ namespace Nektar
{
}
ExpList0D::ExpList0D(const SpatialDomains::VertexComponentSharedPtr &m_geom):
ExpList()
ExpList0D::ExpList0D(const SpatialDomains::VertexComponentSharedPtr &m_geom):
ExpList()
{
m_point = MemoryManager<LocalRegions::PointExp>::AllocateSharedPtr(m_geom);
m_ncoeffs = 1;
m_npoints = 1;
// Set up m_coeffs, m_phys.
m_coeffs = Array<OneD, NekDouble>(m_ncoeffs);
m_phys = Array<OneD, NekDouble>(m_npoints);
m_coeffs[0] = m_point->GetCoeff(0);
m_phys[0] = m_point->GetPhys(0);
m_coeffs = m_point->UpdateCoeffs();
m_phys = m_point->UpdatePhys();
m_point = MemoryManager<LocalRegions::PointExp>::AllocateSharedPtr(m_geom);
m_ncoeffs = 1;
m_npoints = 1;
// Set up m_coeffs, m_phys.
m_coeffs = Array<OneD, NekDouble>(m_ncoeffs);
m_phys = Array<OneD, NekDouble>(m_npoints);
}
/**
/**
* Store expansions for the trace space expansions used in
* DisContField1D.
*
......@@ -92,20 +86,20 @@ namespace Nektar
* instead of just normal segment expansions.
*/
ExpList0D::ExpList0D(
const Array<OneD,const ExpListSharedPtr> &bndConstraint,
const Array<OneD, const SpatialDomains
::BoundaryConditionShPtr> &bndCond,
const StdRegions::StdExpansionVector &locexp,
const SpatialDomains::MeshGraphSharedPtr &graph1D,
const map<int,int> &periodicVertices,
const bool DeclareCoeffPhysArrays):
ExpList()
const Array<OneD, const ExpListSharedPtr> &bndConstraint,
const Array<OneD, const SpatialDomains::BoundaryConditionShPtr>
&bndCond,
const StdRegions::StdExpansionVector &locexp,
const SpatialDomains::MeshGraphSharedPtr &graph1D,
const map<int,int> &periodicVertices,
const bool DeclareCoeffPhysArrays)
: ExpList()
{
int i,j,cnt,id, elmtid=0;
map<int,int> EdgeDone;
map<int,int> NormalSet;
SpatialDomains::VertexComponentSharedPtr PointGeom;
SpatialDomains::VertexComponentSharedPtr PointGeom;
LocalRegions::PointExpSharedPtr Point;
// First loop over boundary conditions to renumber Dirichlet boundaries
......@@ -221,35 +215,35 @@ namespace Nektar
{
}
void ExpList0D::v_GetCoords(NekDouble &x, NekDouble &y, NekDouble &z)
void ExpList0D::v_GetCoords(NekDouble &x, NekDouble &y, NekDouble &z)
{
m_point->GetCoords(x,y,z);
}
m_point->GetCoords(x,y,z);
}
void ExpList0D::v_GetCoord(Array<OneD,NekDouble> &coords)
void ExpList0D::v_GetCoord(Array<OneD,NekDouble> &coords)
{
m_point->GetCoords(coords);
}
m_point->GetCoords(coords);
}
void ExpList0D::v_SetCoeff(NekDouble val)
void ExpList0D::v_SetCoeff(NekDouble val)
{
m_point->SetCoeff(val);
}
m_coeffs[0] = val;
}
void ExpList0D::v_SetPhys(NekDouble val)
void ExpList0D::v_SetPhys(NekDouble val)
{
m_point->SetPhys(val);
}
m_phys[0] = val;
}
const SpatialDomains::VertexComponentSharedPtr &ExpList0D::v_GetGeom(void) const
{
return m_point->GetGeom();
}
const SpatialDomains::VertexComponentSharedPtr &ExpList0D::v_GetGeom(void) const
{
return m_point->GetGeom();
}
const SpatialDomains::VertexComponentSharedPtr &ExpList0D::v_GetVertex(void) const
{
return m_point->GetVertex();
}
const SpatialDomains::VertexComponentSharedPtr &ExpList0D::v_GetVertex(void) const
{
return m_point->GetVertex();
}
/**
* For each local element, copy the normals stored in the element list
......
library/SolverUtils/AdvectionFR.cpp
View file @ 9556cb2a
This diff is collapsed.
library/SolverUtils/EquationSystem.cpp
View file @ 9556cb2a
......@@ -60,15 +60,16 @@ namespace Nektar
/**
* @class EquationSystem
*
* This class is a base class for all solver implementations. It provides
* the underlying generic functionality and interface for solving equations.
* This class is a base class for all solver implementations. It
* provides the underlying generic functionality and interface for
* solving equations.
*
* To solve a steady-state equation, create a derived class from this class
* and reimplement the virtual functions to provide custom implementation
* for the problem.
* To solve a steady-state equation, create a derived class from this
* class and reimplement the virtual functions to provide custom
* implementation for the problem.
*
* To solve unsteady problems, derive from the UnsteadySystem class instead
* which provides general time integration.
* To solve unsteady problems, derive from the UnsteadySystem class
* instead which provides general time integration.
*/
EquationSystemFactory& GetEquationSystemFactory()
{
......@@ -496,12 +497,14 @@ namespace Nektar
}
// Set Default Parameter
m_session->LoadParameter("Time", m_time, 0.0);
m_session->LoadParameter("TimeStep", m_timestep, 0.01);
m_session->LoadParameter("NumSteps", m_steps, 0);
m_session->LoadParameter("IO_CheckSteps", m_checksteps, 0);
m_session->LoadParameter("FinTime", m_fintime, 0);
m_session->LoadParameter("NumQuadPointsError", m_NumQuadPointsError, 0);
m_session->LoadParameter("Time", m_time, 0.0);
m_session->LoadParameter("TimeStep", m_timestep, 0.01);
m_session->LoadParameter("NumSteps", m_steps, 0);
m_session->LoadParameter("IO_CheckSteps", m_checksteps, 0);
m_session->LoadParameter("IO_CheckTime", m_checktime, 0.0);
m_session->LoadParameter("FinTime", m_fintime, 0);
m_session->LoadParameter("NumQuadPointsError",
m_NumQuadPointsError, 0);
// Read in spatial data
int nq = m_fields[0]->GetNpoints();
......@@ -1250,7 +1253,7 @@ namespace Nektar
outname += "_P"+boost::lexical_cast<std::string>(m_comm->GetRank());
}
outname += ".fld";
WriteFld(outname);
WriteFld(outname);
}
/**
......@@ -2046,6 +2049,20 @@ namespace Nektar
out << "\tMax Exp. Order : " << m_fields[0]->EvalBasisNumModesMax()<< endl;
}
if (m_session->DefinesSolverInfo("UpwindType"))
{
std::string UpwindType;
UpwindType = m_session->GetSolverInfo("UpwindType");
if (UpwindType == "Exact")
{
out << "\tRiemann Solver : Exact" <<endl;
}
else if (UpwindType == "Average")
{
out << "\tRiemann Solver : Average" <<endl;
}
}
if (m_session->DefinesSolverInfo("AdvectionType"))
{
std::string AdvectionType;
......@@ -2062,26 +2079,27 @@ namespace Nektar
{
if (AdvectionType == "WeakDG")
{
out << "\tProjection Type : Weak Discontinuous Galerkin" <<endl;
out << "\tProjection Type : Weak Discontinuous Galerkin" <<endl;
}
else if (AdvectionType == "FRDG")
{
out << "\tProjection Type : Flux Reconstruction DG" <<endl;
out << "\tProjection Type : Flux Reconstruction DG" <<endl;
}
else if (AdvectionType == "FRSD")
{
out << "\tProjection Type : Flux Reconstruction SD" <<endl;
out << "\tProjection Type : Flux Reconstruction SD" <<endl;
}
else if (AdvectionType == "FRHU")
{
out << "\tProjection Type : Flux Reconstruction HU" <<endl;
out << "\tProjection Type : Flux Reconstruction HU" <<endl;
}
else if (AdvectionType == "FRc")
{
out << "\tProjection Type : Flux Reconstruction c = c-min" <<endl;
}
else if (AdvectionType == "FRc")
{
out << "\tProjection Type : Flux Reconstruction C" <<endl;
out << "\tProjection Type : Flux Reconstruction c = c-infinity" <<endl;
}
break;
}
......
library/SolverUtils/EquationSystem.h
View file @ 9556cb2a
......@@ -400,33 +400,52 @@ namespace Nektar
SpatialDomains::BoundaryConditionsSharedPtr m_boundaryConditions;
/// Pointer to graph defining mesh.
SpatialDomains::MeshGraphSharedPtr m_graph;
SpatialDomains::SpatialParametersSharedPtr m_spatialParameters;
std::string m_filename; ///< Filename
std::string m_sessionName; ///< Name of the sessions
NekDouble m_time; ///< Continous time
NekDouble m_fintime; ///< Time to be taken during the simulation
NekDouble m_timestep; ///< Time step size
NekDouble m_lambda; ///< Lambda constant in real system if one required
int m_steps; ///< Number of steps to take
int m_checksteps; ///< Number of steps between checkpoints
int m_spacedim; ///< Spatial dimension (> expansion dim)
int m_expdim; ///< Dimension of the expansion
bool m_SingleMode; ///< Flag to determine if use single mode or not.
bool m_HalfMode; ///< Flag to determine if use half mode or not
bool m_MultipleModes; ///< Flag to determine if use multiple mode or not
enum MultiRegions::ProjectionType m_projectionType; ///< Type of projection, i.e. Continuous or Discontinuous
Array<OneD, Array<OneD, NekDouble> > m_traceNormals; ///< Array holding the forward normals
/// Filename.
std::string m_filename;
/// Name of the session.
std::string m_sessionName;
/// Current time of simulation.
NekDouble m_time;
/// Finish time of the simulation.
NekDouble m_fintime;
/// Time step size
NekDouble m_timestep;
/// Lambda constant in real system if one required.
NekDouble m_lambda;
/// Time between checkpoints.
NekDouble m_checktime;
/// Number of steps to take.
int m_steps;
/// Number of steps between checkpoints.
int m_checksteps;
/// Spatial dimension (>= expansion dim).
int m_spacedim;
/// Expansion dimension.
int m_expdim;
/// Flag to determine if single homogeneous mode is used.
bool m_SingleMode;
/// Flag to determine if half homogeneous mode is used.
bool m_HalfMode;
/// Flag to determine if use multiple homogenenous modes are used.
bool m_MultipleModes;
/// Flag to determine if FFT is used for homogeneous transform.
bool m_useFFT;
/// Flag to determine if dealiasing is used for homogeneous
/// simulations.
bool m_dealiasing;
/// Type of projection; e.g continuous or discontinuous.
enum MultiRegions::ProjectionType m_projectionType;
/// Array holding trace normals for DG simulations in the forwards
/// direction.
Array<OneD, Array<OneD, NekDouble> > m_traceNormals;
/// 1 x nvariable x nq
Array<OneD, Array<OneD, Array<OneD,NekDouble> > > m_gradtan;
/// 2 x m_spacedim x nq
Array<OneD, Array<OneD, Array<OneD,NekDouble> > > m_tanbasis;
/// Flag to indicate if the fields should be checked for singularity.
Array<OneD, bool> m_checkIfSystemSingular;
/// Flag to indicate if the fields should be checked for
/// singularity.
Array<OneD, bool> m_checkIfSystemSingular;
/// Number of Quadrature points used to work out the error
int m_NumQuadPointsError;
......@@ -441,8 +460,6 @@ namespace Nektar
eNotHomogeneous
};
bool m_useFFT; ///< flag to determine if use or not the FFT for transformations
bool m_dealiasing; ///< flag to determine if use dealising or not
enum HomogeneousType m_HomogeneousType;
......
library/SolverUtils/UnsteadySystem.cpp
View file @ 9556cb2a
This diff is collapsed.
library/SolverUtils/UnsteadySystem.h
View file @ 9556cb2a
......@@ -50,39 +50,14 @@ namespace Nektar
public:
/// Destructor
SOLVER_UTILS_EXPORT virtual ~UnsteadySystem();
/// Calculate the largest time-step for maintaining a stable
/// solution (using CFL).
SOLVER_UTILS_EXPORT NekDouble GetTimeStep(
const Array<OneD,int> ExpOrder,
const Array<OneD,NekDouble> CFL,
NekDouble timeCFL);
/// Calculate the larger time-step mantaining the problem stable, just for CFLTester equation
/// Calculate the larger time-step mantaining the problem stable.
SOLVER_UTILS_EXPORT NekDouble GetTimeStep(
int ExpOrder,
NekDouble CFL,
NekDouble TimeStability);
const Array<OneD, const Array<OneD, NekDouble> > &inarray);
/// CFL number
NekDouble m_cfl;
// Mapping of the real convective field on the standard element.
// This function gives back the convective filed in the standard
// element to calculate the stability region of the problem in a
// unique way.
SOLVER_UTILS_EXPORT Array<OneD,NekDouble> GetStdVelocity(
const Array<OneD, Array<OneD,NekDouble> > inarray);
/// Function to calculate the stability limit for DG/CG
SOLVER_UTILS_EXPORT NekDouble GetStabilityLimit(int n);
/// Function to calculate the stability limit for DG/CG (a vector
/// of them)
SOLVER_UTILS_EXPORT Array<OneD,NekDouble>
GetStabilityLimitVector(
const Array<OneD,int> &ExpOrder);
/// CFL safety factor (comprise between 0 to 1).
NekDouble m_cflSafetyFactor;
protected:
/// Number of time steps between outputting status information.
int m_infosteps;
......@@ -104,9 +79,14 @@ namespace Nektar
std::vector<FilterSharedPtr> m_filters;
/// Initialises UnsteadySystem class members.
SOLVER_UTILS_EXPORT UnsteadySystem(const LibUtilities::SessionReaderSharedPtr& pSession);
SOLVER_UTILS_EXPORT UnsteadySystem(
const LibUtilities::SessionReaderSharedPtr& pSession);
/// Init object for UnsteadySystem class.
SOLVER_UTILS_EXPORT virtual void v_InitObject();
/// Get the maximum timestep estimator for cfl control.
SOLVER_UTILS_EXPORT NekDouble MaxTimeStepEstimator();
/// Solves an unsteady problem.
SOLVER_UTILS_EXPORT virtual void v_DoSolve();
......@@ -118,7 +98,8 @@ namespace Nektar
SOLVER_UTILS_EXPORT virtual void v_PrintSummary(std::ostream &out);
/// Print the solution at each solution point in a txt file
SOLVER_UTILS_EXPORT virtual void v_AppendOutput1D(Array<OneD, Array<OneD, NekDouble> > &solution1D);
SOLVER_UTILS_EXPORT virtual void v_AppendOutput1D(
Array<OneD, Array<OneD, NekDouble> > &solution1D);
///
SOLVER_UTILS_EXPORT virtual void v_NumericalFlux(
......@@ -149,17 +130,9 @@ namespace Nektar
const int j,
Array<OneD, Array<OneD, NekDouble> > &physfield,
Array<OneD, Array<OneD, NekDouble> > &flux);
/// Virtual function to get the time step
SOLVER_UTILS_EXPORT virtual NekDouble v_GetTimeStep(
const Array<OneD,int> ExpOrder,
const Array<OneD,NekDouble> CFL, NekDouble timeCFL);
/// Virtual function to get the time step
SOLVER_UTILS_EXPORT virtual NekDouble v_GetTimeStep(
int ExpOrder,
NekDouble CFL,
NekDouble TimeStability);
const Array<OneD, const Array<OneD, NekDouble> > &inarray);
private:
///
......
library/SpatialDomains/Conditions.h
View file @ 9556cb2a
......@@ -78,6 +78,8 @@ namespace Nektar
eRterminal,
eCRterminal,
eRCRterminal,
eInflowCFE,
eOutflowCFE,
eNoUserDefined
};
......@@ -104,6 +106,8 @@ namespace Nektar
known_type["Symmetry"] = eSymmetry;
known_type["TimeDependent"] = eTimeDependent;
known_type["IsentropicVortex"] = eIsentropicVortex;
known_type["InflowCFE"] = eInflowCFE;
known_type["OutflowCFE"] = eOutflowCFE;
known_type["NoUserDefined"] = eNoUserDefined;
std::map<const std::string, BndUserDefinedType>::const_iterator it = known_type.find(userDefined);
......
solvers/ADRSolver/EquationSystems/CFLtester.cpp
View file @ 9556cb2a
......@@ -65,7 +65,7 @@ namespace Nektar
if (m_explicitAdvection)
{
m_ode.DefineOdeRhs (&CFLtester::DoOdeRhs, this);
m_ode.DefineOdeRhs (&CFLtester::DoOdeRhs, this);
m_ode.DefineProjection (&CFLtester::DoOdeProjection, this);
}
else
......@@ -78,9 +78,10 @@ namespace Nektar
{
}
void CFLtester::DoOdeRhs(const Array<OneD, const Array<OneD, NekDouble> >&inarray,
Array<OneD, Array<OneD, NekDouble> >&outarray,
const NekDouble time)
void CFLtester::DoOdeRhs(
const Array<OneD, const Array<OneD, NekDouble> >&inarray,
Array<OneD, Array<OneD, NekDouble> >&outarray,
const NekDouble time)
{
int j;
int nvariables = inarray.num_elements();
......@@ -88,7 +89,7 @@ namespace Nektar
switch (m_projectionType)
{
case MultiRegions::eDiscontinuous:
case MultiRegions::eDiscontinuous:
{
int ncoeffs = inarray[0].num_elements();
Array<OneD, Array<OneD, NekDouble> > WeakAdv(nvariables);
......@@ -99,38 +100,37 @@ namespace Nektar
WeakAdv[j] = WeakAdv[j-1] + ncoeffs;
}
WeakDGAdvection(inarray, WeakAdv,true,true);
WeakDGAdvection(inarray, WeakAdv, true, true);
for(j = 0; j < nvariables; ++j)
{
m_fields[j]->MultiplyByElmtInvMass(WeakAdv[j],
WeakAdv[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:
case MultiRegions::eMixed_CG_Discontinuous:
{
// Calculate -V\cdot Grad(u);
for(j = 0; j < nvariables; ++j)
{
AdvectionNonConservativeForm(m_velocity,
AdvectionNonConservativeForm(m_velocity,
inarray[j],
outarray[j]);
Vmath::Neg(npoints,outarray[j],1);
Vmath::Neg(npoints, outarray[j], 1);
}
break;
}
}
}
void CFLtester::DoOdeProjection(const Array<OneD,
const Array<OneD, NekDouble> >&inarray,
Array<OneD, Array<OneD, NekDouble> >&outarray,
const NekDouble time)
void CFLtester::DoOdeProjection(
const Array<OneD, const Array<OneD, NekDouble> >&inarray,
Array<OneD, Array<OneD, NekDouble> >&outarray,
const NekDouble time)
{
int j;
int nvariables = inarray.num_elements();
......@@ -138,7 +138,7 @@ namespace Nektar
switch(m_projectionType)
{
case MultiRegions::eDiscontinuous:
case MultiRegions::eDiscontinuous:
{
// Just copy over array
int npoints = GetNpoints();
......@@ -149,8 +149,8 @@ namespace Nektar
}
}
break;
case MultiRegions::eGalerkin:
case MultiRegions::eMixed_CG_Discontinuous:
case MultiRegions::eGalerkin:
case MultiRegions::eMixed_CG_Discontinuous:
{
Array<OneD, NekDouble> coeffs(m_fields[0]->GetNcoeffs());
......@@ -161,29 +161,35 @@ namespace Nektar
}
break;
}
default:
ASSERTL0(false,"Unknown projection scheme");
break;
default:
ASSERTL0(false,"Unknown projection scheme");
break;
}
}
void CFLtester::v_GetFluxVector(const int i, Array<OneD, Array<OneD, NekDouble> > &physfield,
Array<OneD, Array<OneD, NekDouble> > &flux)
void CFLtester::v_GetFluxVector(
const int i,
Array<OneD, Array<OneD, NekDouble> > &physfield,
Array<OneD, Array<OneD, NekDouble> > &flux)
{
ASSERTL1(flux.num_elements() == m_velocity.num_elements(),"Dimension of flux array and velocity array do not match");
ASSERTL1(flux.num_elements() == m_velocity.num_elements(),
"Dimension of flux array and velocity array do not match");
for(int j = 0; j < flux.num_elements(); ++j)
for (int j = 0; j < flux.num_elements(); ++j)
{
Vmath::Vmul(GetNpoints(),physfield[i],1,
m_velocity[j],1,flux[j],1);
Vmath::Vmul(GetNpoints(),
physfield[i], 1,
m_velocity[j], 1,
flux[j], 1);
}
}
void CFLtester::v_NumericalFlux(Array<OneD, Array<OneD, NekDouble> > &physfield, Array<OneD, Array<OneD, NekDouble> > &numflux)
void CFLtester::v_NumericalFlux(
Array<OneD, Array<OneD, NekDouble> > &physfield,
Array<OneD, Array<OneD, NekDouble> > &numflux)
{
int i;
int nTraceNumPoints = GetTraceNpoints();
int nvel = m_spacedim;
......@@ -192,41 +198,52 @@ namespace Nektar
Array<OneD, NekDouble > Vn (nTraceNumPoints,0.0);
// Get Edge Velocity - Could be stored if time independent
for(i = 0; i < nvel; ++i)
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);
Vmath::Vvtvp(nTraceNumPoints,
m_traceNormals[i], 1,
Fwd, 1, Vn, 1, Vn, 1);
}
for(i = 0; i < numflux.num_elements(); ++i)
for (i = 0; i < numflux.num_elements(); ++i)
{
m_fields[i]->GetFwdBwdTracePhys(physfield[i],Fwd,Bwd);
m_fields[i]->GetFwdBwdTracePhys(physfield[i], Fwd, Bwd);
//evaulate upwinded m_fields[i]
m_fields[i]->GetTrace()->Upwind(Vn,Fwd,Bwd,numflux[i]);
m_fields[i]->GetTrace()->Upwind(Vn,