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

#include <DiffusionLFRNS.h>

Inheritance diagram for Nektar::SolverUtils::DiffusionLFRNS:
Inheritance graph
[legend]
Collaboration diagram for Nektar::SolverUtils::DiffusionLFRNS:
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

 DiffusionLFRNS (std::string diffType)
 DiffusionLFRNS 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 DiffusionLFRNS 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 nConvective, const Array< OneD, MultiRegions::ExpListSharedPtr > &fields, const Array< OneD, Array< OneD, NekDouble > > &inarray, Array< OneD, Array< OneD, NekDouble > > &outarray)
 Calculate FR Diffusion for the Navier-Stokes (NS) equations using an LDG interface flux. More...
 
virtual void v_NumericalFluxO1 (const Array< OneD, MultiRegions::ExpListSharedPtr > &fields, const Array< OneD, Array< OneD, NekDouble > > &inarray, Array< OneD, Array< OneD, Array< OneD, NekDouble > > > &numericalFluxO1)
 Builds the numerical flux for the 1st order derivatives. More...
 
virtual void v_WeakPenaltyO1 (const Array< OneD, MultiRegions::ExpListSharedPtr > &fields, const Array< OneD, Array< OneD, NekDouble > > &inarray, Array< OneD, Array< OneD, NekDouble > > &penaltyfluxO1)
 Imposes appropriate bcs for the 1st order derivatives. More...
 
virtual void v_NumericalFluxO2 (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_WeakPenaltyO2 (const Array< OneD, MultiRegions::ExpListSharedPtr > &fields, const int var, const int dir, const Array< OneD, const NekDouble > &qfield, Array< OneD, NekDouble > &penaltyflux)
 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...
 
virtual void v_FluxVec (Array< OneD, Array< OneD, Array< OneD, NekDouble > > > &fluxvector)
 
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_traceVel
 
Array< OneD, Array< OneD,
NekDouble > > 		m_traceNormals
 
LibUtilities::SessionReaderSharedPtr m_session
 
NekDouble m_gamma
 
NekDouble m_gasConstant
 
NekDouble m_Twall
 
std::string m_ViscosityType
 
NekDouble m_mu
 
NekDouble m_thermalConductivity
 
NekDouble m_rhoInf
 
NekDouble m_pInf
 
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,
Array< OneD, NekDouble > > > 		m_viscTensor
 
Array< OneD, Array< OneD,
NekDouble > > 		m_viscFlux
 
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
 
Array< OneD, Array< OneD,
NekDouble > > 		m_homoDerivs
 
int m_spaceDim
 
int m_diffDim
 
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 DiffusionLFRNS.h.

Constructor & Destructor Documentation

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

DiffusionLFRNS 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 68 of file DiffusionLFRNS.cpp.

Referenced by create().

68  :m_diffType(diffType)
69  {
70  }

Member Function Documentation

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

Definition at line 48 of file DiffusionLFRNS.h.

References DiffusionLFRNS().

49  {
50  return DiffusionSharedPtr(new DiffusionLFRNS(diffType));
51  }
Nektar::SolverUtils::DiffusionLFRNS::DiffusionLFRNS
DiffusionLFRNS(std::string diffType)
DiffusionLFRNS uses the Flux Reconstruction (FR) approach to compute the diffusion term...
Definition: DiffusionLFRNS.cpp:68
Nektar::SolverUtils::DiffusionSharedPtr
boost::shared_ptr< Diffusion > DiffusionSharedPtr
A shared pointer to an EquationSystem object.
Definition: Diffusion.h:164
void Nektar::SolverUtils::DiffusionLFRNS::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 1651 of file DiffusionLFRNS.cpp.

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

Referenced by v_Diffuse().

1657  {
1658  int n;
1659  int nLocalSolutionPts, phys_offset;
1660 
1661  Array<OneD, NekDouble> auxArray1, auxArray2;
1662  Array<TwoD, const NekDouble> gmat;
1663  Array<OneD, const NekDouble> jac;
1664 
1665  LibUtilities::PointsKeyVector ptsKeys;
1666  LibUtilities::BasisSharedPtr Basis;
1667  Basis = fields[0]->GetExp(0)->GetBasis(0);
1668 
1669  int nElements = fields[0]->GetExpSize();
1670  int nPts = fields[0]->GetTotPoints();
1671 
1672  // Arrays to store the derivatives of the correction flux
1673  Array<OneD, NekDouble> DCL(nPts/nElements, 0.0);
1674  Array<OneD, NekDouble> DCR(nPts/nElements, 0.0);
1675 
1676  // Arrays to store the intercell numerical flux jumps
1677  Array<OneD, NekDouble> JumpL(nElements);
1678  Array<OneD, NekDouble> JumpR(nElements);
1679 
1680  for (n = 0; n < nElements; ++n)
1681  {
1682  nLocalSolutionPts = fields[0]->GetExp(n)->GetTotPoints();
1683  phys_offset = fields[0]->GetPhys_Offset(n);
1684 
1685  Array<OneD, NekDouble> tmparrayX1(nLocalSolutionPts, 0.0);
1686  NekDouble tmpFluxVertex = 0;
1687  Vmath::Vcopy(nLocalSolutionPts,
1688  &flux[phys_offset], 1,
1689  &tmparrayX1[0], 1);
1690 
1691  fields[0]->GetExp(n)->GetVertexPhysVals(0, tmparrayX1,
1692  tmpFluxVertex);
1693  JumpL[n] = iFlux[n] - tmpFluxVertex;
1694 
1695  fields[0]->GetExp(n)->GetVertexPhysVals(1, tmparrayX1,
1696  tmpFluxVertex);
1697  JumpR[n] = iFlux[n+1] - tmpFluxVertex;
1698  }
1699 
1700  for (n = 0; n < nElements; ++n)
1701  {
1702  ptsKeys = fields[0]->GetExp(n)->GetPointsKeys();
1703  nLocalSolutionPts = fields[0]->GetExp(n)->GetTotPoints();
1704  phys_offset = fields[0]->GetPhys_Offset(n);
1705  jac = fields[0]->GetExp(n)
1706  ->as<LocalRegions::Expansion1D>()->GetGeom1D()
1707  ->GetMetricInfo()->GetJac(ptsKeys);
1708 
1709  JumpL[n] = JumpL[n] * jac[0];
1710  JumpR[n] = JumpR[n] * jac[0];
1711 
1712  // Left jump multiplied by left derivative of C function
1713  Vmath::Smul(nLocalSolutionPts, JumpL[n],
1714  &m_dGL_xi1[n][0], 1, &DCL[0], 1);
1715 
1716  // Right jump multiplied by right derivative of C function
1717  Vmath::Smul(nLocalSolutionPts, JumpR[n],
1718  &m_dGR_xi1[n][0], 1, &DCR[0], 1);
1719 
1720  // Assembling derivative of the correction flux
1721  Vmath::Vadd(nLocalSolutionPts, &DCL[0], 1, &DCR[0], 1,
1722  &derCFlux[phys_offset], 1);
1723  }
1724  }
Nektar::SolverUtils::DiffusionLFRNS::m_dGR_xi1
Array< OneD, Array< OneD, NekDouble > > m_dGR_xi1
Definition: DiffusionLFRNS.h:64
Nektar::LibUtilities::PointsKeyVector
std::vector< PointsKey > PointsKeyVector
Definition: Points.h:220
Nektar::SolverUtils::DiffusionLFRNS::m_dGL_xi1
Array< OneD, Array< OneD, NekDouble > > m_dGL_xi1
Definition: DiffusionLFRNS.h:63
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::StdRegions::StdExpansion::GetMetricInfo
const boost::shared_ptr< SpatialDomains::GeomFactors > & GetMetricInfo(void) const
Definition: StdExpansion.h:1179
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::DiffusionLFRNS::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 1741 of file DiffusionLFRNS.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().

1748  {
1749  int n, e, i, j, cnt;
1750 
1751  Array<OneD, const NekDouble> jac;
1752 
1753  int nElements = fields[0]->GetExpSize();
1754  int trace_offset, phys_offset;
1755  int nLocalSolutionPts;
1756  int nquad0, nquad1;
1757  int nEdgePts;
1758 
1759  LibUtilities::PointsKeyVector ptsKeys;
1760  Array<OneD, NekDouble> auxArray1, auxArray2;
1761  Array<OneD, LibUtilities::BasisSharedPtr> base;
1762  Array<OneD, Array<OneD, LocalRegions::ExpansionSharedPtr> >
1763  &elmtToTrace = fields[0]->GetTraceMap()->GetElmtToTrace();
1764 
1765  // Loop on the elements
1766  for (n = 0; n < nElements; ++n)
1767  {
1768  // Offset of the element on the global vector
1769  phys_offset = fields[0]->GetPhys_Offset(n);
1770  nLocalSolutionPts = fields[0]->GetExp(n)->GetTotPoints();
1771  ptsKeys = fields[0]->GetExp(n)->GetPointsKeys();
1772 
1773  jac = fields[0]->GetExp(n)->as<LocalRegions::Expansion2D>()
1774  ->GetGeom2D()->GetMetricInfo()->GetJac(ptsKeys);
1775 
1776  base = fields[0]->GetExp(n)->GetBase();
1777  nquad0 = base[0]->GetNumPoints();
1778  nquad1 = base[1]->GetNumPoints();
1779 
1780  Array<OneD, NekDouble> divCFluxE0(nLocalSolutionPts, 0.0);
1781  Array<OneD, NekDouble> divCFluxE1(nLocalSolutionPts, 0.0);
1782  Array<OneD, NekDouble> divCFluxE2(nLocalSolutionPts, 0.0);
1783  Array<OneD, NekDouble> divCFluxE3(nLocalSolutionPts, 0.0);
1784 
1785  // Loop on the edges
1786  for (e = 0; e < fields[0]->GetExp(n)->GetNedges(); ++e)
1787  {
1788  // Number of edge points of edge 'e'
1789  nEdgePts = fields[0]->GetExp(n)->GetEdgeNumPoints(e);
1790 
1791  // Array for storing volumetric fluxes on each edge
1792  Array<OneD, NekDouble> tmparray(nEdgePts, 0.0);
1793 
1794  // Offset of the trace space correspondent to edge 'e'
1795  trace_offset = fields[0]->GetTrace()->GetPhys_Offset(
1796  elmtToTrace[n][e]->GetElmtId());
1797 
1798  // Get the normals of edge 'e'
1799  //const Array<OneD, const Array<OneD, NekDouble> > &normals =
1800  //fields[0]->GetExp(n)->GetEdgeNormal(e);
1801 
1802  // Extract the edge values of the volumetric fluxes
1803  // on edge 'e' and order them accordingly to the order
1804  // of the trace space
1805  fields[0]->GetExp(n)->GetEdgePhysVals(
1806  e, elmtToTrace[n][e],
1807  flux + phys_offset,
1808  auxArray1 = tmparray);
1809 
1810  // Splitting inarray into the 'j' direction
1811  /*Vmath::Vmul(nEdgePts,
1812  &m_traceNormals[direction][trace_offset], 1,
1813  &tmparray[0], 1, &tmparray[0], 1);*/
1814 
1815  // Compute the fluxJumps per each edge 'e' and each
1816  // flux point
1817  Array<OneD, NekDouble> fluxJumps(nEdgePts, 0.0);
1818  Vmath::Vsub(nEdgePts, &iFlux[trace_offset], 1,
1819  &tmparray[0], 1, &fluxJumps[0], 1);
1820 
1821  // Check the ordering of the fluxJumps and reverse
1822  // it in case of backward definition of edge 'e'
1823  if (fields[0]->GetExp(n)->GetEorient(e) ==
1824  StdRegions::eBackwards)
1825  {
1826  Vmath::Reverse(nEdgePts,
1827  &fluxJumps[0], 1,
1828  &fluxJumps[0], 1);
1829  }
1830 
1831  // Deformed elements
1832  if (fields[0]->GetExp(n)->as<LocalRegions::Expansion2D>()
1833  ->GetGeom2D()->GetMetricInfo()->GetGtype()
1834  == SpatialDomains::eDeformed)
1835  {
1836  // Extract the Jacobians along edge 'e'
1837  Array<OneD, NekDouble> jacEdge(nEdgePts, 0.0);
1838  fields[0]->GetExp(n)->GetEdgePhysVals(
1839  e, elmtToTrace[n][e],
1840  jac, auxArray1 = jacEdge);
1841 
1842  // Check the ordering of the fluxJumps and reverse
1843  // it in case of backward definition of edge 'e'
1844  if (fields[0]->GetExp(n)->GetEorient(e) ==
1845  StdRegions::eBackwards)
1846  {
1847  Vmath::Reverse(nEdgePts, &jacEdge[0], 1,
1848  &jacEdge[0], 1);
1849  }
1850 
1851  // Multiply the fluxJumps by the edge 'e' Jacobians
1852  // to bring the fluxJumps into the standard space
1853  for (j = 0; j < nEdgePts; j++)
1854  {
1855  fluxJumps[j] = fluxJumps[j] * jacEdge[j];
1856  }
1857  }
1858  // Non-deformed elements
1859  else
1860  {
1861  // Multiply the fluxJumps by the edge 'e' Jacobians
1862  // to bring the fluxJumps into the standard space
1863  Vmath::Smul(nEdgePts, jac[0], fluxJumps, 1,
1864  fluxJumps, 1);
1865  }
1866 
1867  // Multiply jumps by derivatives of the correction functions
1868  // All the quntities at this level should be defined into
1869  // the standard space
1870  switch (e)
1871  {
1872  case 0:
1873  for (i = 0; i < nquad0; ++i)
1874  {
1875  // Multiply fluxJumps by Q factors
1876  //fluxJumps[i] = -(m_Q2D_e0[n][i]) * fluxJumps[i];
1877 
1878  for (j = 0; j < nquad1; ++j)
1879  {
1880  cnt = i + j*nquad0;
1881  divCFluxE0[cnt] = fluxJumps[i] * m_dGL_xi2[n][j];
1882  }
1883  }
1884  break;
1885  case 1:
1886  for (i = 0; i < nquad1; ++i)
1887  {
1888  // Multiply fluxJumps by Q factors
1889  //fluxJumps[i] = (m_Q2D_e1[n][i]) * fluxJumps[i];
1890 
1891  for (j = 0; j < nquad0; ++j)
1892  {
1893  cnt = (nquad0)*i + j;
1894  divCFluxE1[cnt] = fluxJumps[i] * m_dGR_xi1[n][j];
1895  }
1896  }
1897  break;
1898  case 2:
1899  for (i = 0; i < nquad0; ++i)
1900  {
1901  // Multiply fluxJumps by Q factors
1902  //fluxJumps[i] = (m_Q2D_e2[n][i]) * fluxJumps[i];
1903 
1904  for (j = 0; j < nquad1; ++j)
1905  {
1906  cnt = (nquad0 - 1) + j*nquad0 - i;
1907  divCFluxE2[cnt] = fluxJumps[i] * m_dGR_xi2[n][j];
1908  }
1909  }
1910  break;
1911  case 3:
1912  for (i = 0; i < nquad1; ++i)
1913  {
1914  // Multiply fluxJumps by Q factors
1915  //fluxJumps[i] = -(m_Q2D_e3[n][i]) * fluxJumps[i];
1916  for (j = 0; j < nquad0; ++j)
1917  {
1918  cnt = (nquad0*nquad1 - nquad0) + j - i*nquad0;
1919  divCFluxE3[cnt] = fluxJumps[i] * m_dGL_xi1[n][j];
1920  }
1921  }
1922  break;
1923 
1924  default:
1925  ASSERTL0(false, "edge value (< 3) is out of range");
1926  break;
1927  }
1928  }
1929 
1930 
1931  // Sum all the edge contributions since I am passing the
1932  // component of the flux x and y already. So I should not
1933  // need to sum E0-E2 to get the derivative wrt xi2 and E1-E3
1934  // to get the derivative wrt xi1
1935 
1936  if (direction == 0)
1937  {
1938  Vmath::Vadd(nLocalSolutionPts, &divCFluxE1[0], 1,
1939  &divCFluxE3[0], 1, &derCFlux[phys_offset], 1);
1940  }
1941  else if (direction == 1)
1942  {
1943  Vmath::Vadd(nLocalSolutionPts, &divCFluxE0[0], 1,
1944  &divCFluxE2[0], 1, &derCFlux[phys_offset], 1);
1945  }
1946  }
1947  }
Nektar::SolverUtils::DiffusionLFRNS::m_dGR_xi1
Array< OneD, Array< OneD, NekDouble > > m_dGR_xi1
Definition: DiffusionLFRNS.h:64
ASSERTL0
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:161
Nektar::LibUtilities::PointsKeyVector
std::vector< PointsKey > PointsKeyVector
Definition: Points.h:220
Nektar::SolverUtils::DiffusionLFRNS::m_dGR_xi2
Array< OneD, Array< OneD, NekDouble > > m_dGR_xi2
Definition: DiffusionLFRNS.h:66
Nektar::SolverUtils::DiffusionLFRNS::m_dGL_xi1
Array< OneD, Array< OneD, NekDouble > > m_dGL_xi1
Definition: DiffusionLFRNS.h:63
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
Nektar::SolverUtils::DiffusionLFRNS::m_dGL_xi2
Array< OneD, Array< OneD, NekDouble > > m_dGL_xi2
Definition: DiffusionLFRNS.h:65
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::StdRegions::StdExpansion::GetMetricInfo
const boost::shared_ptr< SpatialDomains::GeomFactors > & GetMetricInfo(void) const
Definition: StdExpansion.h:1179
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::DiffusionLFRNS::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 Navier-Stokes (NS) equations using an LDG interface flux.

The equations that need a diffusion operator are those related with the velocities and with the energy.

Implements Nektar::SolverUtils::Diffusion.

Definition at line 901 of file DiffusionLFRNS.cpp.

References ASSERTL0, Nektar::LibUtilities::eGaussGaussLegendre, m_BD1, m_D1, m_DD1, m_DFC1, m_DFC2, m_diffDim, m_divFC, m_divFD, m_DU1, Nektar::SolverUtils::Diffusion::m_fluxVectorNS, m_gmat, m_homoDerivs, m_IF1, m_jac, m_tmp1, m_tmp2, m_viscFlux, m_viscTensor, v_DerCFlux_1D(), v_DerCFlux_2D(), v_DivCFlux_2D(), v_DivCFlux_2D_Gauss(), v_NumericalFluxO1(), v_NumericalFluxO2(), Vmath::Vadd(), Vmath::Vcopy(), Vmath::Vdiv(), and Vmath::Vmul().

906  {
907  int i, j, n;
908  int phys_offset;
909 
910  Array<TwoD, const NekDouble> gmat;
911  Array<OneD, const NekDouble> jac;
912  Array<OneD, NekDouble> auxArray1, auxArray2;
913 
914  Array<OneD, LibUtilities::BasisSharedPtr> Basis;
915  Basis = fields[0]->GetExp(0)->GetBase();
916 
917  int nElements = fields[0]->GetExpSize();
918  int nDim = fields[0]->GetCoordim(0);
919  int nScalars = inarray.num_elements();
920  int nSolutionPts = fields[0]->GetTotPoints();
921  int nCoeffs = fields[0]->GetNcoeffs();
922 
923  Array<OneD, Array<OneD, NekDouble> > outarrayCoeff(nConvectiveFields);
924  for (i = 0; i < nConvectiveFields; ++i)
925  {
926  outarrayCoeff[i] = Array<OneD, NekDouble>(nCoeffs);
927  }
928 
929  // Compute interface numerical fluxes for inarray in physical space
930  v_NumericalFluxO1(fields, inarray, m_IF1);
931 
932  switch(nDim)
933  {
934  // 1D problems
935  case 1:
936  {
937  for (i = 0; i < nScalars; ++i)
938  {
939  // Computing the physical first-order discountinuous
940  // derivative
941  for (n = 0; n < nElements; n++)
942  {
943  phys_offset = fields[0]->GetPhys_Offset(n);
944 
945  fields[i]->GetExp(n)->PhysDeriv(0,
946  auxArray1 = inarray[i] + phys_offset,
947  auxArray2 = m_DU1[i][0] + phys_offset);
948  }
949 
950  // Computing the standard first-order correction
951  // derivative
952  v_DerCFlux_1D(nConvectiveFields, fields, inarray[i],
953  m_IF1[i][0], m_DFC1[i][0]);
954 
955  // Back to the physical space using global operations
956  Vmath::Vdiv(nSolutionPts, &m_DFC1[i][0][0], 1,
957  &m_jac[0], 1, &m_DFC1[i][0][0], 1);
958  Vmath::Vdiv(nSolutionPts, &m_DFC1[i][0][0], 1,
959  &m_jac[0], 1, &m_DFC1[i][0][0], 1);
960 
961  // Computing total first order derivatives
962  Vmath::Vadd(nSolutionPts, &m_DFC1[i][0][0], 1,
963  &m_DU1[i][0][0], 1, &m_D1[i][0][0], 1);
964 
965  Vmath::Vcopy(nSolutionPts, &m_D1[i][0][0], 1,
966  &m_tmp1[i][0][0], 1);
967  }
968 
969  // Computing the viscous tensor
970  m_fluxVectorNS(inarray, m_D1, m_viscTensor);
971 
972  // Compute u from q_{\eta} and q_{\xi}
973  // Obtain numerical fluxes
974  v_NumericalFluxO2(fields, inarray, m_viscTensor,
975  m_viscFlux);
976 
977  for (i = 0; i < nConvectiveFields; ++i)
978  {
979  // Computing the physical second-order discountinuous
980  // derivative
981  for (n = 0; n < nElements; n++)
982  {
983  phys_offset = fields[0]->GetPhys_Offset(n);
984 
985  fields[i]->GetExp(n)->PhysDeriv(0,
986  auxArray1 = m_viscTensor[0][i] + phys_offset,
987  auxArray2 = m_DD1[i][0] + phys_offset);
988  }
989 
990  // Computing the standard second-order correction
991  // derivative
992  v_DerCFlux_1D(nConvectiveFields, fields,
993  m_viscTensor[0][i], m_viscFlux[i],
994  m_DFC2[i][0]);
995 
996  // Back to the physical space using global operations
997  Vmath::Vdiv(nSolutionPts, &m_DFC2[i][0][0], 1,
998  &m_jac[0], 1, &m_DFC2[i][0][0], 1);
999  Vmath::Vdiv(nSolutionPts, &m_DFC2[i][0][0], 1,
1000  &m_jac[0], 1, &m_DFC2[i][0][0], 1);
1001 
1002  // Adding the total divergence to outarray (RHS)
1003  Vmath::Vadd(nSolutionPts, &m_DFC2[i][0][0], 1,
1004  &m_DD1[i][0][0], 1, &outarray[i][0], 1);
1005 
1006  // Primitive Dealiasing 1D
1007  if(!(Basis[0]->Collocation()))
1008  {
1009  fields[i]->FwdTrans(outarray[i], outarrayCoeff[i]);
1010  fields[i]->BwdTrans(outarrayCoeff[i], outarray[i]);
1011  }
1012  }
1013  break;
1014  }
1015  // 2D problems
1016  case 2:
1017  {
1018  for(i = 0; i < nScalars; ++i)
1019  {
1020  for (j = 0; j < nDim; ++j)
1021  {
1022  // Temporary vectors
1023  Array<OneD, NekDouble> u1_hat(nSolutionPts, 0.0);
1024  Array<OneD, NekDouble> u2_hat(nSolutionPts, 0.0);
1025 
1026  if (j == 0)
1027  {
1028  Vmath::Vmul(nSolutionPts, &inarray[i][0], 1,
1029  &m_gmat[0][0], 1, &u1_hat[0], 1);
1030 
1031  Vmath::Vmul(nSolutionPts, &u1_hat[0], 1,
1032  &m_jac[0], 1, &u1_hat[0], 1);
1033 
1034  Vmath::Vmul(nSolutionPts, &inarray[i][0], 1,
1035  &m_gmat[1][0], 1, &u2_hat[0], 1);
1036 
1037  Vmath::Vmul(nSolutionPts, &u2_hat[0], 1,
1038  &m_jac[0], 1, &u2_hat[0], 1);
1039  }
1040  else if (j == 1)
1041  {
1042  Vmath::Vmul(nSolutionPts, &inarray[i][0], 1,
1043  &m_gmat[2][0], 1, &u1_hat[0], 1);
1044 
1045  Vmath::Vmul(nSolutionPts, &u1_hat[0], 1,
1046  &m_jac[0], 1, &u1_hat[0], 1);
1047 
1048  Vmath::Vmul(nSolutionPts, &inarray[i][0], 1,
1049  &m_gmat[3][0], 1, &u2_hat[0], 1);
1050 
1051  Vmath::Vmul(nSolutionPts, &u2_hat[0], 1,
1052  &m_jac[0], 1, &u2_hat[0], 1);
1053  }
1054 
1055  for (n = 0; n < nElements; n++)
1056  {
1057  phys_offset = fields[0]->GetPhys_Offset(n);
1058 
1059  fields[i]->GetExp(n)->StdPhysDeriv(0,
1060  auxArray1 = u1_hat + phys_offset,
1061  auxArray2 = m_tmp1[i][j] + phys_offset);
1062 
1063  fields[i]->GetExp(n)->StdPhysDeriv(1,
1064  auxArray1 = u2_hat + phys_offset,
1065  auxArray2 = m_tmp2[i][j] + phys_offset);
1066  }
1067 
1068  Vmath::Vadd(nSolutionPts, &m_tmp1[i][j][0], 1,
1069  &m_tmp2[i][j][0], 1,
1070  &m_DU1[i][j][0], 1);
1071 
1072  // Divide by the metric jacobian
1073  Vmath::Vdiv(nSolutionPts, &m_DU1[i][j][0], 1,
1074  &m_jac[0], 1, &m_DU1[i][j][0], 1);
1075 
1076  // Computing the standard first-order correction
1077  // derivatives
1078  v_DerCFlux_2D(nConvectiveFields, j, fields,
1079  inarray[i], m_IF1[i][j],
1080  m_DFC1[i][j]);
1081  }
1082 
1083  // Multiplying derivatives by B matrix to get auxiliary
1084  // variables
1085  for (j = 0; j < nSolutionPts; ++j)
1086  {
1087  // std(q1)
1088  m_BD1[i][0][j] =
1089  (m_gmat[0][j]*m_gmat[0][j] +
1090  m_gmat[2][j]*m_gmat[2][j]) *
1091  m_DFC1[i][0][j] +
1092  (m_gmat[1][j]*m_gmat[0][j] +
1093  m_gmat[3][j]*m_gmat[2][j]) *
1094  m_DFC1[i][1][j];
1095 
1096  // std(q2)
1097  m_BD1[i][1][j] =
1098  (m_gmat[1][j]*m_gmat[0][j] +
1099  m_gmat[3][j]*m_gmat[2][j]) *
1100  m_DFC1[i][0][j] +
1101  (m_gmat[1][j]*m_gmat[1][j] +
1102  m_gmat[3][j]*m_gmat[3][j]) *
1103  m_DFC1[i][1][j];
1104  }
1105 
1106  // Multiplying derivatives by A^(-1) to get back
1107  // into the physical space
1108  for (j = 0; j < nSolutionPts; j++)
1109  {
1110  // q1 = A11^(-1)*std(q1) + A12^(-1)*std(q2)
1111  m_DFC1[i][0][j] =
1112  m_gmat[3][j] * m_BD1[i][0][j] -
1113  m_gmat[2][j] * m_BD1[i][1][j];
1114 
1115  // q2 = A21^(-1)*std(q1) + A22^(-1)*std(q2)
1116  m_DFC1[i][1][j] =
1117  m_gmat[0][j] * m_BD1[i][1][j] -
1118  m_gmat[1][j] * m_BD1[i][0][j];
1119  }
1120 
1121  // Computing the physical first-order derivatives
1122  for (j = 0; j < nDim; ++j)
1123  {
1124  Vmath::Vadd(nSolutionPts, &m_DU1[i][j][0], 1,
1125  &m_DFC1[i][j][0], 1,
1126  &m_D1[j][i][0], 1);
1127  }
1128  }// Close loop on nScalars
1129 
1130  // For 3D Homogeneous 1D only take derivatives in 3rd direction
1131  if (m_diffDim == 1)
1132  {
1133  for (i = 0; i < nScalars; ++i)
1134  {
1135  m_D1[2][i] = m_homoDerivs[i];
1136  }
1137 
1138  }
1139 
1140  // Computing the viscous tensor
1141  m_fluxVectorNS(inarray, m_D1, m_viscTensor);
1142 
1143  // Compute u from q_{\eta} and q_{\xi}
1144  // Obtain numerical fluxes
1145  v_NumericalFluxO2(fields, inarray, m_viscTensor,
1146  m_viscFlux);
1147 
1148  // Computing the standard second-order discontinuous
1149  // derivatives
1150  for (i = 0; i < nConvectiveFields; ++i)
1151  {
1152  // Temporary vectors
1153  Array<OneD, NekDouble> f_hat(nSolutionPts, 0.0);
1154  Array<OneD, NekDouble> g_hat(nSolutionPts, 0.0);
1155 
1156  for (j = 0; j < nSolutionPts; j++)
1157  {
1158  f_hat[j] = (m_viscTensor[0][i][j] * m_gmat[0][j] +
1159  m_viscTensor[1][i][j] * m_gmat[2][j])*m_jac[j];
1160 
1161  g_hat[j] = (m_viscTensor[0][i][j] * m_gmat[1][j] +
1162  m_viscTensor[1][i][j] * m_gmat[3][j])*m_jac[j];
1163  }
1164 
1165  for (n = 0; n < nElements; n++)
1166  {
1167  phys_offset = fields[0]->GetPhys_Offset(n);
1168 
1169  fields[0]->GetExp(n)->StdPhysDeriv(0,
1170  auxArray1 = f_hat + phys_offset,
1171  auxArray2 = m_DD1[i][0] + phys_offset);
1172 
1173  fields[0]->GetExp(n)->StdPhysDeriv(1,
1174  auxArray1 = g_hat + phys_offset,
1175  auxArray2 = m_DD1[i][1] + phys_offset);
1176  }
1177 
1178  // Divergence of the standard discontinuous flux
1179  Vmath::Vadd(nSolutionPts, &m_DD1[i][0][0], 1,
1180  &m_DD1[i][1][0], 1, &m_divFD[i][0], 1);
1181 
1182  // Divergence of the standard correction flux
1183  if (Basis[0]->GetPointsType() ==
1184  LibUtilities::eGaussGaussLegendre &&
1185  Basis[1]->GetPointsType() ==
1186  LibUtilities::eGaussGaussLegendre)
1187  {
1188  v_DivCFlux_2D_Gauss(nConvectiveFields, fields,
1189  f_hat, g_hat,
1190  m_viscFlux[i], m_divFC[i]);
1191  }
1192  else
1193  {
1194  v_DivCFlux_2D(nConvectiveFields, fields,
1195  m_viscTensor[0][i], m_viscTensor[1][i],
1196  m_viscFlux[i], m_divFC[i]);
1197 
1198  }
1199 
1200 
1201  // Divergence of the standard final flux
1202  Vmath::Vadd(nSolutionPts, &m_divFD[i][0], 1,
1203  &m_divFC[i][0], 1, &outarray[i][0], 1);
1204 
1205  // Dividing by the jacobian to get back into
1206  // physical space
1207  Vmath::Vdiv(nSolutionPts, &outarray[i][0], 1,
1208  &m_jac[0], 1, &outarray[i][0], 1);
1209 
1210  // Primitive Dealiasing 2D
1211  if(!(Basis[0]->Collocation()))
1212  {
1213  fields[i]->FwdTrans(outarray[i], outarrayCoeff[i]);
1214  fields[i]->BwdTrans(outarrayCoeff[i], outarray[i]);
1215  }
1216  }
1217  break;
1218  }
1219  // 3D problems
1220  case 3:
1221  {
1222  ASSERTL0(false, "3D FRDG case not implemented yet");
1223  break;
1224  }
1225  }
1226  }
Nektar::SolverUtils::DiffusionLFRNS::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: DiffusionLFRNS.cpp:1741
ASSERTL0
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:161
Nektar::SolverUtils::DiffusionLFRNS::m_DU1
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > m_DU1
Definition: DiffusionLFRNS.h:89
Nektar::SolverUtils::DiffusionLFRNS::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: DiffusionLFRNS.cpp:1962
Nektar::SolverUtils::DiffusionLFRNS::v_NumericalFluxO1
virtual void v_NumericalFluxO1(const Array< OneD, MultiRegions::ExpListSharedPtr > &fields, const Array< OneD, Array< OneD, NekDouble > > &inarray, Array< OneD, Array< OneD, Array< OneD, NekDouble > > > &numericalFluxO1)
Builds the numerical flux for the 1st order derivatives.
Definition: DiffusionLFRNS.cpp:1233
Nektar::SolverUtils::DiffusionLFRNS::m_BD1
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > m_BD1
Definition: DiffusionLFRNS.h:91
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::SolverUtils::DiffusionLFRNS::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: DiffusionLFRNS.cpp:2158
Nektar::SolverUtils::Diffusion::m_fluxVectorNS
DiffusionFluxVecCBNS m_fluxVectorNS
Definition: Diffusion.h:133
Nektar::SolverUtils::DiffusionLFRNS::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: DiffusionLFRNS.cpp:1651
Nektar::SolverUtils::DiffusionLFRNS::m_homoDerivs
Array< OneD, Array< OneD, NekDouble > > m_homoDerivs
Definition: DiffusionLFRNS.h:103
Nektar::SolverUtils::DiffusionLFRNS::m_jac
Array< OneD, NekDouble > m_jac
Definition: DiffusionLFRNS.h:55
Nektar::LibUtilities::eGaussGaussLegendre
1D Gauss-Gauss-Legendre quadrature points
Definition: PointsType.h:47
Nektar::SolverUtils::DiffusionLFRNS::m_tmp2
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > m_tmp2
Definition: DiffusionLFRNS.h:101
Nektar::SolverUtils::DiffusionLFRNS::m_tmp1
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > m_tmp1
Definition: DiffusionLFRNS.h:100
Nektar::SolverUtils::DiffusionLFRNS::m_divFD
Array< OneD, Array< OneD, NekDouble > > m_divFD
Definition: DiffusionLFRNS.h:97
Nektar::SolverUtils::DiffusionLFRNS::m_IF1
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > m_IF1
Definition: DiffusionLFRNS.h:88
Nektar::SolverUtils::DiffusionLFRNS::m_DD1
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > m_DD1
Definition: DiffusionLFRNS.h:93
Nektar::SolverUtils::DiffusionLFRNS::m_gmat
Array< OneD, Array< OneD, NekDouble > > m_gmat
Definition: DiffusionLFRNS.h:56
Nektar::SolverUtils::DiffusionLFRNS::m_DFC1
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > m_DFC1
Definition: DiffusionLFRNS.h:90
Nektar::SolverUtils::DiffusionLFRNS::m_viscTensor
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > m_viscTensor
Definition: DiffusionLFRNS.h:94
Nektar::SolverUtils::DiffusionLFRNS::m_D1
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > m_D1
Definition: DiffusionLFRNS.h:92
Nektar::SolverUtils::DiffusionLFRNS::m_DFC2
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > m_DFC2
Definition: DiffusionLFRNS.h:96
Nektar::SolverUtils::DiffusionLFRNS::m_viscFlux
Array< OneD, Array< OneD, NekDouble > > m_viscFlux
Definition: DiffusionLFRNS.h:95
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
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::DiffusionLFRNS::v_NumericalFluxO2
virtual void v_NumericalFluxO2(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: DiffusionLFRNS.cpp:1498
Nektar::SolverUtils::DiffusionLFRNS::m_divFC
Array< OneD, Array< OneD, NekDouble > > m_divFC
Definition: DiffusionLFRNS.h:98
void Nektar::SolverUtils::DiffusionLFRNS::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 1962 of file DiffusionLFRNS.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().

1969  {
1970  int n, e, i, j, cnt;
1971 
1972  int nElements = fields[0]->GetExpSize();
1973  int nLocalSolutionPts;
1974  int nEdgePts;
1975  int trace_offset;
1976  int phys_offset;
1977  int nquad0;
1978  int nquad1;
1979 
1980  Array<OneD, NekDouble> auxArray1, auxArray2;
1981  Array<OneD, LibUtilities::BasisSharedPtr> base;
1982 
1983  Array<OneD, Array<OneD, LocalRegions::ExpansionSharedPtr> >
1984  &elmtToTrace = fields[0]->GetTraceMap()->GetElmtToTrace();
1985 
1986  // Loop on the elements
1987  for(n = 0; n < nElements; ++n)
1988  {
1989  // Offset of the element on the global vector
1990  phys_offset = fields[0]->GetPhys_Offset(n);
1991  nLocalSolutionPts = fields[0]->GetExp(n)->GetTotPoints();
1992 
1993  base = fields[0]->GetExp(n)->GetBase();
1994  nquad0 = base[0]->GetNumPoints();
1995  nquad1 = base[1]->GetNumPoints();
1996 
1997  Array<OneD, NekDouble> divCFluxE0(nLocalSolutionPts, 0.0);
1998  Array<OneD, NekDouble> divCFluxE1(nLocalSolutionPts, 0.0);
1999  Array<OneD, NekDouble> divCFluxE2(nLocalSolutionPts, 0.0);
2000  Array<OneD, NekDouble> divCFluxE3(nLocalSolutionPts, 0.0);
2001 
2002  // Loop on the edges
2003  for(e = 0; e < fields[0]->GetExp(n)->GetNedges(); ++e)
2004  {
2005  // Number of edge points of edge e
2006  nEdgePts = fields[0]->GetExp(n)->GetEdgeNumPoints(e);
2007 
2008  Array<OneD, NekDouble> tmparrayX1(nEdgePts, 0.0);
2009  Array<OneD, NekDouble> tmparrayX2(nEdgePts, 0.0);
2010  Array<OneD, NekDouble> fluxN (nEdgePts, 0.0);
2011  Array<OneD, NekDouble> fluxT (nEdgePts, 0.0);
2012  Array<OneD, NekDouble> fluxJumps(nEdgePts, 0.0);
2013 
2014  // Offset of the trace space correspondent to edge e
2015  trace_offset = fields[0]->GetTrace()->GetPhys_Offset(
2016  elmtToTrace[n][e]->GetElmtId());
2017 
2018  // Get the normals of edge e
2019  const Array<OneD, const Array<OneD, NekDouble> > &normals =
2020  fields[0]->GetExp(n)->GetEdgeNormal(e);
2021 
2022  // Extract the edge values of flux-x on edge e and order
2023  // them accordingly to the order of the trace space
2024  fields[0]->GetExp(n)->GetEdgePhysVals(
2025  e, elmtToTrace[n][e],
2026  fluxX1 + phys_offset,
2027  auxArray1 = tmparrayX1);
2028 
2029  // Extract the edge values of flux-y on edge e and order
2030  // them accordingly to the order of the trace space
2031  fields[0]->GetExp(n)->GetEdgePhysVals(
2032  e, elmtToTrace[n][e],
2033  fluxX2 + phys_offset,
2034  auxArray1 = tmparrayX2);
2035 
2036  // Multiply the edge components of the flux by the normal
2037  for (i = 0; i < nEdgePts; ++i)
2038  {
2039  fluxN[i] =
2040  tmparrayX1[i]*m_traceNormals[0][trace_offset+i] +
2041  tmparrayX2[i]*m_traceNormals[1][trace_offset+i];
2042  }
2043 
2044  // Subtract to the Riemann flux the discontinuous flux
2045  Vmath::Vsub(nEdgePts,
2046  &numericalFlux[trace_offset], 1,
2047  &fluxN[0], 1, &fluxJumps[0], 1);
2048 
2049  // Check the ordering of the jump vectors
2050  if (fields[0]->GetExp(n)->GetEorient(e) ==
2051  StdRegions::eBackwards)
2052  {
2053  Vmath::Reverse(nEdgePts,
2054  auxArray1 = fluxJumps, 1,
2055  auxArray2 = fluxJumps, 1);
2056  }
2057 
2058  NekDouble fac = fields[0]->GetExp(n)->EdgeNormalNegated(e) ?
2059  -1.0 : 1.0;
2060 
2061  for (i = 0; i < nEdgePts; ++i)
2062  {
2063  if (m_traceNormals[0][trace_offset+i] != fac*normals[0][i]
2064  || m_traceNormals[1][trace_offset+i] != fac*normals[1][i])
2065  {
2066  fluxJumps[i] = -fluxJumps[i];
2067  }
2068  }
2069 
2070  // Multiply jumps by derivatives of the correction functions
2071  switch (e)
2072  {
2073  case 0:
2074  for (i = 0; i < nquad0; ++i)
2075  {
2076  // Multiply fluxJumps by Q factors
2077  fluxJumps[i] = -(m_Q2D_e0[n][i]) * fluxJumps[i];
2078 
2079  for (j = 0; j < nquad1; ++j)
2080  {
2081  cnt = i + j*nquad0;
2082  divCFluxE0[cnt] = fluxJumps[i] * m_dGL_xi2[n][j];
2083  }
2084  }
2085  break;
2086  case 1:
2087  for (i = 0; i < nquad1; ++i)
2088  {
2089  // Multiply fluxJumps by Q factors
2090  fluxJumps[i] = (m_Q2D_e1[n][i]) * fluxJumps[i];
2091 
2092  for (j = 0; j < nquad0; ++j)
2093  {
2094  cnt = (nquad0)*i + j;
2095  divCFluxE1[cnt] = fluxJumps[i] * m_dGR_xi1[n][j];
2096  }
2097  }
2098  break;
2099  case 2:
2100  for (i = 0; i < nquad0; ++i)
2101  {
2102  // Multiply fluxJumps by Q factors
2103  fluxJumps[i] = (m_Q2D_e2[n][i]) * fluxJumps[i];
2104 
2105  for (j = 0; j < nquad1; ++j)
2106  {
2107  cnt = (nquad0 - 1) + j*nquad0 - i;
2108  divCFluxE2[cnt] = fluxJumps[i] * m_dGR_xi2[n][j];
2109  }
2110  }
2111  break;
2112  case 3:
2113  for (i = 0; i < nquad1; ++i)
2114  {
2115  // Multiply fluxJumps by Q factors
2116  fluxJumps[i] = -(m_Q2D_e3[n][i]) * fluxJumps[i];
2117  for (j = 0; j < nquad0; ++j)
2118  {
2119  cnt = (nquad0*nquad1 - nquad0) + j - i*nquad0;
2120  divCFluxE3[cnt] = fluxJumps[i] * m_dGL_xi1[n][j];
2121 
2122  }
2123  }
2124  break;
2125 
2126  default:
2127  ASSERTL0(false,"edge value (< 3) is out of range");
2128  break;
2129  }
2130  }
2131 
2132  // Sum each edge contribution
2133  Vmath::Vadd(nLocalSolutionPts, &divCFluxE0[0], 1,
2134  &divCFluxE1[0], 1, &divCFlux[phys_offset], 1);
2135 
2136  Vmath::Vadd(nLocalSolutionPts, &divCFlux[phys_offset], 1,
2137  &divCFluxE2[0], 1, &divCFlux[phys_offset], 1);
2138 
2139  Vmath::Vadd(nLocalSolutionPts, &divCFlux[phys_offset], 1,
2140  &divCFluxE3[0], 1, &divCFlux[phys_offset], 1);
2141  }
2142  }
Nektar::SolverUtils::DiffusionLFRNS::m_dGR_xi1
Array< OneD, Array< OneD, NekDouble > > m_dGR_xi1
Definition: DiffusionLFRNS.h:64
ASSERTL0
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:161
Nektar::SolverUtils::DiffusionLFRNS::m_dGR_xi2
Array< OneD, Array< OneD, NekDouble > > m_dGR_xi2
Definition: DiffusionLFRNS.h:66
Nektar::SolverUtils::DiffusionLFRNS::m_traceNormals
Array< OneD, Array< OneD, NekDouble > > m_traceNormals
Definition: DiffusionLFRNS.h:76
Nektar::SolverUtils::DiffusionLFRNS::m_Q2D_e0
Array< OneD, Array< OneD, NekDouble > > m_Q2D_e0
Definition: DiffusionLFRNS.h:58
Nektar::SolverUtils::DiffusionLFRNS::m_dGL_xi1
Array< OneD, Array< OneD, NekDouble > > m_dGL_xi1
Definition: DiffusionLFRNS.h:63
Nektar::SolverUtils::DiffusionLFRNS::m_Q2D_e3
Array< OneD, Array< OneD, NekDouble > > m_Q2D_e3
Definition: DiffusionLFRNS.h:61
Vmath::Reverse
void Reverse(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.cpp:1071
Nektar::NekDouble
double NekDouble
Definition: NektarUnivTypeDefs.hpp:44
Nektar::SolverUtils::DiffusionLFRNS::m_dGL_xi2
Array< OneD, Array< OneD, NekDouble > > m_dGL_xi2
Definition: DiffusionLFRNS.h:65
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::SolverUtils::DiffusionLFRNS::m_Q2D_e2
Array< OneD, Array< OneD, NekDouble > > m_Q2D_e2
Definition: DiffusionLFRNS.h:60
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
Nektar::SolverUtils::DiffusionLFRNS::m_Q2D_e1
Array< OneD, Array< OneD, NekDouble > > m_Q2D_e1
Definition: DiffusionLFRNS.h:59
void Nektar::SolverUtils::DiffusionLFRNS::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 2158 of file DiffusionLFRNS.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().

2165  {
2166  int n, e, i, j, cnt;
2167 
2168  int nElements = fields[0]->GetExpSize();
2169  int nLocalSolutionPts;
2170  int nEdgePts;
2171  int trace_offset;
2172  int phys_offset;
2173  int nquad0;
2174  int nquad1;
2175 
2176  Array<OneD, NekDouble> auxArray1, auxArray2;
2177  Array<OneD, LibUtilities::BasisSharedPtr> base;
2178 
2179  Array<OneD, Array<OneD, LocalRegions::ExpansionSharedPtr> >
2180  &elmtToTrace = fields[0]->GetTraceMap()->GetElmtToTrace();
2181 
2182  // Loop on the elements
2183  for(n = 0; n < nElements; ++n)
2184  {
2185  // Offset of the element on the global vector
2186  phys_offset = fields[0]->GetPhys_Offset(n);
2187  nLocalSolutionPts = fields[0]->GetExp(n)->GetTotPoints();
2188 
2189  base = fields[0]->GetExp(n)->GetBase();
2190  nquad0 = base[0]->GetNumPoints();
2191  nquad1 = base[1]->GetNumPoints();
2192 
2193  Array<OneD, NekDouble> divCFluxE0(nLocalSolutionPts, 0.0);
2194  Array<OneD, NekDouble> divCFluxE1(nLocalSolutionPts, 0.0);
2195  Array<OneD, NekDouble> divCFluxE2(nLocalSolutionPts, 0.0);
2196  Array<OneD, NekDouble> divCFluxE3(nLocalSolutionPts, 0.0);
2197 
2198  // Loop on the edges
2199  for(e = 0; e < fields[0]->GetExp(n)->GetNedges(); ++e)
2200  {
2201  // Number of edge points of edge e
2202  nEdgePts = fields[0]->GetExp(n)->GetEdgeNumPoints(e);
2203 
2204  Array<OneD, NekDouble> fluxN (nEdgePts, 0.0);
2205  Array<OneD, NekDouble> fluxT (nEdgePts, 0.0);
2206  Array<OneD, NekDouble> fluxN_R (nEdgePts, 0.0);
2207  Array<OneD, NekDouble> fluxN_D (nEdgePts, 0.0);
2208  Array<OneD, NekDouble> fluxJumps (nEdgePts, 0.0);
2209 
2210  // Offset of the trace space correspondent to edge e
2211  trace_offset = fields[0]->GetTrace()->GetPhys_Offset(
2212  elmtToTrace[n][e]->GetElmtId());
2213 
2214  // Get the normals of edge e
2215  const Array<OneD, const Array<OneD, NekDouble> > &normals =
2216  fields[0]->GetExp(n)->GetEdgeNormal(e);
2217 
2218  // Extract the trasformed normal flux at each edge
2219  switch (e)
2220  {
2221  case 0:
2222  // Extract the edge values of transformed flux-y on
2223  // edge e and order them accordingly to the order of
2224  // the trace space
2225  fields[0]->GetExp(n)->GetEdgePhysVals(
2226  e, elmtToTrace[n][e],
2227  fluxX2 + phys_offset,
2228  auxArray1 = fluxN_D);
2229 
2230  Vmath::Neg (nEdgePts, fluxN_D, 1);
2231 
2232  // Extract the physical Rieamann flux at each edge
2233  Vmath::Vcopy(nEdgePts,
2234  &numericalFlux[trace_offset], 1,
2235  &fluxN[0], 1);
2236 
2237  // Check the ordering of vectors
2238  if (fields[0]->GetExp(n)->GetEorient(e) ==
2239  StdRegions::eBackwards)
2240  {
2241  Vmath::Reverse(nEdgePts,
2242  auxArray1 = fluxN, 1,
2243  auxArray2 = fluxN, 1);
2244 
2245  Vmath::Reverse(nEdgePts,
2246  auxArray1 = fluxN_D, 1,
2247  auxArray2 = fluxN_D, 1);
2248  }
2249 
2250  // Transform Riemann Fluxes in the standard element
2251  for (i = 0; i < nquad0; ++i)
2252  {
2253  // Multiply Riemann Flux by Q factors
2254  fluxN_R[i] = (m_Q2D_e0[n][i]) * fluxN[i];
2255  }
2256 
2257  for (i = 0; i < nEdgePts; ++i)
2258  {
2259  if (m_traceNormals[0][trace_offset+i]
2260  != normals[0][i] ||
2261  m_traceNormals[1][trace_offset+i]
2262  != normals[1][i])
2263  {
2264  fluxN_R[i] = -fluxN_R[i];
2265  }
2266  }
2267 
2268  // Subtract to the Riemann flux the discontinuous
2269  // flux
2270  Vmath::Vsub(nEdgePts,
2271  &fluxN_R[0], 1,
2272  &fluxN_D[0], 1, &fluxJumps[0], 1);
2273 
2274  // Multiplicate the flux jumps for the correction
2275  // function
2276  for (i = 0; i < nquad0; ++i)
2277  {
2278  for (j = 0; j < nquad1; ++j)
2279  {
2280  cnt = i + j*nquad0;
2281  divCFluxE0[cnt] =
2282  -fluxJumps[i] * m_dGL_xi2[n][j];
2283  }
2284  }
2285  break;
2286  case 1:
2287  // Extract the edge values of transformed flux-x on
2288  // edge e and order them accordingly to the order of
2289  // the trace space
2290  fields[0]->GetExp(n)->GetEdgePhysVals(
2291  e, elmtToTrace[n][e],
2292  fluxX1 + phys_offset,
2293  auxArray1 = fluxN_D);
2294 
2295  // Extract the physical Rieamann flux at each edge
2296  Vmath::Vcopy(nEdgePts,
2297  &numericalFlux[trace_offset], 1,
2298  &fluxN[0], 1);
2299 
2300  // Check the ordering of vectors
2301  if (fields[0]->GetExp(n)->GetEorient(e) ==
2302  StdRegions::eBackwards)
2303  {
2304  Vmath::Reverse(nEdgePts,
2305  auxArray1 = fluxN, 1,
2306  auxArray2 = fluxN, 1);
2307 
2308  Vmath::Reverse(nEdgePts,
2309  auxArray1 = fluxN_D, 1,
2310  auxArray2 = fluxN_D, 1);
2311  }
2312 
2313  // Transform Riemann Fluxes in the standard element
2314  for (i = 0; i < nquad1; ++i)
2315  {
2316  // Multiply Riemann Flux by Q factors
2317  fluxN_R[i] = (m_Q2D_e1[n][i]) * fluxN[i];
2318  }
2319 
2320  for (i = 0; i < nEdgePts; ++i)
2321  {
2322  if (m_traceNormals[0][trace_offset+i]
2323  != normals[0][i] ||
2324  m_traceNormals[1][trace_offset+i]
2325  != normals[1][i])
2326  {
2327  fluxN_R[i] = -fluxN_R[i];
2328  }
2329  }
2330 
2331  // Subtract to the Riemann flux the discontinuous
2332  // flux
2333  Vmath::Vsub(nEdgePts,
2334  &fluxN_R[0], 1,
2335  &fluxN_D[0], 1, &fluxJumps[0], 1);
2336 
2337  // Multiplicate the flux jumps for the correction
2338  // function
2339  for (i = 0; i < nquad1; ++i)
2340  {
2341  for (j = 0; j < nquad0; ++j)
2342  {
2343  cnt = (nquad0)*i + j;
2344  divCFluxE1[cnt] =
2345  fluxJumps[i] * m_dGR_xi1[n][j];
2346  }
2347  }
2348  break;
2349  case 2:
2350 
2351  // Extract the edge values of transformed flux-y on
2352  // edge e and order them accordingly to the order of
2353  // the trace space
2354 
2355  fields[0]->GetExp(n)->GetEdgePhysVals(
2356  e, elmtToTrace[n][e],
2357  fluxX2 + phys_offset,
2358  auxArray1 = fluxN_D);
2359 
2360  // Extract the physical Rieamann flux at each edge
2361  Vmath::Vcopy(nEdgePts,
2362  &numericalFlux[trace_offset], 1,
2363  &fluxN[0], 1);
2364 
2365  // Check the ordering of vectors
2366  if (fields[0]->GetExp(n)->GetEorient(e) ==
2367  StdRegions::eBackwards)
2368  {
2369  Vmath::Reverse(nEdgePts,
2370  auxArray1 = fluxN, 1,
2371  auxArray2 = fluxN, 1);
2372 
2373  Vmath::Reverse(nEdgePts,
2374  auxArray1 = fluxN_D, 1,
2375  auxArray2 = fluxN_D, 1);
2376  }
2377 
2378  // Transform Riemann Fluxes in the standard element
2379  for (i = 0; i < nquad0; ++i)
2380  {
2381  // Multiply Riemann Flux by Q factors
2382  fluxN_R[i] = (m_Q2D_e2[n][i]) * fluxN[i];
2383  }
2384 
2385  for (i = 0; i < nEdgePts; ++i)
2386  {
2387  if (m_traceNormals[0][trace_offset+i]
2388  != normals[0][i] ||
2389  m_traceNormals[1][trace_offset+i]
2390  != normals[1][i])
2391  {
2392  fluxN_R[i] = -fluxN_R[i];
2393  }
2394  }
2395 
2396  // Subtract to the Riemann flux the discontinuous
2397  // flux
2398 
2399  Vmath::Vsub(nEdgePts,
2400  &fluxN_R[0], 1,
2401  &fluxN_D[0], 1, &fluxJumps[0], 1);
2402 
2403  // Multiplicate the flux jumps for the correction
2404  // function
2405  for (i = 0; i < nquad0; ++i)
2406  {
2407  for (j = 0; j < nquad1; ++j)
2408  {
2409  cnt = (nquad0 - 1) + j*nquad0 - i;
2410  divCFluxE2[cnt] =
2411  fluxJumps[i] * m_dGR_xi2[n][j];
2412  }
2413  }
2414  break;
2415  case 3:
2416  // Extract the edge values of transformed flux-x on
2417  // edge e and order them accordingly to the order of
2418  // the trace space
2419 
2420  fields[0]->GetExp(n)->GetEdgePhysVals(
2421  e, elmtToTrace[n][e],
2422  fluxX1 + phys_offset,
2423  auxArray1 = fluxN_D);
2424  Vmath::Neg (nEdgePts, fluxN_D, 1);
2425 
2426  // Extract the physical Rieamann flux at each edge
2427  Vmath::Vcopy(nEdgePts,
2428  &numericalFlux[trace_offset], 1,
2429  &fluxN[0], 1);
2430 
2431  // Check the ordering of vectors
2432  if (fields[0]->GetExp(n)->GetEorient(e) ==
2433  StdRegions::eBackwards)
2434  {
2435  Vmath::Reverse(nEdgePts,
2436  auxArray1 = fluxN, 1,
2437  auxArray2 = fluxN, 1);
2438 
2439  Vmath::Reverse(nEdgePts,
2440  auxArray1 = fluxN_D, 1,
2441  auxArray2 = fluxN_D, 1);
2442  }
2443 
2444  // Transform Riemann Fluxes in the standard element
2445  for (i = 0; i < nquad1; ++i)
2446  {
2447  // Multiply Riemann Flux by Q factors
2448  fluxN_R[i] = (m_Q2D_e3[n][i]) * fluxN[i];
2449  }
2450 
2451  for (i = 0; i < nEdgePts; ++i)
2452  {
2453  if (m_traceNormals[0][trace_offset+i]
2454  != normals[0][i] ||
2455  m_traceNormals[1][trace_offset+i]
2456  != normals[1][i])
2457  {
2458  fluxN_R[i] = -fluxN_R[i];
2459  }
2460  }
2461 
2462  // Subtract to the Riemann flux the discontinuous
2463  // flux
2464 
2465  Vmath::Vsub(nEdgePts,
2466  &fluxN_R[0], 1,
2467  &fluxN_D[0], 1, &fluxJumps[0], 1);
2468 
2469  // Multiplicate the flux jumps for the correction
2470  // function
2471  for (i = 0; i < nquad1; ++i)
2472  {
2473  for (j = 0; j < nquad0; ++j)
2474  {
2475  cnt = (nquad0*nquad1 - nquad0) + j
2476  - i*nquad0;
2477  divCFluxE3[cnt] =
2478  -fluxJumps[i] * m_dGL_xi1[n][j];
2479  }
2480  }
2481  break;
2482  default:
2483  ASSERTL0(false,"edge value (< 3) is out of range");
2484  break;
2485  }
2486  }
2487 
2488 
2489  // Sum each edge contribution
2490  Vmath::Vadd(nLocalSolutionPts, &divCFluxE0[0], 1,
2491  &divCFluxE1[0], 1, &divCFlux[phys_offset], 1);
2492 
2493  Vmath::Vadd(nLocalSolutionPts, &divCFlux[phys_offset], 1,
2494  &divCFluxE2[0], 1, &divCFlux[phys_offset], 1);
2495 
2496  Vmath::Vadd(nLocalSolutionPts, &divCFlux[phys_offset], 1,
2497  &divCFluxE3[0], 1, &divCFlux[phys_offset], 1);
2498  }
2499  }
Nektar::SolverUtils::DiffusionLFRNS::m_dGR_xi1
Array< OneD, Array< OneD, NekDouble > > m_dGR_xi1
Definition: DiffusionLFRNS.h:64
ASSERTL0
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:161
Nektar::SolverUtils::DiffusionLFRNS::m_dGR_xi2
Array< OneD, Array< OneD, NekDouble > > m_dGR_xi2
Definition: DiffusionLFRNS.h:66
Nektar::SolverUtils::DiffusionLFRNS::m_traceNormals
Array< OneD, Array< OneD, NekDouble > > m_traceNormals
Definition: DiffusionLFRNS.h:76
Nektar::SolverUtils::DiffusionLFRNS::m_Q2D_e0
Array< OneD, Array< OneD, NekDouble > > m_Q2D_e0
Definition: DiffusionLFRNS.h:58
Nektar::SolverUtils::DiffusionLFRNS::m_dGL_xi1
Array< OneD, Array< OneD, NekDouble > > m_dGL_xi1
Definition: DiffusionLFRNS.h:63
Nektar::SolverUtils::DiffusionLFRNS::m_Q2D_e3
Array< OneD, Array< OneD, NekDouble > > m_Q2D_e3
Definition: DiffusionLFRNS.h:61
Vmath::Reverse
void Reverse(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.cpp:1071
Vmath::Neg
void Neg(int n, T *x, const int incx)
Negate x = -x.
Definition: Vmath.cpp:382
Nektar::SolverUtils::DiffusionLFRNS::m_dGL_xi2
Array< OneD, Array< OneD, NekDouble > > m_dGL_xi2
Definition: DiffusionLFRNS.h:65
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::SolverUtils::DiffusionLFRNS::m_Q2D_e2
Array< OneD, Array< OneD, NekDouble > > m_Q2D_e2
Definition: DiffusionLFRNS.h:60
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
Nektar::SolverUtils::DiffusionLFRNS::m_Q2D_e1
Array< OneD, Array< OneD, NekDouble > > m_Q2D_e1
Definition: DiffusionLFRNS.h:59
virtual void Nektar::SolverUtils::DiffusionLFRNS::v_FluxVec ( Array< OneD, Array< OneD, Array< OneD, NekDouble > > > &  fluxvector)
inlineprotectedvirtual

Definition at line 183 of file DiffusionLFRNS.h.

References m_viscTensor.

185  {
186  fluxvector = m_viscTensor;
187  };
Nektar::SolverUtils::DiffusionLFRNS::m_viscTensor
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > m_viscTensor
Definition: DiffusionLFRNS.h:94
virtual Array<OneD, Array<OneD, Array<OneD, NekDouble> > >& Nektar::SolverUtils::DiffusionLFRNS::v_GetFluxTensor ( )
inlineprotectedvirtual

Reimplemented from Nektar::SolverUtils::Diffusion.

Definition at line 195 of file DiffusionLFRNS.h.

References m_viscTensor.

196  {
197  return m_viscTensor;
198  }
Nektar::SolverUtils::DiffusionLFRNS::m_viscTensor
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > m_viscTensor
Definition: DiffusionLFRNS.h:94
void Nektar::SolverUtils::DiffusionLFRNS::v_InitObject ( LibUtilities::SessionReaderSharedPtr  pSession,
Array< OneD, MultiRegions::ExpListSharedPtr >  pFields 
)
protectedvirtual

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

Parameters
pSessionPointer to session reader.
pFieldsPointer to fields.

Reimplemented from Nektar::SolverUtils::Diffusion.

Definition at line 83 of file DiffusionLFRNS.cpp.

References m_BD1, m_D1, m_DD1, m_DFC1, m_DFC2, m_diffDim, m_divFC, m_divFD, m_DU1, m_gamma, m_gasConstant, m_IF1, m_mu, m_pInf, m_rhoInf, m_session, m_spaceDim, m_thermalConductivity, m_tmp1, m_tmp2, m_traceVel, m_Twall, m_viscFlux, m_ViscosityType, m_viscTensor, v_SetupCFunctions(), and v_SetupMetrics().

86  {
87  m_session = pSession;
88  m_session->LoadParameter ("Gamma", m_gamma, 1.4);
89  m_session->LoadParameter ("GasConstant", m_gasConstant, 287.058);
90  m_session->LoadParameter ("Twall", m_Twall, 300.15);
91  m_session->LoadSolverInfo("ViscosityType", m_ViscosityType,
92  "Constant");
93  m_session->LoadParameter ("mu", m_mu, 1.78e-05);
94  m_session->LoadParameter ("thermalConductivity",
95  m_thermalConductivity, 0.0257);
96  m_session->LoadParameter ("rhoInf", m_rhoInf, 1.225);
97  m_session->LoadParameter ("pInf", m_pInf, 101325);
98  v_SetupMetrics(pSession, pFields);
99  v_SetupCFunctions(pSession, pFields);
100 
101  // Initialising arrays
102  int i, j;
103  int nConvectiveFields = pFields.num_elements();
104  int nScalars = nConvectiveFields - 1;
105  int nDim = pFields[0]->GetCoordim(0);
106  int nSolutionPts = pFields[0]->GetTotPoints();
107  int nTracePts = pFields[0]->GetTrace()->GetTotPoints();
108 
109  m_spaceDim = nDim;
110  if (pSession->DefinesSolverInfo("HOMOGENEOUS"))
111  {
112  m_spaceDim = 3;
113  }
114 
115  m_diffDim = m_spaceDim - nDim;
116 
117  m_traceVel = Array<OneD, Array<OneD, NekDouble> >(m_spaceDim);
118 
119  for (i = 0; i < m_spaceDim; ++i)
120  {
121  m_traceVel[i] = Array<OneD, NekDouble> (nTracePts, 0.0);
122  }
123 
124  m_IF1 = Array<OneD, Array<OneD, Array<OneD, NekDouble> > > (
125  nScalars);
126  m_DU1 = Array<OneD, Array<OneD, Array<OneD, NekDouble> > > (
127  nScalars);
128  m_DFC1 = Array<OneD, Array<OneD, Array<OneD, NekDouble> > > (
129  nScalars);
130  m_tmp1 = Array<OneD, Array<OneD, Array<OneD, NekDouble> > > (
131  nScalars);
132  m_tmp2 = Array<OneD, Array<OneD, Array<OneD, NekDouble> > > (
133  nScalars);
134  m_BD1 = Array<OneD, Array<OneD, Array<OneD, NekDouble> > > (
135  nScalars);
136 
137 
138  m_DFC2 = Array<OneD, Array<OneD, Array<OneD, NekDouble> > > (
139  nConvectiveFields);
140  m_DD1 = Array<OneD, Array<OneD, Array<OneD, NekDouble> > > (
141  nConvectiveFields);
142  m_viscFlux = Array<OneD, Array<OneD, NekDouble> > (
143  nConvectiveFields);
144  m_divFD = Array<OneD, Array<OneD, NekDouble> > (nConvectiveFields);
145  m_divFC = Array<OneD, Array<OneD, NekDouble> > (nConvectiveFields);
146 
147  m_D1 = Array<OneD, Array<OneD, Array<OneD, NekDouble> > >
148  (m_spaceDim);
149  m_viscTensor = Array<OneD, Array<OneD, Array<OneD, NekDouble> > >
150  (m_spaceDim);
151  for (i = 0; i < nScalars; ++i)
152  {
153  m_IF1[i] = Array<OneD, Array<OneD, NekDouble> >(m_spaceDim);
154  m_DU1[i] = Array<OneD, Array<OneD, NekDouble> >(nDim);
155  m_DFC1[i] = Array<OneD, Array<OneD, NekDouble> >(nDim);
156  m_tmp1[i] = Array<OneD, Array<OneD, NekDouble> >(nDim);
157  m_tmp2[i] = Array<OneD, Array<OneD, NekDouble> >(nDim);
158  m_BD1[i] = Array<OneD, Array<OneD, NekDouble> >(nDim);
159 
160  for (j = 0; j < nDim; ++j)
161  {
162  m_IF1[i][j] = Array<OneD, NekDouble>(nTracePts, 0.0);
163  m_DU1[i][j] = Array<OneD, NekDouble>(nSolutionPts, 0.0);
164  m_DFC1[i][j] = Array<OneD, NekDouble>(nSolutionPts, 0.0);
165  m_tmp1[i][j] = Array<OneD, NekDouble>(nSolutionPts, 0.0);
166  m_tmp2[i][j] = Array<OneD, NekDouble>(nSolutionPts, 0.0);
167  m_BD1[i][j] = Array<OneD, NekDouble>(nSolutionPts, 0.0);
168  }
169  }
170 
171  for (i = 0; i < nConvectiveFields; ++i)
172  {
173  m_DFC2[i] = Array<OneD, Array<OneD, NekDouble> >(nDim);
174  m_DD1[i] = Array<OneD, Array<OneD, NekDouble> >(nDim);
175  m_viscFlux[i] = Array<OneD, NekDouble>(nTracePts, 0.0);
176  m_divFD[i] = Array<OneD, NekDouble>(nSolutionPts, 0.0);
177  m_divFC[i] = Array<OneD, NekDouble>(nSolutionPts, 0.0);
178 
179  for (j = 0; j < nDim; ++j)
180  {
181  m_DFC2[i][j] = Array<OneD, NekDouble>(nSolutionPts, 0.0);
182  m_DD1[i][j] = Array<OneD, NekDouble>(nSolutionPts, 0.0);
183  }
184  }
185 
186  for (i = 0; i < m_spaceDim; ++i)
187  {
188  m_D1[i] = Array<OneD, Array<OneD, NekDouble> >(nScalars);
189 
190  for (j = 0; j < nScalars; ++j)
191  {
192  m_D1[i][j] = Array<OneD, NekDouble>(nSolutionPts, 0.0);
193  }
194  }
195 
196  for (i = 0; i < m_spaceDim; ++i)
197  {
198  m_viscTensor[i] = Array<OneD, Array<OneD, NekDouble> >(nScalars+1);
199 
200  for (j = 0; j < nScalars+1; ++j)
201  {
202  m_viscTensor[i][j] = Array<OneD, NekDouble>(nSolutionPts, 0.0);
203  }
204  }
205  }
Nektar::SolverUtils::DiffusionLFRNS::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: DiffusionLFRNS.cpp:223
Nektar::SolverUtils::DiffusionLFRNS::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: DiffusionLFRNS.cpp:373
Nektar::SolverUtils::DiffusionLFRNS::m_DU1
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > m_DU1
Definition: DiffusionLFRNS.h:89
Nektar::SolverUtils::DiffusionLFRNS::m_BD1
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > m_BD1
Definition: DiffusionLFRNS.h:91
Nektar::SolverUtils::DiffusionLFRNS::m_session
LibUtilities::SessionReaderSharedPtr m_session
Definition: DiffusionLFRNS.h:77
Nektar::SolverUtils::DiffusionLFRNS::m_tmp2
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > m_tmp2
Definition: DiffusionLFRNS.h:101
Nektar::SolverUtils::DiffusionLFRNS::m_tmp1
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > m_tmp1
Definition: DiffusionLFRNS.h:100
Nektar::SolverUtils::DiffusionLFRNS::m_divFD
Array< OneD, Array< OneD, NekDouble > > m_divFD
Definition: DiffusionLFRNS.h:97
Nektar::SolverUtils::DiffusionLFRNS::m_IF1
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > m_IF1
Definition: DiffusionLFRNS.h:88
Nektar::SolverUtils::DiffusionLFRNS::m_DD1
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > m_DD1
Definition: DiffusionLFRNS.h:93
Nektar::SolverUtils::DiffusionLFRNS::m_DFC1
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > m_DFC1
Definition: DiffusionLFRNS.h:90
Nektar::SolverUtils::DiffusionLFRNS::m_viscTensor
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > m_viscTensor
Definition: DiffusionLFRNS.h:94
Nektar::SolverUtils::DiffusionLFRNS::m_D1
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > m_D1
Definition: DiffusionLFRNS.h:92
Nektar::SolverUtils::DiffusionLFRNS::m_DFC2
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > m_DFC2
Definition: DiffusionLFRNS.h:96
Nektar::SolverUtils::DiffusionLFRNS::m_traceVel
Array< OneD, Array< OneD, NekDouble > > m_traceVel
Definition: DiffusionLFRNS.h:75
Nektar::SolverUtils::DiffusionLFRNS::m_viscFlux
Array< OneD, Array< OneD, NekDouble > > m_viscFlux
Definition: DiffusionLFRNS.h:95
Nektar::SolverUtils::DiffusionLFRNS::m_divFC
Array< OneD, Array< OneD, NekDouble > > m_divFC
Definition: DiffusionLFRNS.h:98
void Nektar::SolverUtils::DiffusionLFRNS::v_NumericalFluxO1 ( const Array< OneD, MultiRegions::ExpListSharedPtr > &  fields,
const Array< OneD, Array< OneD, NekDouble > > &  inarray,
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > &  numericalFluxO1 
)
protectedvirtual

Builds the numerical flux for the 1st order derivatives.

Definition at line 1233 of file DiffusionLFRNS.cpp.

References m_traceNormals, m_traceVel, v_WeakPenaltyO1(), Vmath::Vcopy(), and Vmath::Vvtvp().

Referenced by v_Diffuse().

1238  {
1239  int i, j;
1240  int nTracePts = fields[0]->GetTrace()->GetTotPoints();
1241  int nScalars = inarray.num_elements();
1242  int nDim = fields[0]->GetCoordim(0);
1243 
1244  Array<OneD, NekDouble > Vn (nTracePts, 0.0);
1245  Array<OneD, NekDouble > fluxtemp(nTracePts, 0.0);
1246 
1247  // Get the normal velocity Vn
1248  for (i = 0; i < nDim; ++i)
1249  {
1250  fields[0]->ExtractTracePhys(inarray[i], m_traceVel[i]);
1251  Vmath::Vvtvp(nTracePts, m_traceNormals[i], 1,
1252  m_traceVel[i], 1, Vn, 1, Vn, 1);
1253  }
1254 
1255  // Store forwards/backwards space along trace space
1256  Array<OneD, Array<OneD, NekDouble> > Fwd (nScalars);
1257  Array<OneD, Array<OneD, NekDouble> > Bwd (nScalars);
1258  Array<OneD, Array<OneD, NekDouble> > numflux(nScalars);
1259 
1260  for (i = 0; i < nScalars; ++i)
1261  {
1262  Fwd[i] = Array<OneD, NekDouble>(nTracePts);
1263  Bwd[i] = Array<OneD, NekDouble>(nTracePts);
1264  numflux[i] = Array<OneD, NekDouble>(nTracePts);
1265  fields[i]->GetFwdBwdTracePhys(inarray[i], Fwd[i], Bwd[i]);
1266  fields[0]->GetTrace()->Upwind(Vn, Fwd[i], Bwd[i], numflux[i]);
1267  }
1268 
1269  // Modify the values in case of boundary interfaces
1270  if (fields[0]->GetBndCondExpansions().num_elements())
1271  {
1272  v_WeakPenaltyO1(fields, inarray, numflux);
1273  }
1274 
1275  // Splitting the numerical flux into the dimensions
1276  for (j = 0; j < nDim; ++j)
1277  {
1278  for (i = 0; i < nScalars; ++i)
1279  {
1280  Vmath::Vcopy(nTracePts,numflux[i], 1,
1281  numericalFluxO1[i][j], 1);
1282  }
1283  }
1284  }
Nektar::SolverUtils::DiffusionLFRNS::m_traceNormals
Array< OneD, Array< OneD, NekDouble > > m_traceNormals
Definition: DiffusionLFRNS.h:76
Nektar::SolverUtils::DiffusionLFRNS::v_WeakPenaltyO1
virtual void v_WeakPenaltyO1(const Array< OneD, MultiRegions::ExpListSharedPtr > &fields, const Array< OneD, Array< OneD, NekDouble > > &inarray, Array< OneD, Array< OneD, NekDouble > > &penaltyfluxO1)
Imposes appropriate bcs for the 1st order derivatives.
Definition: DiffusionLFRNS.cpp:1290
Vmath::Vvtvp
void Vvtvp(int n, const T *w, const int incw, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
vvtvp (vector times vector plus vector): z = w*x + y
Definition: Vmath.cpp:428
Nektar::SolverUtils::DiffusionLFRNS::m_traceVel
Array< OneD, Array< OneD, NekDouble > > m_traceVel
Definition: DiffusionLFRNS.h:75
Vmath::Vcopy
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.cpp:1047
void Nektar::SolverUtils::DiffusionLFRNS::v_NumericalFluxO2 ( 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 1498 of file DiffusionLFRNS.cpp.

References m_traceNormals, m_traceVel, v_WeakPenaltyO2(), Vmath::Vadd(), Vmath::Vmul(), and Vmath::Vvtvp().

Referenced by v_Diffuse().

1503  {
1504  int i, j;
1505  int nTracePts = fields[0]->GetTrace()->GetTotPoints();
1506  int nVariables = fields.num_elements();
1507  int nDim = fields[0]->GetCoordim(0);
1508 
1509  Array<OneD, NekDouble > Fwd(nTracePts);
1510  Array<OneD, NekDouble > Bwd(nTracePts);
1511  Array<OneD, NekDouble > Vn (nTracePts, 0.0);
1512 
1513  Array<OneD, NekDouble > qFwd (nTracePts);
1514  Array<OneD, NekDouble > qBwd (nTracePts);
1515  Array<OneD, NekDouble > qfluxtemp(nTracePts, 0.0);
1516 
1517  // Get the normal velocity Vn
1518  for (i = 0; i < nDim; ++i)
1519  {
1520  fields[0]->ExtractTracePhys(ufield[i], m_traceVel[i]);
1521  Vmath::Vvtvp(nTracePts, m_traceNormals[i], 1,
1522  m_traceVel[i], 1, Vn, 1, Vn, 1);
1523  }
1524 
1525  // Evaulate Riemann flux
1526  // qflux = \hat{q} \cdot u = q \cdot n
1527  // Notice: i = 1 (first row of the viscous tensor is zero)
1528  for (i = 1; i < nVariables; ++i)
1529  {
1530  qflux[i] = Array<OneD, NekDouble> (nTracePts, 0.0);
1531  for (j = 0; j < nDim; ++j)
1532  {
1533  // Compute qFwd and qBwd value of qfield in position 'ji'
1534  fields[i]->GetFwdBwdTracePhys(qfield[j][i], qFwd, qBwd);
1535 
1536  // Get Riemann flux of qflux --> LDG implies upwind
1537  fields[i]->GetTrace()->Upwind(Vn, qBwd, qFwd, qfluxtemp);
1538 
1539  // Multiply the Riemann flux by the trace normals
1540  Vmath::Vmul(nTracePts, m_traceNormals[j], 1, qfluxtemp, 1,
1541  qfluxtemp, 1);
1542 
1543  // Impose weak boundary condition with flux
1544  if (fields[0]->GetBndCondExpansions().num_elements())
1545  {
1546  v_WeakPenaltyO2(fields, i, j, qfield[j][i], qfluxtemp);
1547  }
1548 
1549  // Store the final flux into qflux
1550  Vmath::Vadd(nTracePts, qfluxtemp, 1, qflux[i], 1,
1551  qflux[i], 1);
1552  }
1553  }
1554  }
Nektar::SolverUtils::DiffusionLFRNS::m_traceNormals
Array< OneD, Array< OneD, NekDouble > > m_traceNormals
Definition: DiffusionLFRNS.h:76
Vmath::Vvtvp
void Vvtvp(int n, const T *w, const int incw, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
vvtvp (vector times vector plus vector): z = w*x + y
Definition: Vmath.cpp:428
Nektar::SolverUtils::DiffusionLFRNS::v_WeakPenaltyO2
virtual void v_WeakPenaltyO2(const Array< OneD, MultiRegions::ExpListSharedPtr > &fields, const int var, const int dir, const Array< OneD, const NekDouble > &qfield, Array< OneD, NekDouble > &penaltyflux)
Imposes appropriate bcs for the 2nd order derivatives.
Definition: DiffusionLFRNS.cpp:1561
Nektar::SolverUtils::DiffusionLFRNS::m_traceVel
Array< OneD, Array< OneD, NekDouble > > m_traceVel
Definition: DiffusionLFRNS.h:75
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
virtual void Nektar::SolverUtils::DiffusionLFRNS::v_SetHomoDerivs ( Array< OneD, Array< OneD, NekDouble > > &  deriv)
inlineprotectedvirtual

Reimplemented from Nektar::SolverUtils::Diffusion.

Definition at line 189 of file DiffusionLFRNS.h.

191  {
192  m_homoDerivs = deriv;
193  }
Nektar::SolverUtils::DiffusionLFRNS::m_homoDerivs
Array< OneD, Array< OneD, NekDouble > > m_homoDerivs
Definition: DiffusionLFRNS.h:103
void Nektar::SolverUtils::DiffusionLFRNS::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 to compute the corrective fluxes.

Parameters
pSessionPointer to session reader.
pFieldsPointer to fields.

Definition at line 373 of file DiffusionLFRNS.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().

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

226  {
227  int i, n;
228  int nquad0, nquad1;
229  int phys_offset;
230  int nLocalSolutionPts;
231  int nElements = pFields[0]->GetExpSize();
232  int nDimensions = pFields[0]->GetCoordim(0);
233  int nSolutionPts = pFields[0]->GetTotPoints();
234  int nTracePts = pFields[0]->GetTrace()->GetTotPoints();
235 
236  m_traceNormals = Array<OneD, Array<OneD, NekDouble> >(nDimensions);
237  for (i = 0; i < nDimensions; ++i)
238  {
239  m_traceNormals[i] = Array<OneD, NekDouble> (nTracePts, 0.0);
240  }
241  pFields[0]->GetTrace()->GetNormals(m_traceNormals);
242 
243  Array<TwoD, const NekDouble> gmat;
244  Array<OneD, const NekDouble> jac;
245 
246  m_jac = Array<OneD, NekDouble>(nSolutionPts);
247 
248  Array<OneD, NekDouble> auxArray1;
249  Array<OneD, LibUtilities::BasisSharedPtr> base;
250  LibUtilities::PointsKeyVector ptsKeys;
251 
252  switch (nDimensions)
253  {
254  case 1:
255  {
256  for (n = 0; n < nElements; ++n)
257  {
258  ptsKeys = pFields[0]->GetExp(n)->GetPointsKeys();
259  nLocalSolutionPts = pFields[0]->GetExp(n)->GetTotPoints();
260  phys_offset = pFields[0]->GetPhys_Offset(n);
261  jac = pFields[0]->GetExp(n)
262  ->as<LocalRegions::Expansion1D>()->GetGeom1D()
263  ->GetMetricInfo()->GetJac(ptsKeys);
264  for (i = 0; i < nLocalSolutionPts; ++i)
265  {
266  m_jac[i+phys_offset] = jac[0];
267  }
268  }
269  break;
270  }
271  case 2:
272  {
273  m_gmat = Array<OneD, Array<OneD, NekDouble> >(4);
274  m_gmat[0] = Array<OneD, NekDouble>(nSolutionPts);
275  m_gmat[1] = Array<OneD, NekDouble>(nSolutionPts);
276  m_gmat[2] = Array<OneD, NekDouble>(nSolutionPts);
277  m_gmat[3] = Array<OneD, NekDouble>(nSolutionPts);
278 
279  m_Q2D_e0 = Array<OneD, Array<OneD, NekDouble> >(nElements);
280  m_Q2D_e1 = Array<OneD, Array<OneD, NekDouble> >(nElements);
281  m_Q2D_e2 = Array<OneD, Array<OneD, NekDouble> >(nElements);
282  m_Q2D_e3 = Array<OneD, Array<OneD, NekDouble> >(nElements);
283 
284  for (n = 0; n < nElements; ++n)
285  {
286  base = pFields[0]->GetExp(n)->GetBase();
287  nquad0 = base[0]->GetNumPoints();
288  nquad1 = base[1]->GetNumPoints();
289 
290  m_Q2D_e0[n] = Array<OneD, NekDouble>(nquad0);
291  m_Q2D_e1[n] = Array<OneD, NekDouble>(nquad1);
292  m_Q2D_e2[n] = Array<OneD, NekDouble>(nquad0);
293  m_Q2D_e3[n] = Array<OneD, NekDouble>(nquad1);
294 
295  // Extract the Q factors at each edge point
296  pFields[0]->GetExp(n)->GetEdgeQFactors(
297  0, auxArray1 = m_Q2D_e0[n]);
298  pFields[0]->GetExp(n)->GetEdgeQFactors(
299  1, auxArray1 = m_Q2D_e1[n]);
300  pFields[0]->GetExp(n)->GetEdgeQFactors(
301  2, auxArray1 = m_Q2D_e2[n]);
302  pFields[0]->GetExp(n)->GetEdgeQFactors(
303  3, auxArray1 = m_Q2D_e3[n]);
304 
305  ptsKeys = pFields[0]->GetExp(n)->GetPointsKeys();
306  nLocalSolutionPts = pFields[0]->GetExp(n)->GetTotPoints();
307  phys_offset = pFields[0]->GetPhys_Offset(n);
308 
309  jac = pFields[0]->GetExp(n)
310  ->as<LocalRegions::Expansion2D>()->GetGeom2D()
311  ->GetMetricInfo()->GetJac(ptsKeys);
312  gmat = pFields[0]->GetExp(n)
313  ->as<LocalRegions::Expansion2D>()->GetGeom2D()
314  ->GetMetricInfo()->GetDerivFactors(ptsKeys);
315 
316  if (pFields[0]->GetExp(n)
317  ->as<LocalRegions::Expansion2D>()->GetGeom2D()
318  ->GetMetricInfo()->GetGtype()
319  == SpatialDomains::eDeformed)
320  {
321  for (i = 0; i < nLocalSolutionPts; ++i)
322  {
323  m_jac[i+phys_offset] = jac[i];
324  m_gmat[0][i+phys_offset] = gmat[0][i];
325  m_gmat[1][i+phys_offset] = gmat[1][i];
326  m_gmat[2][i+phys_offset] = gmat[2][i];
327  m_gmat[3][i+phys_offset] = gmat[3][i];
328  }
329  }
330  else
331  {
332  for (i = 0; i < nLocalSolutionPts; ++i)
333  {
334  m_jac[i+phys_offset] = jac[0];
335  m_gmat[0][i+phys_offset] = gmat[0][0];
336  m_gmat[1][i+phys_offset] = gmat[1][0];
337  m_gmat[2][i+phys_offset] = gmat[2][0];
338  m_gmat[3][i+phys_offset] = gmat[3][0];
339  }
340  }
341  }
342  break;
343  }
344  case 3:
345  {
346  ASSERTL0(false,"3DFR Metric terms not implemented yet");
347  break;
348  }
349  default:
350  {
351  ASSERTL0(false, "Expansion dimension not recognised");
352  break;
353  }
354  }
355  }
ASSERTL0
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:161
Nektar::LibUtilities::PointsKeyVector
std::vector< PointsKey > PointsKeyVector
Definition: Points.h:220
Nektar::SolverUtils::DiffusionLFRNS::m_traceNormals
Array< OneD, Array< OneD, NekDouble > > m_traceNormals
Definition: DiffusionLFRNS.h:76
Nektar::SolverUtils::DiffusionLFRNS::m_Q2D_e0
Array< OneD, Array< OneD, NekDouble > > m_Q2D_e0
Definition: DiffusionLFRNS.h:58
Nektar::SolverUtils::DiffusionLFRNS::m_jac
Array< OneD, NekDouble > m_jac
Definition: DiffusionLFRNS.h:55
Nektar::SolverUtils::DiffusionLFRNS::m_Q2D_e3
Array< OneD, Array< OneD, NekDouble > > m_Q2D_e3
Definition: DiffusionLFRNS.h:61
Nektar::SolverUtils::DiffusionLFRNS::m_Q2D_e2
Array< OneD, Array< OneD, NekDouble > > m_Q2D_e2
Definition: DiffusionLFRNS.h:60
Nektar::SolverUtils::DiffusionLFRNS::m_gmat
Array< OneD, Array< OneD, NekDouble > > m_gmat
Definition: DiffusionLFRNS.h:56
Nektar::StdRegions::StdExpansion::GetMetricInfo
const boost::shared_ptr< SpatialDomains::GeomFactors > & GetMetricInfo(void) const
Definition: StdExpansion.h:1179
Nektar::SpatialDomains::eDeformed
Geometry is curved or has non-constant factors.
Definition: SpatialDomains.hpp:57
Nektar::SolverUtils::DiffusionLFRNS::m_Q2D_e1
Array< OneD, Array< OneD, NekDouble > > m_Q2D_e1
Definition: DiffusionLFRNS.h:59
void Nektar::SolverUtils::DiffusionLFRNS::v_WeakPenaltyO1 ( const Array< OneD, MultiRegions::ExpListSharedPtr > &  fields,
const Array< OneD, Array< OneD, NekDouble > > &  inarray,
Array< OneD, Array< OneD, NekDouble > > &  penaltyfluxO1 
)
protectedvirtual

Imposes appropriate bcs for the 1st order derivatives.

Definition at line 1290 of file DiffusionLFRNS.cpp.

References Nektar::SpatialDomains::eDirichlet, Nektar::SpatialDomains::eNeumann, m_gamma, m_gasConstant, m_Twall, Vmath::Smul(), Vmath::Vadd(), Vmath::Vcopy(), Vmath::Vdiv(), Vmath::Vmul(), Vmath::Vsub(), and Vmath::Zero().

Referenced by v_NumericalFluxO1().

1294  {
1295  int cnt;
1296  int i, j, e;
1297  int id1, id2;
1298 
1299  int nBndEdgePts, nBndEdges, nBndRegions;
1300 
1301  int nTracePts = fields[0]->GetTrace()->GetTotPoints();
1302  int nScalars = inarray.num_elements();
1303 
1304  Array<OneD, NekDouble> tmp1(nTracePts, 0.0);
1305  Array<OneD, NekDouble> tmp2(nTracePts, 0.0);
1306  Array<OneD, NekDouble> Tw(nTracePts, m_Twall);
1307 
1308  Array< OneD, Array<OneD, NekDouble > > scalarVariables(nScalars);
1309  Array< OneD, Array<OneD, NekDouble > > uplus(nScalars);
1310 
1311  // Extract internal values of the scalar variables for Neumann bcs
1312  for (i = 0; i < nScalars; ++i)
1313  {
1314  scalarVariables[i] = Array<OneD, NekDouble>(nTracePts, 0.0);
1315 
1316  uplus[i] = Array<OneD, NekDouble>(nTracePts, 0.0);
1317  fields[i]->ExtractTracePhys(inarray[i], uplus[i]);
1318  }
1319 
1320  // Compute boundary conditions for velocities
1321  for (i = 0; i < nScalars-1; ++i)
1322  {
1323  // Note that cnt has to loop on nBndRegions and nBndEdges
1324  // and has to be reset to zero per each equation
1325  cnt = 0;
1326  nBndRegions = fields[i+1]->
1327  GetBndCondExpansions().num_elements();
1328  for (j = 0; j < nBndRegions; ++j)
1329  {
1330  nBndEdges = fields[i+1]->
1331  GetBndCondExpansions()[j]->GetExpSize();
1332  for (e = 0; e < nBndEdges; ++e)
1333  {
1334  nBndEdgePts = fields[i+1]->
1335  GetBndCondExpansions()[j]->GetExp(e)->GetTotPoints();
1336 
1337  id1 = fields[i+1]->
1338  GetBndCondExpansions()[j]->GetPhys_Offset(e);
1339 
1340  id2 = fields[0]->GetTrace()->
1341  GetPhys_Offset(fields[0]->GetTraceMap()->
1342  GetBndCondTraceToGlobalTraceMap(cnt++));
1343 
1344  // Reinforcing bcs for velocity in case of Wall bcs
1345  if (boost::iequals(fields[i]->GetBndConditions()[j]->
1346  GetUserDefined(),"WallViscous") ||
1347  boost::iequals(fields[i]->GetBndConditions()[j]->
1348  GetUserDefined(),"WallAdiabatic"))
1349  {
1350  Vmath::Zero(nBndEdgePts,
1351  &scalarVariables[i][id2], 1);
1352  }
1353 
1354  // Imposing velocity bcs if not Wall
1355  else if (fields[i]->GetBndConditions()[j]->
1356  GetBoundaryConditionType() ==
1357  SpatialDomains::eDirichlet)
1358  {
1359  Vmath::Vdiv(nBndEdgePts,
1360  &(fields[i+1]->
1361  GetBndCondExpansions()[j]->
1362  UpdatePhys())[id1], 1,
1363  &(fields[0]->
1364  GetBndCondExpansions()[j]->
1365  UpdatePhys())[id1], 1,
1366  &scalarVariables[i][id2], 1);
1367  }
1368 
1369  // For Dirichlet boundary condition: uflux = u_bcs
1370  if (fields[i]->GetBndConditions()[j]->
1371  GetBoundaryConditionType() ==
1372  SpatialDomains::eDirichlet)
1373  {
1374  Vmath::Vcopy(nBndEdgePts,
1375  &scalarVariables[i][id2], 1,
1376  &penaltyfluxO1[i][id2], 1);
1377  }
1378 
1379  // For Neumann boundary condition: uflux = u_+
1380  else if ((fields[i]->GetBndConditions()[j])->
1381  GetBoundaryConditionType() ==
1382  SpatialDomains::eNeumann)
1383  {
1384  Vmath::Vcopy(nBndEdgePts,
1385  &uplus[i][id2], 1,
1386  &penaltyfluxO1[i][id2], 1);
1387  }
1388 
1389  // Building kinetic energy to be used for T bcs
1390  Vmath::Vmul(nBndEdgePts,
1391  &scalarVariables[i][id2], 1,
1392  &scalarVariables[i][id2], 1,
1393  &tmp1[id2], 1);
1394 
1395  Vmath::Smul(nBndEdgePts, 0.5,
1396  &tmp1[id2], 1,
1397  &tmp1[id2], 1);
1398 
1399  Vmath::Vadd(nBndEdgePts,
1400  &tmp2[id2], 1,
1401  &tmp1[id2], 1,
1402  &tmp2[id2], 1);
1403  }
1404  }
1405  }
1406 
1407  // Compute boundary conditions for temperature
1408  cnt = 0;
1409  nBndRegions = fields[nScalars]->
1410  GetBndCondExpansions().num_elements();
1411  for (j = 0; j < nBndRegions; ++j)
1412  {
1413  nBndEdges = fields[nScalars]->
1414  GetBndCondExpansions()[j]->GetExpSize();
1415  for (e = 0; e < nBndEdges; ++e)
1416  {
1417  nBndEdgePts = fields[nScalars]->
1418  GetBndCondExpansions()[j]->GetExp(e)->GetTotPoints();
1419 
1420  id1 = fields[nScalars]->
1421  GetBndCondExpansions()[j]->GetPhys_Offset(e);
1422 
1423  id2 = fields[0]->GetTrace()->
1424  GetPhys_Offset(fields[0]->GetTraceMap()->
1425  GetBndCondTraceToGlobalTraceMap(cnt++));
1426 
1427  // Imposing Temperature Twall at the wall
1428  if (boost::iequals(fields[i]->GetBndConditions()[j]->
1429  GetUserDefined(),"WallViscous"))
1430  {
1431  Vmath::Vcopy(nBndEdgePts,
1432  &Tw[0], 1,
1433  &scalarVariables[nScalars-1][id2], 1);
1434  }
1435  // Imposing Temperature through condition on the Energy
1436  // for no wall boundaries (e.g. farfield)
1437  else if (fields[i]->GetBndConditions()[j]->
1438  GetBoundaryConditionType() ==
1439  SpatialDomains::eDirichlet)
1440  {
1441  // Divide E by rho
1442  Vmath::Vdiv(nBndEdgePts,
1443  &(fields[nScalars]->
1444  GetBndCondExpansions()[j]->
1445  GetPhys())[id1], 1,
1446  &(fields[0]->
1447  GetBndCondExpansions()[j]->
1448  GetPhys())[id1], 1,
1449  &scalarVariables[nScalars-1][id2], 1);
1450 
1451  // Subtract kinetic energy to E/rho
1452  Vmath::Vsub(nBndEdgePts,
1453  &scalarVariables[nScalars-1][id2], 1,
1454  &tmp2[id2], 1,
1455  &scalarVariables[nScalars-1][id2], 1);
1456 
1457  // Multiply by constant factor (gamma-1)/R
1458  Vmath::Smul(nBndEdgePts, (m_gamma - 1)/m_gasConstant,
1459  &scalarVariables[nScalars-1][id2], 1,
1460  &scalarVariables[nScalars-1][id2], 1);
1461  }
1462 
1463  // For Dirichlet boundary condition: uflux = u_bcs
1464  if (fields[nScalars]->GetBndConditions()[j]->
1465  GetBoundaryConditionType() ==
1466  SpatialDomains::eDirichlet &&
1467  !boost::iequals(
1468  fields[nScalars]->GetBndConditions()[j]
1469  ->GetUserDefined(), "WallAdiabatic"))
1470  {
1471  Vmath::Vcopy(nBndEdgePts,
1472  &scalarVariables[nScalars-1][id2], 1,
1473  &penaltyfluxO1[nScalars-1][id2], 1);
1474 
1475  }
1476 
1477  // For Neumann boundary condition: uflux = u_+
1478  else if (((fields[nScalars]->GetBndConditions()[j])->
1479  GetBoundaryConditionType() ==
1480  SpatialDomains::eNeumann) ||
1481  boost::iequals(
1482  fields[nScalars]->GetBndConditions()[j]
1483  ->GetUserDefined(), "WallAdiabatic"))
1484  {
1485  Vmath::Vcopy(nBndEdgePts,
1486  &uplus[nScalars-1][id2], 1,
1487  &penaltyfluxO1[nScalars-1][id2], 1);
1488 
1489  }
1490  }
1491  }
1492  }
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
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
Vmath::Zero
void Zero(int n, T *x, const int incx)
Zero vector.
Definition: Vmath.cpp:359
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
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::DiffusionLFRNS::v_WeakPenaltyO2 ( const Array< OneD, MultiRegions::ExpListSharedPtr > &  fields,
const int  var,
const int  dir,
const Array< OneD, const NekDouble > &  qfield,
Array< OneD, NekDouble > &  penaltyflux 
)
protectedvirtual

Imposes appropriate bcs for the 2nd order derivatives.

Definition at line 1561 of file DiffusionLFRNS.cpp.

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

Referenced by v_NumericalFluxO2().

1567  {
1568  int cnt = 0;
1569  int nBndEdges, nBndEdgePts;
1570  int i, e;
1571  int id2;
1572 
1573  int nTracePts = fields[0]->GetTrace()->GetTotPoints();
1574  int nBndRegions = fields[var]->GetBndCondExpansions().num_elements();
1575 
1576  Array<OneD, NekDouble > uterm(nTracePts);
1577  Array<OneD, NekDouble > qtemp(nTracePts);
1578 
1579  // Extract the physical values of the solution at the boundaries
1580  fields[var]->ExtractTracePhys(qfield, qtemp);
1581 
1582  // Loop on the boundary regions to apply appropriate bcs
1583  for (i = 0; i < nBndRegions; ++i)
1584  {
1585  // Number of boundary regions related to region 'i'
1586  nBndEdges = fields[var]->
1587  GetBndCondExpansions()[i]->GetExpSize();
1588 
1589  // Weakly impose bcs by modifying flux values
1590  for (e = 0; e < nBndEdges; ++e)
1591  {
1592  nBndEdgePts = fields[var]->
1593  GetBndCondExpansions()[i]->GetExp(e)->GetTotPoints();
1594 
1595  id2 = fields[0]->GetTrace()->
1596  GetPhys_Offset(fields[0]->GetTraceMap()->
1597  GetBndCondTraceToGlobalTraceMap(cnt++));
1598 
1599 
1600  // In case of Dirichlet bcs:
1601  // uflux = gD
1602  if(fields[var]->GetBndConditions()[i]->
1603  GetBoundaryConditionType() == SpatialDomains::eDirichlet
1604  && !boost::iequals(fields[var]->GetBndConditions()[i]
1605  ->GetUserDefined(), "WallAdiabatic"))
1606  {
1607  Vmath::Vmul(nBndEdgePts,
1608  &m_traceNormals[dir][id2], 1,
1609  &qtemp[id2], 1,
1610  &penaltyflux[id2], 1);
1611  }
1612  // 3.4) In case of Neumann bcs:
1613  // uflux = u+
1614  else if((fields[var]->GetBndConditions()[i])->
1615  GetBoundaryConditionType() == SpatialDomains::eNeumann)
1616  {
1617  ASSERTL0(false,
1618  "Neumann bcs not implemented for LFRNS");
1619  }
1620  else if(boost::iequals(fields[var]->GetBndConditions()[i]
1621  ->GetUserDefined(), "WallAdiabatic"))
1622  {
1623  if ((var == m_spaceDim + 1))
1624  {
1625  Vmath::Zero(nBndEdgePts, &penaltyflux[id2], 1);
1626  }
1627  else
1628  {
1629 
1630  Vmath::Vmul(nBndEdgePts,
1631  &m_traceNormals[dir][id2], 1,
1632  &qtemp[id2], 1,
1633  &penaltyflux[id2], 1);
1634  }
1635 
1636  }
1637  }
1638  }
1639  }
ASSERTL0
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:161
Nektar::SolverUtils::DiffusionLFRNS::m_traceNormals
Array< OneD, Array< OneD, NekDouble > > m_traceNormals
Definition: DiffusionLFRNS.h:76
Vmath::Zero
void Zero(int n, T *x, const int incx)
Zero vector.
Definition: Vmath.cpp:359
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::DiffusionLFRNS::m_BD1
protected

Definition at line 91 of file DiffusionLFRNS.h.

Referenced by v_Diffuse(), and v_InitObject().

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

Definition at line 92 of file DiffusionLFRNS.h.

Referenced by v_Diffuse(), and v_InitObject().

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

Definition at line 93 of file DiffusionLFRNS.h.

Referenced by v_Diffuse(), and v_InitObject().

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

Definition at line 90 of file DiffusionLFRNS.h.

Referenced by v_Diffuse(), and v_InitObject().

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

Definition at line 96 of file DiffusionLFRNS.h.

Referenced by v_Diffuse(), and v_InitObject().

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

Definition at line 63 of file DiffusionLFRNS.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::DiffusionLFRNS::m_dGL_xi2

Definition at line 65 of file DiffusionLFRNS.h.

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

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

Definition at line 67 of file DiffusionLFRNS.h.

Referenced by v_SetupCFunctions().

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

Definition at line 64 of file DiffusionLFRNS.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::DiffusionLFRNS::m_dGR_xi2

Definition at line 66 of file DiffusionLFRNS.h.

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

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

Definition at line 68 of file DiffusionLFRNS.h.

Referenced by v_SetupCFunctions().

int Nektar::SolverUtils::DiffusionLFRNS::m_diffDim
protected

Definition at line 106 of file DiffusionLFRNS.h.

Referenced by v_Diffuse(), and v_InitObject().

std::string Nektar::SolverUtils::DiffusionLFRNS::m_diffType
protected

Definition at line 108 of file DiffusionLFRNS.h.

Referenced by v_SetupCFunctions().

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

Definition at line 98 of file DiffusionLFRNS.h.

Referenced by v_Diffuse(), and v_InitObject().

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

Definition at line 97 of file DiffusionLFRNS.h.

Referenced by v_Diffuse(), and v_InitObject().

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

Definition at line 89 of file DiffusionLFRNS.h.

Referenced by v_Diffuse(), and v_InitObject().

NekDouble Nektar::SolverUtils::DiffusionLFRNS::m_gamma
protected

Definition at line 78 of file DiffusionLFRNS.h.

Referenced by v_InitObject(), and v_WeakPenaltyO1().

NekDouble Nektar::SolverUtils::DiffusionLFRNS::m_gasConstant
protected

Definition at line 79 of file DiffusionLFRNS.h.

Referenced by v_InitObject(), and v_WeakPenaltyO1().

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

Definition at line 56 of file DiffusionLFRNS.h.

Referenced by v_Diffuse(), and v_SetupMetrics().

Array<OneD, Array<OneD, NekDouble> > Nektar::SolverUtils::DiffusionLFRNS::m_homoDerivs
protected

Definition at line 103 of file DiffusionLFRNS.h.

Referenced by v_Diffuse().

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

Definition at line 88 of file DiffusionLFRNS.h.

Referenced by v_Diffuse(), and v_InitObject().

DNekMatSharedPtr Nektar::SolverUtils::DiffusionLFRNS::m_Ixm

Definition at line 69 of file DiffusionLFRNS.h.

DNekMatSharedPtr Nektar::SolverUtils::DiffusionLFRNS::m_Ixp

Definition at line 70 of file DiffusionLFRNS.h.

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

Definition at line 55 of file DiffusionLFRNS.h.

Referenced by v_Diffuse(), and v_SetupMetrics().

NekDouble Nektar::SolverUtils::DiffusionLFRNS::m_mu
protected

Definition at line 82 of file DiffusionLFRNS.h.

Referenced by v_InitObject().

NekDouble Nektar::SolverUtils::DiffusionLFRNS::m_pInf
protected

Definition at line 85 of file DiffusionLFRNS.h.

Referenced by v_InitObject().

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

Definition at line 58 of file DiffusionLFRNS.h.

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

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

Definition at line 59 of file DiffusionLFRNS.h.

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

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

Definition at line 60 of file DiffusionLFRNS.h.

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

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

Definition at line 61 of file DiffusionLFRNS.h.

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

NekDouble Nektar::SolverUtils::DiffusionLFRNS::m_rhoInf
protected

Definition at line 84 of file DiffusionLFRNS.h.

Referenced by v_InitObject().

LibUtilities::SessionReaderSharedPtr Nektar::SolverUtils::DiffusionLFRNS::m_session
protected

Definition at line 77 of file DiffusionLFRNS.h.

Referenced by v_InitObject().

int Nektar::SolverUtils::DiffusionLFRNS::m_spaceDim
protected

Definition at line 105 of file DiffusionLFRNS.h.

Referenced by v_InitObject(), and v_WeakPenaltyO2().

NekDouble Nektar::SolverUtils::DiffusionLFRNS::m_thermalConductivity
protected

Definition at line 83 of file DiffusionLFRNS.h.

Referenced by v_InitObject().

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

Definition at line 100 of file DiffusionLFRNS.h.

Referenced by v_Diffuse(), and v_InitObject().

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

Definition at line 101 of file DiffusionLFRNS.h.

Referenced by v_Diffuse(), and v_InitObject().

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

Definition at line 76 of file DiffusionLFRNS.h.

Referenced by v_DivCFlux_2D(), v_DivCFlux_2D_Gauss(), v_NumericalFluxO1(), v_NumericalFluxO2(), v_SetupMetrics(), and v_WeakPenaltyO2().

Array<OneD, Array<OneD, NekDouble> > Nektar::SolverUtils::DiffusionLFRNS::m_traceVel
protected

Definition at line 75 of file DiffusionLFRNS.h.

Referenced by v_InitObject(), v_NumericalFluxO1(), and v_NumericalFluxO2().

NekDouble Nektar::SolverUtils::DiffusionLFRNS::m_Twall
protected

Definition at line 80 of file DiffusionLFRNS.h.

Referenced by v_InitObject(), and v_WeakPenaltyO1().

Array<OneD, Array<OneD, NekDouble> > Nektar::SolverUtils::DiffusionLFRNS::m_viscFlux
protected

Definition at line 95 of file DiffusionLFRNS.h.

Referenced by v_Diffuse(), and v_InitObject().

std::string Nektar::SolverUtils::DiffusionLFRNS::m_ViscosityType
protected

Definition at line 81 of file DiffusionLFRNS.h.

Referenced by v_InitObject().

Array<OneD, Array<OneD, Array<OneD, NekDouble> > > Nektar::SolverUtils::DiffusionLFRNS::m_viscTensor
protected

Definition at line 94 of file DiffusionLFRNS.h.

Referenced by v_Diffuse(), v_FluxVec(), v_GetFluxTensor(), and v_InitObject().

std::string Nektar::SolverUtils::DiffusionLFRNS::type
static
Initial value:
= {
GetDiffusionFactory().RegisterCreatorFunction(
"LFRDGNS", DiffusionLFRNS::create),
GetDiffusionFactory().RegisterCreatorFunction(
"LFRSDNS", DiffusionLFRNS::create),
GetDiffusionFactory().RegisterCreatorFunction(
"LFRHUNS", DiffusionLFRNS::create),
GetDiffusionFactory().RegisterCreatorFunction(
"LFRcminNS", DiffusionLFRNS::create),
GetDiffusionFactory().RegisterCreatorFunction(
"LFRcinfNS", DiffusionLFRNS::create)}

Definition at line 53 of file DiffusionLFRNS.h.

Generated on Wed Mar 23 2016 20:14:36 for Nektar++ by   doxygen 1.8.8