Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
Static Public Member Functions | Public Attributes | Static Public Attributes | Protected Member Functions | Protected Attributes | List of all members
Nektar::SolverUtils::AdvectionFR Class Reference

#include <AdvectionFR.h>

Inheritance diagram for Nektar::SolverUtils::AdvectionFR:
Inheritance graph
[legend]
Collaboration diagram for Nektar::SolverUtils::AdvectionFR:
Collaboration graph
[legend]

Static Public Member Functions

static AdvectionSharedPtr create (std::string advType)
 

Public Attributes

Array< OneD, NekDoublem_jac
 
Array< OneD, Array< OneD,
NekDouble > > 
m_gmat
 
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
 

Static Public Attributes

static std::string type []
 

Protected Member Functions

 AdvectionFR (std::string advType)
 AdvectionFR uses the Flux Reconstruction (FR) approach to compute the advection term. The implementation is only for segments, quadrilaterals and hexahedra at the moment. More...
 
virtual void v_InitObject (LibUtilities::SessionReaderSharedPtr pSession, Array< OneD, MultiRegions::ExpListSharedPtr > pFields)
 Initiliase AdvectionFR objects and store them before starting the time-stepping. More...
 
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)
 Compute the advection term at each time-step using the Flux Reconstruction approach (FR). More...
 
virtual void v_SetupMetrics (LibUtilities::SessionReaderSharedPtr pSession, Array< OneD, MultiRegions::ExpListSharedPtr > pFields)
 Setup the metric terms to compute the contravariant fluxes. (i.e. this special metric terms transform the fluxes at the interfaces of each element from the physical space to the standard space). More...
 
virtual void v_SetupCFunctions (LibUtilities::SessionReaderSharedPtr pSession, Array< OneD, MultiRegions::ExpListSharedPtr > pFields)
 Setup the derivatives of the correction functions. For more details see J Sci Comput (2011) 47: 50–72. More...
 
virtual void v_DivCFlux_1D (const int nConvectiveFields, const Array< OneD, MultiRegions::ExpListSharedPtr > &fields, const Array< OneD, const NekDouble > &fluxX1, const Array< OneD, const NekDouble > &numericalFlux, Array< OneD, NekDouble > &divCFlux)
 Compute the divergence of the corrective flux for 1D problems. More...
 
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)
 Compute the divergence of the corrective flux for 2D problems. More...
 
virtual void v_DivCFlux_2D_Gauss (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)
 Compute the divergence of the corrective flux for 2D problems where POINTSTYPE="GaussGaussLegendre". More...
 
virtual void v_DivCFlux_3D (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 > &fluxX3, const Array< OneD, const NekDouble > &numericalFlux, Array< OneD, NekDouble > &divCFlux)
 Compute the divergence of the corrective flux for 3D problems. More...
 
- Protected Member Functions inherited from Nektar::SolverUtils::Advection
virtual SOLVER_UTILS_EXPORT void v_SetBaseFlow (const Array< OneD, Array< OneD, NekDouble > > &inarray)
 Overrides the base flow used during linearised advection. More...
 

Protected Attributes

Array< OneD, Array< OneD,
NekDouble > > 
m_traceNormals
 
std::string m_advType
 
- Protected Attributes inherited from Nektar::SolverUtils::Advection
AdvectionFluxVecCB m_fluxVector
 Callback function to the flux vector (set when advection is in conservative form). More...
 
RiemannSolverSharedPtr m_riemann
 Riemann solver for DG-type schemes. More...
 
int m_spaceDim
 Storage for space dimension. Used for homogeneous extension. More...
 

Additional Inherited Members

- Public Member Functions inherited from Nektar::SolverUtils::Advection
SOLVER_UTILS_EXPORT void InitObject (LibUtilities::SessionReaderSharedPtr pSession, Array< OneD, MultiRegions::ExpListSharedPtr > pFields)
 Interface function to initialise the advection object. More...
 
SOLVER_UTILS_EXPORT void 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)
 Interface function to advect the vector field. More...
 
template<typename FuncPointerT , typename ObjectPointerT >
void SetFluxVector (FuncPointerT func, ObjectPointerT obj)
 Set the flux vector callback function. More...
 
void SetRiemannSolver (RiemannSolverSharedPtr riemann)
 Set a Riemann solver object for this advection object. More...
 
void SetFluxVector (AdvectionFluxVecCB fluxVector)
 Set the flux vector callback function. More...
 
void SetBaseFlow (const Array< OneD, Array< OneD, NekDouble > > &inarray)
 Set the base flow used for linearised advection objects. More...
 

Detailed Description

Definition at line 47 of file AdvectionFR.h.

Constructor & Destructor Documentation

Nektar::SolverUtils::AdvectionFR::AdvectionFR ( std::string  advType)
protected

AdvectionFR uses the Flux Reconstruction (FR) approach to compute the advection term. The implementation is only for segments, quadrilaterals and hexahedra at the moment.

Todo:
Extension to triangles, tetrahedra and other shapes. (Long term objective)

Definition at line 77 of file AdvectionFR.cpp.

Referenced by create().

77  :m_advType(advType)
78  {
79  }

Member Function Documentation

static AdvectionSharedPtr Nektar::SolverUtils::AdvectionFR::create ( std::string  advType)
inlinestatic

Definition at line 50 of file AdvectionFR.h.

References AdvectionFR().

51  {
52  return AdvectionSharedPtr(new AdvectionFR(advType));
53  }
AdvectionFR(std::string advType)
AdvectionFR uses the Flux Reconstruction (FR) approach to compute the advection term. The implementation is only for segments, quadrilaterals and hexahedra at the moment.
Definition: AdvectionFR.cpp:77
boost::shared_ptr< Advection > AdvectionSharedPtr
A shared pointer to an Advection object.
Definition: Advection.h:158
void Nektar::SolverUtils::AdvectionFR::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 
)
protectedvirtual

Compute the advection term at each time-step using the Flux Reconstruction approach (FR).

Parameters
nConvectiveFieldsNumber of fields.
fieldsPointer to fields.
advVelAdvection velocities.
inarraySolution at the previous time-step.
outarrayAdvection term to be passed at the time integration class.

Implements Nektar::SolverUtils::Advection.

Definition at line 801 of file AdvectionFR.cpp.

References ASSERTL0, Nektar::LibUtilities::eGaussGaussLegendre, Nektar::SolverUtils::Advection::m_fluxVector, m_gmat, m_jac, Nektar::SolverUtils::Advection::m_riemann, Nektar::SolverUtils::Advection::m_spaceDim, v_DivCFlux_1D(), v_DivCFlux_2D(), v_DivCFlux_2D_Gauss(), Vmath::Vadd(), Vmath::Vdiv(), Vmath::Vmul(), and Vmath::Vvtvvtp().

808  {
809  int i, j, n;
810  int phys_offset;
811 
812  Array<OneD, NekDouble> auxArray1, auxArray2;
813 
814  Array<OneD, LibUtilities::BasisSharedPtr> Basis;
815  Basis = fields[0]->GetExp(0)->GetBase();
816 
817  int nElements = fields[0]->GetExpSize();
818  int nDimensions = fields[0]->GetCoordim(0);
819  int nSolutionPts = fields[0]->GetTotPoints();
820  int nTracePts = fields[0]->GetTrace()->GetTotPoints();
821  int nCoeffs = fields[0]->GetNcoeffs();
822 
823  // Storage for flux vector.
824  Array<OneD, Array<OneD, Array<OneD, NekDouble> > > fluxvector
825  (nConvectiveFields);
826  // Outarray for Galerkin projection in case of primitive dealising
827  Array<OneD, Array<OneD, NekDouble> > outarrayCoeff
828  (nConvectiveFields);
829  // Store forwards/backwards space along trace space
830  Array<OneD, Array<OneD, NekDouble> > Fwd (nConvectiveFields);
831  Array<OneD, Array<OneD, NekDouble> > Bwd (nConvectiveFields);
832  Array<OneD, Array<OneD, NekDouble> > numflux(nConvectiveFields);
833 
834  // Set up storage for flux vector.
835  for (i = 0; i < nConvectiveFields; ++i)
836  {
837  fluxvector[i] =
838  Array<OneD, Array<OneD, NekDouble> >(m_spaceDim);
839 
840  for (j = 0; j < m_spaceDim; ++j)
841  {
842  fluxvector[i][j] = Array<OneD, NekDouble>(nSolutionPts);
843  }
844  }
845 
846  for (i = 0; i < nConvectiveFields; ++i)
847  {
848  outarrayCoeff[i] = Array<OneD, NekDouble>(nCoeffs);
849  Fwd[i] = Array<OneD, NekDouble>(nTracePts);
850  Bwd[i] = Array<OneD, NekDouble>(nTracePts);
851  numflux[i] = Array<OneD, NekDouble>(nTracePts);
852  fields[i]->GetFwdBwdTracePhys(inarray[i], Fwd[i], Bwd[i]);
853  }
854 
855  // Computing the interface flux at each trace point
856  m_riemann->Solve(m_spaceDim, Fwd, Bwd, numflux);
857 
858  // Divergence of the flux (computing the RHS)
859  switch(nDimensions)
860  {
861  // 1D problems
862  case 1:
863  {
864  Array<OneD, NekDouble> DfluxvectorX1(nSolutionPts);
865  Array<OneD, NekDouble> divFC (nSolutionPts);
866 
867  // Get the advection flux vector (solver specific)
868  m_fluxVector(inarray, fluxvector);
869 
870  // Get the discontinuous flux divergence
871  for(i = 0; i < nConvectiveFields; ++i)
872  {
873  for (n = 0; n < nElements; n++)
874  {
875  phys_offset = fields[0]->GetPhys_Offset(n);
876 
877  fields[i]->GetExp(n)->PhysDeriv(
878  0, fluxvector[i][0] + phys_offset,
879  auxArray2 = DfluxvectorX1 + phys_offset);
880  }
881 
882  // Get the correction flux divergence
883  v_DivCFlux_1D(nConvectiveFields, fields,
884  fluxvector[i][0], numflux[i], divFC);
885 
886  // Back to the physical space using global operations
887  Vmath::Vdiv(nSolutionPts, &divFC[0], 1, &m_jac[0], 1,
888  &outarray[i][0], 1);
889 
890  // Adding the total divergence to outarray (RHS)
891  Vmath::Vadd(nSolutionPts, &outarray[i][0], 1,
892  &DfluxvectorX1[0], 1, &outarray[i][0], 1);
893 
894  // Primitive Dealiasing 1D
895  if (!(Basis[0]->Collocation()))
896  {
897  fields[i]->FwdTrans(outarray[i], outarrayCoeff[i]);
898  fields[i]->BwdTrans(outarrayCoeff[i], outarray[i]);
899  }
900  }
901  break;
902  }
903  // 2D problems
904  case 2:
905  {
906  Array<OneD, NekDouble> DfluxvectorX1(nSolutionPts);
907  Array<OneD, NekDouble> DfluxvectorX2(nSolutionPts);
908  Array<OneD, NekDouble> divFD(nSolutionPts);
909  Array<OneD, NekDouble> divFC(nSolutionPts);
910 
911  // Get the advection flux vector (solver specific)
912  m_fluxVector(inarray, fluxvector);
913 
914  for (i = 0; i < nConvectiveFields; ++i)
915  {
916  // Temporary vectors
917  Array<OneD, NekDouble> f_hat(nSolutionPts);
918  Array<OneD, NekDouble> g_hat(nSolutionPts);
919 
920  Vmath::Vvtvvtp(nSolutionPts,
921  &m_gmat[0][0], 1,
922  &fluxvector[i][0][0], 1,
923  &m_gmat[2][0], 1,
924  &fluxvector[i][1][0], 1,
925  &f_hat[0], 1);
926 
927  Vmath::Vmul(nSolutionPts, &m_jac[0], 1, &f_hat[0], 1,
928  &f_hat[0], 1);
929 
930  Vmath::Vvtvvtp(nSolutionPts,
931  &m_gmat[1][0], 1,
932  &fluxvector[i][0][0], 1,
933  &m_gmat[3][0], 1,
934  &fluxvector[i][1][0], 1,
935  &g_hat[0], 1);
936 
937  Vmath::Vmul(nSolutionPts, &m_jac[0], 1, &g_hat[0], 1,
938  &g_hat[0], 1);
939 
940  // Get the discontinuous flux derivatives
941  for (n = 0; n < nElements; n++)
942  {
943  phys_offset = fields[0]->GetPhys_Offset(n);
944  fields[0]->GetExp(n)->StdPhysDeriv(0,
945  auxArray1 = f_hat + phys_offset,
946  auxArray2 = DfluxvectorX1 + phys_offset);
947  fields[0]->GetExp(n)->StdPhysDeriv(1,
948  auxArray1 = g_hat + phys_offset,
949  auxArray2 = DfluxvectorX2 + phys_offset);
950  }
951 
952  // Divergence of the discontinuous flux
953  Vmath::Vadd(nSolutionPts, DfluxvectorX1, 1,
954  DfluxvectorX2, 1, divFD, 1);
955 
956  // Divergence of the correction flux
957  if (Basis[0]->GetPointsType() ==
959  Basis[1]->GetPointsType() ==
961  {
963  nConvectiveFields, fields, f_hat, g_hat,
964  numflux[i], divFC);
965  }
966  else
967  {
969  nConvectiveFields, fields,
970  fluxvector[i][0], fluxvector[i][1],
971  numflux[i], divFC);
972 
973  }
974 
975  // Divergence of the final flux
976  Vmath::Vadd(nSolutionPts, divFD, 1, divFC, 1,
977  outarray[i], 1);
978 
979  // Back to the physical space using a global operation
980  Vmath::Vdiv(nSolutionPts, &outarray[i][0], 1,
981  &m_jac[0], 1, &outarray[i][0], 1);
982 
983  // Primitive Dealiasing 2D
984  if (!(Basis[0]->Collocation()))
985  {
986  fields[i]->FwdTrans(outarray[i], outarrayCoeff[i]);
987  fields[i]->BwdTrans(outarrayCoeff[i], outarray[i]);
988  }
989  } // close nConvectiveFields loop
990  break;
991  }
992  // 3D problems
993  case 3:
994  {
995  ASSERTL0(false,"3D FRDG case not implemented yet");
996  break;
997  }
998  }
999  }
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:161
Array< OneD, Array< OneD, NekDouble > > m_gmat
Definition: AdvectionFR.h:58
RiemannSolverSharedPtr m_riemann
Riemann solver for DG-type schemes.
Definition: Advection.h:134
void Vdiv(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Multiply vector z = x/y.
Definition: Vmath.cpp:227
1D Gauss-Gauss-Legendre quadrature points
Definition: PointsType.h:47
virtual void v_DivCFlux_2D_Gauss(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)
Compute the divergence of the corrective flux for 2D problems where POINTSTYPE="GaussGaussLegendre".
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)
Compute the divergence of the corrective flux for 2D problems.
void Vvtvvtp(int n, const T *v, int incv, const T *w, int incw, const T *x, int incx, const T *y, int incy, T *z, int incz)
vvtvvtp (vector times vector plus vector times vector):
Definition: Vmath.cpp:523
virtual void v_DivCFlux_1D(const int nConvectiveFields, const Array< OneD, MultiRegions::ExpListSharedPtr > &fields, const Array< OneD, const NekDouble > &fluxX1, const Array< OneD, const NekDouble > &numericalFlux, Array< OneD, NekDouble > &divCFlux)
Compute the divergence of the corrective flux for 1D problems.
AdvectionFluxVecCB m_fluxVector
Callback function to the flux vector (set when advection is in conservative form).
Definition: Advection.h:132
int m_spaceDim
Storage for space dimension. Used for homogeneous extension.
Definition: Advection.h:136
Array< OneD, NekDouble > m_jac
Definition: AdvectionFR.h:57
void Vadd(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Add vector z = x+y.
Definition: Vmath.cpp:285
void Vmul(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Multiply vector z = x*y.
Definition: Vmath.cpp:169
void Nektar::SolverUtils::AdvectionFR::v_DivCFlux_1D ( const int  nConvectiveFields,
const Array< OneD, MultiRegions::ExpListSharedPtr > &  fields,
const Array< OneD, const NekDouble > &  fluxX1,
const Array< OneD, const NekDouble > &  numericalFlux,
Array< OneD, NekDouble > &  divCFlux 
)
protectedvirtual

Compute the divergence of the corrective flux for 1D problems.

Parameters
nConvectiveFieldsNumber of fields.
fieldsPointer to fields.
fluxX1X1-volumetric flux in physical space.
numericalFluxInterface flux in physical space.
divCFluxDivergence of the corrective flux.

Definition at line 1011 of file AdvectionFR.cpp.

References m_dGL_xi1, m_dGR_xi1, Vmath::Smul(), Vmath::Vadd(), and Vmath::Vcopy().

Referenced by v_Advect().

1017  {
1018  int n;
1019  int nLocalSolutionPts, phys_offset, t_offset;
1020 
1022  Basis = fields[0]->GetExp(0)->GetBasis(0);
1023 
1024  int nElements = fields[0]->GetExpSize();
1025  int nSolutionPts = fields[0]->GetTotPoints();
1026 
1027 
1028  vector<bool> negatedFluxNormal = (boost::static_pointer_cast<MultiRegions::DisContField1D>(fields[0]))->GetNegatedFluxNormal();
1029 
1030  // Arrays to store the derivatives of the correction flux
1031  Array<OneD, NekDouble> DCL(nSolutionPts/nElements, 0.0);
1032  Array<OneD, NekDouble> DCR(nSolutionPts/nElements, 0.0);
1033 
1034  // Arrays to store the intercell numerical flux jumps
1035  Array<OneD, NekDouble> JumpL(nElements);
1036  Array<OneD, NekDouble> JumpR(nElements);
1037 
1038  Array<OneD, Array<OneD, LocalRegions::ExpansionSharedPtr> >
1039  &elmtToTrace = fields[0]->GetTraceMap()->GetElmtToTrace();
1040 
1041  for (n = 0; n < nElements; ++n)
1042  {
1043  nLocalSolutionPts = fields[0]->GetExp(n)->GetTotPoints();
1044  phys_offset = fields[0]->GetPhys_Offset(n);
1045 
1046  Array<OneD, NekDouble> tmparrayX1(nLocalSolutionPts, 0.0);
1047  NekDouble tmpFluxVertex = 0;
1048  Vmath::Vcopy(nLocalSolutionPts,
1049  &fluxX1[phys_offset], 1,
1050  &tmparrayX1[0], 1);
1051 
1052  fields[0]->GetExp(n)->GetVertexPhysVals(0, tmparrayX1,
1053  tmpFluxVertex);
1054 
1055  t_offset = fields[0]->GetTrace()
1056  ->GetPhys_Offset(elmtToTrace[n][0]->GetElmtId());
1057 
1058  if(negatedFluxNormal[2*n])
1059  {
1060  JumpL[n] = numericalFlux[t_offset] - tmpFluxVertex;
1061  }
1062  else
1063  {
1064  JumpL[n] = -numericalFlux[t_offset] - tmpFluxVertex;
1065  }
1066 
1067  fields[0]->GetExp(n)->GetVertexPhysVals(1, tmparrayX1,
1068  tmpFluxVertex);
1069 
1070  t_offset = fields[0]->GetTrace()
1071  ->GetPhys_Offset(elmtToTrace[n][1]->GetElmtId());
1072 
1073  if(negatedFluxNormal[2*n+1])
1074  {
1075  JumpR[n] = -numericalFlux[t_offset] - tmpFluxVertex;
1076  }
1077  else
1078  {
1079  JumpR[n] = numericalFlux[t_offset] - tmpFluxVertex;
1080  }
1081  }
1082 
1083  for (n = 0; n < nElements; ++n)
1084  {
1085  nLocalSolutionPts = fields[0]->GetExp(n)->GetTotPoints();
1086  phys_offset = fields[0]->GetPhys_Offset(n);
1087 
1088  // Left jump multiplied by left derivative of C function
1089  Vmath::Smul(nLocalSolutionPts, JumpL[n], &m_dGL_xi1[n][0], 1,
1090  &DCL[0], 1);
1091 
1092  // Right jump multiplied by right derivative of C function
1093  Vmath::Smul(nLocalSolutionPts, JumpR[n], &m_dGR_xi1[n][0], 1,
1094  &DCR[0], 1);
1095 
1096  // Assembling divergence of the correction flux
1097  Vmath::Vadd(nLocalSolutionPts, &DCL[0], 1, &DCR[0], 1,
1098  &divCFlux[phys_offset], 1);
1099  }
1100  }
Array< OneD, Array< OneD, NekDouble > > m_dGL_xi1
Definition: AdvectionFR.h:65
Array< OneD, Array< OneD, NekDouble > > m_dGR_xi1
Definition: AdvectionFR.h:66
void Smul(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Scalar multiply y = alpha*y.
Definition: Vmath.cpp:199
double NekDouble
boost::shared_ptr< Basis > BasisSharedPtr
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.cpp:1047
void Vadd(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Add vector z = x+y.
Definition: Vmath.cpp:285
void Nektar::SolverUtils::AdvectionFR::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 
)
protectedvirtual

Compute the divergence of the corrective flux for 2D problems.

Parameters
nConvectiveFieldsNumber of fields.
fieldsPointer to fields.
fluxX1X1-volumetric flux in physical space.
fluxX2X2-volumetric flux in physical space.
numericalFluxInterface flux in physical space.
divCFluxDivergence of the corrective flux.
Todo:
: Switch on shapes eventually here.

Definition at line 1114 of file AdvectionFR.cpp.

References ASSERTL0, Nektar::StdRegions::eBackwards, m_dGL_xi1, m_dGL_xi2, m_dGR_xi1, m_dGR_xi2, m_Q2D_e0, m_Q2D_e1, m_Q2D_e2, m_Q2D_e3, m_traceNormals, Vmath::Reverse(), Vmath::Vadd(), Vmath::Vsub(), and Vmath::Vvtvvtp().

Referenced by v_Advect().

1121  {
1122  int n, e, i, j, cnt;
1123 
1124  int nElements = fields[0]->GetExpSize();
1125 
1126  int nLocalSolutionPts, nEdgePts;
1127  int trace_offset, phys_offset;
1128  int nquad0, nquad1;
1129 
1130  Array<OneD, NekDouble> auxArray1;
1131  Array<OneD, LibUtilities::BasisSharedPtr> base;
1132 
1133  Array<OneD, Array<OneD, LocalRegions::ExpansionSharedPtr> >
1134  &elmtToTrace = fields[0]->GetTraceMap()->GetElmtToTrace();
1135 
1136  // Loop on the elements
1137  for(n = 0; n < nElements; ++n)
1138  {
1139  // Offset of the element on the global vector
1140  phys_offset = fields[0]->GetPhys_Offset(n);
1141  nLocalSolutionPts = fields[0]->GetExp(n)->GetTotPoints();
1142 
1143  base = fields[0]->GetExp(n)->GetBase();
1144  nquad0 = base[0]->GetNumPoints();
1145  nquad1 = base[1]->GetNumPoints();
1146 
1147  Array<OneD, NekDouble> divCFluxE0(nLocalSolutionPts, 0.0);
1148  Array<OneD, NekDouble> divCFluxE1(nLocalSolutionPts, 0.0);
1149  Array<OneD, NekDouble> divCFluxE2(nLocalSolutionPts, 0.0);
1150  Array<OneD, NekDouble> divCFluxE3(nLocalSolutionPts, 0.0);
1151 
1152  // Loop on the edges
1153  for(e = 0; e < fields[0]->GetExp(n)->GetNedges(); ++e)
1154  {
1155  // Number of edge points of edge e
1156  nEdgePts = fields[0]->GetExp(n)->GetEdgeNumPoints(e);
1157 
1158  Array<OneD, NekDouble> tmparrayX1(nEdgePts, 0.0);
1159  Array<OneD, NekDouble> tmparrayX2(nEdgePts, 0.0);
1160  Array<OneD, NekDouble> fluxN (nEdgePts, 0.0);
1161  Array<OneD, NekDouble> fluxT (nEdgePts, 0.0);
1162  Array<OneD, NekDouble> fluxJumps (nEdgePts, 0.0);
1163 
1164  // Offset of the trace space correspondent to edge e
1165  trace_offset = fields[0]->GetTrace()->GetPhys_Offset(
1166  elmtToTrace[n][e]->GetElmtId());
1167 
1168  // Get the normals of edge e
1169  const Array<OneD, const Array<OneD, NekDouble> > &normals =
1170  fields[0]->GetExp(n)->GetEdgeNormal(e);
1171 
1172  // Extract the edge values of flux-x on edge e and order
1173  // them accordingly to the order of the trace space
1174  fields[0]->GetExp(n)->GetEdgePhysVals(
1175  e, elmtToTrace[n][e],
1176  fluxX1 + phys_offset,
1177  auxArray1 = tmparrayX1);
1178 
1179  // Extract the edge values of flux-y on edge e and order
1180  // them accordingly to the order of the trace space
1181  fields[0]->GetExp(n)->GetEdgePhysVals(
1182  e, elmtToTrace[n][e],
1183  fluxX2 + phys_offset,
1184  auxArray1 = tmparrayX2);
1185 
1186  // Multiply the edge components of the flux by the normal
1187  Vmath::Vvtvvtp(nEdgePts, &tmparrayX1[0], 1,
1188  &m_traceNormals[0][trace_offset], 1,
1189  &tmparrayX2[0], 1,
1190  &m_traceNormals[1][trace_offset], 1,
1191  &fluxN[0], 1);
1192 
1193  // Subtract to the Riemann flux the discontinuous flux
1194  Vmath::Vsub(nEdgePts, &numericalFlux[trace_offset], 1,
1195  &fluxN[0], 1, &fluxJumps[0], 1);
1196 
1197  // Check the ordering of the jump vectors
1198  if (fields[0]->GetExp(n)->GetEorient(e) ==
1200  {
1201  Vmath::Reverse(nEdgePts, &fluxJumps[0], 1,
1202  &fluxJumps[0], 1);
1203  }
1204 
1205  NekDouble fac = fields[0]->GetExp(n)->EdgeNormalNegated(e) ?
1206  -1.0 : 1.0;
1207 
1208  for (i = 0; i < nEdgePts; ++i)
1209  {
1210  if (m_traceNormals[0][trace_offset+i] != fac*normals[0][i]
1211  || m_traceNormals[1][trace_offset+i] != fac*normals[1][i])
1212  {
1213  fluxJumps[i] = -fluxJumps[i];
1214  }
1215  }
1216 
1217  // Multiply jumps by derivatives of the correction functions
1218  switch (e)
1219  {
1220  case 0:
1221  for (i = 0; i < nquad0; ++i)
1222  {
1223  // Multiply fluxJumps by Q factors
1224  fluxJumps[i] = -(m_Q2D_e0[n][i]) * fluxJumps[i];
1225 
1226  for (j = 0; j < nquad1; ++j)
1227  {
1228  cnt = i + j*nquad0;
1229  divCFluxE0[cnt] = fluxJumps[i] * m_dGL_xi2[n][j];
1230  }
1231  }
1232  break;
1233  case 1:
1234  for (i = 0; i < nquad1; ++i)
1235  {
1236  // Multiply fluxJumps by Q factors
1237  fluxJumps[i] = (m_Q2D_e1[n][i]) * fluxJumps[i];
1238 
1239  for (j = 0; j < nquad0; ++j)
1240  {
1241  cnt = (nquad0)*i + j;
1242  divCFluxE1[cnt] = fluxJumps[i] * m_dGR_xi1[n][j];
1243  }
1244  }
1245  break;
1246  case 2:
1247  for (i = 0; i < nquad0; ++i)
1248  {
1249  // Multiply fluxJumps by Q factors
1250  fluxJumps[i] = (m_Q2D_e2[n][i]) * fluxJumps[i];
1251 
1252  for (j = 0; j < nquad1; ++j)
1253  {
1254  cnt = (nquad0 - 1) + j*nquad0 - i;
1255  divCFluxE2[cnt] = fluxJumps[i] * m_dGR_xi2[n][j];
1256  }
1257  }
1258  break;
1259  case 3:
1260  for (i = 0; i < nquad1; ++i)
1261  {
1262  // Multiply fluxJumps by Q factors
1263  fluxJumps[i] = -(m_Q2D_e3[n][i]) * fluxJumps[i];
1264  for (j = 0; j < nquad0; ++j)
1265  {
1266  cnt = (nquad0*nquad1 - nquad0) + j - i*nquad0;
1267  divCFluxE3[cnt] = fluxJumps[i] * m_dGL_xi1[n][j];
1268  }
1269  }
1270  break;
1271  default:
1272  ASSERTL0(false,"edge value (< 3) is out of range");
1273  break;
1274  }
1275  }
1276 
1277  // Sum each edge contribution
1278  Vmath::Vadd(nLocalSolutionPts, &divCFluxE0[0], 1,
1279  &divCFluxE1[0], 1, &divCFlux[phys_offset], 1);
1280 
1281  Vmath::Vadd(nLocalSolutionPts, &divCFlux[phys_offset], 1,
1282  &divCFluxE2[0], 1, &divCFlux[phys_offset], 1);
1283 
1284  Vmath::Vadd(nLocalSolutionPts, &divCFlux[phys_offset], 1,
1285  &divCFluxE3[0], 1, &divCFlux[phys_offset], 1);
1286  }
1287  }
Array< OneD, Array< OneD, NekDouble > > m_Q2D_e1
Definition: AdvectionFR.h:61
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:161
Array< OneD, Array< OneD, NekDouble > > m_dGL_xi1
Definition: AdvectionFR.h:65
Array< OneD, Array< OneD, NekDouble > > m_dGR_xi1
Definition: AdvectionFR.h:66
Array< OneD, Array< OneD, NekDouble > > m_Q2D_e0
Definition: AdvectionFR.h:60
Array< OneD, Array< OneD, NekDouble > > m_dGR_xi2
Definition: AdvectionFR.h:68
Array< OneD, Array< OneD, NekDouble > > m_Q2D_e3
Definition: AdvectionFR.h:63
Array< OneD, Array< OneD, NekDouble > > m_Q2D_e2
Definition: AdvectionFR.h:62
void Reverse(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.cpp:1071
double NekDouble
Array< OneD, Array< OneD, NekDouble > > m_dGL_xi2
Definition: AdvectionFR.h:67
void Vsub(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Subtract vector z = x-y.
Definition: Vmath.cpp:329
void Vvtvvtp(int n, const T *v, int incv, const T *w, int incw, const T *x, int incx, const T *y, int incy, T *z, int incz)
vvtvvtp (vector times vector plus vector times vector):
Definition: Vmath.cpp:523
Array< OneD, Array< OneD, NekDouble > > m_traceNormals
Definition: AdvectionFR.h:77
void Vadd(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Add vector z = x+y.
Definition: Vmath.cpp:285
void Nektar::SolverUtils::AdvectionFR::v_DivCFlux_2D_Gauss ( 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 
)
protectedvirtual

Compute the divergence of the corrective flux for 2D problems where POINTSTYPE="GaussGaussLegendre".

Parameters
nConvectiveFieldsNumber of fields.
fieldsPointer to fields.
fluxX1X1-volumetric flux in physical space.
fluxX2X2-volumetric flux in physical space.
numericalFluxInterface flux in physical space.
divCFluxDivergence of the corrective flux.
Todo:
: Switch on shapes eventually here.

Definition at line 1303 of file AdvectionFR.cpp.

References ASSERTL0, Nektar::StdRegions::eBackwards, m_dGL_xi1, m_dGL_xi2, m_dGR_xi1, m_dGR_xi2, m_Q2D_e0, m_Q2D_e1, m_Q2D_e2, m_Q2D_e3, m_traceNormals, Vmath::Neg(), Vmath::Reverse(), Vmath::Vadd(), Vmath::Vcopy(), and Vmath::Vsub().

Referenced by v_Advect().

1310  {
1311  int n, e, i, j, cnt;
1312 
1313  int nElements = fields[0]->GetExpSize();
1314  int nLocalSolutionPts;
1315  int nEdgePts;
1316  int trace_offset;
1317  int phys_offset;
1318  int nquad0;
1319  int nquad1;
1320 
1321  NekDouble fac;
1322 
1323  Array<OneD, NekDouble> auxArray1, auxArray2;
1324  Array<OneD, LibUtilities::BasisSharedPtr> base;
1325 
1326  Array<OneD, Array<OneD, LocalRegions::ExpansionSharedPtr> >
1327  &elmtToTrace = fields[0]->GetTraceMap()->GetElmtToTrace();
1328 
1329  // Loop on the elements
1330  for(n = 0; n < nElements; ++n)
1331  {
1332  // Offset of the element on the global vector
1333  phys_offset = fields[0]->GetPhys_Offset(n);
1334  nLocalSolutionPts = fields[0]->GetExp(n)->GetTotPoints();
1335 
1336  base = fields[0]->GetExp(n)->GetBase();
1337  nquad0 = base[0]->GetNumPoints();
1338  nquad1 = base[1]->GetNumPoints();
1339 
1340  Array<OneD, NekDouble> divCFluxE0(nLocalSolutionPts, 0.0);
1341  Array<OneD, NekDouble> divCFluxE1(nLocalSolutionPts, 0.0);
1342  Array<OneD, NekDouble> divCFluxE2(nLocalSolutionPts, 0.0);
1343  Array<OneD, NekDouble> divCFluxE3(nLocalSolutionPts, 0.0);
1344 
1345  // Loop on the edges
1346  for(e = 0; e < fields[0]->GetExp(n)->GetNedges(); ++e)
1347  {
1348  // Number of edge points of edge e
1349  nEdgePts = fields[0]->GetExp(n)->GetEdgeNumPoints(e);
1350 
1351  Array<OneD, NekDouble> fluxN (nEdgePts, 0.0);
1352  Array<OneD, NekDouble> fluxT (nEdgePts, 0.0);
1353  Array<OneD, NekDouble> fluxN_R (nEdgePts, 0.0);
1354  Array<OneD, NekDouble> fluxN_D (nEdgePts, 0.0);
1355  Array<OneD, NekDouble> fluxJumps (nEdgePts, 0.0);
1356 
1357  // Offset of the trace space correspondent to edge e
1358  trace_offset = fields[0]->GetTrace()->GetPhys_Offset(
1359  elmtToTrace[n][e]->GetElmtId());
1360 
1361  // Get the normals of edge e
1362  const Array<OneD, const Array<OneD, NekDouble> > &normals =
1363  fields[0]->GetExp(n)->GetEdgeNormal(e);
1364 
1365  // Extract the trasformed normal flux at each edge
1366  switch (e)
1367  {
1368  case 0:
1369  // Extract the edge values of transformed flux-y on
1370  // edge e and order them accordingly to the order of
1371  // the trace space
1372  fields[0]->GetExp(n)->GetEdgePhysVals(
1373  e, elmtToTrace[n][e],
1374  fluxX2 + phys_offset,
1375  auxArray1 = fluxN_D);
1376 
1377  Vmath::Neg (nEdgePts, fluxN_D, 1);
1378 
1379  // Extract the physical Rieamann flux at each edge
1380  Vmath::Vcopy(nEdgePts,
1381  &numericalFlux[trace_offset], 1,
1382  &fluxN[0], 1);
1383 
1384  // Check the ordering of vectors
1385  if (fields[0]->GetExp(n)->GetEorient(e) ==
1387  {
1388  Vmath::Reverse(nEdgePts,
1389  auxArray1 = fluxN, 1,
1390  auxArray2 = fluxN, 1);
1391 
1392  Vmath::Reverse(nEdgePts,
1393  auxArray1 = fluxN_D, 1,
1394  auxArray2 = fluxN_D, 1);
1395  }
1396 
1397  // Transform Riemann Fluxes in the standard element
1398  for (i = 0; i < nquad0; ++i)
1399  {
1400  // Multiply Riemann Flux by Q factors
1401  fluxN_R[i] = (m_Q2D_e0[n][i]) * fluxN[i];
1402  }
1403 
1404  fac = fields[0]->GetExp(n)->EdgeNormalNegated(e) ?
1405  -1.0 : 1.0;
1406 
1407  for (i = 0; i < nEdgePts; ++i)
1408  {
1409  if (m_traceNormals[0][trace_offset+i]
1410  != fac*normals[0][i] ||
1411  m_traceNormals[1][trace_offset+i]
1412  != fac*normals[1][i])
1413  {
1414  fluxN_R[i] = -fluxN_R[i];
1415  }
1416  }
1417 
1418  // Subtract to the Riemann flux the discontinuous
1419  // flux
1420  Vmath::Vsub(nEdgePts,
1421  &fluxN_R[0], 1,
1422  &fluxN_D[0], 1, &fluxJumps[0], 1);
1423 
1424  // Multiplicate the flux jumps for the correction
1425  // function
1426  for (i = 0; i < nquad0; ++i)
1427  {
1428  for (j = 0; j < nquad1; ++j)
1429  {
1430  cnt = i + j*nquad0;
1431  divCFluxE0[cnt] =
1432  -fluxJumps[i] * m_dGL_xi2[n][j];
1433  }
1434  }
1435  break;
1436  case 1:
1437  // Extract the edge values of transformed flux-x on
1438  // edge e and order them accordingly to the order of
1439  // the trace space
1440  fields[0]->GetExp(n)->GetEdgePhysVals(
1441  e, elmtToTrace[n][e],
1442  fluxX1 + phys_offset,
1443  auxArray1 = fluxN_D);
1444 
1445  // Extract the physical Rieamann flux at each edge
1446  Vmath::Vcopy(nEdgePts,
1447  &numericalFlux[trace_offset], 1,
1448  &fluxN[0], 1);
1449 
1450  // Check the ordering of vectors
1451  if (fields[0]->GetExp(n)->GetEorient(e) ==
1453  {
1454  Vmath::Reverse(nEdgePts,
1455  auxArray1 = fluxN, 1,
1456  auxArray2 = fluxN, 1);
1457 
1458  Vmath::Reverse(nEdgePts,
1459  auxArray1 = fluxN_D, 1,
1460  auxArray2 = fluxN_D, 1);
1461  }
1462 
1463  // Transform Riemann Fluxes in the standard element
1464  for (i = 0; i < nquad1; ++i)
1465  {
1466  // Multiply Riemann Flux by Q factors
1467  fluxN_R[i] = (m_Q2D_e1[n][i]) * fluxN[i];
1468  }
1469 
1470  fac = fields[0]->GetExp(n)->EdgeNormalNegated(e) ?
1471  -1.0 : 1.0;
1472 
1473  for (i = 0; i < nEdgePts; ++i)
1474  {
1475  if (m_traceNormals[0][trace_offset+i]
1476  != fac*normals[0][i] ||
1477  m_traceNormals[1][trace_offset+i]
1478  != fac*normals[1][i])
1479  {
1480  fluxN_R[i] = -fluxN_R[i];
1481  }
1482  }
1483 
1484  // Subtract to the Riemann flux the discontinuous
1485  // flux
1486  Vmath::Vsub(nEdgePts,
1487  &fluxN_R[0], 1,
1488  &fluxN_D[0], 1, &fluxJumps[0], 1);
1489 
1490  // Multiplicate the flux jumps for the correction
1491  // function
1492  for (i = 0; i < nquad1; ++i)
1493  {
1494  for (j = 0; j < nquad0; ++j)
1495  {
1496  cnt = (nquad0)*i + j;
1497  divCFluxE1[cnt] =
1498  fluxJumps[i] * m_dGR_xi1[n][j];
1499  }
1500  }
1501  break;
1502  case 2:
1503 
1504  // Extract the edge values of transformed flux-y on
1505  // edge e and order them accordingly to the order of
1506  // the trace space
1507 
1508  fields[0]->GetExp(n)->GetEdgePhysVals(
1509  e, elmtToTrace[n][e],
1510  fluxX2 + phys_offset,
1511  auxArray1 = fluxN_D);
1512 
1513  // Extract the physical Rieamann flux at each edge
1514  Vmath::Vcopy(nEdgePts,
1515  &numericalFlux[trace_offset], 1,
1516  &fluxN[0], 1);
1517 
1518  // Check the ordering of vectors
1519  if (fields[0]->GetExp(n)->GetEorient(e) ==
1521  {
1522  Vmath::Reverse(nEdgePts,
1523  auxArray1 = fluxN, 1,
1524  auxArray2 = fluxN, 1);
1525 
1526  Vmath::Reverse(nEdgePts,
1527  auxArray1 = fluxN_D, 1,
1528  auxArray2 = fluxN_D, 1);
1529  }
1530 
1531  // Transform Riemann Fluxes in the standard element
1532  for (i = 0; i < nquad0; ++i)
1533  {
1534  // Multiply Riemann Flux by Q factors
1535  fluxN_R[i] = (m_Q2D_e2[n][i]) * fluxN[i];
1536  }
1537 
1538  fac = fields[0]->GetExp(n)->EdgeNormalNegated(e) ?
1539  -1.0 : 1.0;
1540 
1541  for (i = 0; i < nEdgePts; ++i)
1542  {
1543  if (m_traceNormals[0][trace_offset+i]
1544  != fac*normals[0][i] ||
1545  m_traceNormals[1][trace_offset+i]
1546  != fac*normals[1][i])
1547  {
1548  fluxN_R[i] = -fluxN_R[i];
1549  }
1550  }
1551 
1552  // Subtract to the Riemann flux the discontinuous
1553  // flux
1554 
1555  Vmath::Vsub(nEdgePts,
1556  &fluxN_R[0], 1,
1557  &fluxN_D[0], 1, &fluxJumps[0], 1);
1558 
1559  // Multiplicate the flux jumps for the correction
1560  // function
1561  for (i = 0; i < nquad0; ++i)
1562  {
1563  for (j = 0; j < nquad1; ++j)
1564  {
1565  cnt = (nquad0 - 1) + j*nquad0 - i;
1566  divCFluxE2[cnt] =
1567  fluxJumps[i] * m_dGR_xi2[n][j];
1568  }
1569  }
1570  break;
1571  case 3:
1572  // Extract the edge values of transformed flux-x on
1573  // edge e and order them accordingly to the order of
1574  // the trace space
1575 
1576  fields[0]->GetExp(n)->GetEdgePhysVals(
1577  e, elmtToTrace[n][e],
1578  fluxX1 + phys_offset,
1579  auxArray1 = fluxN_D);
1580  Vmath::Neg (nEdgePts, fluxN_D, 1);
1581 
1582  // Extract the physical Rieamann flux at each edge
1583  Vmath::Vcopy(nEdgePts,
1584  &numericalFlux[trace_offset], 1,
1585  &fluxN[0], 1);
1586 
1587  // Check the ordering of vectors
1588  if (fields[0]->GetExp(n)->GetEorient(e) ==
1590  {
1591  Vmath::Reverse(nEdgePts,
1592  auxArray1 = fluxN, 1,
1593  auxArray2 = fluxN, 1);
1594 
1595  Vmath::Reverse(nEdgePts,
1596  auxArray1 = fluxN_D, 1,
1597  auxArray2 = fluxN_D, 1);
1598  }
1599 
1600  // Transform Riemann Fluxes in the standard element
1601  for (i = 0; i < nquad1; ++i)
1602  {
1603  // Multiply Riemann Flux by Q factors
1604  fluxN_R[i] = (m_Q2D_e3[n][i]) * fluxN[i];
1605  }
1606 
1607  fac = fields[0]->GetExp(n)->EdgeNormalNegated(e) ?
1608  -1.0 : 1.0;
1609 
1610  for (i = 0; i < nEdgePts; ++i)
1611  {
1612  if (m_traceNormals[0][trace_offset+i]
1613  != fac*normals[0][i] ||
1614  m_traceNormals[1][trace_offset+i]
1615  != fac*normals[1][i])
1616  {
1617  fluxN_R[i] = -fluxN_R[i];
1618  }
1619  }
1620 
1621  // Subtract to the Riemann flux the discontinuous
1622  // flux
1623 
1624  Vmath::Vsub(nEdgePts,
1625  &fluxN_R[0], 1,
1626  &fluxN_D[0], 1, &fluxJumps[0], 1);
1627 
1628  // Multiplicate the flux jumps for the correction
1629  // function
1630  for (i = 0; i < nquad1; ++i)
1631  {
1632  for (j = 0; j < nquad0; ++j)
1633  {
1634  cnt = (nquad0*nquad1 - nquad0) + j
1635  - i*nquad0;
1636  divCFluxE3[cnt] =
1637  -fluxJumps[i] * m_dGL_xi1[n][j];
1638  }
1639  }
1640  break;
1641  default:
1642  ASSERTL0(false,"edge value (< 3) is out of range");
1643  break;
1644  }
1645  }
1646 
1647 
1648  // Sum each edge contribution
1649  Vmath::Vadd(nLocalSolutionPts, &divCFluxE0[0], 1,
1650  &divCFluxE1[0], 1, &divCFlux[phys_offset], 1);
1651 
1652  Vmath::Vadd(nLocalSolutionPts, &divCFlux[phys_offset], 1,
1653  &divCFluxE2[0], 1, &divCFlux[phys_offset], 1);
1654 
1655  Vmath::Vadd(nLocalSolutionPts, &divCFlux[phys_offset], 1,
1656  &divCFluxE3[0], 1, &divCFlux[phys_offset], 1);
1657  }
1658  }
Array< OneD, Array< OneD, NekDouble > > m_Q2D_e1
Definition: AdvectionFR.h:61
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:161
Array< OneD, Array< OneD, NekDouble > > m_dGL_xi1
Definition: AdvectionFR.h:65
Array< OneD, Array< OneD, NekDouble > > m_dGR_xi1
Definition: AdvectionFR.h:66
Array< OneD, Array< OneD, NekDouble > > m_Q2D_e0
Definition: AdvectionFR.h:60
Array< OneD, Array< OneD, NekDouble > > m_dGR_xi2
Definition: AdvectionFR.h:68
Array< OneD, Array< OneD, NekDouble > > m_Q2D_e3
Definition: AdvectionFR.h:63
Array< OneD, Array< OneD, NekDouble > > m_Q2D_e2
Definition: AdvectionFR.h:62
void Reverse(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.cpp:1071
void Neg(int n, T *x, const int incx)
Negate x = -x.
Definition: Vmath.cpp:382
double NekDouble
Array< OneD, Array< OneD, NekDouble > > m_dGL_xi2
Definition: AdvectionFR.h:67
void Vsub(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Subtract vector z = x-y.
Definition: Vmath.cpp:329
Array< OneD, Array< OneD, NekDouble > > m_traceNormals
Definition: AdvectionFR.h:77
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.cpp:1047
void Vadd(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Add vector z = x+y.
Definition: Vmath.cpp:285
void Nektar::SolverUtils::AdvectionFR::v_DivCFlux_3D ( 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 > &  fluxX3,
const Array< OneD, const NekDouble > &  numericalFlux,
Array< OneD, NekDouble > &  divCFlux 
)
protectedvirtual

Compute the divergence of the corrective flux for 3D problems.

Parameters
nConvectiveFieldsNumber of fields.
fieldsPointer to fields.
fluxX1X1-volumetric flux in physical space.
fluxX2X2-volumetric flux in physical space.
fluxX3X3-volumetric flux in physical space.
numericalFluxInterface flux in physical space.
divCFluxDivergence of the corrective flux.
Todo:
: To be implemented. Switch on shapes eventually here.

Definition at line 1674 of file AdvectionFR.cpp.

1682  {
1683 
1684  }
void Nektar::SolverUtils::AdvectionFR::v_InitObject ( LibUtilities::SessionReaderSharedPtr  pSession,
Array< OneD, MultiRegions::ExpListSharedPtr pFields 
)
protectedvirtual

Initiliase AdvectionFR objects and store them before starting the time-stepping.

This routine calls the virtual functions v_SetupMetrics and v_SetupCFunctions to initialise the objects needed by AdvectionFR.

Parameters
pSessionPointer to session reader.
pFieldsPointer to fields.

Reimplemented from Nektar::SolverUtils::Advection.

Definition at line 91 of file AdvectionFR.cpp.

References Nektar::SolverUtils::Advection::v_InitObject(), v_SetupCFunctions(), and v_SetupMetrics().

94  {
95  Advection::v_InitObject(pSession, pFields);
96  v_SetupMetrics (pSession, pFields);
97  v_SetupCFunctions (pSession, pFields);
98  }
virtual void v_SetupMetrics(LibUtilities::SessionReaderSharedPtr pSession, Array< OneD, MultiRegions::ExpListSharedPtr > pFields)
Setup the metric terms to compute the contravariant fluxes. (i.e. this special metric terms transform...
virtual SOLVER_UTILS_EXPORT void v_InitObject(LibUtilities::SessionReaderSharedPtr pSession, Array< OneD, MultiRegions::ExpListSharedPtr > pFields)
Initialises the advection object.
Definition: Advection.cpp:97
virtual void v_SetupCFunctions(LibUtilities::SessionReaderSharedPtr pSession, Array< OneD, MultiRegions::ExpListSharedPtr > pFields)
Setup the derivatives of the correction functions. For more details see J Sci Comput (2011) 47: 50–7...
void Nektar::SolverUtils::AdvectionFR::v_SetupCFunctions ( LibUtilities::SessionReaderSharedPtr  pSession,
Array< OneD, MultiRegions::ExpListSharedPtr pFields 
)
protectedvirtual

Setup the derivatives of the correction functions. For more details see J Sci Comput (2011) 47: 50–72.

This routine calls 3 different bases: #eDG_DG_Left - #eDG_DG_Left which recovers a nodal DG scheme, #eDG_SD_Left - #eDG_SD_Left which recovers the SD scheme, #eDG_HU_Left - #eDG_HU_Left which recovers the Huynh scheme. The values of the derivatives of the correction function are then stored into global variables and reused into the virtual functions v_DivCFlux_1D, v_DivCFlux_2D, v_DivCFlux_3D to compute the the divergence of the correction flux for 1D, 2D or 3D problems.

Parameters
pSessionPointer to session reader.
pFieldsPointer to fields.

Definition at line 270 of file AdvectionFR.cpp.

References ASSERTL0, Polylib::jacobd(), m_advType, m_dGL_xi1, m_dGL_xi2, m_dGL_xi3, m_dGR_xi1, m_dGR_xi2, and m_dGR_xi3.

Referenced by v_InitObject().

273  {
274  int i, n;
275  NekDouble c0 = 0.0;
276  NekDouble c1 = 0.0;
277  NekDouble c2 = 0.0;
278  int nquad0, nquad1, nquad2;
279  int nmodes0, nmodes1, nmodes2;
280  Array<OneD, LibUtilities::BasisSharedPtr> base;
281 
282  int nElements = pFields[0]->GetExpSize();
283  int nDimensions = pFields[0]->GetCoordim(0);
284 
285  switch (nDimensions)
286  {
287  case 1:
288  {
289  m_dGL_xi1 = Array<OneD, Array<OneD, NekDouble> >(nElements);
290  m_dGR_xi1 = Array<OneD, Array<OneD, NekDouble> >(nElements);
291 
292  for (n = 0; n < nElements; ++n)
293  {
294  base = pFields[0]->GetExp(n)->GetBase();
295  nquad0 = base[0]->GetNumPoints();
296  nmodes0 = base[0]->GetNumModes();
297  Array<OneD, const NekDouble> z0;
298  Array<OneD, const NekDouble> w0;
299 
300  base[0]->GetZW(z0, w0);
301 
302  m_dGL_xi1[n] = Array<OneD, NekDouble>(nquad0);
303  m_dGR_xi1[n] = Array<OneD, NekDouble>(nquad0);
304 
305  // Auxiliary vectors to build up the auxiliary
306  // derivatives of the Legendre polynomials
307  Array<OneD,NekDouble> dLp0 (nquad0, 0.0);
308  Array<OneD,NekDouble> dLpp0(nquad0, 0.0);
309  Array<OneD,NekDouble> dLpm0(nquad0, 0.0);
310 
311  // Degree of the correction functions
312  int p0 = nmodes0 - 1;
313 
314  // Function sign to compute left correction function
315  NekDouble sign0 = pow(-1.0, p0);
316 
317  // Factorial factor to build the scheme
318  NekDouble ap0 = boost::math::tgamma(2 * p0 + 1)
319  / (pow(2.0, p0)
320  * boost::math::tgamma(p0 + 1)
321  * boost::math::tgamma(p0 + 1));
322 
323  // Scalar parameter which recovers the FR schemes
324  if (m_advType == "FRDG")
325  {
326  c0 = 0.0;
327  }
328  else if (m_advType == "FRSD")
329  {
330  c0 = 2.0 * p0 / ((2.0 * p0 + 1.0) * (p0 + 1.0)
331  * (ap0 * boost::math::tgamma(p0 + 1))
332  * (ap0 * boost::math::tgamma(p0 + 1)));
333  }
334  else if (m_advType == "FRHU")
335  {
336  c0 = 2.0 * (p0 + 1.0) / ((2.0 * p0 + 1.0) * p0
337  * (ap0 * boost::math::tgamma(p0 + 1))
338  * (ap0 * boost::math::tgamma(p0 + 1)));
339  }
340  else if (m_advType == "FRcmin")
341  {
342  c0 = -2.0 / ((2.0 * p0 + 1.0)
343  * (ap0 * boost::math::tgamma(p0 + 1))
344  * (ap0 * boost::math::tgamma(p0 + 1)));
345  }
346  else if (m_advType == "FRcinf")
347  {
348  c0 = 10000000000000000.0;
349  }
350 
351  NekDouble etap0 = 0.5 * c0 * (2.0 * p0 + 1.0)
352  * (ap0 * boost::math::tgamma(p0 + 1))
353  * (ap0 * boost::math::tgamma(p0 + 1));
354 
355  NekDouble overeta0 = 1.0 / (1.0 + etap0);
356 
357  // Derivative of the Legendre polynomials
358  // dLp = derivative of the Legendre polynomial of order p
359  // dLpp = derivative of the Legendre polynomial of order p+1
360  // dLpm = derivative of the Legendre polynomial of order p-1
361  Polylib::jacobd(nquad0, z0.data(), &(dLp0[0]), p0, 0.0, 0.0);
362  Polylib::jacobd(nquad0, z0.data(), &(dLpp0[0]), p0+1, 0.0, 0.0);
363  Polylib::jacobd(nquad0, z0.data(), &(dLpm0[0]), p0-1, 0.0, 0.0);
364 
365  // Building the DG_c_Left
366  for(i = 0; i < nquad0; ++i)
367  {
368  m_dGL_xi1[n][i] = etap0 * dLpm0[i];
369  m_dGL_xi1[n][i] += dLpp0[i];
370  m_dGL_xi1[n][i] *= overeta0;
371  m_dGL_xi1[n][i] = dLp0[i] - m_dGL_xi1[n][i];
372  m_dGL_xi1[n][i] = 0.5 * sign0 * m_dGL_xi1[n][i];
373  }
374 
375  // Building the DG_c_Right
376  for(i = 0; i < nquad0; ++i)
377  {
378  m_dGR_xi1[n][i] = etap0 * dLpm0[i];
379  m_dGR_xi1[n][i] += dLpp0[i];
380  m_dGR_xi1[n][i] *= overeta0;
381  m_dGR_xi1[n][i] += dLp0[i];
382  m_dGR_xi1[n][i] = 0.5 * m_dGR_xi1[n][i];
383  }
384  }
385  break;
386  }
387  case 2:
388  {
389  LibUtilities::BasisSharedPtr BasisFR_Left0;
390  LibUtilities::BasisSharedPtr BasisFR_Right0;
391  LibUtilities::BasisSharedPtr BasisFR_Left1;
392  LibUtilities::BasisSharedPtr BasisFR_Right1;
393 
394  m_dGL_xi1 = Array<OneD, Array<OneD, NekDouble> >(nElements);
395  m_dGR_xi1 = Array<OneD, Array<OneD, NekDouble> >(nElements);
396  m_dGL_xi2 = Array<OneD, Array<OneD, NekDouble> >(nElements);
397  m_dGR_xi2 = Array<OneD, Array<OneD, NekDouble> >(nElements);
398 
399  for (n = 0; n < nElements; ++n)
400  {
401  base = pFields[0]->GetExp(n)->GetBase();
402  nquad0 = base[0]->GetNumPoints();
403  nquad1 = base[1]->GetNumPoints();
404  nmodes0 = base[0]->GetNumModes();
405  nmodes1 = base[1]->GetNumModes();
406 
407  Array<OneD, const NekDouble> z0;
408  Array<OneD, const NekDouble> w0;
409  Array<OneD, const NekDouble> z1;
410  Array<OneD, const NekDouble> w1;
411 
412  base[0]->GetZW(z0, w0);
413  base[1]->GetZW(z1, w1);
414 
415  m_dGL_xi1[n] = Array<OneD, NekDouble>(nquad0);
416  m_dGR_xi1[n] = Array<OneD, NekDouble>(nquad0);
417  m_dGL_xi2[n] = Array<OneD, NekDouble>(nquad1);
418  m_dGR_xi2[n] = Array<OneD, NekDouble>(nquad1);
419 
420  // Auxiliary vectors to build up the auxiliary
421  // derivatives of the Legendre polynomials
422  Array<OneD,NekDouble> dLp0 (nquad0, 0.0);
423  Array<OneD,NekDouble> dLpp0 (nquad0, 0.0);
424  Array<OneD,NekDouble> dLpm0 (nquad0, 0.0);
425  Array<OneD,NekDouble> dLp1 (nquad1, 0.0);
426  Array<OneD,NekDouble> dLpp1 (nquad1, 0.0);
427  Array<OneD,NekDouble> dLpm1 (nquad1, 0.0);
428 
429  // Degree of the correction functions
430  int p0 = nmodes0 - 1;
431  int p1 = nmodes1 - 1;
432 
433  // Function sign to compute left correction function
434  NekDouble sign0 = pow(-1.0, p0);
435  NekDouble sign1 = pow(-1.0, p1);
436 
437  // Factorial factor to build the scheme
438  NekDouble ap0 = boost::math::tgamma(2 * p0 + 1)
439  / (pow(2.0, p0)
440  * boost::math::tgamma(p0 + 1)
441  * boost::math::tgamma(p0 + 1));
442 
443  NekDouble ap1 = boost::math::tgamma(2 * p1 + 1)
444  / (pow(2.0, p1)
445  * boost::math::tgamma(p1 + 1)
446  * boost::math::tgamma(p1 + 1));
447 
448  // Scalar parameter which recovers the FR schemes
449  if (m_advType == "FRDG")
450  {
451  c0 = 0.0;
452  c1 = 0.0;
453  }
454  else if (m_advType == "FRSD")
455  {
456  c0 = 2.0 * p0 / ((2.0 * p0 + 1.0) * (p0 + 1.0)
457  * (ap0 * boost::math::tgamma(p0 + 1))
458  * (ap0 * boost::math::tgamma(p0 + 1)));
459 
460  c1 = 2.0 * p1 / ((2.0 * p1 + 1.0) * (p1 + 1.0)
461  * (ap1 * boost::math::tgamma(p1 + 1))
462  * (ap1 * boost::math::tgamma(p1 + 1)));
463  }
464  else if (m_advType == "FRHU")
465  {
466  c0 = 2.0 * (p0 + 1.0) / ((2.0 * p0 + 1.0) * p0
467  * (ap0 * boost::math::tgamma(p0 + 1))
468  * (ap0 * boost::math::tgamma(p0 + 1)));
469 
470  c1 = 2.0 * (p1 + 1.0) / ((2.0 * p1 + 1.0) * p1
471  * (ap1 * boost::math::tgamma(p1 + 1))
472  * (ap1 * boost::math::tgamma(p1 + 1)));
473  }
474  else if (m_advType == "FRcmin")
475  {
476  c0 = -2.0 / ((2.0 * p0 + 1.0)
477  * (ap0 * boost::math::tgamma(p0 + 1))
478  * (ap0 * boost::math::tgamma(p0 + 1)));
479 
480  c1 = -2.0 / ((2.0 * p1 + 1.0)
481  * (ap1 * boost::math::tgamma(p1 + 1))
482  * (ap1 * boost::math::tgamma(p1 + 1)));
483  }
484  else if (m_advType == "FRcinf")
485  {
486  c0 = 10000000000000000.0;
487  c1 = 10000000000000000.0;
488  }
489 
490  NekDouble etap0 = 0.5 * c0 * (2.0 * p0 + 1.0)
491  * (ap0 * boost::math::tgamma(p0 + 1))
492  * (ap0 * boost::math::tgamma(p0 + 1));
493 
494  NekDouble etap1 = 0.5 * c1 * (2.0 * p1 + 1.0)
495  * (ap1 * boost::math::tgamma(p1 + 1))
496  * (ap1 * boost::math::tgamma(p1 + 1));
497 
498  NekDouble overeta0 = 1.0 / (1.0 + etap0);
499  NekDouble overeta1 = 1.0 / (1.0 + etap1);
500 
501  // Derivative of the Legendre polynomials
502  // dLp = derivative of the Legendre polynomial of order p
503  // dLpp = derivative of the Legendre polynomial of order p+1
504  // dLpm = derivative of the Legendre polynomial of order p-1
505  Polylib::jacobd(nquad0, z0.data(), &(dLp0[0]), p0, 0.0, 0.0);
506  Polylib::jacobd(nquad0, z0.data(), &(dLpp0[0]), p0+1, 0.0, 0.0);
507  Polylib::jacobd(nquad0, z0.data(), &(dLpm0[0]), p0-1, 0.0, 0.0);
508 
509  Polylib::jacobd(nquad1, z1.data(), &(dLp1[0]), p1, 0.0, 0.0);
510  Polylib::jacobd(nquad1, z1.data(), &(dLpp1[0]), p1+1, 0.0, 0.0);
511  Polylib::jacobd(nquad1, z1.data(), &(dLpm1[0]), p1-1, 0.0, 0.0);
512 
513  // Building the DG_c_Left
514  for(i = 0; i < nquad0; ++i)
515  {
516  m_dGL_xi1[n][i] = etap0 * dLpm0[i];
517  m_dGL_xi1[n][i] += dLpp0[i];
518  m_dGL_xi1[n][i] *= overeta0;
519  m_dGL_xi1[n][i] = dLp0[i] - m_dGL_xi1[n][i];
520  m_dGL_xi1[n][i] = 0.5 * sign0 * m_dGL_xi1[n][i];
521  }
522 
523  // Building the DG_c_Left
524  for(i = 0; i < nquad1; ++i)
525  {
526  m_dGL_xi2[n][i] = etap1 * dLpm1[i];
527  m_dGL_xi2[n][i] += dLpp1[i];
528  m_dGL_xi2[n][i] *= overeta1;
529  m_dGL_xi2[n][i] = dLp1[i] - m_dGL_xi2[n][i];
530  m_dGL_xi2[n][i] = 0.5 * sign1 * m_dGL_xi2[n][i];
531  }
532 
533  // Building the DG_c_Right
534  for(i = 0; i < nquad0; ++i)
535  {
536  m_dGR_xi1[n][i] = etap0 * dLpm0[i];
537  m_dGR_xi1[n][i] += dLpp0[i];
538  m_dGR_xi1[n][i] *= overeta0;
539  m_dGR_xi1[n][i] += dLp0[i];
540  m_dGR_xi1[n][i] = 0.5 * m_dGR_xi1[n][i];
541  }
542 
543  // Building the DG_c_Right
544  for(i = 0; i < nquad1; ++i)
545  {
546  m_dGR_xi2[n][i] = etap1 * dLpm1[i];
547  m_dGR_xi2[n][i] += dLpp1[i];
548  m_dGR_xi2[n][i] *= overeta1;
549  m_dGR_xi2[n][i] += dLp1[i];
550  m_dGR_xi2[n][i] = 0.5 * m_dGR_xi2[n][i];
551  }
552  }
553  break;
554  }
555  case 3:
556  {
557  m_dGL_xi1 = Array<OneD, Array<OneD, NekDouble> >(nElements);
558  m_dGR_xi1 = Array<OneD, Array<OneD, NekDouble> >(nElements);
559  m_dGL_xi2 = Array<OneD, Array<OneD, NekDouble> >(nElements);
560  m_dGR_xi2 = Array<OneD, Array<OneD, NekDouble> >(nElements);
561  m_dGL_xi3 = Array<OneD, Array<OneD, NekDouble> >(nElements);
562  m_dGR_xi3 = Array<OneD, Array<OneD, NekDouble> >(nElements);
563 
564  for (n = 0; n < nElements; ++n)
565  {
566  base = pFields[0]->GetExp(n)->GetBase();
567  nquad0 = base[0]->GetNumPoints();
568  nquad1 = base[1]->GetNumPoints();
569  nquad2 = base[2]->GetNumPoints();
570  nmodes0 = base[0]->GetNumModes();
571  nmodes1 = base[1]->GetNumModes();
572  nmodes2 = base[2]->GetNumModes();
573 
574  Array<OneD, const NekDouble> z0;
575  Array<OneD, const NekDouble> w0;
576  Array<OneD, const NekDouble> z1;
577  Array<OneD, const NekDouble> w1;
578  Array<OneD, const NekDouble> z2;
579  Array<OneD, const NekDouble> w2;
580 
581  base[0]->GetZW(z0, w0);
582  base[1]->GetZW(z1, w1);
583  base[1]->GetZW(z2, w2);
584 
585  m_dGL_xi1[n] = Array<OneD, NekDouble>(nquad0);
586  m_dGR_xi1[n] = Array<OneD, NekDouble>(nquad0);
587  m_dGL_xi2[n] = Array<OneD, NekDouble>(nquad1);
588  m_dGR_xi2[n] = Array<OneD, NekDouble>(nquad1);
589  m_dGL_xi3[n] = Array<OneD, NekDouble>(nquad2);
590  m_dGR_xi3[n] = Array<OneD, NekDouble>(nquad2);
591 
592  // Auxiliary vectors to build up the auxiliary
593  // derivatives of the Legendre polynomials
594  Array<OneD,NekDouble> dLp0 (nquad0, 0.0);
595  Array<OneD,NekDouble> dLpp0 (nquad0, 0.0);
596  Array<OneD,NekDouble> dLpm0 (nquad0, 0.0);
597  Array<OneD,NekDouble> dLp1 (nquad1, 0.0);
598  Array<OneD,NekDouble> dLpp1 (nquad1, 0.0);
599  Array<OneD,NekDouble> dLpm1 (nquad1, 0.0);
600  Array<OneD,NekDouble> dLp2 (nquad2, 0.0);
601  Array<OneD,NekDouble> dLpp2 (nquad2, 0.0);
602  Array<OneD,NekDouble> dLpm2 (nquad2, 0.0);
603 
604  // Degree of the correction functions
605  int p0 = nmodes0 - 1;
606  int p1 = nmodes1 - 1;
607  int p2 = nmodes2 - 1;
608 
609  // Function sign to compute left correction function
610  NekDouble sign0 = pow(-1.0, p0);
611  NekDouble sign1 = pow(-1.0, p1);
612 
613  // Factorial factor to build the scheme
614  NekDouble ap0 = boost::math::tgamma(2 * p0 + 1)
615  / (pow(2.0, p0)
616  * boost::math::tgamma(p0 + 1)
617  * boost::math::tgamma(p0 + 1));
618 
619  // Factorial factor to build the scheme
620  NekDouble ap1 = boost::math::tgamma(2 * p1 + 1)
621  / (pow(2.0, p1)
622  * boost::math::tgamma(p1 + 1)
623  * boost::math::tgamma(p1 + 1));
624 
625  // Factorial factor to build the scheme
626  NekDouble ap2 = boost::math::tgamma(2 * p2 + 1)
627  / (pow(2.0, p2)
628  * boost::math::tgamma(p2 + 1)
629  * boost::math::tgamma(p2 + 1));
630 
631  // Scalar parameter which recovers the FR schemes
632  if (m_advType == "FRDG")
633  {
634  c0 = 0.0;
635  c1 = 0.0;
636  c2 = 0.0;
637  }
638  else if (m_advType == "FRSD")
639  {
640  c0 = 2.0 * p0 / ((2.0 * p0 + 1.0) * (p0 + 1.0)
641  * (ap0 * boost::math::tgamma(p0 + 1))
642  * (ap0 * boost::math::tgamma(p0 + 1)));
643 
644  c1 = 2.0 * p1 / ((2.0 * p1 + 1.0) * (p1 + 1.0)
645  * (ap1 * boost::math::tgamma(p1 + 1))
646  * (ap1 * boost::math::tgamma(p1 + 1)));
647 
648  c2 = 2.0 * p2 / ((2.0 * p2 + 1.0) * (p2 + 1.0)
649  * (ap2 * boost::math::tgamma(p2 + 1))
650  * (ap2 * boost::math::tgamma(p2 + 1)));
651  }
652  else if (m_advType == "FRHU")
653  {
654  c0 = 2.0 * (p0 + 1.0) / ((2.0 * p0 + 1.0) * p0
655  * (ap0 * boost::math::tgamma(p0 + 1))
656  * (ap0 * boost::math::tgamma(p0 + 1)));
657 
658  c1 = 2.0 * (p1 + 1.0) / ((2.0 * p1 + 1.0) * p1
659  * (ap1 * boost::math::tgamma(p1 + 1))
660  * (ap1 * boost::math::tgamma(p1 + 1)));
661 
662  c2 = 2.0 * (p2 + 1.0) / ((2.0 * p2 + 1.0) * p2
663  * (ap2 * boost::math::tgamma(p2 + 1))
664  * (ap2 * boost::math::tgamma(p2 + 1)));
665  }
666  else if (m_advType == "FRcmin")
667  {
668  c0 = -2.0 / ((2.0 * p0 + 1.0)
669  * (ap0 * boost::math::tgamma(p0 + 1))
670  * (ap0 * boost::math::tgamma(p0 + 1)));
671 
672  c1 = -2.0 / ((2.0 * p1 + 1.0)
673  * (ap1 * boost::math::tgamma(p1 + 1))
674  * (ap1 * boost::math::tgamma(p1 + 1)));
675 
676  c2 = -2.0 / ((2.0 * p2 + 1.0)
677  * (ap2 * boost::math::tgamma(p2 + 1))
678  * (ap2 * boost::math::tgamma(p2 + 1)));
679  }
680  else if (m_advType == "FRcinf")
681  {
682  c0 = 10000000000000000.0;
683  c1 = 10000000000000000.0;
684  c2 = 10000000000000000.0;
685  }
686 
687  NekDouble etap0 = 0.5 * c0 * (2.0 * p0 + 1.0)
688  * (ap0 * boost::math::tgamma(p0 + 1))
689  * (ap0 * boost::math::tgamma(p0 + 1));
690 
691  NekDouble etap1 = 0.5 * c1 * (2.0 * p1 + 1.0)
692  * (ap1 * boost::math::tgamma(p1 + 1))
693  * (ap1 * boost::math::tgamma(p1 + 1));
694 
695  NekDouble etap2 = 0.5 * c2 * (2.0 * p2 + 1.0)
696  * (ap2 * boost::math::tgamma(p2 + 1))
697  * (ap2 * boost::math::tgamma(p2 + 1));
698 
699  NekDouble overeta0 = 1.0 / (1.0 + etap0);
700  NekDouble overeta1 = 1.0 / (1.0 + etap1);
701  NekDouble overeta2 = 1.0 / (1.0 + etap2);
702 
703  // Derivative of the Legendre polynomials
704  // dLp = derivative of the Legendre polynomial of order p
705  // dLpp = derivative of the Legendre polynomial of order p+1
706  // dLpm = derivative of the Legendre polynomial of order p-1
707  Polylib::jacobd(nquad0, z0.data(), &(dLp0[0]), p0, 0.0, 0.0);
708  Polylib::jacobd(nquad0, z0.data(), &(dLpp0[0]), p0+1, 0.0, 0.0);
709  Polylib::jacobd(nquad0, z0.data(), &(dLpm0[0]), p0-1, 0.0, 0.0);
710 
711  Polylib::jacobd(nquad1, z1.data(), &(dLp1[0]), p1, 0.0, 0.0);
712  Polylib::jacobd(nquad1, z1.data(), &(dLpp1[0]), p1+1, 0.0, 0.0);
713  Polylib::jacobd(nquad1, z1.data(), &(dLpm1[0]), p1-1, 0.0, 0.0);
714 
715  Polylib::jacobd(nquad2, z2.data(), &(dLp2[0]), p2, 0.0, 0.0);
716  Polylib::jacobd(nquad2, z2.data(), &(dLpp2[0]), p2+1, 0.0, 0.0);
717  Polylib::jacobd(nquad2, z2.data(), &(dLpm2[0]), p2-1, 0.0, 0.0);
718 
719  // Building the DG_c_Left
720  for(i = 0; i < nquad0; ++i)
721  {
722  m_dGL_xi1[n][i] = etap0 * dLpm0[i];
723  m_dGL_xi1[n][i] += dLpp0[i];
724  m_dGL_xi1[n][i] *= overeta0;
725  m_dGL_xi1[n][i] = dLp0[i] - m_dGL_xi1[n][i];
726  m_dGL_xi1[n][i] = 0.5 * sign0 * m_dGL_xi1[n][i];
727  }
728 
729  // Building the DG_c_Left
730  for(i = 0; i < nquad1; ++i)
731  {
732  m_dGL_xi2[n][i] = etap1 * dLpm1[i];
733  m_dGL_xi2[n][i] += dLpp1[i];
734  m_dGL_xi2[n][i] *= overeta1;
735  m_dGL_xi2[n][i] = dLp1[i] - m_dGL_xi2[n][i];
736  m_dGL_xi2[n][i] = 0.5 * sign1 * m_dGL_xi2[n][i];
737  }
738 
739  // Building the DG_c_Left
740  for(i = 0; i < nquad2; ++i)
741  {
742  m_dGL_xi3[n][i] = etap2 * dLpm2[i];
743  m_dGL_xi3[n][i] += dLpp2[i];
744  m_dGL_xi3[n][i] *= overeta2;
745  m_dGL_xi3[n][i] = dLp2[i] - m_dGL_xi3[n][i];
746  m_dGL_xi3[n][i] = 0.5 * sign1 * m_dGL_xi3[n][i];
747  }
748 
749  // Building the DG_c_Right
750  for(i = 0; i < nquad0; ++i)
751  {
752  m_dGR_xi1[n][i] = etap0 * dLpm0[i];
753  m_dGR_xi1[n][i] += dLpp0[i];
754  m_dGR_xi1[n][i] *= overeta0;
755  m_dGR_xi1[n][i] += dLp0[i];
756  m_dGR_xi1[n][i] = 0.5 * m_dGR_xi1[n][i];
757  }
758 
759  // Building the DG_c_Right
760  for(i = 0; i < nquad1; ++i)
761  {
762  m_dGR_xi2[n][i] = etap1 * dLpm1[i];
763  m_dGR_xi2[n][i] += dLpp1[i];
764  m_dGR_xi2[n][i] *= overeta1;
765  m_dGR_xi2[n][i] += dLp1[i];
766  m_dGR_xi2[n][i] = 0.5 * m_dGR_xi2[n][i];
767  }
768 
769  // Building the DG_c_Right
770  for(i = 0; i < nquad2; ++i)
771  {
772  m_dGR_xi3[n][i] = etap2 * dLpm2[i];
773  m_dGR_xi3[n][i] += dLpp2[i];
774  m_dGR_xi3[n][i] *= overeta2;
775  m_dGR_xi3[n][i] += dLp2[i];
776  m_dGR_xi3[n][i] = 0.5 * m_dGR_xi3[n][i];
777  }
778  }
779  break;
780  }
781  default:
782  {
783  ASSERTL0(false,"Expansion dimension not recognised");
784  break;
785  }
786  }
787  }
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:161
Array< OneD, Array< OneD, NekDouble > > m_dGL_xi1
Definition: AdvectionFR.h:65
Array< OneD, Array< OneD, NekDouble > > m_dGL_xi3
Definition: AdvectionFR.h:69
Array< OneD, Array< OneD, NekDouble > > m_dGR_xi1
Definition: AdvectionFR.h:66
Array< OneD, Array< OneD, NekDouble > > m_dGR_xi3
Definition: AdvectionFR.h:70
Array< OneD, Array< OneD, NekDouble > > m_dGR_xi2
Definition: AdvectionFR.h:68
void jacobd(const int np, const double *z, double *polyd, const int n, const double alpha, const double beta)
Calculate the derivative of Jacobi polynomials.
Definition: Polylib.cpp:2146
double NekDouble
Array< OneD, Array< OneD, NekDouble > > m_dGL_xi2
Definition: AdvectionFR.h:67
boost::shared_ptr< Basis > BasisSharedPtr
void Nektar::SolverUtils::AdvectionFR::v_SetupMetrics ( LibUtilities::SessionReaderSharedPtr  pSession,
Array< OneD, MultiRegions::ExpListSharedPtr pFields 
)
protectedvirtual

Setup the metric terms to compute the contravariant fluxes. (i.e. this special metric terms transform the fluxes at the interfaces of each element from the physical space to the standard space).

This routine calls the function #GetEdgeQFactors to compute and store the metric factors following an anticlockwise conventions along the edges/faces of each element. Note: for 1D problem the transformation is not needed.

Parameters
pSessionPointer to session reader.
pFieldsPointer to fields.
Todo:
Add the metric terms for 3D Hexahedra.

Definition at line 116 of file AdvectionFR.cpp.

References ASSERTL0, Nektar::SpatialDomains::eDeformed, Nektar::StdRegions::StdExpansion::GetMetricInfo(), m_gmat, m_jac, m_Q2D_e0, m_Q2D_e1, m_Q2D_e2, m_Q2D_e3, and m_traceNormals.

Referenced by v_InitObject().

119  {
120  int i, n;
121  int nquad0, nquad1;
122  int phys_offset;
123  int nLocalSolutionPts;
124  int nElements = pFields[0]->GetExpSize();
125  int nDimensions = pFields[0]->GetCoordim(0);
126  int nSolutionPts = pFields[0]->GetTotPoints();
127  int nTracePts = pFields[0]->GetTrace()->GetTotPoints();
128 
129  Array<TwoD, const NekDouble> gmat;
130  Array<OneD, const NekDouble> jac;
131 
132  m_jac = Array<OneD, NekDouble>(nSolutionPts);
133 
134  Array<OneD, NekDouble> auxArray1;
135  Array<OneD, LibUtilities::BasisSharedPtr> base;
137 
138  switch (nDimensions)
139  {
140  case 1:
141  {
142  for (n = 0; n < nElements; ++n)
143  {
144  ptsKeys = pFields[0]->GetExp(n)->GetPointsKeys();
145  nLocalSolutionPts = pFields[0]->GetExp(n)->GetTotPoints();
146  phys_offset = pFields[0]->GetPhys_Offset(n);
147  jac = pFields[0]->GetExp(n)
148  ->as<LocalRegions::Expansion1D>()->GetGeom1D()
149  ->GetMetricInfo()->GetJac(ptsKeys);
150  for (i = 0; i < nLocalSolutionPts; ++i)
151  {
152  m_jac[i+phys_offset] = jac[0];
153  }
154  }
155  break;
156  }
157  case 2:
158  {
159  m_gmat = Array<OneD, Array<OneD, NekDouble> >(4);
160  m_gmat[0] = Array<OneD, NekDouble>(nSolutionPts);
161  m_gmat[1] = Array<OneD, NekDouble>(nSolutionPts);
162  m_gmat[2] = Array<OneD, NekDouble>(nSolutionPts);
163  m_gmat[3] = Array<OneD, NekDouble>(nSolutionPts);
164 
165  m_Q2D_e0 = Array<OneD, Array<OneD, NekDouble> >(nElements);
166  m_Q2D_e1 = Array<OneD, Array<OneD, NekDouble> >(nElements);
167  m_Q2D_e2 = Array<OneD, Array<OneD, NekDouble> >(nElements);
168  m_Q2D_e3 = Array<OneD, Array<OneD, NekDouble> >(nElements);
169 
171 
172  for (n = 0; n < nElements; ++n)
173  {
174  base = pFields[0]->GetExp(n)->GetBase();
175  nquad0 = base[0]->GetNumPoints();
176  nquad1 = base[1]->GetNumPoints();
177 
178  m_Q2D_e0[n] = Array<OneD, NekDouble>(nquad0);
179  m_Q2D_e1[n] = Array<OneD, NekDouble>(nquad1);
180  m_Q2D_e2[n] = Array<OneD, NekDouble>(nquad0);
181  m_Q2D_e3[n] = Array<OneD, NekDouble>(nquad1);
182 
183  // Extract the Q factors at each edge point
184  pFields[0]->GetExp(n)->GetEdgeQFactors(
185  0, auxArray1 = m_Q2D_e0[n]);
186  pFields[0]->GetExp(n)->GetEdgeQFactors(
187  1, auxArray1 = m_Q2D_e1[n]);
188  pFields[0]->GetExp(n)->GetEdgeQFactors(
189  2, auxArray1 = m_Q2D_e2[n]);
190  pFields[0]->GetExp(n)->GetEdgeQFactors(
191  3, auxArray1 = m_Q2D_e3[n]);
192 
193  ptsKeys = pFields[0]->GetExp(n)->GetPointsKeys();
194  nLocalSolutionPts = pFields[0]->GetExp(n)->GetTotPoints();
195  phys_offset = pFields[0]->GetPhys_Offset(n);
196 
197  jac = pFields[0]->GetExp(n)
198  ->as<LocalRegions::Expansion2D>()->GetGeom2D()
199  ->GetMetricInfo()->GetJac(ptsKeys);
200  gmat = pFields[0]->GetExp(n)
201  ->as<LocalRegions::Expansion2D>()->GetGeom2D()
202  ->GetMetricInfo()->GetDerivFactors(ptsKeys);
203 
204  if (pFields[0]->GetExp(n)
205  ->as<LocalRegions::Expansion2D>()->GetGeom2D()
206  ->GetMetricInfo()->GetGtype()
208  {
209  for (i = 0; i < nLocalSolutionPts; ++i)
210  {
211  m_jac[i+phys_offset] = jac[i];
212  m_gmat[0][i+phys_offset] = gmat[0][i];
213  m_gmat[1][i+phys_offset] = gmat[1][i];
214  m_gmat[2][i+phys_offset] = gmat[2][i];
215  m_gmat[3][i+phys_offset] = gmat[3][i];
216  }
217  }
218  else
219  {
220  for (i = 0; i < nLocalSolutionPts; ++i)
221  {
222  m_jac[i+phys_offset] = jac[0];
223  m_gmat[0][i+phys_offset] = gmat[0][0];
224  m_gmat[1][i+phys_offset] = gmat[1][0];
225  m_gmat[2][i+phys_offset] = gmat[2][0];
226  m_gmat[3][i+phys_offset] = gmat[3][0];
227  }
228  }
229  }
230 
231  m_traceNormals = Array<OneD, Array<OneD, NekDouble> >(
232  nDimensions);
233  for(i = 0; i < nDimensions; ++i)
234  {
235  m_traceNormals[i] = Array<OneD, NekDouble> (nTracePts);
236  }
237  pFields[0]->GetTrace()->GetNormals(m_traceNormals);
238 
239  break;
240  }
241  case 3:
242  {
243  ASSERTL0(false,"3DFR Metric terms not implemented yet");
244  break;
245  }
246  default:
247  {
248  ASSERTL0(false, "Expansion dimension not recognised");
249  break;
250  }
251  }
252  }
Array< OneD, Array< OneD, NekDouble > > m_Q2D_e1
Definition: AdvectionFR.h:61
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:161
std::vector< PointsKey > PointsKeyVector
Definition: Points.h:220
Array< OneD, Array< OneD, NekDouble > > m_Q2D_e0
Definition: AdvectionFR.h:60
Array< OneD, Array< OneD, NekDouble > > m_gmat
Definition: AdvectionFR.h:58
Array< OneD, Array< OneD, NekDouble > > m_Q2D_e3
Definition: AdvectionFR.h:63
Array< OneD, Array< OneD, NekDouble > > m_Q2D_e2
Definition: AdvectionFR.h:62
Array< OneD, Array< OneD, NekDouble > > m_traceNormals
Definition: AdvectionFR.h:77
Geometry is curved or has non-constant factors.
Array< OneD, NekDouble > m_jac
Definition: AdvectionFR.h:57

Member Data Documentation

std::string Nektar::SolverUtils::AdvectionFR::m_advType
protected

Definition at line 79 of file AdvectionFR.h.

Referenced by v_SetupCFunctions().

Array<OneD, Array<OneD, NekDouble> > Nektar::SolverUtils::AdvectionFR::m_dGL_xi1

Definition at line 65 of file AdvectionFR.h.

Referenced by v_DivCFlux_1D(), v_DivCFlux_2D(), v_DivCFlux_2D_Gauss(), and v_SetupCFunctions().

Array<OneD, Array<OneD, NekDouble> > Nektar::SolverUtils::AdvectionFR::m_dGL_xi2

Definition at line 67 of file AdvectionFR.h.

Referenced by v_DivCFlux_2D(), v_DivCFlux_2D_Gauss(), and v_SetupCFunctions().

Array<OneD, Array<OneD, NekDouble> > Nektar::SolverUtils::AdvectionFR::m_dGL_xi3

Definition at line 69 of file AdvectionFR.h.

Referenced by v_SetupCFunctions().

Array<OneD, Array<OneD, NekDouble> > Nektar::SolverUtils::AdvectionFR::m_dGR_xi1

Definition at line 66 of file AdvectionFR.h.

Referenced by v_DivCFlux_1D(), v_DivCFlux_2D(), v_DivCFlux_2D_Gauss(), and v_SetupCFunctions().

Array<OneD, Array<OneD, NekDouble> > Nektar::SolverUtils::AdvectionFR::m_dGR_xi2

Definition at line 68 of file AdvectionFR.h.

Referenced by v_DivCFlux_2D(), v_DivCFlux_2D_Gauss(), and v_SetupCFunctions().

Array<OneD, Array<OneD, NekDouble> > Nektar::SolverUtils::AdvectionFR::m_dGR_xi3

Definition at line 70 of file AdvectionFR.h.

Referenced by v_SetupCFunctions().

Array<OneD, Array<OneD, NekDouble> > Nektar::SolverUtils::AdvectionFR::m_gmat

Definition at line 58 of file AdvectionFR.h.

Referenced by v_Advect(), and v_SetupMetrics().

DNekMatSharedPtr Nektar::SolverUtils::AdvectionFR::m_Ixm

Definition at line 71 of file AdvectionFR.h.

DNekMatSharedPtr Nektar::SolverUtils::AdvectionFR::m_Ixp

Definition at line 72 of file AdvectionFR.h.

Array<OneD, NekDouble> Nektar::SolverUtils::AdvectionFR::m_jac

Definition at line 57 of file AdvectionFR.h.

Referenced by v_Advect(), and v_SetupMetrics().

Array<OneD, Array<OneD, NekDouble> > Nektar::SolverUtils::AdvectionFR::m_Q2D_e0

Definition at line 60 of file AdvectionFR.h.

Referenced by v_DivCFlux_2D(), v_DivCFlux_2D_Gauss(), and v_SetupMetrics().

Array<OneD, Array<OneD, NekDouble> > Nektar::SolverUtils::AdvectionFR::m_Q2D_e1

Definition at line 61 of file AdvectionFR.h.

Referenced by v_DivCFlux_2D(), v_DivCFlux_2D_Gauss(), and v_SetupMetrics().

Array<OneD, Array<OneD, NekDouble> > Nektar::SolverUtils::AdvectionFR::m_Q2D_e2

Definition at line 62 of file AdvectionFR.h.

Referenced by v_DivCFlux_2D(), v_DivCFlux_2D_Gauss(), and v_SetupMetrics().

Array<OneD, Array<OneD, NekDouble> > Nektar::SolverUtils::AdvectionFR::m_Q2D_e3

Definition at line 63 of file AdvectionFR.h.

Referenced by v_DivCFlux_2D(), v_DivCFlux_2D_Gauss(), and v_SetupMetrics().

Array<OneD, Array<OneD, NekDouble> > Nektar::SolverUtils::AdvectionFR::m_traceNormals
protected

Definition at line 77 of file AdvectionFR.h.

Referenced by v_DivCFlux_2D(), v_DivCFlux_2D_Gauss(), and v_SetupMetrics().

std::string Nektar::SolverUtils::AdvectionFR::type
static