Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
Static Public Member Functions | Public Attributes | Static Public Attributes | Protected Member Functions | Protected Attributes | List of all members
Nektar::SolverUtils::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, const Array< OneD, Array< OneD, NekDouble > > &pFwd=NullNekDoubleArrayofArray, const Array< OneD, Array< OneD, NekDouble > > &pBwd=NullNekDoubleArrayofArray)
 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
 
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, const Array< OneD, Array< OneD, NekDouble > > &pFwd=NullNekDoubleArrayofArray, const Array< OneD, Array< OneD, NekDouble > > &pBwd=NullNekDoubleArrayofArray)
 
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 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:162
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 1653 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().

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

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

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

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

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

References m_viscTensor.

187  {
188  fluxvector = m_viscTensor;
189  };
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 197 of file DiffusionLFRNS.h.

References m_viscTensor.

198  {
199  return m_viscTensor;
200  }
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 1235 of file DiffusionLFRNS.cpp.

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

Referenced by v_Diffuse().

1240  {
1241  int i, j;
1242  int nTracePts = fields[0]->GetTrace()->GetTotPoints();
1243  int nScalars = inarray.num_elements();
1244  int nDim = fields[0]->GetCoordim(0);
1245 
1246  Array<OneD, NekDouble > Vn (nTracePts, 0.0);
1247  Array<OneD, NekDouble > fluxtemp(nTracePts, 0.0);
1248 
1249  // Get the normal velocity Vn
1250  for (i = 0; i < nDim; ++i)
1251  {
1252  fields[0]->ExtractTracePhys(inarray[i], m_traceVel[i]);
1253  Vmath::Vvtvp(nTracePts, m_traceNormals[i], 1,
1254  m_traceVel[i], 1, Vn, 1, Vn, 1);
1255  }
1256 
1257  // Store forwards/backwards space along trace space
1258  Array<OneD, Array<OneD, NekDouble> > Fwd (nScalars);
1259  Array<OneD, Array<OneD, NekDouble> > Bwd (nScalars);
1260  Array<OneD, Array<OneD, NekDouble> > numflux(nScalars);
1261 
1262  for (i = 0; i < nScalars; ++i)
1263  {
1264  Fwd[i] = Array<OneD, NekDouble>(nTracePts);
1265  Bwd[i] = Array<OneD, NekDouble>(nTracePts);
1266  numflux[i] = Array<OneD, NekDouble>(nTracePts);
1267  fields[i]->GetFwdBwdTracePhys(inarray[i], Fwd[i], Bwd[i]);
1268  fields[0]->GetTrace()->Upwind(Vn, Fwd[i], Bwd[i], numflux[i]);
1269  }
1270 
1271  // Modify the values in case of boundary interfaces
1272  if (fields[0]->GetBndCondExpansions().num_elements())
1273  {
1274  v_WeakPenaltyO1(fields, inarray, numflux);
1275  }
1276 
1277  // Splitting the numerical flux into the dimensions
1278  for (j = 0; j < nDim; ++j)
1279  {
1280  for (i = 0; i < nScalars; ++i)
1281  {
1282  Vmath::Vcopy(nTracePts,numflux[i], 1,
1283  numericalFluxO1[i][j], 1);
1284  }
1285  }
1286  }
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:1292
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:442
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:1061
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 1500 of file DiffusionLFRNS.cpp.

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

Referenced by v_Diffuse().

1505  {
1506  int i, j;
1507  int nTracePts = fields[0]->GetTrace()->GetTotPoints();
1508  int nVariables = fields.num_elements();
1509  int nDim = fields[0]->GetCoordim(0);
1510 
1511  Array<OneD, NekDouble > Fwd(nTracePts);
1512  Array<OneD, NekDouble > Bwd(nTracePts);
1513  Array<OneD, NekDouble > Vn (nTracePts, 0.0);
1514 
1515  Array<OneD, NekDouble > qFwd (nTracePts);
1516  Array<OneD, NekDouble > qBwd (nTracePts);
1517  Array<OneD, NekDouble > qfluxtemp(nTracePts, 0.0);
1518 
1519  // Get the normal velocity Vn
1520  for (i = 0; i < nDim; ++i)
1521  {
1522  fields[0]->ExtractTracePhys(ufield[i], m_traceVel[i]);
1523  Vmath::Vvtvp(nTracePts, m_traceNormals[i], 1,
1524  m_traceVel[i], 1, Vn, 1, Vn, 1);
1525  }
1526 
1527  // Evaulate Riemann flux
1528  // qflux = \hat{q} \cdot u = q \cdot n
1529  // Notice: i = 1 (first row of the viscous tensor is zero)
1530  for (i = 1; i < nVariables; ++i)
1531  {
1532  qflux[i] = Array<OneD, NekDouble> (nTracePts, 0.0);
1533  for (j = 0; j < nDim; ++j)
1534  {
1535  // Compute qFwd and qBwd value of qfield in position 'ji'
1536  fields[i]->GetFwdBwdTracePhys(qfield[j][i], qFwd, qBwd);
1537 
1538  // Get Riemann flux of qflux --> LDG implies upwind
1539  fields[i]->GetTrace()->Upwind(Vn, qBwd, qFwd, qfluxtemp);
1540 
1541  // Multiply the Riemann flux by the trace normals
1542  Vmath::Vmul(nTracePts, m_traceNormals[j], 1, qfluxtemp, 1,
1543  qfluxtemp, 1);
1544 
1545  // Impose weak boundary condition with flux
1546  if (fields[0]->GetBndCondExpansions().num_elements())
1547  {
1548  v_WeakPenaltyO2(fields, i, j, qfield[j][i], qfluxtemp);
1549  }
1550 
1551  // Store the final flux into qflux
1552  Vmath::Vadd(nTracePts, qfluxtemp, 1, qflux[i], 1,
1553  qflux[i], 1);
1554  }
1555  }
1556  }
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:442
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:1563
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:299
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:183
virtual void Nektar::SolverUtils::DiffusionLFRNS::v_SetHomoDerivs ( Array< OneD, Array< OneD, NekDouble > > &  deriv)
inlineprotectedvirtual

Reimplemented from Nektar::SolverUtils::Diffusion.

Definition at line 191 of file DiffusionLFRNS.h.

193  {
194  m_homoDerivs = deriv;
195  }
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:198
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:2151
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:80
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:198
Nektar::LibUtilities::PointsKeyVector
std::vector< PointsKey > PointsKeyVector
Definition: Points.h:242
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::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 1292 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().

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

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

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 8 2017 12:30:02 for Nektar++ by   doxygen 1.8.8