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::DiffusionLFR Class Reference

#include <DiffusionLFR.h>

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

Static Public Member Functions

static DiffusionSharedPtr create (std::string diffType)
 

Public Attributes

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

 DiffusionLFR (std::string diffType)
 DiffusionLFR uses the Flux Reconstruction (FR) approach to compute the diffusion 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 DiffusionLFR objects and store them before starting the time-stepping. 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_Diffuse (const int nConvectiveFields, const Array< OneD, MultiRegions::ExpListSharedPtr > &fields, const Array< OneD, Array< OneD, NekDouble > > &inarray, Array< OneD, Array< OneD, NekDouble > > &outarray)
 Calculate FR Diffusion for the linear problems using an LDG interface flux. More...
 
virtual void v_NumFluxforScalar (const Array< OneD, MultiRegions::ExpListSharedPtr > &fields, const Array< OneD, Array< OneD, NekDouble > > &ufield, Array< OneD, Array< OneD, Array< OneD, NekDouble > > > &uflux)
 Builds the numerical flux for the 1st order derivatives. More...
 
virtual void v_WeakPenaltyforScalar (const Array< OneD, MultiRegions::ExpListSharedPtr > &fields, const int var, const Array< OneD, const NekDouble > &ufield, Array< OneD, NekDouble > &penaltyflux)
 Imposes appropriate bcs for the 1st order derivatives. More...
 
virtual void v_NumFluxforVector (const Array< OneD, MultiRegions::ExpListSharedPtr > &fields, const Array< OneD, Array< OneD, NekDouble > > &ufield, Array< OneD, Array< OneD, Array< OneD, NekDouble > > > &qfield, Array< OneD, Array< OneD, NekDouble > > &qflux)
 Build the numerical flux for the 2nd order derivatives. More...
 
virtual void v_WeakPenaltyforVector (const Array< OneD, MultiRegions::ExpListSharedPtr > &fields, const int var, const int dir, const Array< OneD, const NekDouble > &qfield, Array< OneD, NekDouble > &penaltyflux, NekDouble C11)
 Imposes appropriate bcs for the 2nd order derivatives. More...
 
virtual void v_DerCFlux_1D (const int nConvectiveFields, const Array< OneD, MultiRegions::ExpListSharedPtr > &fields, const Array< OneD, const NekDouble > &flux, const Array< OneD, const NekDouble > &iFlux, Array< OneD, NekDouble > &derCFlux)
 Compute the derivative of the corrective flux for 1D problems. More...
 
virtual void v_DerCFlux_2D (const int nConvectiveFields, const int direction, const Array< OneD, MultiRegions::ExpListSharedPtr > &fields, const Array< OneD, const NekDouble > &flux, const Array< OneD, NekDouble > &iFlux, Array< OneD, NekDouble > &derCFlux)
 Compute the derivative of the corrective flux wrt a given coordinate for 2D 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...
 
- Protected Member Functions inherited from Nektar::SolverUtils::Diffusion
virtual void v_SetHomoDerivs (Array< OneD, Array< OneD, NekDouble > > &deriv)
 
virtual Array< OneD, Array
< OneD, Array< OneD, NekDouble > > > & 		v_GetFluxTensor ()
 

Protected Attributes

Array< OneD, Array< OneD,
NekDouble > > 		m_traceNormals
 
LibUtilities::SessionReaderSharedPtr m_session
 
Array< OneD, Array< OneD,
Array< OneD, NekDouble > > > 		m_IF1
 
Array< OneD, Array< OneD,
Array< OneD, NekDouble > > > 		m_DU1
 
Array< OneD, Array< OneD,
Array< OneD, NekDouble > > > 		m_DFC1
 
Array< OneD, Array< OneD,
Array< OneD, NekDouble > > > 		m_BD1
 
Array< OneD, Array< OneD,
Array< OneD, NekDouble > > > 		m_D1
 
Array< OneD, Array< OneD,
Array< OneD, NekDouble > > > 		m_DD1
 
Array< OneD, Array< OneD,
NekDouble > > 		m_IF2
 
Array< OneD, Array< OneD,
Array< OneD, NekDouble > > > 		m_DFC2
 
Array< OneD, Array< OneD,
NekDouble > > 		m_divFD
 
Array< OneD, Array< OneD,
NekDouble > > 		m_divFC
 
Array< OneD, Array< OneD,
Array< OneD, NekDouble > > > 		m_tmp1
 
Array< OneD, Array< OneD,
Array< OneD, NekDouble > > > 		m_tmp2
 
std::string m_diffType
 
- Protected Attributes inherited from Nektar::SolverUtils::Diffusion
DiffusionFluxVecCB m_fluxVector
 
DiffusionFluxVecCBNS m_fluxVectorNS
 
RiemannSolverSharedPtr m_riemann
 
DiffusionArtificialDiffusion m_ArtificialDiffusionVector
 

Additional Inherited Members

- Public Member Functions inherited from Nektar::SolverUtils::Diffusion
SOLVER_UTILS_EXPORT void InitObject (LibUtilities::SessionReaderSharedPtr pSession, Array< OneD, MultiRegions::ExpListSharedPtr > pFields)
 
SOLVER_UTILS_EXPORT void Diffuse (const int nConvectiveFields, const Array< OneD, MultiRegions::ExpListSharedPtr > &fields, const Array< OneD, Array< OneD, NekDouble > > &inarray, Array< OneD, Array< OneD, NekDouble > > &outarray)
 
SOLVER_UTILS_EXPORT void FluxVec (Array< OneD, Array< OneD, Array< OneD, NekDouble > > > &fluxvector)
 
template<typename FuncPointerT , typename ObjectPointerT >
void SetFluxVector (FuncPointerT func, ObjectPointerT obj)
 
void SetFluxVectorVec (DiffusionFluxVecCB fluxVector)
 
template<typename FuncPointerT , typename ObjectPointerT >
void SetFluxVectorNS (FuncPointerT func, ObjectPointerT obj)
 
template<typename FuncPointerT , typename ObjectPointerT >
void SetArtificialDiffusionVector (FuncPointerT func, ObjectPointerT obj)
 
void SetFluxVectorNS (DiffusionFluxVecCBNS fluxVector)
 
void SetRiemannSolver (RiemannSolverSharedPtr riemann)
 
void SetHomoDerivs (Array< OneD, Array< OneD, NekDouble > > &deriv)
 
virtual Array< OneD, Array
< OneD, Array< OneD, NekDouble > > > & 		GetFluxTensor ()
 

Detailed Description

Definition at line 45 of file DiffusionLFR.h.

Constructor & Destructor Documentation

Nektar::SolverUtils::DiffusionLFR::DiffusionLFR ( std::string  diffType)
protected

DiffusionLFR uses the Flux Reconstruction (FR) approach to compute the diffusion 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 71 of file DiffusionLFR.cpp.

Referenced by create().

71  :m_diffType(diffType)
72  {
73  }

Member Function Documentation

static DiffusionSharedPtr Nektar::SolverUtils::DiffusionLFR::create ( std::string  diffType)
inlinestatic

Definition at line 48 of file DiffusionLFR.h.

References DiffusionLFR().

49  {
50  return DiffusionSharedPtr(new DiffusionLFR(diffType));
51  }
Nektar::SolverUtils::DiffusionSharedPtr
boost::shared_ptr< Diffusion > DiffusionSharedPtr
A shared pointer to an EquationSystem object.
Definition: Diffusion.h:164
Nektar::SolverUtils::DiffusionLFR::DiffusionLFR
DiffusionLFR(std::string diffType)
DiffusionLFR uses the Flux Reconstruction (FR) approach to compute the diffusion term. The implementation is only for segments, quadrilaterals and hexahedra at the moment.
Definition: DiffusionLFR.cpp:71
void Nektar::SolverUtils::DiffusionLFR::v_DerCFlux_1D ( const int  nConvectiveFields,
const Array< OneD, MultiRegions::ExpListSharedPtr > &  fields,
const Array< OneD, const NekDouble > &  flux,
const Array< OneD, const NekDouble > &  iFlux,
Array< OneD, NekDouble > &  derCFlux 
)
protectedvirtual

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

Parameters
nConvectiveFieldsNumber of fields.
fieldsPointer to fields.
fluxVolumetric flux in the physical space.
iFluxNumerical interface flux in physical space.
derCFluxDerivative of the corrective flux.

Definition at line 1462 of file DiffusionLFR.cpp.

References Nektar::StdRegions::StdExpansion::GetMetricInfo(), m_dGL_xi1, m_dGR_xi1, Vmath::Smul(), Vmath::Vadd(), and Vmath::Vcopy().

Referenced by v_Diffuse().

1468  {
1469  int n;
1470  int nLocalSolutionPts, phys_offset, t_offset;
1471 
1472  Array<OneD, NekDouble> auxArray1, auxArray2;
1473  Array<TwoD, const NekDouble> gmat;
1474  Array<OneD, const NekDouble> jac;
1475 
1476  LibUtilities::PointsKeyVector ptsKeys;
1477  LibUtilities::BasisSharedPtr Basis;
1478  Basis = fields[0]->GetExp(0)->GetBasis(0);
1479 
1480  int nElements = fields[0]->GetExpSize();
1481  int nSolutionPts = fields[0]->GetTotPoints();
1482 
1483  vector<bool> negatedFluxNormal = (boost::static_pointer_cast<MultiRegions::DisContField1D>(fields[0]))->GetNegatedFluxNormal();
1484 
1485  // Arrays to store the derivatives of the correction flux
1486  Array<OneD, NekDouble> DCL(nSolutionPts/nElements, 0.0);
1487  Array<OneD, NekDouble> DCR(nSolutionPts/nElements, 0.0);
1488 
1489  // Arrays to store the intercell numerical flux jumps
1490  Array<OneD, NekDouble> JumpL(nElements);
1491  Array<OneD, NekDouble> JumpR(nElements);
1492 
1493  Array<OneD, Array<OneD, LocalRegions::ExpansionSharedPtr> >
1494  &elmtToTrace = fields[0]->GetTraceMap()->GetElmtToTrace();
1495 
1496  for (n = 0; n < nElements; ++n)
1497  {
1498  nLocalSolutionPts = fields[0]->GetExp(n)->GetTotPoints();
1499  phys_offset = fields[0]->GetPhys_Offset(n);
1500 
1501  Array<OneD, NekDouble> tmparrayX1(nLocalSolutionPts, 0.0);
1502  NekDouble tmpFluxVertex = 0;
1503  Vmath::Vcopy(nLocalSolutionPts,
1504  &flux[phys_offset], 1,
1505  &tmparrayX1[0], 1);
1506 
1507  fields[0]->GetExp(n)->GetVertexPhysVals(0, tmparrayX1,
1508  tmpFluxVertex);
1509 
1510  t_offset = fields[0]->GetTrace()
1511  ->GetPhys_Offset(elmtToTrace[n][0]->GetElmtId());
1512 
1513  if(negatedFluxNormal[2*n])
1514  {
1515  JumpL[n] = iFlux[t_offset] - tmpFluxVertex;
1516  }
1517  else
1518  {
1519  JumpL[n] = -iFlux[t_offset] - tmpFluxVertex;
1520  }
1521 
1522  t_offset = fields[0]->GetTrace()
1523  ->GetPhys_Offset(elmtToTrace[n][1]->GetElmtId());
1524 
1525  fields[0]->GetExp(n)->GetVertexPhysVals(1, tmparrayX1,
1526  tmpFluxVertex);
1527  if(negatedFluxNormal[2*n+1])
1528  {
1529  JumpR[n] = -iFlux[t_offset] - tmpFluxVertex;
1530  }
1531  else
1532  {
1533  JumpR[n] = iFlux[t_offset] - tmpFluxVertex;
1534  }
1535  }
1536 
1537  for (n = 0; n < nElements; ++n)
1538  {
1539  ptsKeys = fields[0]->GetExp(n)->GetPointsKeys();
1540  nLocalSolutionPts = fields[0]->GetExp(n)->GetTotPoints();
1541  phys_offset = fields[0]->GetPhys_Offset(n);
1542  jac = fields[0]->GetExp(n)
1543  ->as<LocalRegions::Expansion1D>()->GetGeom1D()
1544  ->GetMetricInfo()->GetJac(ptsKeys);
1545 
1546  JumpL[n] = JumpL[n] * jac[0];
1547  JumpR[n] = JumpR[n] * jac[0];
1548 
1549  // Left jump multiplied by left derivative of C function
1550  Vmath::Smul(nLocalSolutionPts, JumpL[n],
1551  &m_dGL_xi1[n][0], 1, &DCL[0], 1);
1552 
1553  // Right jump multiplied by right derivative of C function
1554  Vmath::Smul(nLocalSolutionPts, JumpR[n],
1555  &m_dGR_xi1[n][0], 1, &DCR[0], 1);
1556 
1557  // Assembling derivative of the correction flux
1558  Vmath::Vadd(nLocalSolutionPts, &DCL[0], 1, &DCR[0], 1,
1559  &derCFlux[phys_offset], 1);
1560  }
1561  }
Nektar::SolverUtils::DiffusionLFR::m_dGL_xi1
Array< OneD, Array< OneD, NekDouble > > m_dGL_xi1
Definition: DiffusionLFR.h:63
Nektar::LibUtilities::PointsKeyVector
std::vector< PointsKey > PointsKeyVector
Definition: Points.h:220
Nektar::SolverUtils::DiffusionLFR::m_dGR_xi1
Array< OneD, Array< OneD, NekDouble > > m_dGR_xi1
Definition: DiffusionLFR.h:64
Vmath::Smul
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
Nektar::NekDouble
double NekDouble
Definition: NektarUnivTypeDefs.hpp:44
Nektar::LibUtilities::BasisSharedPtr
boost::shared_ptr< Basis > BasisSharedPtr
Definition: FoundationsFwd.hpp:79
Vmath::Vcopy
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.cpp:1047
Vmath::Vadd
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::DiffusionLFR::v_DerCFlux_2D ( const int  nConvectiveFields,
const int  direction,
const Array< OneD, MultiRegions::ExpListSharedPtr > &  fields,
const Array< OneD, const NekDouble > &  flux,
const Array< OneD, NekDouble > &  iFlux,
Array< OneD, NekDouble > &  derCFlux 
)
protectedvirtual

Compute the derivative of the corrective flux wrt a given coordinate for 2D problems.

There could be a bug for deformed elements since first derivatives are not exactly the same results of DiffusionLDG as expected

Parameters
nConvectiveFieldsNumber of fields.
fieldsPointer to fields.
fluxVolumetric flux in the physical space.
numericalFluxNumerical interface flux in physical space.
derCFluxDerivative of the corrective flux.
Todo:
: Switch on shapes eventually here.

Definition at line 1578 of file DiffusionLFR.cpp.

References ASSERTL0, Nektar::StdRegions::eBackwards, Nektar::SpatialDomains::eDeformed, Nektar::StdRegions::StdExpansion::GetMetricInfo(), m_dGL_xi1, m_dGL_xi2, m_dGR_xi1, m_dGR_xi2, Vmath::Reverse(), Vmath::Smul(), Vmath::Vadd(), and Vmath::Vsub().

Referenced by v_Diffuse().

1585  {
1586  int n, e, i, j, cnt;
1587 
1588  Array<TwoD, const NekDouble> gmat;
1589  Array<OneD, const NekDouble> jac;
1590 
1591  int nElements = fields[0]->GetExpSize();
1592 
1593  int trace_offset, phys_offset;
1594  int nLocalSolutionPts;
1595  int nquad0, nquad1;
1596  int nEdgePts;
1597 
1598  LibUtilities::PointsKeyVector ptsKeys;
1599  Array<OneD, NekDouble> auxArray1, auxArray2;
1600  Array<OneD, LibUtilities::BasisSharedPtr> base;
1601  Array<OneD, Array<OneD, LocalRegions::ExpansionSharedPtr> >
1602  &elmtToTrace = fields[0]->GetTraceMap()->GetElmtToTrace();
1603 
1604  // Loop on the elements
1605  for (n = 0; n < nElements; ++n)
1606  {
1607  // Offset of the element on the global vector
1608  phys_offset = fields[0]->GetPhys_Offset(n);
1609  nLocalSolutionPts = fields[0]->GetExp(n)->GetTotPoints();
1610  ptsKeys = fields[0]->GetExp(n)->GetPointsKeys();
1611 
1612  jac = fields[0]->GetExp(n)->as<LocalRegions::Expansion2D>()
1613  ->GetGeom2D()->GetMetricInfo()->GetJac(ptsKeys);
1614 
1615  base = fields[0]->GetExp(n)->GetBase();
1616  nquad0 = base[0]->GetNumPoints();
1617  nquad1 = base[1]->GetNumPoints();
1618 
1619  Array<OneD, NekDouble> divCFluxE0(nLocalSolutionPts, 0.0);
1620  Array<OneD, NekDouble> divCFluxE1(nLocalSolutionPts, 0.0);
1621  Array<OneD, NekDouble> divCFluxE2(nLocalSolutionPts, 0.0);
1622  Array<OneD, NekDouble> divCFluxE3(nLocalSolutionPts, 0.0);
1623 
1624  // Loop on the edges
1625  for (e = 0; e < fields[0]->GetExp(n)->GetNedges(); ++e)
1626  {
1627  // Number of edge points of edge 'e'
1628  nEdgePts = fields[0]->GetExp(n)->GetEdgeNumPoints(e);
1629 
1630  // Array for storing volumetric fluxes on each edge
1631  Array<OneD, NekDouble> tmparray(nEdgePts, 0.0);
1632 
1633  // Offset of the trace space correspondent to edge 'e'
1634  trace_offset = fields[0]->GetTrace()->GetPhys_Offset(
1635  elmtToTrace[n][e]->GetElmtId());
1636 
1637  // Get the normals of edge 'e'
1638  //const Array<OneD, const Array<OneD, NekDouble> > &normals =
1639  //fields[0]->GetExp(n)->GetEdgeNormal(e);
1640 
1641  // Extract the edge values of the volumetric fluxes
1642  // on edge 'e' and order them accordingly to the order
1643  // of the trace space
1644  fields[0]->GetExp(n)->GetEdgePhysVals(e, elmtToTrace[n][e],
1645  flux + phys_offset,
1646  auxArray1 = tmparray);
1647 
1648  // Compute the fluxJumps per each edge 'e' and each
1649  // flux point
1650  Array<OneD, NekDouble> fluxJumps(nEdgePts, 0.0);
1651  Vmath::Vsub(nEdgePts, &iFlux[trace_offset], 1,
1652  &tmparray[0], 1, &fluxJumps[0], 1);
1653 
1654  // Check the ordering of the fluxJumps and reverse
1655  // it in case of backward definition of edge 'e'
1656  if (fields[0]->GetExp(n)->GetEorient(e) ==
1657  StdRegions::eBackwards)
1658  {
1659  Vmath::Reverse(nEdgePts,
1660  &fluxJumps[0], 1,
1661  &fluxJumps[0], 1);
1662  }
1663 
1664  // Deformed elements
1665  if (fields[0]->GetExp(n)->as<LocalRegions::Expansion2D>()
1666  ->GetGeom2D()->GetMetricInfo()->GetGtype()
1667  == SpatialDomains::eDeformed)
1668  {
1669  // Extract the Jacobians along edge 'e'
1670  Array<OneD, NekDouble> jacEdge(nEdgePts, 0.0);
1671 
1672  fields[0]->GetExp(n)->GetEdgePhysVals(
1673  e, elmtToTrace[n][e],
1674  jac, auxArray1 = jacEdge);
1675 
1676  // Check the ordering of the fluxJumps and reverse
1677  // it in case of backward definition of edge 'e'
1678  if (fields[0]->GetExp(n)->GetEorient(e) ==
1679  StdRegions::eBackwards)
1680  {
1681  Vmath::Reverse(nEdgePts,
1682  &jacEdge[0], 1,
1683  &jacEdge[0], 1);
1684  }
1685 
1686  // Multiply the fluxJumps by the edge 'e' Jacobians
1687  // to bring the fluxJumps into the standard space
1688  for (j = 0; j < nEdgePts; j++)
1689  {
1690  fluxJumps[j] = fluxJumps[j] * jacEdge[j];
1691  }
1692  }
1693  // Non-deformed elements
1694  else
1695  {
1696  // Multiply the fluxJumps by the edge 'e' Jacobians
1697  // to bring the fluxJumps into the standard space
1698  Vmath::Smul(nEdgePts, jac[0], fluxJumps, 1,
1699  fluxJumps, 1);
1700  }
1701 
1702  // Multiply jumps by derivatives of the correction functions
1703  // All the quntities at this level should be defined into
1704  // the standard space
1705  switch (e)
1706  {
1707  case 0:
1708  for (i = 0; i < nquad0; ++i)
1709  {
1710  // Multiply fluxJumps by Q factors
1711  //fluxJumps[i] = -(m_Q2D_e0[n][i]) * fluxJumps[i];
1712 
1713  for (j = 0; j < nquad1; ++j)
1714  {
1715  cnt = i + j*nquad0;
1716  divCFluxE0[cnt] = fluxJumps[i] * m_dGL_xi2[n][j];
1717  }
1718  }
1719  break;
1720  case 1:
1721  for (i = 0; i < nquad1; ++i)
1722  {
1723  // Multiply fluxJumps by Q factors
1724  //fluxJumps[i] = (m_Q2D_e1[n][i]) * fluxJumps[i];
1725 
1726  for (j = 0; j < nquad0; ++j)
1727  {
1728  cnt = (nquad0)*i + j;
1729  divCFluxE1[cnt] = fluxJumps[i] * m_dGR_xi1[n][j];
1730  }
1731  }
1732  break;
1733  case 2:
1734  for (i = 0; i < nquad0; ++i)
1735  {
1736  // Multiply fluxJumps by Q factors
1737  //fluxJumps[i] = (m_Q2D_e2[n][i]) * fluxJumps[i];
1738 
1739  for (j = 0; j < nquad1; ++j)
1740  {
1741  cnt = (nquad0 - 1) + j*nquad0 - i;
1742  divCFluxE2[cnt] = fluxJumps[i] * m_dGR_xi2[n][j];
1743  }
1744  }
1745  break;
1746  case 3:
1747  for (i = 0; i < nquad1; ++i)
1748  {
1749  // Multiply fluxJumps by Q factors
1750  //fluxJumps[i] = -(m_Q2D_e3[n][i]) * fluxJumps[i];
1751  for (j = 0; j < nquad0; ++j)
1752  {
1753  cnt = (nquad0*nquad1 - nquad0) + j - i*nquad0;
1754  divCFluxE3[cnt] = fluxJumps[i] * m_dGL_xi1[n][j];
1755  }
1756  }
1757  break;
1758 
1759  default:
1760  ASSERTL0(false, "edge value (< 3) is out of range");
1761  break;
1762  }
1763  }
1764 
1765  // Summing the component of the flux for calculating the
1766  // derivatives per each direction
1767  if (direction == 0)
1768  {
1769  Vmath::Vadd(nLocalSolutionPts, &divCFluxE1[0], 1,
1770  &divCFluxE3[0], 1, &derCFlux[phys_offset], 1);
1771  }
1772  else if (direction == 1)
1773  {
1774  Vmath::Vadd(nLocalSolutionPts, &divCFluxE0[0], 1,
1775  &divCFluxE2[0], 1, &derCFlux[phys_offset], 1);
1776  }
1777  }
1778  }
Nektar::SolverUtils::DiffusionLFR::m_dGL_xi1
Array< OneD, Array< OneD, NekDouble > > m_dGL_xi1
Definition: DiffusionLFR.h:63
ASSERTL0
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:161
Nektar::LibUtilities::PointsKeyVector
std::vector< PointsKey > PointsKeyVector
Definition: Points.h:220
Nektar::SolverUtils::DiffusionLFR::m_dGR_xi1
Array< OneD, Array< OneD, NekDouble > > m_dGR_xi1
Definition: DiffusionLFR.h:64
Nektar::SolverUtils::DiffusionLFR::m_dGL_xi2
Array< OneD, Array< OneD, NekDouble > > m_dGL_xi2
Definition: DiffusionLFR.h:65
Nektar::SolverUtils::DiffusionLFR::m_dGR_xi2
Array< OneD, Array< OneD, NekDouble > > m_dGR_xi2
Definition: DiffusionLFR.h:66
Vmath::Reverse
void Reverse(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.cpp:1071
Vmath::Smul
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
Vmath::Vsub
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
Nektar::SpatialDomains::eDeformed
Geometry is curved or has non-constant factors.
Definition: SpatialDomains.hpp:57
Vmath::Vadd
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::DiffusionLFR::v_Diffuse ( const int  nConvectiveFields,
const Array< OneD, MultiRegions::ExpListSharedPtr > &  fields,
const Array< OneD, Array< OneD, NekDouble > > &  inarray,
Array< OneD, Array< OneD, NekDouble > > &  outarray 
)
protectedvirtual

Calculate FR Diffusion for the linear problems using an LDG interface flux.

Implements Nektar::SolverUtils::Diffusion.

Definition at line 849 of file DiffusionLFR.cpp.

References ASSERTL0, Nektar::LibUtilities::eGaussGaussLegendre, m_BD1, m_D1, m_DD1, m_DFC1, m_DFC2, m_divFC, m_divFD, m_DU1, m_gmat, m_IF1, m_IF2, m_jac, m_tmp1, m_tmp2, v_DerCFlux_1D(), v_DerCFlux_2D(), v_DivCFlux_2D(), v_DivCFlux_2D_Gauss(), v_NumFluxforScalar(), v_NumFluxforVector(), Vmath::Vadd(), Vmath::Vcopy(), Vmath::Vdiv(), and Vmath::Vmul().

854  {
855  int i, j, n;
856  int phys_offset;
857  //Array<TwoD, const NekDouble> gmat;
858  //Array<OneD, const NekDouble> jac;
859  Array<OneD, NekDouble> auxArray1, auxArray2;
860 
861  Array<OneD, LibUtilities::BasisSharedPtr> Basis;
862  Basis = fields[0]->GetExp(0)->GetBase();
863 
864  int nElements = fields[0]->GetExpSize();
865  int nDim = fields[0]->GetCoordim(0);
866  int nSolutionPts = fields[0]->GetTotPoints();
867  int nCoeffs = fields[0]->GetNcoeffs();
868 
869  Array<OneD, Array<OneD, NekDouble> > outarrayCoeff(nConvectiveFields);
870  for (i = 0; i < nConvectiveFields; ++i)
871  {
872  outarrayCoeff[i] = Array<OneD, NekDouble>(nCoeffs);
873  }
874 
875  // Compute interface numerical fluxes for inarray in physical space
876  v_NumFluxforScalar(fields, inarray, m_IF1);
877 
878  switch(nDim)
879  {
880  // 1D problems
881  case 1:
882  {
883  for (i = 0; i < nConvectiveFields; ++i)
884  {
885  // Computing the physical first-order discountinuous
886  // derivative
887  for (n = 0; n < nElements; n++)
888  {
889  phys_offset = fields[0]->GetPhys_Offset(n);
890 
891  fields[i]->GetExp(n)->PhysDeriv(0,
892  auxArray1 = inarray[i] + phys_offset,
893  auxArray2 = m_DU1[i][0] + phys_offset);
894  }
895 
896  // Computing the standard first-order correction
897  // derivative
898  v_DerCFlux_1D(nConvectiveFields, fields, inarray[i],
899  m_IF1[i][0], m_DFC1[i][0]);
900 
901  // Back to the physical space using global operations
902  Vmath::Vdiv(nSolutionPts, &m_DFC1[i][0][0], 1,
903  &m_jac[0], 1, &m_DFC1[i][0][0], 1);
904  Vmath::Vdiv(nSolutionPts, &m_DFC1[i][0][0], 1,
905  &m_jac[0], 1, &m_DFC1[i][0][0], 1);
906 
907  // Computing total first order derivatives
908  Vmath::Vadd(nSolutionPts, &m_DFC1[i][0][0], 1,
909  &m_DU1[i][0][0], 1, &m_D1[i][0][0], 1);
910 
911  Vmath::Vcopy(nSolutionPts, &m_D1[i][0][0], 1,
912  &m_tmp1[i][0][0], 1);
913  }
914 
915  // Computing interface numerical fluxes for m_D1
916  // in physical space
917  v_NumFluxforVector(fields, inarray, m_tmp1, m_IF2);
918 
919  for (i = 0; i < nConvectiveFields; ++i)
920  {
921  // Computing the physical second-order discountinuous
922  // derivative
923  for (n = 0; n < nElements; n++)
924  {
925  phys_offset = fields[0]->GetPhys_Offset(n);
926 
927  fields[i]->GetExp(n)->PhysDeriv(0,
928  auxArray1 = m_D1[i][0] + phys_offset,
929  auxArray2 = m_DD1[i][0] + phys_offset);
930  }
931 
932  // Computing the standard second-order correction
933  // derivative
934  v_DerCFlux_1D(nConvectiveFields, fields, m_D1[i][0],
935  m_IF2[i], m_DFC2[i][0]);
936 
937  // Back to the physical space using global operations
938  Vmath::Vdiv(nSolutionPts, &m_DFC2[i][0][0], 1,
939  &m_jac[0], 1, &m_DFC2[i][0][0], 1);
940  Vmath::Vdiv(nSolutionPts, &m_DFC2[i][0][0], 1,
941  &m_jac[0], 1, &m_DFC2[i][0][0], 1);
942 
943  // Adding the total divergence to outarray (RHS)
944  Vmath::Vadd(nSolutionPts, &m_DFC2[i][0][0], 1,
945  &m_DD1[i][0][0], 1, &outarray[i][0], 1);
946 
947  // Primitive Dealiasing 1D
948  if (!(Basis[0]->Collocation()))
949  {
950  fields[i]->FwdTrans(outarray[i], outarrayCoeff[i]);
951  fields[i]->BwdTrans(outarrayCoeff[i], outarray[i]);
952  }
953  }
954  break;
955  }
956  // 2D problems
957  case 2:
958  {
959  for(i = 0; i < nConvectiveFields; ++i)
960  {
961  for (j = 0; j < nDim; ++j)
962  {
963  // Temporary vectors
964  Array<OneD, NekDouble> u1_hat(nSolutionPts, 0.0);
965  Array<OneD, NekDouble> u2_hat(nSolutionPts, 0.0);
966 
967  if (j == 0)
968  {
969  Vmath::Vmul(nSolutionPts, &inarray[i][0], 1,
970  &m_gmat[0][0], 1, &u1_hat[0], 1);
971 
972  Vmath::Vmul(nSolutionPts, &u1_hat[0], 1,
973  &m_jac[0], 1, &u1_hat[0], 1);
974 
975  Vmath::Vmul(nSolutionPts, &inarray[i][0], 1,
976  &m_gmat[1][0], 1, &u2_hat[0], 1);
977 
978  Vmath::Vmul(nSolutionPts, &u2_hat[0], 1,
979  &m_jac[0], 1, &u2_hat[0], 1);
980  }
981  else if (j == 1)
982  {
983  Vmath::Vmul(nSolutionPts, &inarray[i][0], 1,
984  &m_gmat[2][0], 1, &u1_hat[0], 1);
985 
986  Vmath::Vmul(nSolutionPts, &u1_hat[0], 1,
987  &m_jac[0], 1, &u1_hat[0], 1);
988 
989  Vmath::Vmul(nSolutionPts, &inarray[i][0], 1,
990  &m_gmat[3][0], 1, &u2_hat[0], 1);
991 
992  Vmath::Vmul(nSolutionPts, &u2_hat[0], 1,
993  &m_jac[0], 1, &u2_hat[0], 1);
994  }
995 
996  for (n = 0; n < nElements; n++)
997  {
998  phys_offset = fields[0]->GetPhys_Offset(n);
999 
1000  fields[i]->GetExp(n)->StdPhysDeriv(0,
1001  auxArray1 = u1_hat + phys_offset,
1002  auxArray2 = m_tmp1[i][j] + phys_offset);
1003 
1004  fields[i]->GetExp(n)->StdPhysDeriv(1,
1005  auxArray1 = u2_hat + phys_offset,
1006  auxArray2 = m_tmp2[i][j] + phys_offset);
1007  }
1008 
1009  Vmath::Vadd(nSolutionPts, &m_tmp1[i][j][0], 1,
1010  &m_tmp2[i][j][0], 1,
1011  &m_DU1[i][j][0], 1);
1012 
1013  // Divide by the metric jacobian
1014  Vmath::Vdiv(nSolutionPts, &m_DU1[i][j][0], 1,
1015  &m_jac[0], 1, &m_DU1[i][j][0], 1);
1016 
1017  // Computing the standard first-order correction
1018  // derivatives
1019  v_DerCFlux_2D(nConvectiveFields, j, fields,
1020  inarray[i], m_IF1[i][j],
1021  m_DFC1[i][j]);
1022  }
1023 
1024  // Multiplying derivatives by B matrix to get auxiliary
1025  // variables
1026  for (j = 0; j < nSolutionPts; ++j)
1027  {
1028  // std(q1)
1029  m_BD1[i][0][j] =
1030  (m_gmat[0][j]*m_gmat[0][j] +
1031  m_gmat[2][j]*m_gmat[2][j]) *
1032  m_DFC1[i][0][j] +
1033  (m_gmat[1][j]*m_gmat[0][j] +
1034  m_gmat[3][j]*m_gmat[2][j]) *
1035  m_DFC1[i][1][j];
1036 
1037  // std(q2)
1038  m_BD1[i][1][j] =
1039  (m_gmat[1][j]*m_gmat[0][j] +
1040  m_gmat[3][j]*m_gmat[2][j]) *
1041  m_DFC1[i][0][j] +
1042  (m_gmat[1][j]*m_gmat[1][j] +
1043  m_gmat[3][j]*m_gmat[3][j]) *
1044  m_DFC1[i][1][j];
1045  }
1046 
1047  // Multiplying derivatives by A^(-1) to get back
1048  // into the physical space
1049  for (j = 0; j < nSolutionPts; j++)
1050  {
1051  // q1 = A11^(-1)*std(q1) + A12^(-1)*std(q2)
1052  m_DFC1[i][0][j] =
1053  m_gmat[3][j] * m_BD1[i][0][j] -
1054  m_gmat[2][j] * m_BD1[i][1][j];
1055 
1056  // q2 = A21^(-1)*std(q1) + A22^(-1)*std(q2)
1057  m_DFC1[i][1][j] =
1058  m_gmat[0][j] * m_BD1[i][1][j] -
1059  m_gmat[1][j] * m_BD1[i][0][j];
1060  }
1061 
1062  // Computing the physical first-order derivatives
1063  for (j = 0; j < nDim; ++j)
1064  {
1065  Vmath::Vadd(nSolutionPts, &m_DU1[i][j][0], 1,
1066  &m_DFC1[i][j][0], 1,
1067  &m_D1[i][j][0], 1);
1068  }
1069  }
1070 
1071  // Computing interface numerical fluxes for m_D1
1072  // in physical space
1073  v_NumFluxforVector(fields, inarray, m_D1, m_IF2);
1074 
1075  // Computing the standard second-order discontinuous
1076  // derivatives
1077  for (i = 0; i < nConvectiveFields; ++i)
1078  {
1079  Array<OneD, NekDouble> f_hat(nSolutionPts, 0.0);
1080  Array<OneD, NekDouble> g_hat(nSolutionPts, 0.0);
1081 
1082  for (j = 0; j < nSolutionPts; j++)
1083  {
1084  f_hat[j] = (m_D1[i][0][j] * m_gmat[0][j] +
1085  m_D1[i][1][j] * m_gmat[2][j])*m_jac[j];
1086 
1087  g_hat[j] = (m_D1[i][0][j] * m_gmat[1][j] +
1088  m_D1[i][1][j] * m_gmat[3][j])*m_jac[j];
1089  }
1090 
1091  for (n = 0; n < nElements; n++)
1092  {
1093  phys_offset = fields[0]->GetPhys_Offset(n);
1094 
1095  fields[0]->GetExp(n)->StdPhysDeriv(0,
1096  auxArray1 = f_hat + phys_offset,
1097  auxArray2 = m_DD1[i][0] + phys_offset);
1098 
1099  fields[0]->GetExp(n)->StdPhysDeriv(1,
1100  auxArray1 = g_hat + phys_offset,
1101  auxArray2 = m_DD1[i][1] + phys_offset);
1102  }
1103 
1104  // Divergence of the standard discontinuous flux
1105  Vmath::Vadd(nSolutionPts,
1106  &m_DD1[i][0][0], 1,
1107  &m_DD1[i][1][0], 1,
1108  &m_divFD[i][0], 1);
1109 
1110  // Divergence of the standard correction flux
1111  if (Basis[0]->GetPointsType() ==
1112  LibUtilities::eGaussGaussLegendre &&
1113  Basis[1]->GetPointsType() ==
1114  LibUtilities::eGaussGaussLegendre)
1115  {
1116  v_DivCFlux_2D_Gauss(nConvectiveFields, fields,
1117  f_hat, g_hat,
1118  m_IF2[i], m_divFC[i]);
1119  }
1120  else
1121  {
1122  v_DivCFlux_2D(nConvectiveFields, fields,
1123  m_D1[i][0], m_D1[i][1],
1124  m_IF2[i], m_divFC[i]);
1125  }
1126 
1127  // Divergence of the standard final flux
1128  Vmath::Vadd(nSolutionPts, &m_divFD[i][0], 1,
1129  &m_divFC[i][0], 1, &outarray[i][0], 1);
1130 
1131  // Dividing by the jacobian to get back into
1132  // physical space
1133  Vmath::Vdiv(nSolutionPts, &outarray[i][0], 1,
1134  &m_jac[0], 1, &outarray[i][0], 1);
1135 
1136  // Primitive Dealiasing 2D
1137  if (!(Basis[0]->Collocation()))
1138  {
1139  fields[i]->FwdTrans(outarray[i], outarrayCoeff[i]);
1140  fields[i]->BwdTrans(outarrayCoeff[i], outarray[i]);
1141  }
1142  }//Close loop on nConvectiveFields
1143  break;
1144  }
1145  // 3D problems
1146  case 3:
1147  {
1148  ASSERTL0(false, "3D FRDG case not implemented yet");
1149  break;
1150  }
1151  }
1152  }
Nektar::SolverUtils::DiffusionLFR::v_DivCFlux_2D
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.
Definition: DiffusionLFR.cpp:1793
ASSERTL0
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:161
Nektar::SolverUtils::DiffusionLFR::v_NumFluxforVector
virtual void v_NumFluxforVector(const Array< OneD, MultiRegions::ExpListSharedPtr > &fields, const Array< OneD, Array< OneD, NekDouble > > &ufield, Array< OneD, Array< OneD, Array< OneD, NekDouble > > > &qfield, Array< OneD, Array< OneD, NekDouble > > &qflux)
Build the numerical flux for the 2nd order derivatives.
Definition: DiffusionLFR.cpp:1299
Nektar::SolverUtils::DiffusionLFR::m_DU1
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > m_DU1
Definition: DiffusionLFR.h:79
Vmath::Vdiv
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
Nektar::LibUtilities::eGaussGaussLegendre
1D Gauss-Gauss-Legendre quadrature points
Definition: PointsType.h:47
Nektar::SolverUtils::DiffusionLFR::m_jac
Array< OneD, NekDouble > m_jac
Definition: DiffusionLFR.h:55
Nektar::SolverUtils::DiffusionLFR::m_divFC
Array< OneD, Array< OneD, NekDouble > > m_divFC
Definition: DiffusionLFR.h:87
Nektar::SolverUtils::DiffusionLFR::v_NumFluxforScalar
virtual void v_NumFluxforScalar(const Array< OneD, MultiRegions::ExpListSharedPtr > &fields, const Array< OneD, Array< OneD, NekDouble > > &ufield, Array< OneD, Array< OneD, Array< OneD, NekDouble > > > &uflux)
Builds the numerical flux for the 1st order derivatives.
Definition: DiffusionLFR.cpp:1158
Nektar::SolverUtils::DiffusionLFR::m_DFC1
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > m_DFC1
Definition: DiffusionLFR.h:80
Nektar::SolverUtils::DiffusionLFR::v_DerCFlux_1D
virtual void v_DerCFlux_1D(const int nConvectiveFields, const Array< OneD, MultiRegions::ExpListSharedPtr > &fields, const Array< OneD, const NekDouble > &flux, const Array< OneD, const NekDouble > &iFlux, Array< OneD, NekDouble > &derCFlux)
Compute the derivative of the corrective flux for 1D problems.
Definition: DiffusionLFR.cpp:1462
Nektar::SolverUtils::DiffusionLFR::v_DerCFlux_2D
virtual void v_DerCFlux_2D(const int nConvectiveFields, const int direction, const Array< OneD, MultiRegions::ExpListSharedPtr > &fields, const Array< OneD, const NekDouble > &flux, const Array< OneD, NekDouble > &iFlux, Array< OneD, NekDouble > &derCFlux)
Compute the derivative of the corrective flux wrt a given coordinate for 2D problems.
Definition: DiffusionLFR.cpp:1578
Nektar::SolverUtils::DiffusionLFR::m_gmat
Array< OneD, Array< OneD, NekDouble > > m_gmat
Definition: DiffusionLFR.h:56
Nektar::SolverUtils::DiffusionLFR::m_DFC2
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > m_DFC2
Definition: DiffusionLFR.h:85
Nektar::SolverUtils::DiffusionLFR::m_tmp1
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > m_tmp1
Definition: DiffusionLFR.h:89
Nektar::SolverUtils::DiffusionLFR::m_tmp2
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > m_tmp2
Definition: DiffusionLFR.h:90
Nektar::SolverUtils::DiffusionLFR::m_IF1
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > m_IF1
Definition: DiffusionLFR.h:78
Nektar::SolverUtils::DiffusionLFR::m_DD1
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > m_DD1
Definition: DiffusionLFR.h:83
Nektar::SolverUtils::DiffusionLFR::m_divFD
Array< OneD, Array< OneD, NekDouble > > m_divFD
Definition: DiffusionLFR.h:86
Nektar::SolverUtils::DiffusionLFR::m_D1
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > m_D1
Definition: DiffusionLFR.h:82
Vmath::Vcopy
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.cpp:1047
Nektar::SolverUtils::DiffusionLFR::m_BD1
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > m_BD1
Definition: DiffusionLFR.h:81
Vmath::Vadd
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
Vmath::Vmul
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
Nektar::SolverUtils::DiffusionLFR::v_DivCFlux_2D_Gauss
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".
Definition: DiffusionLFR.cpp:1989
Nektar::SolverUtils::DiffusionLFR::m_IF2
Array< OneD, Array< OneD, NekDouble > > m_IF2
Definition: DiffusionLFR.h:84
void Nektar::SolverUtils::DiffusionLFR::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 for 2D problems.
Todo:
: Switch on shapes eventually here.

Definition at line 1793 of file DiffusionLFR.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(), and Vmath::Vsub().

Referenced by v_Diffuse().

1800  {
1801  int n, e, i, j, cnt;
1802 
1803  int nElements = fields[0]->GetExpSize();
1804 
1805  int nLocalSolutionPts;
1806  int nEdgePts;
1807  int trace_offset;
1808  int phys_offset;
1809  int nquad0;
1810  int nquad1;
1811 
1812  Array<OneD, NekDouble> auxArray1, auxArray2;
1813  Array<OneD, LibUtilities::BasisSharedPtr> base;
1814 
1815  Array<OneD, Array<OneD, LocalRegions::ExpansionSharedPtr> >
1816  &elmtToTrace = fields[0]->GetTraceMap()->GetElmtToTrace();
1817 
1818  // Loop on the elements
1819  for(n = 0; n < nElements; ++n)
1820  {
1821  // Offset of the element on the global vector
1822  phys_offset = fields[0]->GetPhys_Offset(n);
1823  nLocalSolutionPts = fields[0]->GetExp(n)->GetTotPoints();
1824 
1825  base = fields[0]->GetExp(n)->GetBase();
1826  nquad0 = base[0]->GetNumPoints();
1827  nquad1 = base[1]->GetNumPoints();
1828 
1829  Array<OneD, NekDouble> divCFluxE0(nLocalSolutionPts, 0.0);
1830  Array<OneD, NekDouble> divCFluxE1(nLocalSolutionPts, 0.0);
1831  Array<OneD, NekDouble> divCFluxE2(nLocalSolutionPts, 0.0);
1832  Array<OneD, NekDouble> divCFluxE3(nLocalSolutionPts, 0.0);
1833 
1834  // Loop on the edges
1835  for(e = 0; e < fields[0]->GetExp(n)->GetNedges(); ++e)
1836  {
1837  // Number of edge points of edge e
1838  nEdgePts = fields[0]->GetExp(n)->GetEdgeNumPoints(e);
1839 
1840  Array<OneD, NekDouble> tmparrayX1(nEdgePts, 0.0);
1841  Array<OneD, NekDouble> tmparrayX2(nEdgePts, 0.0);
1842  Array<OneD, NekDouble> fluxN (nEdgePts, 0.0);
1843  Array<OneD, NekDouble> fluxT (nEdgePts, 0.0);
1844  Array<OneD, NekDouble> fluxJumps(nEdgePts, 0.0);
1845 
1846  // Offset of the trace space correspondent to edge e
1847  trace_offset = fields[0]->GetTrace()->GetPhys_Offset(
1848  elmtToTrace[n][e]->GetElmtId());
1849 
1850  // Get the normals of edge e
1851  const Array<OneD, const Array<OneD, NekDouble> > &normals =
1852  fields[0]->GetExp(n)->GetEdgeNormal(e);
1853 
1854  // Extract the edge values of flux-x on edge e and order
1855  // them accordingly to the order of the trace space
1856  fields[0]->GetExp(n)->GetEdgePhysVals(
1857  e, elmtToTrace[n][e],
1858  fluxX1 + phys_offset,
1859  auxArray1 = tmparrayX1);
1860 
1861  // Extract the edge values of flux-y on edge e and order
1862  // them accordingly to the order of the trace space
1863  fields[0]->GetExp(n)->GetEdgePhysVals(
1864  e, elmtToTrace[n][e],
1865  fluxX2 + phys_offset,
1866  auxArray1 = tmparrayX2);
1867 
1868  // Multiply the edge components of the flux by the normal
1869  for (i = 0; i < nEdgePts; ++i)
1870  {
1871  fluxN[i] =
1872  tmparrayX1[i]*m_traceNormals[0][trace_offset+i] +
1873  tmparrayX2[i]*m_traceNormals[1][trace_offset+i];
1874  }
1875 
1876  // Subtract to the Riemann flux the discontinuous flux
1877  Vmath::Vsub(nEdgePts,
1878  &numericalFlux[trace_offset], 1,
1879  &fluxN[0], 1, &fluxJumps[0], 1);
1880 
1881  // Check the ordering of the jump vectors
1882  if (fields[0]->GetExp(n)->GetEorient(e) ==
1883  StdRegions::eBackwards)
1884  {
1885  Vmath::Reverse(nEdgePts,
1886  auxArray1 = fluxJumps, 1,
1887  auxArray2 = fluxJumps, 1);
1888  }
1889 
1890  NekDouble fac = fields[0]->GetExp(n)->EdgeNormalNegated(e) ?
1891  -1.0 : 1.0;
1892 
1893  for (i = 0; i < nEdgePts; ++i)
1894  {
1895  if (m_traceNormals[0][trace_offset+i] != fac*normals[0][i]
1896  || m_traceNormals[1][trace_offset+i] != fac*normals[1][i])
1897  {
1898  fluxJumps[i] = -fluxJumps[i];
1899  }
1900  }
1901 
1902  // Multiply jumps by derivatives of the correction functions
1903  switch (e)
1904  {
1905  case 0:
1906  for (i = 0; i < nquad0; ++i)
1907  {
1908  // Multiply fluxJumps by Q factors
1909  fluxJumps[i] = -(m_Q2D_e0[n][i]) * fluxJumps[i];
1910 
1911  for (j = 0; j < nquad1; ++j)
1912  {
1913  cnt = i + j*nquad0;
1914  divCFluxE0[cnt] = fluxJumps[i] * m_dGL_xi2[n][j];
1915  }
1916  }
1917  break;
1918  case 1:
1919  for (i = 0; i < nquad1; ++i)
1920  {
1921  // Multiply fluxJumps by Q factors
1922  fluxJumps[i] = (m_Q2D_e1[n][i]) * fluxJumps[i];
1923 
1924  for (j = 0; j < nquad0; ++j)
1925  {
1926  cnt = (nquad0)*i + j;
1927  divCFluxE1[cnt] = fluxJumps[i] * m_dGR_xi1[n][j];
1928  }
1929  }
1930  break;
1931  case 2:
1932  for (i = 0; i < nquad0; ++i)
1933  {
1934  // Multiply fluxJumps by Q factors
1935  fluxJumps[i] = (m_Q2D_e2[n][i]) * fluxJumps[i];
1936 
1937  for (j = 0; j < nquad1; ++j)
1938  {
1939  cnt = (nquad0 - 1) + j*nquad0 - i;
1940  divCFluxE2[cnt] = fluxJumps[i] * m_dGR_xi2[n][j];
1941  }
1942  }
1943  break;
1944  case 3:
1945  for (i = 0; i < nquad1; ++i)
1946  {
1947  // Multiply fluxJumps by Q factors
1948  fluxJumps[i] = -(m_Q2D_e3[n][i]) * fluxJumps[i];
1949  for (j = 0; j < nquad0; ++j)
1950  {
1951  cnt = (nquad0*nquad1 - nquad0) + j - i*nquad0;
1952  divCFluxE3[cnt] = fluxJumps[i] * m_dGL_xi1[n][j];
1953  }
1954  }
1955  break;
1956 
1957  default:
1958  ASSERTL0(false,"edge value (< 3) is out of range");
1959  break;
1960  }
1961  }
1962 
1963  // Sum each edge contribution
1964  Vmath::Vadd(nLocalSolutionPts, &divCFluxE0[0], 1,
1965  &divCFluxE1[0], 1, &divCFlux[phys_offset], 1);
1966 
1967  Vmath::Vadd(nLocalSolutionPts, &divCFlux[phys_offset], 1,
1968  &divCFluxE2[0], 1, &divCFlux[phys_offset], 1);
1969 
1970  Vmath::Vadd(nLocalSolutionPts, &divCFlux[phys_offset], 1,
1971  &divCFluxE3[0], 1, &divCFlux[phys_offset], 1);
1972  }
1973  }
Nektar::SolverUtils::DiffusionLFR::m_dGL_xi1
Array< OneD, Array< OneD, NekDouble > > m_dGL_xi1
Definition: DiffusionLFR.h:63
ASSERTL0
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:161
Nektar::SolverUtils::DiffusionLFR::m_Q2D_e2
Array< OneD, Array< OneD, NekDouble > > m_Q2D_e2
Definition: DiffusionLFR.h:60
Nektar::SolverUtils::DiffusionLFR::m_Q2D_e3
Array< OneD, Array< OneD, NekDouble > > m_Q2D_e3
Definition: DiffusionLFR.h:61
Nektar::SolverUtils::DiffusionLFR::m_dGR_xi1
Array< OneD, Array< OneD, NekDouble > > m_dGR_xi1
Definition: DiffusionLFR.h:64
Nektar::SolverUtils::DiffusionLFR::m_dGL_xi2
Array< OneD, Array< OneD, NekDouble > > m_dGL_xi2
Definition: DiffusionLFR.h:65
Nektar::SolverUtils::DiffusionLFR::m_traceNormals
Array< OneD, Array< OneD, NekDouble > > m_traceNormals
Definition: DiffusionLFR.h:75
Nektar::SolverUtils::DiffusionLFR::m_dGR_xi2
Array< OneD, Array< OneD, NekDouble > > m_dGR_xi2
Definition: DiffusionLFR.h:66
Vmath::Reverse
void Reverse(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.cpp:1071
Nektar::SolverUtils::DiffusionLFR::m_Q2D_e1
Array< OneD, Array< OneD, NekDouble > > m_Q2D_e1
Definition: DiffusionLFR.h:59
Nektar::NekDouble
double NekDouble
Definition: NektarUnivTypeDefs.hpp:44
Nektar::SolverUtils::DiffusionLFR::m_Q2D_e0
Array< OneD, Array< OneD, NekDouble > > m_Q2D_e0
Definition: DiffusionLFR.h:58
Vmath::Vsub
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
Vmath::Vadd
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::DiffusionLFR::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 1989 of file DiffusionLFR.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_Diffuse().

1996  {
1997  int n, e, i, j, cnt;
1998 
1999  int nElements = fields[0]->GetExpSize();
2000  int nLocalSolutionPts;
2001  int nEdgePts;
2002  int trace_offset;
2003  int phys_offset;
2004  int nquad0;
2005  int nquad1;
2006 
2007  Array<OneD, NekDouble> auxArray1, auxArray2;
2008  Array<OneD, LibUtilities::BasisSharedPtr> base;
2009 
2010  Array<OneD, Array<OneD, LocalRegions::ExpansionSharedPtr> >
2011  &elmtToTrace = fields[0]->GetTraceMap()->GetElmtToTrace();
2012 
2013  // Loop on the elements
2014  for(n = 0; n < nElements; ++n)
2015  {
2016  // Offset of the element on the global vector
2017  phys_offset = fields[0]->GetPhys_Offset(n);
2018  nLocalSolutionPts = fields[0]->GetExp(n)->GetTotPoints();
2019 
2020  base = fields[0]->GetExp(n)->GetBase();
2021  nquad0 = base[0]->GetNumPoints();
2022  nquad1 = base[1]->GetNumPoints();
2023 
2024  Array<OneD, NekDouble> divCFluxE0(nLocalSolutionPts, 0.0);
2025  Array<OneD, NekDouble> divCFluxE1(nLocalSolutionPts, 0.0);
2026  Array<OneD, NekDouble> divCFluxE2(nLocalSolutionPts, 0.0);
2027  Array<OneD, NekDouble> divCFluxE3(nLocalSolutionPts, 0.0);
2028 
2029  // Loop on the edges
2030  for(e = 0; e < fields[0]->GetExp(n)->GetNedges(); ++e)
2031  {
2032  // Number of edge points of edge e
2033  nEdgePts = fields[0]->GetExp(n)->GetEdgeNumPoints(e);
2034 
2035  Array<OneD, NekDouble> fluxN (nEdgePts, 0.0);
2036  Array<OneD, NekDouble> fluxT (nEdgePts, 0.0);
2037  Array<OneD, NekDouble> fluxN_R (nEdgePts, 0.0);
2038  Array<OneD, NekDouble> fluxN_D (nEdgePts, 0.0);
2039  Array<OneD, NekDouble> fluxJumps (nEdgePts, 0.0);
2040 
2041  // Offset of the trace space correspondent to edge e
2042  trace_offset = fields[0]->GetTrace()->GetPhys_Offset(
2043  elmtToTrace[n][e]->GetElmtId());
2044 
2045  // Get the normals of edge e
2046  const Array<OneD, const Array<OneD, NekDouble> > &normals =
2047  fields[0]->GetExp(n)->GetEdgeNormal(e);
2048 
2049  // Extract the trasformed normal flux at each edge
2050  switch (e)
2051  {
2052  case 0:
2053  // Extract the edge values of transformed flux-y on
2054  // edge e and order them accordingly to the order of
2055  // the trace space
2056  fields[0]->GetExp(n)->GetEdgePhysVals(
2057  e, elmtToTrace[n][e],
2058  fluxX2 + phys_offset,
2059  auxArray1 = fluxN_D);
2060 
2061  Vmath::Neg (nEdgePts, fluxN_D, 1);
2062 
2063  // Extract the physical Rieamann flux at each edge
2064  Vmath::Vcopy(nEdgePts,
2065  &numericalFlux[trace_offset], 1,
2066  &fluxN[0], 1);
2067 
2068  // Check the ordering of vectors
2069  if (fields[0]->GetExp(n)->GetEorient(e) ==
2070  StdRegions::eBackwards)
2071  {
2072  Vmath::Reverse(nEdgePts,
2073  auxArray1 = fluxN, 1,
2074  auxArray2 = fluxN, 1);
2075 
2076  Vmath::Reverse(nEdgePts,
2077  auxArray1 = fluxN_D, 1,
2078  auxArray2 = fluxN_D, 1);
2079  }
2080 
2081  // Transform Riemann Fluxes in the standard element
2082  for (i = 0; i < nquad0; ++i)
2083  {
2084  // Multiply Riemann Flux by Q factors
2085  fluxN_R[i] = (m_Q2D_e0[n][i]) * fluxN[i];
2086  }
2087 
2088  for (i = 0; i < nEdgePts; ++i)
2089  {
2090  if (m_traceNormals[0][trace_offset+i] !=
2091  normals[0][i] ||
2092  m_traceNormals[1][trace_offset+i] !=
2093  normals[1][i])
2094  {
2095  fluxN_R[i] = -fluxN_R[i];
2096  }
2097  }
2098 
2099  // Subtract to the Riemann flux the discontinuous
2100  // flux
2101  Vmath::Vsub(nEdgePts,
2102  &fluxN_R[0], 1,
2103  &fluxN_D[0], 1, &fluxJumps[0], 1);
2104 
2105  // Multiplicate the flux jumps for the correction
2106  // function
2107  for (i = 0; i < nquad0; ++i)
2108  {
2109  for (j = 0; j < nquad1; ++j)
2110  {
2111  cnt = i + j*nquad0;
2112  divCFluxE0[cnt] =
2113  -fluxJumps[i] * m_dGL_xi2[n][j];
2114  }
2115  }
2116  break;
2117  case 1:
2118  // Extract the edge values of transformed flux-x on
2119  // edge e and order them accordingly to the order of
2120  // the trace space
2121  fields[0]->GetExp(n)->GetEdgePhysVals(
2122  e, elmtToTrace[n][e],
2123  fluxX1 + phys_offset,
2124  auxArray1 = fluxN_D);
2125 
2126  // Extract the physical Rieamann flux at each edge
2127  Vmath::Vcopy(nEdgePts,
2128  &numericalFlux[trace_offset], 1,
2129  &fluxN[0], 1);
2130 
2131  // Check the ordering of vectors
2132  if (fields[0]->GetExp(n)->GetEorient(e) ==
2133  StdRegions::eBackwards)
2134  {
2135  Vmath::Reverse(nEdgePts,
2136  auxArray1 = fluxN, 1,
2137  auxArray2 = fluxN, 1);
2138 
2139  Vmath::Reverse(nEdgePts,
2140  auxArray1 = fluxN_D, 1,
2141  auxArray2 = fluxN_D, 1);
2142  }
2143 
2144  // Transform Riemann Fluxes in the standard element
2145  for (i = 0; i < nquad1; ++i)
2146  {
2147  // Multiply Riemann Flux by Q factors
2148  fluxN_R[i] = (m_Q2D_e1[n][i]) * fluxN[i];
2149  }
2150 
2151  for (i = 0; i < nEdgePts; ++i)
2152  {
2153  if (m_traceNormals[0][trace_offset+i] !=
2154  normals[0][i] ||
2155  m_traceNormals[1][trace_offset+i] !=
2156  normals[1][i])
2157  {
2158  fluxN_R[i] = -fluxN_R[i];
2159  }
2160  }
2161 
2162  // Subtract to the Riemann flux the discontinuous
2163  // flux
2164  Vmath::Vsub(nEdgePts,
2165  &fluxN_R[0], 1,
2166  &fluxN_D[0], 1, &fluxJumps[0], 1);
2167 
2168  // Multiplicate the flux jumps for the correction
2169  // function
2170  for (i = 0; i < nquad1; ++i)
2171  {
2172  for (j = 0; j < nquad0; ++j)
2173  {
2174  cnt = (nquad0)*i + j;
2175  divCFluxE1[cnt] =
2176  fluxJumps[i] * m_dGR_xi1[n][j];
2177  }
2178  }
2179  break;
2180  case 2:
2181 
2182  // Extract the edge values of transformed flux-y on
2183  // edge e and order them accordingly to the order of
2184  // the trace space
2185 
2186  fields[0]->GetExp(n)->GetEdgePhysVals(
2187  e, elmtToTrace[n][e],
2188  fluxX2 + phys_offset,
2189  auxArray1 = fluxN_D);
2190 
2191  // Extract the physical Rieamann flux at each edge
2192  Vmath::Vcopy(nEdgePts,
2193  &numericalFlux[trace_offset], 1,
2194  &fluxN[0], 1);
2195 
2196  // Check the ordering of vectors
2197  if (fields[0]->GetExp(n)->GetEorient(e) ==
2198  StdRegions::eBackwards)
2199  {
2200  Vmath::Reverse(nEdgePts,
2201  auxArray1 = fluxN, 1,
2202  auxArray2 = fluxN, 1);
2203 
2204  Vmath::Reverse(nEdgePts,
2205  auxArray1 = fluxN_D, 1,
2206  auxArray2 = fluxN_D, 1);
2207  }
2208 
2209  // Transform Riemann Fluxes in the standard element
2210  for (i = 0; i < nquad0; ++i)
2211  {
2212  // Multiply Riemann Flux by Q factors
2213  fluxN_R[i] = (m_Q2D_e2[n][i]) * fluxN[i];
2214  }
2215 
2216  for (i = 0; i < nEdgePts; ++i)
2217  {
2218  if (m_traceNormals[0][trace_offset+i] !=
2219  normals[0][i] ||
2220  m_traceNormals[1][trace_offset+i] !=
2221  normals[1][i])
2222  {
2223  fluxN_R[i] = -fluxN_R[i];
2224  }
2225  }
2226 
2227  // Subtract to the Riemann flux the discontinuous
2228  // flux
2229 
2230  Vmath::Vsub(nEdgePts,
2231  &fluxN_R[0], 1,
2232  &fluxN_D[0], 1, &fluxJumps[0], 1);
2233 
2234  // Multiplicate the flux jumps for the correction
2235  // function
2236  for (i = 0; i < nquad0; ++i)
2237  {
2238  for (j = 0; j < nquad1; ++j)
2239  {
2240  cnt = (nquad0 - 1) + j*nquad0 - i;
2241  divCFluxE2[cnt] =
2242  fluxJumps[i] * m_dGR_xi2[n][j];
2243  }
2244  }
2245  break;
2246  case 3:
2247  // Extract the edge values of transformed flux-x on
2248  // edge e and order them accordingly to the order of
2249  // the trace space
2250 
2251  fields[0]->GetExp(n)->GetEdgePhysVals(
2252  e, elmtToTrace[n][e],
2253  fluxX1 + phys_offset,
2254  auxArray1 = fluxN_D);
2255  Vmath::Neg (nEdgePts, fluxN_D, 1);
2256 
2257  // Extract the physical Rieamann flux at each edge
2258  Vmath::Vcopy(nEdgePts,
2259  &numericalFlux[trace_offset], 1,
2260  &fluxN[0], 1);
2261 
2262  // Check the ordering of vectors
2263  if (fields[0]->GetExp(n)->GetEorient(e) ==
2264  StdRegions::eBackwards)
2265  {
2266  Vmath::Reverse(nEdgePts,
2267  auxArray1 = fluxN, 1,
2268  auxArray2 = fluxN, 1);
2269 
2270  Vmath::Reverse(nEdgePts,
2271  auxArray1 = fluxN_D, 1,
2272  auxArray2 = fluxN_D, 1);
2273  }
2274 
2275  // Transform Riemann Fluxes in the standard element
2276  for (i = 0; i < nquad1; ++i)
2277  {
2278  // Multiply Riemann Flux by Q factors
2279  fluxN_R[i] = (m_Q2D_e3[n][i]) * fluxN[i];
2280  }
2281 
2282  for (i = 0; i < nEdgePts; ++i)
2283  {
2284  if (m_traceNormals[0][trace_offset+i] !=
2285  normals[0][i] ||
2286  m_traceNormals[1][trace_offset+i] !=
2287  normals[1][i])
2288  {
2289  fluxN_R[i] = -fluxN_R[i];
2290  }
2291  }
2292 
2293  // Subtract to the Riemann flux the discontinuous
2294  // flux
2295 
2296  Vmath::Vsub(nEdgePts,
2297  &fluxN_R[0], 1,
2298  &fluxN_D[0], 1, &fluxJumps[0], 1);
2299 
2300  // Multiplicate the flux jumps for the correction
2301  // function
2302  for (i = 0; i < nquad1; ++i)
2303  {
2304  for (j = 0; j < nquad0; ++j)
2305  {
2306  cnt = (nquad0*nquad1 - nquad0) + j
2307  - i*nquad0;
2308  divCFluxE3[cnt] =
2309  -fluxJumps[i] * m_dGL_xi1[n][j];
2310  }
2311  }
2312  break;
2313  default:
2314  ASSERTL0(false,"edge value (< 3) is out of range");
2315  break;
2316  }
2317  }
2318 
2319 
2320  // Sum each edge contribution
2321  Vmath::Vadd(nLocalSolutionPts, &divCFluxE0[0], 1,
2322  &divCFluxE1[0], 1, &divCFlux[phys_offset], 1);
2323 
2324  Vmath::Vadd(nLocalSolutionPts, &divCFlux[phys_offset], 1,
2325  &divCFluxE2[0], 1, &divCFlux[phys_offset], 1);
2326 
2327  Vmath::Vadd(nLocalSolutionPts, &divCFlux[phys_offset], 1,
2328  &divCFluxE3[0], 1, &divCFlux[phys_offset], 1);
2329  }
2330  }
Nektar::SolverUtils::DiffusionLFR::m_dGL_xi1
Array< OneD, Array< OneD, NekDouble > > m_dGL_xi1
Definition: DiffusionLFR.h:63
ASSERTL0
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:161
Nektar::SolverUtils::DiffusionLFR::m_Q2D_e2
Array< OneD, Array< OneD, NekDouble > > m_Q2D_e2
Definition: DiffusionLFR.h:60
Nektar::SolverUtils::DiffusionLFR::m_Q2D_e3
Array< OneD, Array< OneD, NekDouble > > m_Q2D_e3
Definition: DiffusionLFR.h:61
Nektar::SolverUtils::DiffusionLFR::m_dGR_xi1
Array< OneD, Array< OneD, NekDouble > > m_dGR_xi1
Definition: DiffusionLFR.h:64
Nektar::SolverUtils::DiffusionLFR::m_dGL_xi2
Array< OneD, Array< OneD, NekDouble > > m_dGL_xi2
Definition: DiffusionLFR.h:65
Nektar::SolverUtils::DiffusionLFR::m_traceNormals
Array< OneD, Array< OneD, NekDouble > > m_traceNormals
Definition: DiffusionLFR.h:75
Nektar::SolverUtils::DiffusionLFR::m_dGR_xi2
Array< OneD, Array< OneD, NekDouble > > m_dGR_xi2
Definition: DiffusionLFR.h:66
Vmath::Reverse
void Reverse(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.cpp:1071
Nektar::SolverUtils::DiffusionLFR::m_Q2D_e1
Array< OneD, Array< OneD, NekDouble > > m_Q2D_e1
Definition: DiffusionLFR.h:59
Vmath::Neg
void Neg(int n, T *x, const int incx)
Negate x = -x.
Definition: Vmath.cpp:382
Nektar::SolverUtils::DiffusionLFR::m_Q2D_e0
Array< OneD, Array< OneD, NekDouble > > m_Q2D_e0
Definition: DiffusionLFR.h:58
Vmath::Vsub
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
Vmath::Vcopy
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.cpp:1047
Vmath::Vadd
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::DiffusionLFR::v_InitObject ( LibUtilities::SessionReaderSharedPtr  pSession,
Array< OneD, MultiRegions::ExpListSharedPtr >  pFields 
)
protectedvirtual

Initiliase DiffusionLFR 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 DiffusionLFR.

Parameters
pSessionPointer to session reader.
pFieldsPointer to fields.

Reimplemented from Nektar::SolverUtils::Diffusion.

Definition at line 85 of file DiffusionLFR.cpp.

References m_BD1, m_D1, m_DD1, m_DFC1, m_DFC2, m_divFC, m_divFD, m_DU1, m_IF1, m_IF2, m_session, m_tmp1, m_tmp2, v_SetupCFunctions(), and v_SetupMetrics().

88  {
89  m_session = pSession;
90  v_SetupMetrics(pSession, pFields);
91  v_SetupCFunctions(pSession, pFields);
92 
93  // Initialising arrays
94  int i, j;
95  int nConvectiveFields = pFields.num_elements();
96  int nDim = pFields[0]->GetCoordim(0);
97  int nSolutionPts = pFields[0]->GetTotPoints();
98  int nTracePts = pFields[0]->GetTrace()->GetTotPoints();
99 
100  m_IF1 = Array<OneD, Array<OneD, Array<OneD, NekDouble> > > (
101  nConvectiveFields);
102  m_DU1 = Array<OneD, Array<OneD, Array<OneD, NekDouble> > > (
103  nConvectiveFields);
104  m_DFC1 = Array<OneD, Array<OneD, Array<OneD, NekDouble> > > (
105  nConvectiveFields);
106  m_BD1 = Array<OneD, Array<OneD, Array<OneD, NekDouble> > > (
107  nConvectiveFields);
108  m_D1 = Array<OneD, Array<OneD, Array<OneD, NekDouble> > > (
109  nConvectiveFields);
110  m_DD1 = Array<OneD, Array<OneD, Array<OneD, NekDouble> > > (
111  nConvectiveFields);
112  m_IF2 = Array<OneD, Array<OneD, NekDouble> > (nConvectiveFields);
113  m_DFC2 = Array<OneD, Array<OneD, Array<OneD, NekDouble> > > (
114  nConvectiveFields);
115  m_divFD = Array<OneD, Array<OneD, NekDouble> > (nConvectiveFields);
116  m_divFC = Array<OneD, Array<OneD, NekDouble> > (nConvectiveFields);
117 
118  m_tmp1 = Array<OneD, Array<OneD, Array<OneD, NekDouble> > > (
119  nConvectiveFields);
120  m_tmp2 = Array<OneD, Array<OneD, Array<OneD, NekDouble> > > (
121  nConvectiveFields);
122 
123 
124 
125  for (i = 0; i < nConvectiveFields; ++i)
126  {
127  m_IF1[i] = Array<OneD, Array<OneD, NekDouble> >(nDim);
128  m_DU1[i] = Array<OneD, Array<OneD, NekDouble> >(nDim);
129  m_DFC1[i] = Array<OneD, Array<OneD, NekDouble> >(nDim);
130  m_BD1[i] = Array<OneD, Array<OneD, NekDouble> >(nDim);
131  m_D1[i] = Array<OneD, Array<OneD, NekDouble> >(nDim);
132  m_DD1[i] = Array<OneD, Array<OneD, NekDouble> >(nDim);
133  m_IF2[i] = Array<OneD, NekDouble>(nTracePts, 0.0);
134  m_DFC2[i] = Array<OneD, Array<OneD, NekDouble> >(nDim);
135  m_divFD[i] = Array<OneD, NekDouble>(nSolutionPts, 0.0);
136  m_divFC[i] = Array<OneD, NekDouble>(nSolutionPts, 0.0);
137 
138  m_tmp1[i] = Array<OneD, Array<OneD, NekDouble> >(nDim);
139  m_tmp2[i] = Array<OneD, Array<OneD, NekDouble> >(nDim);
140 
141  for (j = 0; j < nDim; ++j)
142  {
143  m_IF1[i][j] = Array<OneD, NekDouble>(nTracePts, 0.0);
144  m_DU1[i][j] = Array<OneD, NekDouble>(nSolutionPts, 0.0);
145  m_DFC1[i][j] = Array<OneD, NekDouble>(nSolutionPts, 0.0);
146  m_BD1[i][j] = Array<OneD, NekDouble>(nSolutionPts, 0.0);
147  m_D1[i][j] = Array<OneD, NekDouble>(nSolutionPts, 0.0);
148  m_DD1[i][j] = Array<OneD, NekDouble>(nSolutionPts, 0.0);
149  m_DFC2[i][j] = Array<OneD, NekDouble>(nSolutionPts, 0.0);
150 
151  m_tmp1[i][j] = Array<OneD, NekDouble>(nSolutionPts, 0.0);
152  m_tmp2[i][j] = Array<OneD, NekDouble>(nSolutionPts, 0.0);
153  }
154  }
155 
156  }
Nektar::SolverUtils::DiffusionLFR::m_DU1
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > m_DU1
Definition: DiffusionLFR.h:79
Nektar::SolverUtils::DiffusionLFR::m_divFC
Array< OneD, Array< OneD, NekDouble > > m_divFC
Definition: DiffusionLFR.h:87
Nektar::SolverUtils::DiffusionLFR::v_SetupCFunctions
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...
Definition: DiffusionLFR.cpp:324
Nektar::SolverUtils::DiffusionLFR::m_DFC1
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > m_DFC1
Definition: DiffusionLFR.h:80
Nektar::SolverUtils::DiffusionLFR::v_SetupMetrics
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...
Definition: DiffusionLFR.cpp:174
Nektar::SolverUtils::DiffusionLFR::m_DFC2
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > m_DFC2
Definition: DiffusionLFR.h:85
Nektar::SolverUtils::DiffusionLFR::m_tmp1
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > m_tmp1
Definition: DiffusionLFR.h:89
Nektar::SolverUtils::DiffusionLFR::m_tmp2
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > m_tmp2
Definition: DiffusionLFR.h:90
Nektar::SolverUtils::DiffusionLFR::m_IF1
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > m_IF1
Definition: DiffusionLFR.h:78
Nektar::SolverUtils::DiffusionLFR::m_DD1
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > m_DD1
Definition: DiffusionLFR.h:83
Nektar::SolverUtils::DiffusionLFR::m_divFD
Array< OneD, Array< OneD, NekDouble > > m_divFD
Definition: DiffusionLFR.h:86
Nektar::SolverUtils::DiffusionLFR::m_D1
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > m_D1
Definition: DiffusionLFR.h:82
Nektar::SolverUtils::DiffusionLFR::m_BD1
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > m_BD1
Definition: DiffusionLFR.h:81
Nektar::SolverUtils::DiffusionLFR::m_session
LibUtilities::SessionReaderSharedPtr m_session
Definition: DiffusionLFR.h:76
Nektar::SolverUtils::DiffusionLFR::m_IF2
Array< OneD, Array< OneD, NekDouble > > m_IF2
Definition: DiffusionLFR.h:84
void Nektar::SolverUtils::DiffusionLFR::v_NumFluxforScalar ( const Array< OneD, MultiRegions::ExpListSharedPtr > &  fields,
const Array< OneD, Array< OneD, NekDouble > > &  ufield,
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > &  uflux 
)
protectedvirtual

Builds the numerical flux for the 1st order derivatives.

Definition at line 1158 of file DiffusionLFR.cpp.

References m_traceNormals, Vmath::Svtvp(), v_WeakPenaltyforScalar(), and Vmath::Vcopy().

Referenced by v_Diffuse().

1162  {
1163  int i, j;
1164  int nTracePts = fields[0]->GetTrace()->GetTotPoints();
1165  int nvariables = fields.num_elements();
1166  int nDim = fields[0]->GetCoordim(0);
1167 
1168  Array<OneD, NekDouble > Fwd (nTracePts);
1169  Array<OneD, NekDouble > Bwd (nTracePts);
1170  Array<OneD, NekDouble > Vn (nTracePts, 0.0);
1171  Array<OneD, NekDouble > fluxtemp(nTracePts, 0.0);
1172 
1173  // Get the normal velocity Vn
1174  for (i = 0; i < nDim; ++i)
1175  {
1176  Vmath::Svtvp(nTracePts, 1.0, m_traceNormals[i], 1,
1177  Vn, 1, Vn, 1);
1178  }
1179 
1180  // Get the sign of (v \cdot n), v = an arbitrary vector
1181  // Evaluate upwind flux:
1182  // uflux = \hat{u} \phi \cdot u = u^{(+,-)} n
1183  for (j = 0; j < nDim; ++j)
1184  {
1185  for (i = 0; i < nvariables ; ++i)
1186  {
1187  // Compute Fwd and Bwd value of ufield of i direction
1188  fields[i]->GetFwdBwdTracePhys(ufield[i], Fwd, Bwd);
1189 
1190  // if Vn >= 0, flux = uFwd, i.e.,
1191  // edge::eForward, if V*n>=0 <=> V*n_F>=0, pick uflux = uFwd
1192  // edge::eBackward, if V*n>=0 <=> V*n_B<0, pick uflux = uFwd
1193 
1194  // else if Vn < 0, flux = uBwd, i.e.,
1195  // edge::eForward, if V*n<0 <=> V*n_F<0, pick uflux = uBwd
1196  // edge::eBackward, if V*n<0 <=> V*n_B>=0, pick uflux = uBwd
1197 
1198  fields[i]->GetTrace()->Upwind(Vn, Fwd, Bwd, fluxtemp);
1199 
1200  // Imposing weak boundary condition with flux
1201  // if Vn >= 0, uflux = uBwd at Neumann, i.e.,
1202  // edge::eForward, if V*n>=0 <=> V*n_F>=0, pick uflux = uBwd
1203  // edge::eBackward, if V*n>=0 <=> V*n_B<0, pick uflux = uBwd
1204 
1205  // if Vn >= 0, uflux = uFwd at Neumann, i.e.,
1206  // edge::eForward, if V*n<0 <=> V*n_F<0, pick uflux = uFwd
1207  // edge::eBackward, if V*n<0 <=> V*n_B>=0, pick uflux = uFwd
1208 
1209  if(fields[0]->GetBndCondExpansions().num_elements())
1210  {
1211  v_WeakPenaltyforScalar(fields, i, ufield[i], fluxtemp);
1212  }
1213 
1214  // if Vn >= 0, flux = uFwd*(tan_{\xi}^- \cdot \vec{n}),
1215  // i.e,
1216  // edge::eForward, uFwd \(\tan_{\xi}^Fwd \cdot \vec{n})
1217  // edge::eBackward, uFwd \(\tan_{\xi}^Bwd \cdot \vec{n})
1218 
1219  // else if Vn < 0, flux = uBwd*(tan_{\xi}^- \cdot \vec{n}),
1220  // i.e,
1221  // edge::eForward, uBwd \(\tan_{\xi}^Fwd \cdot \vec{n})
1222  // edge::eBackward, uBwd \(\tan_{\xi}^Bwd \cdot \vec{n})
1223 
1224  Vmath::Vcopy(nTracePts, fluxtemp, 1, uflux[i][j], 1);
1225  }
1226  }
1227  }
Vmath::Svtvp
void Svtvp(int n, const T alpha, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
svtvp (scalar times vector plus vector): z = alpha*x + y
Definition: Vmath.cpp:471
Nektar::SolverUtils::DiffusionLFR::m_traceNormals
Array< OneD, Array< OneD, NekDouble > > m_traceNormals
Definition: DiffusionLFR.h:75
Nektar::SolverUtils::DiffusionLFR::v_WeakPenaltyforScalar
virtual void v_WeakPenaltyforScalar(const Array< OneD, MultiRegions::ExpListSharedPtr > &fields, const int var, const Array< OneD, const NekDouble > &ufield, Array< OneD, NekDouble > &penaltyflux)
Imposes appropriate bcs for the 1st order derivatives.
Definition: DiffusionLFR.cpp:1233
Vmath::Vcopy
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.cpp:1047
void Nektar::SolverUtils::DiffusionLFR::v_NumFluxforVector ( const Array< OneD, MultiRegions::ExpListSharedPtr > &  fields,
const Array< OneD, Array< OneD, NekDouble > > &  ufield,
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > &  qfield,
Array< OneD, Array< OneD, NekDouble > > &  qflux 
)
protectedvirtual

Build the numerical flux for the 2nd order derivatives.

Definition at line 1299 of file DiffusionLFR.cpp.

References m_traceNormals, Vmath::Svtvp(), v_WeakPenaltyforVector(), Vmath::Vadd(), and Vmath::Vmul().

Referenced by v_Diffuse().

1304  {
1305  int i, j;
1306  int nTracePts = fields[0]->GetTrace()->GetTotPoints();
1307  int nvariables = fields.num_elements();
1308  int nDim = fields[0]->GetCoordim(0);
1309 
1310  NekDouble C11 = 0.0;
1311  Array<OneD, NekDouble > Fwd(nTracePts);
1312  Array<OneD, NekDouble > Bwd(nTracePts);
1313  Array<OneD, NekDouble > Vn (nTracePts, 0.0);
1314 
1315  Array<OneD, NekDouble > qFwd (nTracePts);
1316  Array<OneD, NekDouble > qBwd (nTracePts);
1317  Array<OneD, NekDouble > qfluxtemp(nTracePts, 0.0);
1318  Array<OneD, NekDouble > uterm(nTracePts);
1319 
1320  // Get the normal velocity Vn
1321  for(i = 0; i < nDim; ++i)
1322  {
1323  Vmath::Svtvp(nTracePts, 1.0, m_traceNormals[i], 1,
1324  Vn, 1, Vn, 1);
1325  }
1326 
1327  // Evaulate upwind flux:
1328  // qflux = \hat{q} \cdot u = q \cdot n - C_(11)*(u^+ - u^-)
1329  for (i = 0; i < nvariables; ++i)
1330  {
1331  qflux[i] = Array<OneD, NekDouble> (nTracePts, 0.0);
1332  for (j = 0; j < nDim; ++j)
1333  {
1334  // Compute Fwd and Bwd value of ufield of jth direction
1335  fields[i]->GetFwdBwdTracePhys(qfield[i][j], qFwd, qBwd);
1336 
1337  // if Vn >= 0, flux = uFwd, i.e.,
1338  // edge::eForward, if V*n>=0 <=> V*n_F>=0, pick
1339  // qflux = qBwd = q+
1340  // edge::eBackward, if V*n>=0 <=> V*n_B<0, pick
1341  // qflux = qBwd = q-
1342 
1343  // else if Vn < 0, flux = uBwd, i.e.,
1344  // edge::eForward, if V*n<0 <=> V*n_F<0, pick
1345  // qflux = qFwd = q-
1346  // edge::eBackward, if V*n<0 <=> V*n_B>=0, pick
1347  // qflux = qFwd = q+
1348 
1349  fields[i]->GetTrace()->Upwind(Vn, qBwd, qFwd, qfluxtemp);
1350 
1351  Vmath::Vmul(nTracePts, m_traceNormals[j], 1, qfluxtemp, 1,
1352  qfluxtemp, 1);
1353 
1354  /*
1355  // Generate Stability term = - C11 ( u- - u+ )
1356  fields[i]->GetFwdBwdTracePhys(ufield[i], Fwd, Bwd);
1357 
1358  Vmath::Vsub(nTracePts,
1359  Fwd, 1, Bwd, 1,
1360  uterm, 1);
1361 
1362  Vmath::Smul(nTracePts,
1363  -1.0 * C11, uterm, 1,
1364  uterm, 1);
1365 
1366  // Flux = {Fwd, Bwd} * (nx, ny, nz) + uterm * (nx, ny)
1367  Vmath::Vadd(nTracePts,
1368  uterm, 1,
1369  qfluxtemp, 1,
1370  qfluxtemp, 1);
1371  */
1372 
1373  // Imposing weak boundary condition with flux
1374  if (fields[0]->GetBndCondExpansions().num_elements())
1375  {
1376  v_WeakPenaltyforVector(fields, i, j, qfield[i][j],
1377  qfluxtemp, C11);
1378  }
1379 
1380  // q_hat \cdot n = (q_xi \cdot n_xi) or (q_eta \cdot n_eta)
1381  // n_xi = n_x * tan_xi_x + n_y * tan_xi_y + n_z * tan_xi_z
1382  // n_xi = n_x * tan_eta_x + n_y * tan_eta_y + n_z*tan_eta_z
1383  Vmath::Vadd(nTracePts, qfluxtemp, 1, qflux[i], 1,
1384  qflux[i], 1);
1385  }
1386  }
1387  }
Vmath::Svtvp
void Svtvp(int n, const T alpha, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
svtvp (scalar times vector plus vector): z = alpha*x + y
Definition: Vmath.cpp:471
Nektar::SolverUtils::DiffusionLFR::m_traceNormals
Array< OneD, Array< OneD, NekDouble > > m_traceNormals
Definition: DiffusionLFR.h:75
Nektar::NekDouble
double NekDouble
Definition: NektarUnivTypeDefs.hpp:44
Nektar::SolverUtils::DiffusionLFR::v_WeakPenaltyforVector
virtual void v_WeakPenaltyforVector(const Array< OneD, MultiRegions::ExpListSharedPtr > &fields, const int var, const int dir, const Array< OneD, const NekDouble > &qfield, Array< OneD, NekDouble > &penaltyflux, NekDouble C11)
Imposes appropriate bcs for the 2nd order derivatives.
Definition: DiffusionLFR.cpp:1395
Vmath::Vadd
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
Vmath::Vmul
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::DiffusionLFR::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 324 of file DiffusionLFR.cpp.

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

Referenced by v_InitObject().

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

177  {
178  int i, n;
179  int nquad0, nquad1;
180  int phys_offset;
181  int nLocalSolutionPts;
182  int nElements = pFields[0]->GetExpSize();
183  int nDimensions = pFields[0]->GetCoordim(0);
184  int nSolutionPts = pFields[0]->GetTotPoints();
185  int nTracePts = pFields[0]->GetTrace()->GetTotPoints();
186 
187  m_traceNormals = Array<OneD, Array<OneD, NekDouble> >(nDimensions);
188  for (i = 0; i < nDimensions; ++i)
189  {
190  m_traceNormals[i] = Array<OneD, NekDouble> (nTracePts);
191  }
192  pFields[0]->GetTrace()->GetNormals(m_traceNormals);
193 
194  Array<TwoD, const NekDouble> gmat;
195  Array<OneD, const NekDouble> jac;
196 
197  m_jac = Array<OneD, NekDouble>(nSolutionPts);
198 
199  Array<OneD, NekDouble> auxArray1;
200  Array<OneD, LibUtilities::BasisSharedPtr> base;
201  LibUtilities::PointsKeyVector ptsKeys;
202 
203  switch (nDimensions)
204  {
205  case 1:
206  {
207  for (n = 0; n < nElements; ++n)
208  {
209  ptsKeys = pFields[0]->GetExp(n)->GetPointsKeys();
210  nLocalSolutionPts = pFields[0]->GetExp(n)->GetTotPoints();
211  phys_offset = pFields[0]->GetPhys_Offset(n);
212  jac = pFields[0]->GetExp(n)
213  ->as<LocalRegions::Expansion1D>()->GetGeom1D()
214  ->GetMetricInfo()->GetJac(ptsKeys);
215  for (i = 0; i < nLocalSolutionPts; ++i)
216  {
217  m_jac[i+phys_offset] = jac[0];
218  }
219  }
220  break;
221  }
222  case 2:
223  {
224  m_gmat = Array<OneD, Array<OneD, NekDouble> >(4);
225  m_gmat[0] = Array<OneD, NekDouble>(nSolutionPts);
226  m_gmat[1] = Array<OneD, NekDouble>(nSolutionPts);
227  m_gmat[2] = Array<OneD, NekDouble>(nSolutionPts);
228  m_gmat[3] = Array<OneD, NekDouble>(nSolutionPts);
229 
230  m_Q2D_e0 = Array<OneD, Array<OneD, NekDouble> >(nElements);
231  m_Q2D_e1 = Array<OneD, Array<OneD, NekDouble> >(nElements);
232  m_Q2D_e2 = Array<OneD, Array<OneD, NekDouble> >(nElements);
233  m_Q2D_e3 = Array<OneD, Array<OneD, NekDouble> >(nElements);
234 
235  for (n = 0; n < nElements; ++n)
236  {
237  base = pFields[0]->GetExp(n)->GetBase();
238  nquad0 = base[0]->GetNumPoints();
239  nquad1 = base[1]->GetNumPoints();
240 
241  m_Q2D_e0[n] = Array<OneD, NekDouble>(nquad0);
242  m_Q2D_e1[n] = Array<OneD, NekDouble>(nquad1);
243  m_Q2D_e2[n] = Array<OneD, NekDouble>(nquad0);
244  m_Q2D_e3[n] = Array<OneD, NekDouble>(nquad1);
245 
246  // Extract the Q factors at each edge point
247  pFields[0]->GetExp(n)->GetEdgeQFactors(
248  0, auxArray1 = m_Q2D_e0[n]);
249  pFields[0]->GetExp(n)->GetEdgeQFactors(
250  1, auxArray1 = m_Q2D_e1[n]);
251  pFields[0]->GetExp(n)->GetEdgeQFactors(
252  2, auxArray1 = m_Q2D_e2[n]);
253  pFields[0]->GetExp(n)->GetEdgeQFactors(
254  3, auxArray1 = m_Q2D_e3[n]);
255 
256  ptsKeys = pFields[0]->GetExp(n)->GetPointsKeys();
257  nLocalSolutionPts = pFields[0]->GetExp(n)->GetTotPoints();
258  phys_offset = pFields[0]->GetPhys_Offset(n);
259 
260  jac = pFields[0]->GetExp(n)
261  ->as<LocalRegions::Expansion2D>()->GetGeom2D()
262  ->GetMetricInfo()->GetJac(ptsKeys);
263  gmat = pFields[0]->GetExp(n)
264  ->as<LocalRegions::Expansion2D>()->GetGeom2D()
265  ->GetMetricInfo()->GetDerivFactors(ptsKeys);
266 
267  if (pFields[0]->GetExp(n)
268  ->as<LocalRegions::Expansion2D>()->GetGeom2D()
269  ->GetMetricInfo()->GetGtype()
270  == SpatialDomains::eDeformed)
271  {
272  for (i = 0; i < nLocalSolutionPts; ++i)
273  {
274  m_jac[i+phys_offset] = jac[i];
275  m_gmat[0][i+phys_offset] = gmat[0][i];
276  m_gmat[1][i+phys_offset] = gmat[1][i];
277  m_gmat[2][i+phys_offset] = gmat[2][i];
278  m_gmat[3][i+phys_offset] = gmat[3][i];
279  }
280  }
281  else
282  {
283  for (i = 0; i < nLocalSolutionPts; ++i)
284  {
285  m_jac[i+phys_offset] = jac[0];
286  m_gmat[0][i+phys_offset] = gmat[0][0];
287  m_gmat[1][i+phys_offset] = gmat[1][0];
288  m_gmat[2][i+phys_offset] = gmat[2][0];
289  m_gmat[3][i+phys_offset] = gmat[3][0];
290  }
291  }
292  }
293  break;
294  }
295  case 3:
296  {
297  ASSERTL0(false,"3DFR Metric terms not implemented yet");
298  break;
299  }
300  default:
301  {
302  ASSERTL0(false, "Expansion dimension not recognised");
303  break;
304  }
305  }
306  }
ASSERTL0
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:161
Nektar::SolverUtils::DiffusionLFR::m_Q2D_e2
Array< OneD, Array< OneD, NekDouble > > m_Q2D_e2
Definition: DiffusionLFR.h:60
Nektar::LibUtilities::PointsKeyVector
std::vector< PointsKey > PointsKeyVector
Definition: Points.h:220
Nektar::SolverUtils::DiffusionLFR::m_Q2D_e3
Array< OneD, Array< OneD, NekDouble > > m_Q2D_e3
Definition: DiffusionLFR.h:61
Nektar::SolverUtils::DiffusionLFR::m_jac
Array< OneD, NekDouble > m_jac
Definition: DiffusionLFR.h:55
Nektar::SolverUtils::DiffusionLFR::m_traceNormals
Array< OneD, Array< OneD, NekDouble > > m_traceNormals
Definition: DiffusionLFR.h:75
Nektar::SolverUtils::DiffusionLFR::m_Q2D_e1
Array< OneD, Array< OneD, NekDouble > > m_Q2D_e1
Definition: DiffusionLFR.h:59
Nektar::SolverUtils::DiffusionLFR::m_gmat
Array< OneD, Array< OneD, NekDouble > > m_gmat
Definition: DiffusionLFR.h:56
Nektar::SolverUtils::DiffusionLFR::m_Q2D_e0
Array< OneD, Array< OneD, NekDouble > > m_Q2D_e0
Definition: DiffusionLFR.h:58
Nektar::SpatialDomains::eDeformed
Geometry is curved or has non-constant factors.
Definition: SpatialDomains.hpp:57
void Nektar::SolverUtils::DiffusionLFR::v_WeakPenaltyforScalar ( const Array< OneD, MultiRegions::ExpListSharedPtr > &  fields,
const int  var,
const Array< OneD, const NekDouble > &  ufield,
Array< OneD, NekDouble > &  penaltyflux 
)
protectedvirtual

Imposes appropriate bcs for the 1st order derivatives.

Definition at line 1233 of file DiffusionLFR.cpp.

References Nektar::SpatialDomains::eDirichlet, Nektar::SpatialDomains::eNeumann, and Vmath::Vcopy().

Referenced by v_NumFluxforScalar().

1238  {
1239  // Variable initialisation
1240  int i, e, id1, id2;
1241  int nBndEdgePts, nBndEdges;
1242  int cnt = 0;
1243  int nBndRegions = fields[var]->GetBndCondExpansions().num_elements();
1244  int nTracePts = fields[0]->GetTrace()->GetTotPoints();
1245 
1246  // Extract physical values of the fields at the trace space
1247  Array<OneD, NekDouble > uplus(nTracePts);
1248  fields[var]->ExtractTracePhys(ufield, uplus);
1249 
1250  // Impose boundary conditions
1251  for (i = 0; i < nBndRegions; ++i)
1252  {
1253  // Number of boundary edges related to bcRegion 'i'
1254  nBndEdges = fields[var]->
1255  GetBndCondExpansions()[i]->GetExpSize();
1256 
1257  // Weakly impose boundary conditions by modifying flux values
1258  for (e = 0; e < nBndEdges ; ++e)
1259  {
1260  // Number of points on the expansion
1261  nBndEdgePts = fields[var]->
1262  GetBndCondExpansions()[i]->GetExp(e)->GetTotPoints();
1263 
1264  // Offset of the boundary expansion
1265  id1 = fields[var]->
1266  GetBndCondExpansions()[i]->GetPhys_Offset(e);
1267 
1268  // Offset of the trace space related to boundary expansion
1269  id2 = fields[0]->GetTrace()->
1270  GetPhys_Offset(fields[0]->GetTraceMap()->
1271  GetBndCondTraceToGlobalTraceMap(cnt++));
1272 
1273  // Dirichlet bcs ==> uflux = gD
1274  if (fields[var]->GetBndConditions()[i]->
1275  GetBoundaryConditionType() == SpatialDomains::eDirichlet)
1276  {
1277  Vmath::Vcopy(nBndEdgePts,
1278  &(fields[var]->
1279  GetBndCondExpansions()[i]->
1280  GetPhys())[id1], 1,
1281  &penaltyflux[id2], 1);
1282  }
1283  // Neumann bcs ==> uflux = u+
1284  else if ((fields[var]->GetBndConditions()[i])->
1285  GetBoundaryConditionType() == SpatialDomains::eNeumann)
1286  {
1287  Vmath::Vcopy(nBndEdgePts,
1288  &uplus[id2], 1,
1289  &penaltyflux[id2], 1);
1290  }
1291  }
1292  }
1293  }
Vmath::Vcopy
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.cpp:1047
void Nektar::SolverUtils::DiffusionLFR::v_WeakPenaltyforVector ( const Array< OneD, MultiRegions::ExpListSharedPtr > &  fields,
const int  var,
const int  dir,
const Array< OneD, const NekDouble > &  qfield,
Array< OneD, NekDouble > &  penaltyflux,
NekDouble  C11 
)
protectedvirtual

Imposes appropriate bcs for the 2nd order derivatives.

Definition at line 1395 of file DiffusionLFR.cpp.

References Nektar::SpatialDomains::eDirichlet, Nektar::SpatialDomains::eNeumann, m_traceNormals, and Vmath::Vmul().

Referenced by v_NumFluxforVector().

1402  {
1403  int i, e, id1, id2;
1404  int nBndEdges, nBndEdgePts;
1405  int nBndRegions = fields[var]->GetBndCondExpansions().num_elements();
1406  int nTracePts = fields[0]->GetTrace()->GetTotPoints();
1407 
1408  Array<OneD, NekDouble > uterm(nTracePts);
1409  Array<OneD, NekDouble > qtemp(nTracePts);
1410  int cnt = 0;
1411 
1412  fields[var]->ExtractTracePhys(qfield, qtemp);
1413 
1414  for (i = 0; i < nBndRegions; ++i)
1415  {
1416  nBndEdges = fields[var]->
1417  GetBndCondExpansions()[i]->GetExpSize();
1418 
1419  // Weakly impose boundary conditions by modifying flux values
1420  for (e = 0; e < nBndEdges ; ++e)
1421  {
1422  nBndEdgePts = fields[var]->
1423  GetBndCondExpansions()[i]->GetExp(e)->GetTotPoints();
1424 
1425  id1 = fields[var]->
1426  GetBndCondExpansions()[i]->GetPhys_Offset(e);
1427 
1428  id2 = fields[0]->GetTrace()->
1429  GetPhys_Offset(fields[0]->GetTraceMap()->
1430  GetBndCondTraceToGlobalTraceMap(cnt++));
1431 
1432  // For Dirichlet boundary condition:
1433  //qflux = q+ - C_11 (u+ - g_D) (nx, ny)
1434  if (fields[var]->GetBndConditions()[i]->
1435  GetBoundaryConditionType() == SpatialDomains::eDirichlet)
1436  {
1437  Vmath::Vmul(nBndEdgePts, &m_traceNormals[dir][id2], 1,
1438  &qtemp[id2], 1, &penaltyflux[id2], 1);
1439  }
1440  // For Neumann boundary condition: qflux = g_N
1441  else if ((fields[var]->GetBndConditions()[i])->
1442  GetBoundaryConditionType() == SpatialDomains::eNeumann)
1443  {
1444  Vmath::Vmul(nBndEdgePts, &m_traceNormals[dir][id2], 1,
1445  &(fields[var]->GetBndCondExpansions()[i]->
1446  GetPhys())[id1], 1, &penaltyflux[id2], 1);
1447  }
1448  }
1449  }
1450  }
Nektar::SolverUtils::DiffusionLFR::m_traceNormals
Array< OneD, Array< OneD, NekDouble > > m_traceNormals
Definition: DiffusionLFR.h:75
Vmath::Vmul
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

Member Data Documentation

Array<OneD, Array<OneD, Array<OneD, NekDouble> > > Nektar::SolverUtils::DiffusionLFR::m_BD1
protected

Definition at line 81 of file DiffusionLFR.h.

Referenced by v_Diffuse(), and v_InitObject().

Array<OneD, Array<OneD, Array<OneD, NekDouble> > > Nektar::SolverUtils::DiffusionLFR::m_D1
protected

Definition at line 82 of file DiffusionLFR.h.

Referenced by v_Diffuse(), and v_InitObject().

Array<OneD, Array<OneD, Array<OneD, NekDouble> > > Nektar::SolverUtils::DiffusionLFR::m_DD1
protected

Definition at line 83 of file DiffusionLFR.h.

Referenced by v_Diffuse(), and v_InitObject().

Array<OneD, Array<OneD, Array<OneD, NekDouble> > > Nektar::SolverUtils::DiffusionLFR::m_DFC1
protected

Definition at line 80 of file DiffusionLFR.h.

Referenced by v_Diffuse(), and v_InitObject().

Array<OneD, Array<OneD, Array<OneD, NekDouble> > > Nektar::SolverUtils::DiffusionLFR::m_DFC2
protected

Definition at line 85 of file DiffusionLFR.h.

Referenced by v_Diffuse(), and v_InitObject().

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

Definition at line 63 of file DiffusionLFR.h.

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

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

Definition at line 65 of file DiffusionLFR.h.

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

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

Definition at line 67 of file DiffusionLFR.h.

Referenced by v_SetupCFunctions().

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

Definition at line 64 of file DiffusionLFR.h.

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

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

Definition at line 66 of file DiffusionLFR.h.

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

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

Definition at line 68 of file DiffusionLFR.h.

Referenced by v_SetupCFunctions().

std::string Nektar::SolverUtils::DiffusionLFR::m_diffType
protected

Definition at line 92 of file DiffusionLFR.h.

Referenced by v_SetupCFunctions().

Array<OneD, Array<OneD, NekDouble> > Nektar::SolverUtils::DiffusionLFR::m_divFC
protected

Definition at line 87 of file DiffusionLFR.h.

Referenced by v_Diffuse(), and v_InitObject().

Array<OneD, Array<OneD, NekDouble> > Nektar::SolverUtils::DiffusionLFR::m_divFD
protected

Definition at line 86 of file DiffusionLFR.h.

Referenced by v_Diffuse(), and v_InitObject().

Array<OneD, Array<OneD, Array<OneD, NekDouble> > > Nektar::SolverUtils::DiffusionLFR::m_DU1
protected

Definition at line 79 of file DiffusionLFR.h.

Referenced by v_Diffuse(), and v_InitObject().

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

Definition at line 56 of file DiffusionLFR.h.

Referenced by v_Diffuse(), and v_SetupMetrics().

Array<OneD, Array<OneD, Array<OneD, NekDouble> > > Nektar::SolverUtils::DiffusionLFR::m_IF1
protected

Definition at line 78 of file DiffusionLFR.h.

Referenced by v_Diffuse(), and v_InitObject().

Array<OneD, Array<OneD, NekDouble> > Nektar::SolverUtils::DiffusionLFR::m_IF2
protected

Definition at line 84 of file DiffusionLFR.h.

Referenced by v_Diffuse(), and v_InitObject().

DNekMatSharedPtr Nektar::SolverUtils::DiffusionLFR::m_Ixm

Definition at line 69 of file DiffusionLFR.h.

DNekMatSharedPtr Nektar::SolverUtils::DiffusionLFR::m_Ixp

Definition at line 70 of file DiffusionLFR.h.

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

Definition at line 55 of file DiffusionLFR.h.

Referenced by v_Diffuse(), and v_SetupMetrics().

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

Definition at line 58 of file DiffusionLFR.h.

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

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

Definition at line 59 of file DiffusionLFR.h.

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

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

Definition at line 60 of file DiffusionLFR.h.

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

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

Definition at line 61 of file DiffusionLFR.h.

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

LibUtilities::SessionReaderSharedPtr Nektar::SolverUtils::DiffusionLFR::m_session
protected

Definition at line 76 of file DiffusionLFR.h.

Referenced by v_InitObject().

Array<OneD, Array<OneD, Array<OneD, NekDouble> > > Nektar::SolverUtils::DiffusionLFR::m_tmp1
protected

Definition at line 89 of file DiffusionLFR.h.

Referenced by v_Diffuse(), and v_InitObject().

Array<OneD, Array<OneD, Array<OneD, NekDouble> > > Nektar::SolverUtils::DiffusionLFR::m_tmp2
protected

Definition at line 90 of file DiffusionLFR.h.

Referenced by v_Diffuse(), and v_InitObject().

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

Definition at line 75 of file DiffusionLFR.h.

Referenced by v_DivCFlux_2D(), v_DivCFlux_2D_Gauss(), v_NumFluxforScalar(), v_NumFluxforVector(), v_SetupMetrics(), and v_WeakPenaltyforVector().

std::string Nektar::SolverUtils::DiffusionLFR::type
static
Initial value:
= {
GetDiffusionFactory().RegisterCreatorFunction(
"LFRDG", DiffusionLFR::create),
GetDiffusionFactory().RegisterCreatorFunction(
"LFRSD", DiffusionLFR::create),
GetDiffusionFactory().RegisterCreatorFunction(
"LFRHU", DiffusionLFR::create),
GetDiffusionFactory().RegisterCreatorFunction(
"LFRcmin", DiffusionLFR::create),
GetDiffusionFactory().RegisterCreatorFunction(
"LFRcinf", DiffusionLFR::create)}

Definition at line 53 of file DiffusionLFR.h.

Generated on Wed May 4 2016 17:54:35 for Nektar++ by   doxygen 1.8.8