Commit 19664e9f authored by Gianmarco Mengaldo's avatar Gianmarco Mengaldo
Browse files

Minor modifications to meet the coding-standard requirements.


git-svn-id: https://gforge.sci.utah.edu/svn/nektar/trunk@4000 305cdda6-5ce1-45b3-a98d-dfc68c8b3305
parent 7063c50e
......@@ -50,18 +50,18 @@ namespace Nektar
LibUtilities::SessionReaderSharedPtr pSession,
Array<OneD, MultiRegions::ExpListSharedPtr> pFields)
{
/// Definition of the FR scheme recovered
// Definition of the FR scheme recovered
if(pSession->DefinesSolverInfo("FRSchemeRecovered"))
{
/// Check the FR scheme to be used
// Check the FR scheme to be used
m_FRSchemeRecovered = pSession->GetSolverInfo("FRSchemeRecovered");
/// Check if the FR scheme recovered is set up properly
// Check if the FR scheme recovered is set up properly
if((m_FRSchemeRecovered!="DG") && (m_FRSchemeRecovered!="SD")
&& (m_FRSchemeRecovered!="HU"))
{
fprintf(stderr,"\n ERROR: You need to specify the FRSchemeRecovered in SOLVERINFO\n");
fprintf(stderr," 3 valid choices: 'DG', 'SD', 'HU'. \n");
cerr << "\n ERROR: You must specify FRSchemeRecovered in\n";
cerr << "SOLVERINFO. 3 valid choices: 'DG', 'SD', 'HU'.\n";
exit(1);
}
}
......@@ -70,37 +70,42 @@ namespace Nektar
m_FRSchemeRecovered = "DG";
}
/// Computation of the derivatives of the correction functions in case of FR
// Computation of the derivatives of the correction functions (DG)
if(m_FRSchemeRecovered == "DG")
{
/// Bases initialisation
// Bases initialisation
LibUtilities::BasisSharedPtr Basis;
LibUtilities::BasisSharedPtr BasisFR_Left;
LibUtilities::BasisSharedPtr BasisFR_Right;
Basis = pFields[0]->GetExp(0)->GetBasis(0);
/// Number of modes
// Number of modes
int nModes = Basis->GetNumModes();
/// Total number of quadrature points
// Total number of quadrature points
int nSolutionPts = Basis->GetNumPoints();
/// Type of points
// Type of points
const LibUtilities::PointsKey FRpoints = Basis->GetPointsKey();
/// Construction of the derivatives
const LibUtilities::BasisKey FRBase_Left (LibUtilities::eDG_DG_Left, nSolutionPts, FRpoints);
const LibUtilities::BasisKey FRBase_Right(LibUtilities::eDG_DG_Right, nSolutionPts, FRpoints);
// Construction of the derivatives
const LibUtilities::BasisKey FRBase_Left (LibUtilities::eDG_DG_Left,
nSolutionPts,
FRpoints);
const LibUtilities::BasisKey FRBase_Right(LibUtilities::eDG_DG_Right,
nSolutionPts,
FRpoints);
BasisFR_Left = LibUtilities::BasisManager()[FRBase_Left];
BasisFR_Right = LibUtilities::BasisManager()[FRBase_Right];
/// Storing the derivatives into two global variables
// Storing the derivatives into two global variables
m_dGL = BasisFR_Left ->GetBdata();
m_dGR = BasisFR_Right->GetBdata();
}
/// Computation of the derivatives of the correction functions in case of FR
// Computation of the derivatives of the correction functions (SD)
else if(m_FRSchemeRecovered == "SD")
{
/// Bases initialisation
......@@ -119,8 +124,13 @@ namespace Nektar
const LibUtilities::PointsKey FRpoints = Basis->GetPointsKey();
/// Construction of the derivatives
const LibUtilities::BasisKey FRBase_Left (LibUtilities::eDG_SD_Left, nSolutionPts, FRpoints);
const LibUtilities::BasisKey FRBase_Right(LibUtilities::eDG_SD_Right, nSolutionPts, FRpoints);
const LibUtilities::BasisKey FRBase_Left (LibUtilities::eDG_SD_Left,
nSolutionPts,
FRpoints);
const LibUtilities::BasisKey FRBase_Right(LibUtilities::eDG_SD_Right,
nSolutionPts,
FRpoints);
BasisFR_Left = LibUtilities::BasisManager()[FRBase_Left];
BasisFR_Right = LibUtilities::BasisManager()[FRBase_Right];
......@@ -130,32 +140,37 @@ namespace Nektar
m_dGR = BasisFR_Right->GetBdata();
}
/// Computation of the derivatives of the correction functions in case of FR
// Computation of the derivatives of the correction functions (HU)
else if(m_FRSchemeRecovered == "HU")
{
/// Bases initialisation
// Bases initialisation
LibUtilities::BasisSharedPtr Basis;
LibUtilities::BasisSharedPtr BasisFR_Left;
LibUtilities::BasisSharedPtr BasisFR_Right;
Basis = pFields[0]->GetExp(0)->GetBasis(0);
/// Number of modes
// Number of modes
int nModes = Basis->GetNumModes();
/// Total number of quadrature points
// Total number of quadrature points
int nSolutionPts = Basis->GetNumPoints();
/// Type of points
// Type of points
const LibUtilities::PointsKey FRpoints = Basis->GetPointsKey();
/// Construction of the derivatives
const LibUtilities::BasisKey FRBase_Left (LibUtilities::eDG_HU_Left, nSolutionPts, FRpoints);
const LibUtilities::BasisKey FRBase_Right(LibUtilities::eDG_HU_Right, nSolutionPts, FRpoints);
// Construction of the derivatives
const LibUtilities::BasisKey FRBase_Left (LibUtilities::eDG_HU_Left,
nSolutionPts,
FRpoints);
const LibUtilities::BasisKey FRBase_Right(LibUtilities::eDG_HU_Right,
nSolutionPts,
FRpoints);
BasisFR_Left = LibUtilities::BasisManager()[FRBase_Left];
BasisFR_Right = LibUtilities::BasisManager()[FRBase_Right];
/// Storing the derivatives into two global variables
// Storing the derivatives into two global variables
m_dGL = BasisFR_Left ->GetBdata();
m_dGR = BasisFR_Right->GetBdata();
}
......@@ -168,34 +183,35 @@ namespace Nektar
const Array<OneD, Array<OneD, NekDouble> > &inarray,
Array<OneD, Array<OneD, NekDouble> > &outarray)
{
/// Counter variables
// Counter variables
int i, j;
/// Number of elements
// Number of elements
int nElements = fields[0]->GetExpSize();
/// Number of spatial dimensions
// Number of spatial dimensions
int nDimensions = fields[0]->GetCoordim(0);
/// Number of solution points
// Number of solution points
int nSolutionPts = fields[0]->GetTotPoints();
/// Number of coefficients
// Number of coefficients
int nCoeffs = fields[0]->GetNcoeffs();
/// Number of trace points
// Number of trace points
int nTracePts = fields[0]->GetTrace()->GetTotPoints();
/// Vector to store the discontinuos flux
// Vector to store the discontinuos flux
Array<OneD, Array<OneD, NekDouble> > fluxvector(nDimensions);
/// Vector to store the derivative of the discontinuous flux
// Vector to store the derivative of the discontinuous flux
Array<OneD, Array<OneD, NekDouble> > divfluxvector(nDimensions);
/// Vector to store the solution in physical space
// Vector to store the solution in physical space
Array<OneD, Array<OneD, NekDouble> > physfield (nConvectiveFields);
/// Resize each column of the flux vector to the number of quadrature points
// Resize each column of the flux vector to the number of
// solution points
for(i = 0; i < nDimensions; ++i)
{
fluxvector[i] = Array<OneD, NekDouble>(nSolutionPts);
......@@ -207,22 +223,17 @@ namespace Nektar
physfield[i] = inarray[i];
}
//for(i = 0; i < nConvectiveFields; ++i)
//{
//physfield[i] = Array<OneD, NekDouble>(nSolutionPts);
//fields[i]->BwdTrans(inarray[i], physfield[i]);
//}
/// Get the discontinuous flux FD
// Get the discontinuous flux FD
for(i = 0; i < nConvectiveFields; ++i)
{
/// Get the ith component of the flux vector in physical space
// Get the ith component of the flux vector in physical space
m_fluxVector(i, physfield, fluxvector);
}
Array<OneD,NekDouble> tmpFD, tmpDFD;
/// Computation of the divergence of the discontinuous flux at each quadrature point
// Computation of the divergence of the discontinuous flux at each
// solution point
LibUtilities::BasisSharedPtr Basis;
Basis = fields[0]->GetExp(0)->GetBasis(0);
......@@ -249,34 +260,37 @@ namespace Nektar
fields[i]->GetFwdBwdTracePhys(inarray[i], Fwd[i], Bwd[i]);
}
/// Computing the Riemann flux at each flux (interface) point
// Computing the Riemann flux at each flux (interface) point
m_riemann->Solve(Fwd, Bwd, numflux);
/// Arrays to store the intercell numerical flux jumps
// Arrays to store the intercell numerical flux jumps
Array<OneD, Array<OneD, NekDouble> > numfluxjumpsLeft(nConvectiveFields);
Array<OneD, Array<OneD, NekDouble> > numfluxjumpsRight(nConvectiveFields);
int offsetStart, offsetEnd;
/// Dimension of each column of the jumps array has to be equal to the number of flux points
// Dimension of each column of the jumps array has to be equal to
// the number of flux points
switch(nDimensions)
{
case 1:
{
/// Temporay array for computing the jumps multiply by the derivatives of the correction functions
// Temporay array for computing the jumps multiply by the
// derivatives of the correction functions
Array<OneD,NekDouble> dercorrfluxLeft(nSolutionPts/nElements, 0.0);
Array<OneD,NekDouble> dercorrfluxRight(nSolutionPts/nElements, 0.0);
Array<OneD,NekDouble> tmp, tmparray;
/// The dimension of each column of the jumps arrays is equal to number of trace points minus one
// The dimension of each column of the jumps arrays is equal
// to number of trace points minus one
for(i = 0; i < nConvectiveFields; ++i)
{
numfluxjumpsLeft[i] = Array<OneD, NekDouble>(nTracePts - 1);
numfluxjumpsRight[i] = Array<OneD, NekDouble>(nTracePts - 1);
}
/// Loop to compute the left and the right jump of the flux
// Loop to compute the left and the right jump of the flux
for(i = 0; i < nElements; i++)
{
offsetStart = fields[0]->GetPhys_Offset(i);
......@@ -311,31 +325,31 @@ namespace Nektar
}
case 2:
{
// HOW TO GET THE CORRECT DIMENSION OF THE FLUXJUMPS ARRAY IN 2D
// HOW TO GET CORRECT DIMENSION OF FLUXJUMPS ARRAY IN 2D
ASSERTL0(false,"2D FR case not implemented yet");
break;
}
case 3:
{
// HOW TO GET THE CORRECT DIMENSION OF THE FLUXJUMPS ARRAY IN 3D
// HOW TO GET CORRECT DIMENSION OF FLUXJUMPS ARRAY IN 3D
ASSERTL0(false,"3D FR case not implemented yet");
break;
}
}
/// Array to store the Jacobian and its inverse
// Array to store the Jacobian and its inverse
Array<OneD, const NekDouble>jac(nElements);
Array<OneD, NekDouble> jacobian(nElements);
Array<OneD, NekDouble> tmparray;
/// Evaluation of the jacobian of each element
// Evaluation of the jacobian of each element
for(i = 0; i < nElements; i++)
{
jac = fields[0]->GetExp(i)->GetGeom1D()->GetJac();
jacobian[i] = jac[0];
}
/// Operations to compute the RHS
// Operations to compute the RHS
for(i = 0; i < nConvectiveFields; ++i)
{
for(j = 0; j < nElements; j++)
......
......@@ -58,9 +58,9 @@ namespace Nektar
static std::string type;
std::string m_FRSchemeRecovered; ///< FR scheme recovered (DG, SD, HU)
Array<OneD, NekDouble> m_dGL; ///< Left derivatives of the correction functions
Array<OneD, NekDouble> m_dGR; ///< Right derivatives of the correction functions
std::string m_FRSchemeRecovered;
Array<OneD, NekDouble> m_dGL;
Array<OneD, NekDouble> m_dGR;
protected:
AdvectionFR();
......
......@@ -40,7 +40,8 @@ namespace Nektar
namespace SolverUtils
{
std::string AdvectionNonConservative::type = GetAdvectionFactory().
RegisterCreatorFunction("NonConservative", AdvectionNonConservative::create);
RegisterCreatorFunction("NonConservative",
AdvectionNonConservative::create);
AdvectionNonConservative::AdvectionNonConservative()
{
......@@ -48,7 +49,7 @@ namespace Nektar
}
void AdvectionNonConservative::v_Advect(
const int nConvectiveFields,
const int nConvectiveFields,
const Array<OneD, MultiRegions::ExpListSharedPtr> &fields,
const Array<OneD, Array<OneD, NekDouble> > &advVel,
const Array<OneD, Array<OneD, NekDouble> > &inarray,
......@@ -72,23 +73,45 @@ namespace Nektar
for (int i = 0; i < nConvectiveFields; ++i)
{
// Evaluate V\cdot Grad(u)
// Evaluate V \cdot Grad(u)
switch(nDim)
{
case 1:
fields[0]->PhysDeriv(inarray[i], grad0);
Vmath::Vmul(nPointsTot,grad0,1,advVel[0],1,outarray[i],1);
Vmath::Vmul(nPointsTot,
grad0, 1,
advVel[0], 1,
outarray[i], 1);
break;
case 2:
fields[0]->PhysDeriv(inarray[i],grad0,grad1);
Vmath::Vmul (nPointsTot,grad0,1,advVel[0],1,outarray[i],1);
Vmath::Vvtvp(nPointsTot,grad1,1,advVel[1],1,outarray[i],1,outarray[i],1);
fields[0]->PhysDeriv(inarray[i], grad0, grad1);
Vmath::Vmul (nPointsTot,
grad0, 1,
advVel[0], 1,
outarray[i], 1);
Vmath::Vvtvp(nPointsTot,
grad1, 1,
advVel[1], 1,
outarray[i], 1,
outarray[i], 1);
break;
case 3:
fields[0]->PhysDeriv(inarray[i],grad0,grad1,grad2);
Vmath::Vmul (nPointsTot,grad0,1,advVel[0],1,outarray[i],1);
Vmath::Vvtvvtp(nPointsTot,grad1,1,advVel[1],1,grad2,1,advVel[2],1,
outarray[i],1);
fields[0]->PhysDeriv(inarray[i], grad0, grad1, grad2);
Vmath::Vmul (nPointsTot,
grad0, 1,
advVel[0], 1,
outarray[i], 1);
Vmath::Vvtvvtp(nPointsTot,
grad1, 1,
advVel[1], 1,
grad2, 1,
advVel[2], 1,
outarray[i], 1);
break;
default:
ASSERTL0(false,"dimension unknown");
......
......@@ -48,7 +48,7 @@ namespace Nektar
}
void AdvectionWeakDG::v_Advect(
const int nConvectiveFields,
const int nConvectiveFields,
const Array<OneD, MultiRegions::ExpListSharedPtr> &fields,
const Array<OneD, Array<OneD, NekDouble> > &advVel,
const Array<OneD, Array<OneD, NekDouble> > &inarray,
......@@ -63,7 +63,8 @@ namespace Nektar
Array<OneD, Array<OneD, NekDouble> > fluxvector(nVelDim);
Array<OneD, Array<OneD, NekDouble> > tmp (nConvectiveFields);
ASSERTL1(m_riemann, "Riemann solver must be provided for AdvectionWeakDG.");
ASSERTL1(m_riemann,
"Riemann solver must be provided for AdvectionWeakDG.");
for(i = 0; i < nVelDim; ++i)
{
......@@ -80,7 +81,8 @@ namespace Nektar
for (j = 0; j < nVelDim; ++j)
{
fields[i]->IProductWRTDerivBase(j, fluxvector[j], outarray[i]);
fields[i]->IProductWRTDerivBase(j, fluxvector[j],
outarray[i]);
Vmath::Vadd(nCoeffs, outarray[i], 1, tmp[i], 1, tmp[i], 1);
}
}
......
......@@ -34,7 +34,6 @@
///////////////////////////////////////////////////////////////////////////////
#include <iostream>
#include <ADRSolver/EquationSystems/UnsteadyAdvection.h>
namespace Nektar
......@@ -52,40 +51,41 @@ namespace Nektar
*/
void UnsteadyAdvection::v_InitObject()
{
/// Call to the initialisation object of UnsteadySystem
// Call to the initialisation object of UnsteadySystem
UnsteadySystem::v_InitObject();
/// Check the AdvectionType to be used
// Check the AdvectionType to be used
m_advectionType = m_session->GetSolverInfo("AdvectionType");
/// Check if the AdvectionType is set up properly
if((m_advectionType!="NonConservative")&&(m_advectionType!="WeakDG")&&(m_advectionType!="FR"))
// Check if the AdvectionType is set up properly
if((m_advectionType!="NonConservative") && (m_advectionType!="WeakDG")
&& (m_advectionType!="FR"))
{
fprintf(stderr,"\n ERROR: You need to specify the AdvectionType in SOLVERINFO\n");
fprintf(stderr," Three valid choices: NonConservative', WeakDG, FR. \n");
cerr << "\n ERROR: You must specify AdvectionType in SOLVERINFO\n";
cerr << " Three valid choices: NonConservative', WeakDG, FR.\n";
exit(1);
}
/// Define the normal velocity fields
// Define the normal velocity fields
if (m_fields[0]->GetTrace())
{
m_traceVn = Array<OneD, NekDouble>(GetTraceNpoints());
}
/// Read the advection velocities from session file
// Read the advection velocities from session file
std::vector<std::string> vel;
vel.push_back("Vx");
vel.push_back("Vy");
vel.push_back("Vz");
/// Resize the advection velocitites vector to the dimension of the problem
// Resize the advection velocities vector to dimension of the problem
vel.resize(m_spacedim);
/// Store in the global variable m_velocity the advection velocities
// Store in the global variable m_velocity the advection velocities
m_velocity = Array<OneD, Array<OneD, NekDouble> >(m_spacedim);
EvaluateFunction(vel, m_velocity, "AdvectionVelocity");
/// Create Riemann solver instance, depending on the UPWINDTYPE specified in the session file
// Create Riemann solver instance, depending on UPWINDTYPE tag
if((m_advectionType=="WeakDG")||(m_advectionType=="FR"))
{
m_riemannType = m_session->GetSolverInfo("UpwindType");
......@@ -93,7 +93,7 @@ namespace Nektar
m_riemannSolver->AddScalar("Vn", &UnsteadyAdvection::GetNormalVelocity, this);
}
/// Create an advection object
// Create an advection object
m_advection = SolverUtils::GetAdvectionFactory().CreateInstance(m_advectionType);
m_advection->SetFluxVector (&UnsteadyAdvection::GetFluxVector, this);
if((m_advectionType=="WeakDG")||(m_advectionType=="FR"))
......@@ -102,13 +102,13 @@ namespace Nektar
m_advection->InitObject (m_session, m_fields);
}
/// If explicit it computes the RHS and the PROJECTION for the time integration
// If explicit it computes RHS and PROJECTION for the time integration
if (m_explicitAdvection)
{
m_ode.DefineOdeRhs (&UnsteadyAdvection::DoOdeRhs, this);
m_ode.DefineProjection (&UnsteadyAdvection::DoOdeProjection, this);
}
/// Otherwise it gives an error because there is no implicit integration at the moment
// Otherwise it gives an error (no implicit integration)
else
{
ASSERTL0(false, "Implicit unsteady Advection not set up.");
......@@ -127,20 +127,25 @@ namespace Nektar
*/
Array<OneD, NekDouble> &UnsteadyAdvection::GetNormalVelocity()
{
/// Number of trace (interface) points
// Number of trace (interface) points
int nTracePts = GetTraceNpoints();
/// Auxiliary variable to compute the normal velocity
// Auxiliary variable to compute the normal velocity
Array<OneD, NekDouble> tmp(nTracePts);
/// Reset the normal velocity
// Reset the normal velocity
Vmath::Zero(nTracePts, m_traceVn, 1);
/// Compute the normal velocity
// Compute the normal velocity
for (int i = 0; i < m_velocity.num_elements(); ++i)
{
m_fields[0]->ExtractTracePhys(m_velocity[i], tmp);
Vmath::Vvtvp(nTracePts, m_traceNormals[i], 1, tmp, 1, m_traceVn, 1, m_traceVn, 1);
Vmath::Vvtvp(nTracePts,
m_traceNormals[i], 1,
tmp, 1,
m_traceVn, 1,
m_traceVn, 1);
}
return m_traceVn;
......@@ -153,23 +158,28 @@ namespace Nektar
* @param outarray Calculated solution.
* @param time Time.
*/
void UnsteadyAdvection::DoOdeRhs(const Array<OneD, const Array<OneD, NekDouble> >&inarray,
Array<OneD, Array<OneD, NekDouble> >&outarray,
const NekDouble time)
void UnsteadyAdvection::DoOdeRhs(
const Array<OneD, const Array<OneD, NekDouble> >&inarray,
Array<OneD, Array<OneD, NekDouble> >&outarray,
const NekDouble time)
{
/// Counter variable
// Counter variable
int i;
/// Number of fields (variables of the problem)
// Number of fields (variables of the problem)
int nVariables = inarray.num_elements();
/// Number of solution points
// Number of solution points
int nSolutionPts = GetNpoints();
/// RHS computation using the new advection base class
m_advection->Advect(nVariables, m_fields, m_velocity, inarray, outarray);
// RHS computation using the new advection base class
m_advection->Advect(nVariables,
m_fields,
m_velocity,
inarray,
outarray);
/// Negate the RHS
// Negate the RHS
for (i = 0; i < nVariables; ++i)
{
Vmath::Neg(nSolutionPts, outarray[i], 1);
......@@ -183,29 +193,30 @@ namespace Nektar
* @param outarray Calculated solution.
* @param time Time.
*/
void UnsteadyAdvection::DoOdeProjection(const Array<OneD, const Array<OneD, NekDouble> >&inarray,
Array<OneD, Array<OneD, NekDouble> >&outarray,
const NekDouble time)
void UnsteadyAdvection::DoOdeProjection(
const Array<OneD, const Array<OneD, NekDouble> >&inarray,
Array<OneD, Array<OneD, NekDouble> >&outarray,
const NekDouble time)
{
/// Counter variable
// Counter variable
int i;
/// Number of fields (variables of the problem)
// Number of fields (variables of the problem)
int nVariables = inarray.num_elements();
/// Set the boundary conditions
// Set the boundary conditions
SetBoundaryConditions(time);
/// Switch on the projection type (Discontinuous or Continuous)
// Switch on the projection type (Discontinuous or Continuous)
switch(m_projectionType)
{
/// Discontinuous projection
// Discontinuous projection
case MultiRegions::eDiscontinuousGalerkin:
{
/// Number of quadrature points
// Number of quadrature points
int nQuadraturePts = GetNpoints();
/// Just copy over array
// Just copy over array
for(i = 0; i < nVariables; ++i)
{
Vmath::Vcopy(nQuadraturePts,inarray[i],1,outarray[i],1);
......@@ -213,7 +224,7 @@ namespace Nektar