Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties 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, const Array< OneD, Array< OneD, NekDouble > > &pFwd=NullNekDoubleArrayofArray, const Array< OneD, Array< OneD, NekDouble > > &pBwd=NullNekDoubleArrayofArray)
 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, const Array< OneD, MultiRegions::ExpListSharedPtr > &fields)
 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, const Array< OneD, Array< OneD, NekDouble > > &pFwd=NullNekDoubleArrayofArray, const Array< OneD, Array< OneD, NekDouble > > &pBwd=NullNekDoubleArrayofArray)
 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, const Array< OneD, MultiRegions::ExpListSharedPtr > &fields)
 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:165
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,
const Array< OneD, Array< OneD, NekDouble > > &  pFwd = NullNekDoubleArrayofArray,
const Array< OneD, Array< OneD, NekDouble > > &  pBwd = NullNekDoubleArrayofArray 
)
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().

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

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

Referenced by v_Advect().

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

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

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

1684  {
1685 
1686  }
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:100
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:198
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:2151
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:198
std::vector< PointsKey > PointsKeyVector
Definition: Points.h:242
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