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 75 of file AdvectionFR.cpp.

Referenced by create().

75  :m_advType(advType)
76  {
77  }

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:75
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 799 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().

806  {
807  int i, j, n;
808  int phys_offset;
809 
810  Array<OneD, NekDouble> auxArray1, auxArray2;
811 
813  Basis = fields[0]->GetExp(0)->GetBase();
814 
815  int nElements = fields[0]->GetExpSize();
816  int nDimensions = fields[0]->GetCoordim(0);
817  int nSolutionPts = fields[0]->GetTotPoints();
818  int nTracePts = fields[0]->GetTrace()->GetTotPoints();
819  int nCoeffs = fields[0]->GetNcoeffs();
820 
821  // Storage for flux vector.
823  (nConvectiveFields);
824  // Outarray for Galerkin projection in case of primitive dealising
825  Array<OneD, Array<OneD, NekDouble> > outarrayCoeff
826  (nConvectiveFields);
827  // Store forwards/backwards space along trace space
828  Array<OneD, Array<OneD, NekDouble> > Fwd (nConvectiveFields);
829  Array<OneD, Array<OneD, NekDouble> > Bwd (nConvectiveFields);
830  Array<OneD, Array<OneD, NekDouble> > numflux(nConvectiveFields);
831 
832  // Set up storage for flux vector.
833  for (i = 0; i < nConvectiveFields; ++i)
834  {
835  fluxvector[i] =
837 
838  for (j = 0; j < m_spaceDim; ++j)
839  {
840  fluxvector[i][j] = Array<OneD, NekDouble>(nSolutionPts);
841  }
842  }
843 
844  for (i = 0; i < nConvectiveFields; ++i)
845  {
846  outarrayCoeff[i] = Array<OneD, NekDouble>(nCoeffs);
847  Fwd[i] = Array<OneD, NekDouble>(nTracePts);
848  Bwd[i] = Array<OneD, NekDouble>(nTracePts);
849  numflux[i] = Array<OneD, NekDouble>(nTracePts);
850  fields[i]->GetFwdBwdTracePhys(inarray[i], Fwd[i], Bwd[i]);
851  }
852 
853  // Computing the interface flux at each trace point
854  m_riemann->Solve(m_spaceDim, Fwd, Bwd, numflux);
855 
856  // Divergence of the flux (computing the RHS)
857  switch(nDimensions)
858  {
859  // 1D problems
860  case 1:
861  {
862  Array<OneD, NekDouble> DfluxvectorX1(nSolutionPts);
863  Array<OneD, NekDouble> divFC (nSolutionPts);
864 
865  // Get the advection flux vector (solver specific)
866  m_fluxVector(inarray, fluxvector);
867 
868  // Get the discontinuous flux divergence
869  for(i = 0; i < nConvectiveFields; ++i)
870  {
871  for (n = 0; n < nElements; n++)
872  {
873  phys_offset = fields[0]->GetPhys_Offset(n);
874 
875  fields[i]->GetExp(n)->PhysDeriv(
876  0, fluxvector[i][0] + phys_offset,
877  auxArray2 = DfluxvectorX1 + phys_offset);
878  }
879 
880  // Get the correction flux divergence
881  v_DivCFlux_1D(nConvectiveFields, fields,
882  fluxvector[i][0], numflux[i], divFC);
883 
884  // Back to the physical space using global operations
885  Vmath::Vdiv(nSolutionPts, &divFC[0], 1, &m_jac[0], 1,
886  &outarray[i][0], 1);
887 
888  // Adding the total divergence to outarray (RHS)
889  Vmath::Vadd(nSolutionPts, &outarray[i][0], 1,
890  &DfluxvectorX1[0], 1, &outarray[i][0], 1);
891 
892  // Primitive Dealiasing 1D
893  if (!(Basis[0]->Collocation()))
894  {
895  fields[i]->FwdTrans(outarray[i], outarrayCoeff[i]);
896  fields[i]->BwdTrans(outarrayCoeff[i], outarray[i]);
897  }
898  }
899  break;
900  }
901  // 2D problems
902  case 2:
903  {
904  Array<OneD, NekDouble> DfluxvectorX1(nSolutionPts);
905  Array<OneD, NekDouble> DfluxvectorX2(nSolutionPts);
906  Array<OneD, NekDouble> divFD(nSolutionPts);
907  Array<OneD, NekDouble> divFC(nSolutionPts);
908 
909  // Get the advection flux vector (solver specific)
910  m_fluxVector(inarray, fluxvector);
911 
912  for (i = 0; i < nConvectiveFields; ++i)
913  {
914  // Temporary vectors
915  Array<OneD, NekDouble> f_hat(nSolutionPts);
916  Array<OneD, NekDouble> g_hat(nSolutionPts);
917 
918  Vmath::Vvtvvtp(nSolutionPts,
919  &m_gmat[0][0], 1,
920  &fluxvector[i][0][0], 1,
921  &m_gmat[2][0], 1,
922  &fluxvector[i][1][0], 1,
923  &f_hat[0], 1);
924 
925  Vmath::Vmul(nSolutionPts, &m_jac[0], 1, &f_hat[0], 1,
926  &f_hat[0], 1);
927 
928  Vmath::Vvtvvtp(nSolutionPts,
929  &m_gmat[1][0], 1,
930  &fluxvector[i][0][0], 1,
931  &m_gmat[3][0], 1,
932  &fluxvector[i][1][0], 1,
933  &g_hat[0], 1);
934 
935  Vmath::Vmul(nSolutionPts, &m_jac[0], 1, &g_hat[0], 1,
936  &g_hat[0], 1);
937 
938  // Get the discontinuous flux derivatives
939  for (n = 0; n < nElements; n++)
940  {
941  phys_offset = fields[0]->GetPhys_Offset(n);
942  fields[0]->GetExp(n)->StdPhysDeriv(0,
943  auxArray1 = f_hat + phys_offset,
944  auxArray2 = DfluxvectorX1 + phys_offset);
945  fields[0]->GetExp(n)->StdPhysDeriv(1,
946  auxArray1 = g_hat + phys_offset,
947  auxArray2 = DfluxvectorX2 + phys_offset);
948  }
949 
950  // Divergence of the discontinuous flux
951  Vmath::Vadd(nSolutionPts, DfluxvectorX1, 1,
952  DfluxvectorX2, 1, divFD, 1);
953 
954  // Divergence of the correction flux
955  if (Basis[0]->GetPointsType() ==
957  Basis[1]->GetPointsType() ==
959  {
961  nConvectiveFields, fields, f_hat, g_hat,
962  numflux[i], divFC);
963  }
964  else
965  {
967  nConvectiveFields, fields,
968  fluxvector[i][0], fluxvector[i][1],
969  numflux[i], divFC);
970 
971  }
972 
973  // Divergence of the final flux
974  Vmath::Vadd(nSolutionPts, divFD, 1, divFC, 1,
975  outarray[i], 1);
976 
977  // Back to the physical space using a global operation
978  Vmath::Vdiv(nSolutionPts, &outarray[i][0], 1,
979  &m_jac[0], 1, &outarray[i][0], 1);
980 
981  // Primitive Dealiasing 2D
982  if (!(Basis[0]->Collocation()))
983  {
984  fields[i]->FwdTrans(outarray[i], outarrayCoeff[i]);
985  fields[i]->BwdTrans(outarrayCoeff[i], outarray[i]);
986  }
987  } // close nConvectiveFields loop
988  break;
989  }
990  // 3D problems
991  case 3:
992  {
993  ASSERTL0(false,"3D FRDG case not implemented yet");
994  break;
995  }
996  }
997  }
#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 1009 of file AdvectionFR.cpp.

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

Referenced by v_Advect().

1015  {
1016  int n;
1017  int nLocalSolutionPts, phys_offset, t_offset;
1018 
1020  Basis = fields[0]->GetExp(0)->GetBasis(0);
1021 
1022  int nElements = fields[0]->GetExpSize();
1023  int nSolutionPts = fields[0]->GetTotPoints();
1024 
1025 
1026  vector<bool> negatedFluxNormal = (boost::static_pointer_cast<MultiRegions::DisContField1D>(fields[0]))->GetNegatedFluxNormal();
1027 
1028  // Arrays to store the derivatives of the correction flux
1029  Array<OneD, NekDouble> DCL(nSolutionPts/nElements, 0.0);
1030  Array<OneD, NekDouble> DCR(nSolutionPts/nElements, 0.0);
1031 
1032  // Arrays to store the intercell numerical flux jumps
1033  Array<OneD, NekDouble> JumpL(nElements);
1034  Array<OneD, NekDouble> JumpR(nElements);
1035 
1037  &elmtToTrace = fields[0]->GetTraceMap()->GetElmtToTrace();
1038 
1039  for (n = 0; n < nElements; ++n)
1040  {
1041  nLocalSolutionPts = fields[0]->GetExp(n)->GetTotPoints();
1042  phys_offset = fields[0]->GetPhys_Offset(n);
1043 
1044  Array<OneD, NekDouble> tmparrayX1(nLocalSolutionPts, 0.0);
1045  NekDouble tmpFluxVertex = 0;
1046  Vmath::Vcopy(nLocalSolutionPts,
1047  &fluxX1[phys_offset], 1,
1048  &tmparrayX1[0], 1);
1049 
1050  fields[0]->GetExp(n)->GetVertexPhysVals(0, tmparrayX1,
1051  tmpFluxVertex);
1052 
1053  t_offset = fields[0]->GetTrace()
1054  ->GetPhys_Offset(elmtToTrace[n][0]->GetElmtId());
1055 
1056  if(negatedFluxNormal[2*n])
1057  {
1058  JumpL[n] = numericalFlux[t_offset] - tmpFluxVertex;
1059  }
1060  else
1061  {
1062  JumpL[n] = -numericalFlux[t_offset] - tmpFluxVertex;
1063  }
1064 
1065  fields[0]->GetExp(n)->GetVertexPhysVals(1, tmparrayX1,
1066  tmpFluxVertex);
1067 
1068  t_offset = fields[0]->GetTrace()
1069  ->GetPhys_Offset(elmtToTrace[n][1]->GetElmtId());
1070 
1071  if(negatedFluxNormal[2*n+1])
1072  {
1073  JumpR[n] = -numericalFlux[t_offset] - tmpFluxVertex;
1074  }
1075  else
1076  {
1077  JumpR[n] = numericalFlux[t_offset] - tmpFluxVertex;
1078  }
1079  }
1080 
1081  for (n = 0; n < nElements; ++n)
1082  {
1083  nLocalSolutionPts = fields[0]->GetExp(n)->GetTotPoints();
1084  phys_offset = fields[0]->GetPhys_Offset(n);
1085 
1086  // Left jump multiplied by left derivative of C function
1087  Vmath::Smul(nLocalSolutionPts, JumpL[n], &m_dGL_xi1[n][0], 1,
1088  &DCL[0], 1);
1089 
1090  // Right jump multiplied by right derivative of C function
1091  Vmath::Smul(nLocalSolutionPts, JumpR[n], &m_dGR_xi1[n][0], 1,
1092  &DCR[0], 1);
1093 
1094  // Assembling divergence of the correction flux
1095  Vmath::Vadd(nLocalSolutionPts, &DCL[0], 1, &DCR[0], 1,
1096  &divCFlux[phys_offset], 1);
1097  }
1098  }
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
This class is the abstraction of a global discontinuous two- dimensional spectral/hp element expansio...
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 1112 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().

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

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

1680  {
1681 
1682  }
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 89 of file AdvectionFR.cpp.

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

92  {
93  Advection::v_InitObject(pSession, pFields);
94  v_SetupMetrics (pSession, pFields);
95  v_SetupCFunctions (pSession, pFields);
96  }
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 268 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().

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

117  {
118  int i, n;
119  int nquad0, nquad1;
120  int phys_offset;
121  int nLocalSolutionPts;
122  int nElements = pFields[0]->GetExpSize();
123  int nDimensions = pFields[0]->GetCoordim(0);
124  int nSolutionPts = pFields[0]->GetTotPoints();
125  int nTracePts = pFields[0]->GetTrace()->GetTotPoints();
126 
129 
130  m_jac = Array<OneD, NekDouble>(nSolutionPts);
131 
132  Array<OneD, NekDouble> auxArray1;
135 
136  switch (nDimensions)
137  {
138  case 1:
139  {
140  for (n = 0; n < nElements; ++n)
141  {
142  ptsKeys = pFields[0]->GetExp(n)->GetPointsKeys();
143  nLocalSolutionPts = pFields[0]->GetExp(n)->GetTotPoints();
144  phys_offset = pFields[0]->GetPhys_Offset(n);
145  jac = pFields[0]->GetExp(n)
146  ->as<LocalRegions::Expansion1D>()->GetGeom1D()
147  ->GetMetricInfo()->GetJac(ptsKeys);
148  for (i = 0; i < nLocalSolutionPts; ++i)
149  {
150  m_jac[i+phys_offset] = jac[0];
151  }
152  }
153  break;
154  }
155  case 2:
156  {
158  m_gmat[0] = Array<OneD, NekDouble>(nSolutionPts);
159  m_gmat[1] = Array<OneD, NekDouble>(nSolutionPts);
160  m_gmat[2] = Array<OneD, NekDouble>(nSolutionPts);
161  m_gmat[3] = Array<OneD, NekDouble>(nSolutionPts);
162 
167 
169 
170  for (n = 0; n < nElements; ++n)
171  {
172  base = pFields[0]->GetExp(n)->GetBase();
173  nquad0 = base[0]->GetNumPoints();
174  nquad1 = base[1]->GetNumPoints();
175 
176  m_Q2D_e0[n] = Array<OneD, NekDouble>(nquad0);
177  m_Q2D_e1[n] = Array<OneD, NekDouble>(nquad1);
178  m_Q2D_e2[n] = Array<OneD, NekDouble>(nquad0);
179  m_Q2D_e3[n] = Array<OneD, NekDouble>(nquad1);
180 
181  // Extract the Q factors at each edge point
182  pFields[0]->GetExp(n)->GetEdgeQFactors(
183  0, auxArray1 = m_Q2D_e0[n]);
184  pFields[0]->GetExp(n)->GetEdgeQFactors(
185  1, auxArray1 = m_Q2D_e1[n]);
186  pFields[0]->GetExp(n)->GetEdgeQFactors(
187  2, auxArray1 = m_Q2D_e2[n]);
188  pFields[0]->GetExp(n)->GetEdgeQFactors(
189  3, auxArray1 = m_Q2D_e3[n]);
190 
191  ptsKeys = pFields[0]->GetExp(n)->GetPointsKeys();
192  nLocalSolutionPts = pFields[0]->GetExp(n)->GetTotPoints();
193  phys_offset = pFields[0]->GetPhys_Offset(n);
194 
195  jac = pFields[0]->GetExp(n)
196  ->as<LocalRegions::Expansion2D>()->GetGeom2D()
197  ->GetMetricInfo()->GetJac(ptsKeys);
198  gmat = pFields[0]->GetExp(n)
199  ->as<LocalRegions::Expansion2D>()->GetGeom2D()
200  ->GetMetricInfo()->GetDerivFactors(ptsKeys);
201 
202  if (pFields[0]->GetExp(n)
203  ->as<LocalRegions::Expansion2D>()->GetGeom2D()
204  ->GetMetricInfo()->GetGtype()
206  {
207  for (i = 0; i < nLocalSolutionPts; ++i)
208  {
209  m_jac[i+phys_offset] = jac[i];
210  m_gmat[0][i+phys_offset] = gmat[0][i];
211  m_gmat[1][i+phys_offset] = gmat[1][i];
212  m_gmat[2][i+phys_offset] = gmat[2][i];
213  m_gmat[3][i+phys_offset] = gmat[3][i];
214  }
215  }
216  else
217  {
218  for (i = 0; i < nLocalSolutionPts; ++i)
219  {
220  m_jac[i+phys_offset] = jac[0];
221  m_gmat[0][i+phys_offset] = gmat[0][0];
222  m_gmat[1][i+phys_offset] = gmat[1][0];
223  m_gmat[2][i+phys_offset] = gmat[2][0];
224  m_gmat[3][i+phys_offset] = gmat[3][0];
225  }
226  }
227  }
228 
230  nDimensions);
231  for(i = 0; i < nDimensions; ++i)
232  {
233  m_traceNormals[i] = Array<OneD, NekDouble> (nTracePts);
234  }
235  pFields[0]->GetTrace()->GetNormals(m_traceNormals);
236 
237  break;
238  }
239  case 3:
240  {
241  ASSERTL0(false,"3DFR Metric terms not implemented yet");
242  break;
243  }
244  default:
245  {
246  ASSERTL0(false, "Expansion dimension not recognised");
247  break;
248  }
249  }
250  }
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
const boost::shared_ptr< SpatialDomains::GeomFactors > & GetMetricInfo(void) const
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