#include <DiffusionLFRNS.h>
Inheritance diagram for Nektar::SolverUtils::DiffusionLFRNS:
Static Public Member Functions | |
| static DiffusionSharedPtr | create (std::string diffType) |
Static Public Attributes | |
| static std::string | type [] |
Protected Attributes | |
| std::string | m_diffType |
 Protected Attributes inherited from Nektar::SolverUtils::Diffusion | |
| Array< OneD, NekDouble > | m_divVel |
| Params for Ducros sensor. More... | |
| Array< OneD, NekDouble > | m_divVelSquare |
| Array< OneD, NekDouble > | m_curlVelSquare |
| DiffusionFluxVecCB | m_fluxVector |
| DiffusionFluxVecCBNS | m_fluxVectorNS |
| DiffusionFluxPenaltyNS | m_fluxPenaltyNS |
| DiffusionFluxCons | m_FunctorDiffusionfluxCons |
| DiffusionFluxCons | m_FunctorDiffusionfluxConsTrace |
| SpecialBndTreat | m_SpecialBndTreat |
| DiffusionSymmFluxCons | m_FunctorSymmetricfluxCons |
| NekDouble | m_time = 0.0 |
Private Member Functions | |
| void | 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... | |
| void | 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... | |
| void | 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... | |
| void | 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... | |
| void | 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... | |
| void | 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... | |
| void | 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... | |
| void | 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... | |
| void | 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... | |
| void | 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... | |
Detailed Description
Definition at line 42 of file DiffusionLFRNS.h.
Constructor & Destructor Documentation
◆ DiffusionLFRNS()
|
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 65 of file DiffusionLFRNS.cpp.
Referenced by create().
Member Function Documentation
◆ create()
|
inlinestatic |
Definition at line 45 of file DiffusionLFRNS.h.
46 {
48 }
DiffusionLFRNS(std::string diffType)
DiffusionLFRNS uses the Flux Reconstruction (FR) approach to compute the diffusion term....
Definition: DiffusionLFRNS.cpp:65
std::shared_ptr< Diffusion > DiffusionSharedPtr
A shared pointer to an EquationSystem object.
Definition: Diffusion.h:54
References DiffusionLFRNS().
◆ DerCFlux_1D()
|
private |
Compute the derivative of the corrective flux for 1D problems.
- Parameters
-
nConvectiveFields Number of fields. fields Pointer to fields. flux Volumetric flux in the physical space. iFlux Numerical interface flux in physical space. derCFlux Derivative of the corrective flux.
Definition at line 1650 of file DiffusionLFRNS.cpp.
1655{
1656 int n;
1657 int nLocalSolutionPts, phys_offset;
1658
1659 Array<OneD, NekDouble> auxArray1, auxArray2;
1660 Array<TwoD, const NekDouble> gmat;
1661 Array<OneD, const NekDouble> jac;
1662
1663 LibUtilities::PointsKeyVector ptsKeys;
1664 LibUtilities::BasisSharedPtr Basis;
1665 Basis = fields[0]->GetExp(0)->GetBasis(0);
1666
1667 int nElements = fields[0]->GetExpSize();
1668 int nPts = fields[0]->GetTotPoints();
1669
1670 // Arrays to store the derivatives of the correction flux
1671 Array<OneD, NekDouble> DCL(nPts / nElements, 0.0);
1672 Array<OneD, NekDouble> DCR(nPts / nElements, 0.0);
1673
1674 // Arrays to store the intercell numerical flux jumps
1675 Array<OneD, NekDouble> JumpL(nElements);
1676 Array<OneD, NekDouble> JumpR(nElements);
1677
1678 Array<OneD, Array<OneD, LocalRegions::ExpansionSharedPtr>> &elmtToTrace =
1679 fields[0]->GetTraceMap()->GetElmtToTrace();
1680
1681 for (n = 0; n < nElements; ++n)
1682 {
1683 nLocalSolutionPts = fields[0]->GetExp(n)->GetTotPoints();
1684 phys_offset = fields[0]->GetPhys_Offset(n);
1685
1686 Array<OneD, NekDouble> tmparrayX1(nLocalSolutionPts, 0.0);
1687 Array<OneD, NekDouble> tmpFluxVertex(1, 0.0);
1688 Vmath::Vcopy(nLocalSolutionPts, &flux[phys_offset], 1, &tmparrayX1[0],
1689 1);
1690
1691 fields[0]->GetExp(n)->GetTracePhysVals(0, elmtToTrace[n][0], tmparrayX1,
1692 tmpFluxVertex);
1693 JumpL[n] = iFlux[n] - tmpFluxVertex[0];
1694
1695 fields[0]->GetExp(n)->GetTracePhysVals(1, elmtToTrace[n][1], tmparrayX1,
1696 tmpFluxVertex);
1697 JumpR[n] = iFlux[n + 1] - tmpFluxVertex[0];
1698 }
1699
1700 for (n = 0; n < nElements; ++n)
1701 {
1702 ptsKeys = fields[0]->GetExp(n)->GetPointsKeys();
1703 nLocalSolutionPts = fields[0]->GetExp(n)->GetTotPoints();
1704 phys_offset = fields[0]->GetPhys_Offset(n);
1705 jac = fields[0]
1706 ->GetExp(n)
1707 ->as<LocalRegions::Expansion1D>()
1708 ->GetGeom1D()
1709 ->GetMetricInfo()
1710 ->GetJac(ptsKeys);
1711
1712 JumpL[n] = JumpL[n] * jac[0];
1713 JumpR[n] = JumpR[n] * jac[0];
1714
1715 // Left jump multiplied by left derivative of C function
1717 1);
1718
1719 // Right jump multiplied by right derivative of C function
1721 1);
1722
1723 // Assembling derivative of the correction flux
1724 Vmath::Vadd(nLocalSolutionPts, &DCL[0], 1, &DCR[0], 1,
1725 &derCFlux[phys_offset], 1);
1726 }
1727}
Array< OneD, Array< OneD, NekDouble > > m_dGR_xi1
Definition: DiffusionLFRNS.h:92
Array< OneD, Array< OneD, NekDouble > > m_dGL_xi1
Definition: DiffusionLFRNS.h:91
std::shared_ptr< Basis > BasisSharedPtr
Definition: FoundationsFwd.hpp:69
std::vector< PointsKey > PointsKeyVector
Definition: Points.h:231
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.hpp:180
void Smul(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Scalar multiply y = alpha*x.
Definition: Vmath.hpp:100
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.hpp:825
References Nektar::LocalRegions::Expansion::GetMetricInfo(), m_dGL_xi1, m_dGR_xi1, Vmath::Smul(), Vmath::Vadd(), and Vmath::Vcopy().
Referenced by v_Diffuse().
◆ DerCFlux_2D()
|
private |
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
-
nConvectiveFields Number of fields. fields Pointer to fields. flux Volumetric flux in the physical space. numericalFlux Numerical interface flux in physical space. derCFlux Derivative of the corrective flux.
- Todo:
- : Switch on shapes eventually here.
Definition at line 1744 of file DiffusionLFRNS.cpp.
1749{
1750 int n, e, i, j, cnt;
1751
1752 Array<OneD, const NekDouble> jac;
1753
1754 int nElements = fields[0]->GetExpSize();
1755 int trace_offset, phys_offset;
1756 int nLocalSolutionPts;
1757 int nquad0, nquad1;
1758 int nEdgePts;
1759
1760 LibUtilities::PointsKeyVector ptsKeys;
1761 Array<OneD, NekDouble> auxArray1, auxArray2;
1762 Array<OneD, LibUtilities::BasisSharedPtr> base;
1763 Array<OneD, Array<OneD, LocalRegions::ExpansionSharedPtr>> &elmtToTrace =
1764 fields[0]->GetTraceMap()->GetElmtToTrace();
1765
1766 // Loop on the elements
1767 for (n = 0; n < nElements; ++n)
1768 {
1769 // Offset of the element on the global vector
1770 phys_offset = fields[0]->GetPhys_Offset(n);
1771 nLocalSolutionPts = fields[0]->GetExp(n)->GetTotPoints();
1772 ptsKeys = fields[0]->GetExp(n)->GetPointsKeys();
1773
1774 jac = fields[0]
1775 ->GetExp(n)
1776 ->as<LocalRegions::Expansion2D>()
1777 ->GetGeom2D()
1778 ->GetMetricInfo()
1779 ->GetJac(ptsKeys);
1780
1781 base = fields[0]->GetExp(n)->GetBase();
1782 nquad0 = base[0]->GetNumPoints();
1783 nquad1 = base[1]->GetNumPoints();
1784
1785 Array<OneD, NekDouble> divCFluxE0(nLocalSolutionPts, 0.0);
1786 Array<OneD, NekDouble> divCFluxE1(nLocalSolutionPts, 0.0);
1787 Array<OneD, NekDouble> divCFluxE2(nLocalSolutionPts, 0.0);
1788 Array<OneD, NekDouble> divCFluxE3(nLocalSolutionPts, 0.0);
1789
1790 // Loop on the edges
1791 for (e = 0; e < fields[0]->GetExp(n)->GetNtraces(); ++e)
1792 {
1793 // Number of edge points of edge 'e'
1794 nEdgePts = fields[0]->GetExp(n)->GetTraceNumPoints(e);
1795
1796 // Array for storing volumetric fluxes on each edge
1797 Array<OneD, NekDouble> tmparray(nEdgePts, 0.0);
1798
1799 // Offset of the trace space correspondent to edge 'e'
1800 trace_offset = fields[0]->GetTrace()->GetPhys_Offset(
1801 elmtToTrace[n][e]->GetElmtId());
1802
1803 // Extract the edge values of the volumetric fluxes
1804 // on edge 'e' and order them accordingly to the order
1805 // of the trace space
1806 fields[0]->GetExp(n)->GetTracePhysVals(
1807 e, elmtToTrace[n][e], flux + phys_offset, auxArray1 = tmparray);
1808
1809 // Compute the fluxJumps per each edge 'e' and each
1810 // flux point
1811 Array<OneD, NekDouble> fluxJumps(nEdgePts, 0.0);
1812 Vmath::Vsub(nEdgePts, &iFlux[trace_offset], 1, &tmparray[0], 1,
1813 &fluxJumps[0], 1);
1814
1815 // Check the ordering of the fluxJumps and reverse
1816 // it in case of backward definition of edge 'e'
1817 if (fields[0]->GetExp(n)->GetTraceOrient(e) ==
1818 StdRegions::eBackwards)
1819 {
1820 Vmath::Reverse(nEdgePts, &fluxJumps[0], 1, &fluxJumps[0], 1);
1821 }
1822
1823 // Deformed elements
1824 if (fields[0]
1825 ->GetExp(n)
1826 ->as<LocalRegions::Expansion2D>()
1827 ->GetGeom2D()
1828 ->GetMetricInfo()
1829 ->GetGtype() == SpatialDomains::eDeformed)
1830 {
1831 // Extract the Jacobians along edge 'e'
1832 Array<OneD, NekDouble> jacEdge(nEdgePts, 0.0);
1833
1834 fields[0]->GetExp(n)->GetTracePhysVals(
1835 e, elmtToTrace[n][e], jac, auxArray1 = jacEdge);
1836
1837 // Check the ordering of the fluxJumps and reverse
1838 // it in case of backward definition of edge 'e'
1839 if (fields[0]->GetExp(n)->GetTraceOrient(e) ==
1840 StdRegions::eBackwards)
1841 {
1842 Vmath::Reverse(nEdgePts, &jacEdge[0], 1, &jacEdge[0], 1);
1843 }
1844
1845 // Multiply the fluxJumps by the edge 'e' Jacobians
1846 // to bring the fluxJumps into the standard space
1847 for (j = 0; j < nEdgePts; j++)
1848 {
1849 fluxJumps[j] = fluxJumps[j] * jacEdge[j];
1850 }
1851 }
1852 // Non-deformed elements
1853 else
1854 {
1855 // Multiply the fluxJumps by the edge 'e' Jacobians
1856 // to bring the fluxJumps into the standard space
1857 Vmath::Smul(nEdgePts, jac[0], fluxJumps, 1, fluxJumps, 1);
1858 }
1859
1860 // Multiply jumps by derivatives of the correction functions
1861 // All the quntities at this level should be defined into
1862 // the standard space
1863 switch (e)
1864 {
1865 case 0:
1866 for (i = 0; i < nquad0; ++i)
1867 {
1868 for (j = 0; j < nquad1; ++j)
1869 {
1870 cnt = i + j * nquad0;
1871 divCFluxE0[cnt] = fluxJumps[i] * m_dGL_xi2[n][j];
1872 }
1873 }
1874 break;
1875 case 1:
1876 for (i = 0; i < nquad1; ++i)
1877 {
1878 for (j = 0; j < nquad0; ++j)
1879 {
1880 cnt = (nquad0)*i + j;
1881 divCFluxE1[cnt] = fluxJumps[i] * m_dGR_xi1[n][j];
1882 }
1883 }
1884 break;
1885 case 2:
1886 for (i = 0; i < nquad0; ++i)
1887 {
1888 for (j = 0; j < nquad1; ++j)
1889 {
1890 cnt = j * nquad0 + i;
1891 divCFluxE2[cnt] = fluxJumps[i] * m_dGR_xi2[n][j];
1892 }
1893 }
1894 break;
1895 case 3:
1896 for (i = 0; i < nquad1; ++i)
1897 {
1898 for (j = 0; j < nquad0; ++j)
1899 {
1900 cnt = j + i * nquad0;
1901 divCFluxE3[cnt] = fluxJumps[i] * m_dGL_xi1[n][j];
1902 }
1903 }
1904 break;
1905
1906 default:
1908 break;
1909 }
1910 }
1911
1912 // Summing the component of the flux for calculating the
1913 // derivatives per each direction
1914 if (direction == 0)
1915 {
1916 Vmath::Vadd(nLocalSolutionPts, &divCFluxE1[0], 1, &divCFluxE3[0], 1,
1917 &derCFlux[phys_offset], 1);
1918 }
1919 else if (direction == 1)
1920 {
1921 Vmath::Vadd(nLocalSolutionPts, &divCFluxE0[0], 1, &divCFluxE2[0], 1,
1922 &derCFlux[phys_offset], 1);
1923 }
1924 }
1925}
Array< OneD, Array< OneD, NekDouble > > m_dGR_xi2
Definition: DiffusionLFRNS.h:94
Array< OneD, Array< OneD, NekDouble > > m_dGL_xi2
Definition: DiffusionLFRNS.h:93
void Reverse(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.hpp:844
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.hpp:220
References ASSERTL0, Nektar::StdRegions::eBackwards, Nektar::SpatialDomains::eDeformed, Nektar::LocalRegions::Expansion::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().
◆ DivCFlux_2D()
|
private |
Compute the divergence of the corrective flux for 2D problems.
- Parameters
-
nConvectiveFields Number of fields. fields Pointer to fields. fluxX1 X1 - volumetric flux in physical space. fluxX2 X2 - volumetric flux in physical space. numericalFlux Interface flux in physical space. divCFlux Divergence of the corrective flux for 2D problems.
- Todo:
- : Switch on shapes eventually here.
Definition at line 1940 of file DiffusionLFRNS.cpp.
1947{
1948 int n, e, i, j, cnt;
1949
1950 int nElements = fields[0]->GetExpSize();
1951
1952 int nLocalSolutionPts;
1953 int nEdgePts;
1954 int trace_offset;
1955 int phys_offset;
1956 int nquad0;
1957 int nquad1;
1958
1959 Array<OneD, NekDouble> auxArray1, auxArray2;
1960 Array<OneD, LibUtilities::BasisSharedPtr> base;
1961
1962 Array<OneD, Array<OneD, LocalRegions::ExpansionSharedPtr>> &elmtToTrace =
1963 fields[0]->GetTraceMap()->GetElmtToTrace();
1964
1965 // Loop on the elements
1966 for (n = 0; n < nElements; ++n)
1967 {
1968 // Offset of the element on the global vector
1969 phys_offset = fields[0]->GetPhys_Offset(n);
1970 nLocalSolutionPts = fields[0]->GetExp(n)->GetTotPoints();
1971
1972 base = fields[0]->GetExp(n)->GetBase();
1973 nquad0 = base[0]->GetNumPoints();
1974 nquad1 = base[1]->GetNumPoints();
1975
1976 Array<OneD, NekDouble> divCFluxE0(nLocalSolutionPts, 0.0);
1977 Array<OneD, NekDouble> divCFluxE1(nLocalSolutionPts, 0.0);
1978 Array<OneD, NekDouble> divCFluxE2(nLocalSolutionPts, 0.0);
1979 Array<OneD, NekDouble> divCFluxE3(nLocalSolutionPts, 0.0);
1980
1981 // Loop on the edges
1982 for (e = 0; e < fields[0]->GetExp(n)->GetNtraces(); ++e)
1983 {
1984 // Number of edge points of edge e
1985 nEdgePts = fields[0]->GetExp(n)->GetTraceNumPoints(e);
1986
1987 Array<OneD, NekDouble> tmparrayX1(nEdgePts, 0.0);
1988 Array<OneD, NekDouble> tmparrayX2(nEdgePts, 0.0);
1989 Array<OneD, NekDouble> fluxN(nEdgePts, 0.0);
1990 Array<OneD, NekDouble> fluxT(nEdgePts, 0.0);
1991 Array<OneD, NekDouble> fluxJumps(nEdgePts, 0.0);
1992
1993 // Offset of the trace space correspondent to edge e
1994 trace_offset = fields[0]->GetTrace()->GetPhys_Offset(
1995 elmtToTrace[n][e]->GetElmtId());
1996
1997 // Get the normals of edge e
1998 const Array<OneD, const Array<OneD, NekDouble>> &normals =
1999 fields[0]->GetExp(n)->GetTraceNormal(e);
2000
2001 // Extract the edge values of flux-x on edge e and order
2002 // them accordingly to the order of the trace space
2003 fields[0]->GetExp(n)->GetTracePhysVals(e, elmtToTrace[n][e],
2004 fluxX1 + phys_offset,
2005 auxArray1 = tmparrayX1);
2006
2007 // Extract the edge values of flux-y on edge e and order
2008 // them accordingly to the order of the trace space
2009 fields[0]->GetExp(n)->GetTracePhysVals(e, elmtToTrace[n][e],
2010 fluxX2 + phys_offset,
2011 auxArray1 = tmparrayX2);
2012
2013 // Multiply the edge components of the flux by the normal
2014 for (i = 0; i < nEdgePts; ++i)
2015 {
2016 fluxN[i] = tmparrayX1[i] * m_traceNormals[0][trace_offset + i] +
2017 tmparrayX2[i] * m_traceNormals[1][trace_offset + i];
2018 }
2019
2020 // Subtract to the Riemann flux the discontinuous flux
2021 Vmath::Vsub(nEdgePts, &numericalFlux[trace_offset], 1, &fluxN[0], 1,
2022 &fluxJumps[0], 1);
2023
2024 // Check the ordering of the jump vectors
2025 if (fields[0]->GetExp(n)->GetTraceOrient(e) ==
2026 StdRegions::eBackwards)
2027 {
2028 Vmath::Reverse(nEdgePts, auxArray1 = fluxJumps, 1,
2029 auxArray2 = fluxJumps, 1);
2030 }
2031
2032 for (i = 0; i < nEdgePts; ++i)
2033 {
2035 m_traceNormals[1][trace_offset + i] != normals[1][i])
2036 {
2037 fluxJumps[i] = -fluxJumps[i];
2038 }
2039 }
2040
2041 // Multiply jumps by derivatives of the correction functions
2042 switch (e)
2043 {
2044 case 0:
2045 for (i = 0; i < nquad0; ++i)
2046 {
2047 // Multiply fluxJumps by Q factors
2048 fluxJumps[i] = -(m_Q2D_e0[n][i]) * fluxJumps[i];
2049
2050 for (j = 0; j < nquad1; ++j)
2051 {
2052 cnt = i + j * nquad0;
2053 divCFluxE0[cnt] = fluxJumps[i] * m_dGL_xi2[n][j];
2054 }
2055 }
2056 break;
2057 case 1:
2058 for (i = 0; i < nquad1; ++i)
2059 {
2060 // Multiply fluxJumps by Q factors
2061 fluxJumps[i] = (m_Q2D_e1[n][i]) * fluxJumps[i];
2062
2063 for (j = 0; j < nquad0; ++j)
2064 {
2065 cnt = (nquad0)*i + j;
2066 divCFluxE1[cnt] = fluxJumps[i] * m_dGR_xi1[n][j];
2067 }
2068 }
2069 break;
2070 case 2:
2071 for (i = 0; i < nquad0; ++i)
2072 {
2073 // Multiply fluxJumps by Q factors
2074 fluxJumps[i] = (m_Q2D_e2[n][i]) * fluxJumps[i];
2075
2076 for (j = 0; j < nquad1; ++j)
2077 {
2078 cnt = j * nquad0 + i;
2079 divCFluxE2[cnt] = fluxJumps[i] * m_dGR_xi2[n][j];
2080 }
2081 }
2082 break;
2083 case 3:
2084 for (i = 0; i < nquad1; ++i)
2085 {
2086 // Multiply fluxJumps by Q factors
2087 fluxJumps[i] = -(m_Q2D_e3[n][i]) * fluxJumps[i];
2088 for (j = 0; j < nquad0; ++j)
2089 {
2090 cnt = j + i * nquad0;
2091 divCFluxE3[cnt] = fluxJumps[i] * m_dGL_xi1[n][j];
2092 }
2093 }
2094 break;
2095
2096 default:
2098 break;
2099 }
2100 }
2101
2102 // Sum each edge contribution
2103 Vmath::Vadd(nLocalSolutionPts, &divCFluxE0[0], 1, &divCFluxE1[0], 1,
2104 &divCFlux[phys_offset], 1);
2105
2106 Vmath::Vadd(nLocalSolutionPts, &divCFlux[phys_offset], 1,
2107 &divCFluxE2[0], 1, &divCFlux[phys_offset], 1);
2108
2109 Vmath::Vadd(nLocalSolutionPts, &divCFlux[phys_offset], 1,
2110 &divCFluxE3[0], 1, &divCFlux[phys_offset], 1);
2111 }
2112}
Array< OneD, Array< OneD, NekDouble > > m_Q2D_e1
Definition: DiffusionLFRNS.h:87
Array< OneD, Array< OneD, NekDouble > > m_Q2D_e3
Definition: DiffusionLFRNS.h:89
Array< OneD, Array< OneD, NekDouble > > m_traceNormals
Definition: DiffusionLFRNS.h:80
Array< OneD, Array< OneD, NekDouble > > m_Q2D_e2
Definition: DiffusionLFRNS.h:88
Array< OneD, Array< OneD, NekDouble > > m_Q2D_e0
Definition: DiffusionLFRNS.h:86
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().
◆ DivCFlux_2D_Gauss()
|
private |
Compute the divergence of the corrective flux for 2D problems where POINTSTYPE="GaussGaussLegendre".
- Parameters
-
nConvectiveFields Number of fields. fields Pointer to fields. fluxX1 X1-volumetric flux in physical space. fluxX2 X2-volumetric flux in physical space. numericalFlux Interface flux in physical space. divCFlux Divergence of the corrective flux.
- Todo:
- : Switch on shapes eventually here.
Definition at line 2128 of file DiffusionLFRNS.cpp.
2135{
2136 int n, e, i, j, cnt;
2137
2138 int nElements = fields[0]->GetExpSize();
2139 int nLocalSolutionPts;
2140 int nEdgePts;
2141 int trace_offset;
2142 int phys_offset;
2143 int nquad0;
2144 int nquad1;
2145
2146 Array<OneD, NekDouble> auxArray1, auxArray2;
2147 Array<OneD, LibUtilities::BasisSharedPtr> base;
2148
2149 Array<OneD, Array<OneD, LocalRegions::ExpansionSharedPtr>> &elmtToTrace =
2150 fields[0]->GetTraceMap()->GetElmtToTrace();
2151
2152 // Loop on the elements
2153 for (n = 0; n < nElements; ++n)
2154 {
2155 // Offset of the element on the global vector
2156 phys_offset = fields[0]->GetPhys_Offset(n);
2157 nLocalSolutionPts = fields[0]->GetExp(n)->GetTotPoints();
2158
2159 base = fields[0]->GetExp(n)->GetBase();
2160 nquad0 = base[0]->GetNumPoints();
2161 nquad1 = base[1]->GetNumPoints();
2162
2163 Array<OneD, NekDouble> divCFluxE0(nLocalSolutionPts, 0.0);
2164 Array<OneD, NekDouble> divCFluxE1(nLocalSolutionPts, 0.0);
2165 Array<OneD, NekDouble> divCFluxE2(nLocalSolutionPts, 0.0);
2166 Array<OneD, NekDouble> divCFluxE3(nLocalSolutionPts, 0.0);
2167
2168 // Loop on the edges
2169 for (e = 0; e < fields[0]->GetExp(n)->GetNtraces(); ++e)
2170 {
2171 // Number of edge points of edge e
2172 nEdgePts = fields[0]->GetExp(n)->GetTraceNumPoints(e);
2173
2174 Array<OneD, NekDouble> fluxN(nEdgePts, 0.0);
2175 Array<OneD, NekDouble> fluxT(nEdgePts, 0.0);
2176 Array<OneD, NekDouble> fluxN_R(nEdgePts, 0.0);
2177 Array<OneD, NekDouble> fluxN_D(nEdgePts, 0.0);
2178 Array<OneD, NekDouble> fluxJumps(nEdgePts, 0.0);
2179
2180 // Offset of the trace space correspondent to edge e
2181 trace_offset = fields[0]->GetTrace()->GetPhys_Offset(
2182 elmtToTrace[n][e]->GetElmtId());
2183
2184 // Get the normals of edge e
2185 const Array<OneD, const Array<OneD, NekDouble>> &normals =
2186 fields[0]->GetExp(n)->GetTraceNormal(e);
2187
2188 // Extract the trasformed normal flux at each edge
2189 switch (e)
2190 {
2191 case 0:
2192 // Extract the edge values of transformed flux-y on
2193 // edge e and order them accordingly to the order of
2194 // the trace space
2195 fields[0]->GetExp(n)->GetTracePhysVals(e, elmtToTrace[n][e],
2196 fluxX2 + phys_offset,
2197 auxArray1 = fluxN_D);
2198
2199 Vmath::Neg(nEdgePts, fluxN_D, 1);
2200
2201 // Extract the physical Rieamann flux at each edge
2202 Vmath::Vcopy(nEdgePts, &numericalFlux[trace_offset], 1,
2203 &fluxN[0], 1);
2204
2205 // Check the ordering of vectors
2206 if (fields[0]->GetExp(n)->GetTraceOrient(e) ==
2207 StdRegions::eBackwards)
2208 {
2209 Vmath::Reverse(nEdgePts, auxArray1 = fluxN, 1,
2210 auxArray2 = fluxN, 1);
2211
2212 Vmath::Reverse(nEdgePts, auxArray1 = fluxN_D, 1,
2213 auxArray2 = fluxN_D, 1);
2214 }
2215
2216 // Transform Riemann Fluxes in the standard element
2217 for (i = 0; i < nquad0; ++i)
2218 {
2219 // Multiply Riemann Flux by Q factors
2220 fluxN_R[i] = (m_Q2D_e0[n][i]) * fluxN[i];
2221 }
2222
2223 for (i = 0; i < nEdgePts; ++i)
2224 {
2226 normals[0][i] ||
2227 m_traceNormals[1][trace_offset + i] !=
2228 normals[1][i])
2229 {
2230 fluxN_R[i] = -fluxN_R[i];
2231 }
2232 }
2233
2234 // Subtract to the Riemann flux the discontinuous
2235 // flux
2236 Vmath::Vsub(nEdgePts, &fluxN_R[0], 1, &fluxN_D[0], 1,
2237 &fluxJumps[0], 1);
2238
2239 // Multiplicate the flux jumps for the correction
2240 // function
2241 for (i = 0; i < nquad0; ++i)
2242 {
2243 for (j = 0; j < nquad1; ++j)
2244 {
2245 cnt = i + j * nquad0;
2246 divCFluxE0[cnt] = -fluxJumps[i] * m_dGL_xi2[n][j];
2247 }
2248 }
2249 break;
2250 case 1:
2251 // Extract the edge values of transformed flux-x on
2252 // edge e and order them accordingly to the order of
2253 // the trace space
2254 fields[0]->GetExp(n)->GetTracePhysVals(e, elmtToTrace[n][e],
2255 fluxX1 + phys_offset,
2256 auxArray1 = fluxN_D);
2257
2258 // Extract the physical Rieamann flux at each edge
2259 Vmath::Vcopy(nEdgePts, &numericalFlux[trace_offset], 1,
2260 &fluxN[0], 1);
2261
2262 // Check the ordering of vectors
2263 if (fields[0]->GetExp(n)->GetTraceOrient(e) ==
2264 StdRegions::eBackwards)
2265 {
2266 Vmath::Reverse(nEdgePts, auxArray1 = fluxN, 1,
2267 auxArray2 = fluxN, 1);
2268
2269 Vmath::Reverse(nEdgePts, auxArray1 = fluxN_D, 1,
2270 auxArray2 = fluxN_D, 1);
2271 }
2272
2273 // Transform Riemann Fluxes in the standard element
2274 for (i = 0; i < nquad1; ++i)
2275 {
2276 // Multiply Riemann Flux by Q factors
2277 fluxN_R[i] = (m_Q2D_e1[n][i]) * fluxN[i];
2278 }
2279
2280 for (i = 0; i < nEdgePts; ++i)
2281 {
2283 normals[0][i] ||
2284 m_traceNormals[1][trace_offset + i] !=
2285 normals[1][i])
2286 {
2287 fluxN_R[i] = -fluxN_R[i];
2288 }
2289 }
2290
2291 // Subtract to the Riemann flux the discontinuous
2292 // flux
2293 Vmath::Vsub(nEdgePts, &fluxN_R[0], 1, &fluxN_D[0], 1,
2294 &fluxJumps[0], 1);
2295
2296 // Multiplicate the flux jumps for the correction
2297 // function
2298 for (i = 0; i < nquad1; ++i)
2299 {
2300 for (j = 0; j < nquad0; ++j)
2301 {
2302 cnt = (nquad0)*i + j;
2303 divCFluxE1[cnt] = fluxJumps[i] * m_dGR_xi1[n][j];
2304 }
2305 }
2306 break;
2307 case 2:
2308
2309 // Extract the edge values of transformed flux-y on
2310 // edge e and order them accordingly to the order of
2311 // the trace space
2312
2313 fields[0]->GetExp(n)->GetTracePhysVals(e, elmtToTrace[n][e],
2314 fluxX2 + phys_offset,
2315 auxArray1 = fluxN_D);
2316
2317 // Extract the physical Rieamann flux at each edge
2318 Vmath::Vcopy(nEdgePts, &numericalFlux[trace_offset], 1,
2319 &fluxN[0], 1);
2320
2321 // Check the ordering of vectors
2322 if (fields[0]->GetExp(n)->GetTraceOrient(e) ==
2323 StdRegions::eBackwards)
2324 {
2325 Vmath::Reverse(nEdgePts, auxArray1 = fluxN, 1,
2326 auxArray2 = fluxN, 1);
2327
2328 Vmath::Reverse(nEdgePts, auxArray1 = fluxN_D, 1,
2329 auxArray2 = fluxN_D, 1);
2330 }
2331
2332 // Transform Riemann Fluxes in the standard element
2333 for (i = 0; i < nquad0; ++i)
2334 {
2335 // Multiply Riemann Flux by Q factors
2336 fluxN_R[i] = (m_Q2D_e2[n][i]) * fluxN[i];
2337 }
2338
2339 for (i = 0; i < nEdgePts; ++i)
2340 {
2342 normals[0][i] ||
2343 m_traceNormals[1][trace_offset + i] !=
2344 normals[1][i])
2345 {
2346 fluxN_R[i] = -fluxN_R[i];
2347 }
2348 }
2349
2350 // Subtract to the Riemann flux the discontinuous
2351 // flux
2352
2353 Vmath::Vsub(nEdgePts, &fluxN_R[0], 1, &fluxN_D[0], 1,
2354 &fluxJumps[0], 1);
2355
2356 // Multiplicate the flux jumps for the correction
2357 // function
2358 for (i = 0; i < nquad0; ++i)
2359 {
2360 for (j = 0; j < nquad1; ++j)
2361 {
2362 cnt = j * nquad0 + i;
2363 divCFluxE2[cnt] = fluxJumps[i] * m_dGR_xi2[n][j];
2364 }
2365 }
2366 break;
2367 case 3:
2368 // Extract the edge values of transformed flux-x on
2369 // edge e and order them accordingly to the order of
2370 // the trace space
2371
2372 fields[0]->GetExp(n)->GetTracePhysVals(e, elmtToTrace[n][e],
2373 fluxX1 + phys_offset,
2374 auxArray1 = fluxN_D);
2375 Vmath::Neg(nEdgePts, fluxN_D, 1);
2376
2377 // Extract the physical Rieamann flux at each edge
2378 Vmath::Vcopy(nEdgePts, &numericalFlux[trace_offset], 1,
2379 &fluxN[0], 1);
2380
2381 // Check the ordering of vectors
2382 if (fields[0]->GetExp(n)->GetTraceOrient(e) ==
2383 StdRegions::eBackwards)
2384 {
2385 Vmath::Reverse(nEdgePts, auxArray1 = fluxN, 1,
2386 auxArray2 = fluxN, 1);
2387
2388 Vmath::Reverse(nEdgePts, auxArray1 = fluxN_D, 1,
2389 auxArray2 = fluxN_D, 1);
2390 }
2391
2392 // Transform Riemann Fluxes in the standard element
2393 for (i = 0; i < nquad1; ++i)
2394 {
2395 // Multiply Riemann Flux by Q factors
2396 fluxN_R[i] = (m_Q2D_e3[n][i]) * fluxN[i];
2397 }
2398
2399 for (i = 0; i < nEdgePts; ++i)
2400 {
2402 normals[0][i] ||
2403 m_traceNormals[1][trace_offset + i] !=
2404 normals[1][i])
2405 {
2406 fluxN_R[i] = -fluxN_R[i];
2407 }
2408 }
2409
2410 // Subtract to the Riemann flux the discontinuous
2411 // flux
2412
2413 Vmath::Vsub(nEdgePts, &fluxN_R[0], 1, &fluxN_D[0], 1,
2414 &fluxJumps[0], 1);
2415
2416 // Multiplicate the flux jumps for the correction
2417 // function
2418 for (i = 0; i < nquad1; ++i)
2419 {
2420 for (j = 0; j < nquad0; ++j)
2421 {
2422 cnt = j + i * nquad0;
2423 divCFluxE3[cnt] = -fluxJumps[i] * m_dGL_xi1[n][j];
2424 }
2425 }
2426 break;
2427 default:
2429 break;
2430 }
2431 }
2432
2433 // Sum each edge contribution
2434 Vmath::Vadd(nLocalSolutionPts, &divCFluxE0[0], 1, &divCFluxE1[0], 1,
2435 &divCFlux[phys_offset], 1);
2436
2437 Vmath::Vadd(nLocalSolutionPts, &divCFlux[phys_offset], 1,
2438 &divCFluxE2[0], 1, &divCFlux[phys_offset], 1);
2439
2440 Vmath::Vadd(nLocalSolutionPts, &divCFlux[phys_offset], 1,
2441 &divCFluxE3[0], 1, &divCFlux[phys_offset], 1);
2442 }
2443}
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().
◆ NumericalFluxO1()
|
private |
Builds the numerical flux for the 1st order derivatives.
Definition at line 1227 of file DiffusionLFRNS.cpp.
1231{
1232 int i, j;
1233 int nTracePts = fields[0]->GetTrace()->GetTotPoints();
1234 int nScalars = inarray.size();
1235 int nDim = fields[0]->GetCoordim(0);
1236
1237 Array<OneD, NekDouble> Vn(nTracePts, 0.0);
1238 Array<OneD, NekDouble> fluxtemp(nTracePts, 0.0);
1239
1240 // Get the normal velocity Vn
1241 for (i = 0; i < nDim; ++i)
1242 {
1243 fields[0]->ExtractTracePhys(inarray[i], m_traceVel[i]);
1245 Vn, 1);
1246 }
1247
1248 // Store forwards/backwards space along trace space
1249 Array<OneD, Array<OneD, NekDouble>> Fwd(nScalars);
1250 Array<OneD, Array<OneD, NekDouble>> Bwd(nScalars);
1251 Array<OneD, Array<OneD, NekDouble>> numflux(nScalars);
1252
1253 for (i = 0; i < nScalars; ++i)
1254 {
1255 Fwd[i] = Array<OneD, NekDouble>(nTracePts);
1256 Bwd[i] = Array<OneD, NekDouble>(nTracePts);
1257 numflux[i] = Array<OneD, NekDouble>(nTracePts);
1258 fields[i]->GetFwdBwdTracePhys(inarray[i], Fwd[i], Bwd[i]);
1259 fields[0]->GetTrace()->Upwind(Vn, Fwd[i], Bwd[i], numflux[i]);
1260 }
1261
1262 // Modify the values in case of boundary interfaces
1263 if (fields[0]->GetBndCondExpansions().size())
1264 {
1265 WeakPenaltyO1(fields, inarray, numflux);
1266 }
1267
1268 // Splitting the numerical flux into the dimensions
1269 for (j = 0; j < nDim; ++j)
1270 {
1271 for (i = 0; i < nScalars; ++i)
1272 {
1273 Vmath::Vcopy(nTracePts, numflux[i], 1, numericalFluxO1[i][j], 1);
1274 }
1275 }
1276}
Array< OneD, Array< OneD, NekDouble > > m_traceVel
Definition: DiffusionLFRNS.h:79
void 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:1282
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.hpp:366
References m_traceNormals, m_traceVel, Vmath::Vcopy(), Vmath::Vvtvp(), and WeakPenaltyO1().
Referenced by v_Diffuse().
◆ NumericalFluxO2()
|
private |
Build the numerical flux for the 2nd order derivatives.
Definition at line 1497 of file DiffusionLFRNS.cpp.
1502{
1503 int i, j;
1504 int nTracePts = fields[0]->GetTrace()->GetTotPoints();
1505 int nVariables = fields.size();
1506 int nDim = fields[0]->GetCoordim(0);
1507
1508 Array<OneD, NekDouble> Fwd(nTracePts);
1509 Array<OneD, NekDouble> Bwd(nTracePts);
1510 Array<OneD, NekDouble> Vn(nTracePts, 0.0);
1511
1512 Array<OneD, NekDouble> qFwd(nTracePts);
1513 Array<OneD, NekDouble> qBwd(nTracePts);
1514 Array<OneD, NekDouble> qfluxtemp(nTracePts, 0.0);
1515
1516 // Get the normal velocity Vn
1517 for (i = 0; i < nDim; ++i)
1518 {
1519 fields[0]->ExtractTracePhys(ufield[i], m_traceVel[i]);
1521 Vn, 1);
1522 }
1523
1524 // Evaulate Riemann flux
1525 // qflux = \hat{q} \cdot u = q \cdot n
1526 // Notice: i = 1 (first row of the viscous tensor is zero)
1527 for (i = 1; i < nVariables; ++i)
1528 {
1529 qflux[i] = Array<OneD, NekDouble>(nTracePts, 0.0);
1530 for (j = 0; j < nDim; ++j)
1531 {
1532 // Compute qFwd and qBwd value of qfield in position 'ji'
1533 fields[i]->GetFwdBwdTracePhys(qfield[j][i], qFwd, qBwd);
1534
1535 // Get Riemann flux of qflux --> LDG implies upwind
1536 fields[i]->GetTrace()->Upwind(Vn, qBwd, qFwd, qfluxtemp);
1537
1538 // Multiply the Riemann flux by the trace normals
1540 qfluxtemp, 1);
1541
1542 // Impose weak boundary condition with flux
1543 if (fields[0]->GetBndCondExpansions().size())
1544 {
1545 WeakPenaltyO2(fields, i, j, qfield[j][i], qfluxtemp);
1546 }
1547
1548 // Store the final flux into qflux
1549 Vmath::Vadd(nTracePts, qfluxtemp, 1, qflux[i], 1, qflux[i], 1);
1550 }
1551 }
1552}
void 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:1558
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.hpp:72
References m_traceNormals, m_traceVel, Vmath::Vadd(), Vmath::Vmul(), Vmath::Vvtvp(), and WeakPenaltyO2().
Referenced by v_Diffuse().
◆ SetupCFunctions()
|
private |
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
-
pSession Pointer to session reader. pFields Pointer to fields.
Definition at line 369 of file DiffusionLFRNS.cpp.
372{
373 int i, n;
374 NekDouble c0 = 0.0;
375 NekDouble c1 = 0.0;
376 NekDouble c2 = 0.0;
377 int nquad0, nquad1, nquad2;
378 int nmodes0, nmodes1, nmodes2;
379 Array<OneD, LibUtilities::BasisSharedPtr> base;
380
381 int nElements = pFields[0]->GetExpSize();
382 int nDim = pFields[0]->GetCoordim(0);
383
384 switch (nDim)
385 {
386 case 1:
387 {
388 m_dGL_xi1 = Array<OneD, Array<OneD, NekDouble>>(nElements);
389 m_dGR_xi1 = Array<OneD, Array<OneD, NekDouble>>(nElements);
390
391 for (n = 0; n < nElements; ++n)
392 {
393 base = pFields[0]->GetExp(n)->GetBase();
394 nquad0 = base[0]->GetNumPoints();
395 nmodes0 = base[0]->GetNumModes();
396 Array<OneD, const NekDouble> z0;
397 Array<OneD, const NekDouble> w0;
398
399 base[0]->GetZW(z0, w0);
400
401 m_dGL_xi1[n] = Array<OneD, NekDouble>(nquad0);
402 m_dGR_xi1[n] = Array<OneD, NekDouble>(nquad0);
403
404 // Auxiliary vectors to build up the auxiliary
405 // derivatives of the Legendre polynomials
406 Array<OneD, NekDouble> dLp0(nquad0, 0.0);
407 Array<OneD, NekDouble> dLpp0(nquad0, 0.0);
408 Array<OneD, NekDouble> dLpm0(nquad0, 0.0);
409
410 // Degree of the correction functions
411 int p0 = nmodes0 - 1;
412
413 // Function sign to compute left correction function
414 NekDouble sign0 = pow(-1.0, p0);
415
416 // Factorial factor to build the scheme
417 NekDouble ap0 =
418 std::tgamma(2 * p0 + 1) /
419 (pow(2.0, p0) * std::tgamma(p0 + 1) * std::tgamma(p0 + 1));
420
421 // Scalar parameter which recovers the FR schemes
423 {
424 c0 = 0.0;
425 }
427 {
428 c0 = 2.0 * p0 /
429 ((2.0 * p0 + 1.0) * (p0 + 1.0) *
430 (ap0 * std::tgamma(p0 + 1)) *
431 (ap0 * std::tgamma(p0 + 1)));
432 }
434 {
435 c0 = 2.0 * (p0 + 1.0) /
436 ((2.0 * p0 + 1.0) * p0 * (ap0 * std::tgamma(p0 + 1)) *
437 (ap0 * std::tgamma(p0 + 1)));
438 }
440 {
441 c0 =
442 -2.0 / ((2.0 * p0 + 1.0) * (ap0 * std::tgamma(p0 + 1)) *
443 (ap0 * std::tgamma(p0 + 1)));
444 }
446 {
447 c0 = 10000000000000000.0;
448 }
449
450 NekDouble etap0 = 0.5 * c0 * (2.0 * p0 + 1.0) *
451 (ap0 * std::tgamma(p0 + 1)) *
452 (ap0 * std::tgamma(p0 + 1));
453
454 NekDouble overeta0 = 1.0 / (1.0 + etap0);
455
456 // Derivative of the Legendre polynomials
457 // dLp = derivative of the Legendre polynomial of order p
458 // dLpp = derivative of the Legendre polynomial of order p+1
459 // dLpm = derivative of the Legendre polynomial of order p-1
460 Polylib::jacobd(nquad0, z0.data(), &(dLp0[0]), p0, 0.0, 0.0);
461 Polylib::jacobd(nquad0, z0.data(), &(dLpp0[0]), p0 + 1, 0.0,
462 0.0);
463 Polylib::jacobd(nquad0, z0.data(), &(dLpm0[0]), p0 - 1, 0.0,
464 0.0);
465
466 // Building the DG_c_Left
467 for (i = 0; i < nquad0; ++i)
468 {
469 m_dGL_xi1[n][i] = etap0 * dLpm0[i];
470 m_dGL_xi1[n][i] += dLpp0[i];
471 m_dGL_xi1[n][i] *= overeta0;
474 }
475
476 // Building the DG_c_Right
477 for (i = 0; i < nquad0; ++i)
478 {
479 m_dGR_xi1[n][i] = etap0 * dLpm0[i];
480 m_dGR_xi1[n][i] += dLpp0[i];
481 m_dGR_xi1[n][i] *= overeta0;
482 m_dGR_xi1[n][i] += dLp0[i];
484 }
485 }
486 break;
487 }
488 case 2:
489 {
490 LibUtilities::BasisSharedPtr BasisFR_Left0;
491 LibUtilities::BasisSharedPtr BasisFR_Right0;
492 LibUtilities::BasisSharedPtr BasisFR_Left1;
493 LibUtilities::BasisSharedPtr BasisFR_Right1;
494
495 m_dGL_xi1 = Array<OneD, Array<OneD, NekDouble>>(nElements);
496 m_dGR_xi1 = Array<OneD, Array<OneD, NekDouble>>(nElements);
497 m_dGL_xi2 = Array<OneD, Array<OneD, NekDouble>>(nElements);
498 m_dGR_xi2 = Array<OneD, Array<OneD, NekDouble>>(nElements);
499
500 for (n = 0; n < nElements; ++n)
501 {
502 base = pFields[0]->GetExp(n)->GetBase();
503 nquad0 = base[0]->GetNumPoints();
504 nquad1 = base[1]->GetNumPoints();
505 nmodes0 = base[0]->GetNumModes();
506 nmodes1 = base[1]->GetNumModes();
507
508 Array<OneD, const NekDouble> z0;
509 Array<OneD, const NekDouble> w0;
510 Array<OneD, const NekDouble> z1;
511 Array<OneD, const NekDouble> w1;
512
513 base[0]->GetZW(z0, w0);
514 base[1]->GetZW(z1, w1);
515
516 m_dGL_xi1[n] = Array<OneD, NekDouble>(nquad0);
517 m_dGR_xi1[n] = Array<OneD, NekDouble>(nquad0);
518 m_dGL_xi2[n] = Array<OneD, NekDouble>(nquad1);
519 m_dGR_xi2[n] = Array<OneD, NekDouble>(nquad1);
520
521 // Auxiliary vectors to build up the auxiliary
522 // derivatives of the Legendre polynomials
523 Array<OneD, NekDouble> dLp0(nquad0, 0.0);
524 Array<OneD, NekDouble> dLpp0(nquad0, 0.0);
525 Array<OneD, NekDouble> dLpm0(nquad0, 0.0);
526 Array<OneD, NekDouble> dLp1(nquad1, 0.0);
527 Array<OneD, NekDouble> dLpp1(nquad1, 0.0);
528 Array<OneD, NekDouble> dLpm1(nquad1, 0.0);
529
530 // Degree of the correction functions
531 int p0 = nmodes0 - 1;
532 int p1 = nmodes1 - 1;
533
534 // Function sign to compute left correction function
535 NekDouble sign0 = pow(-1.0, p0);
536 NekDouble sign1 = pow(-1.0, p1);
537
538 // Factorial factor to build the scheme
539 NekDouble ap0 =
540 std::tgamma(2 * p0 + 1) /
541 (pow(2.0, p0) * std::tgamma(p0 + 1) * std::tgamma(p0 + 1));
542
543 NekDouble ap1 =
544 std::tgamma(2 * p1 + 1) /
545 (pow(2.0, p1) * std::tgamma(p1 + 1) * std::tgamma(p1 + 1));
546
547 // Scalar parameter which recovers the FR schemes
549 {
550 c0 = 0.0;
551 c1 = 0.0;
552 }
554 {
555 c0 = 2.0 * p0 /
556 ((2.0 * p0 + 1.0) * (p0 + 1.0) *
557 (ap0 * std::tgamma(p0 + 1)) *
558 (ap0 * std::tgamma(p0 + 1)));
559
560 c1 = 2.0 * p1 /
561 ((2.0 * p1 + 1.0) * (p1 + 1.0) *
562 (ap1 * std::tgamma(p1 + 1)) *
563 (ap1 * std::tgamma(p1 + 1)));
564 }
566 {
567 c0 = 2.0 * (p0 + 1.0) /
568 ((2.0 * p0 + 1.0) * p0 * (ap0 * std::tgamma(p0 + 1)) *
569 (ap0 * std::tgamma(p0 + 1)));
570
571 c1 = 2.0 * (p1 + 1.0) /
572 ((2.0 * p1 + 1.0) * p1 * (ap1 * std::tgamma(p1 + 1)) *
573 (ap1 * std::tgamma(p1 + 1)));
574 }
576 {
577 c0 =
578 -2.0 / ((2.0 * p0 + 1.0) * (ap0 * std::tgamma(p0 + 1)) *
579 (ap0 * std::tgamma(p0 + 1)));
580
581 c1 =
582 -2.0 / ((2.0 * p1 + 1.0) * (ap1 * std::tgamma(p1 + 1)) *
583 (ap1 * std::tgamma(p1 + 1)));
584 }
586 {
587 c0 = 10000000000000000.0;
588 c1 = 10000000000000000.0;
589 }
590
591 NekDouble etap0 = 0.5 * c0 * (2.0 * p0 + 1.0) *
592 (ap0 * std::tgamma(p0 + 1)) *
593 (ap0 * std::tgamma(p0 + 1));
594
595 NekDouble etap1 = 0.5 * c1 * (2.0 * p1 + 1.0) *
596 (ap1 * std::tgamma(p1 + 1)) *
597 (ap1 * std::tgamma(p1 + 1));
598
599 NekDouble overeta0 = 1.0 / (1.0 + etap0);
600 NekDouble overeta1 = 1.0 / (1.0 + etap1);
601
602 // Derivative of the Legendre polynomials
603 // dLp = derivative of the Legendre polynomial of order p
604 // dLpp = derivative of the Legendre polynomial of order p+1
605 // dLpm = derivative of the Legendre polynomial of order p-1
606 Polylib::jacobd(nquad0, z0.data(), &(dLp0[0]), p0, 0.0, 0.0);
607 Polylib::jacobd(nquad0, z0.data(), &(dLpp0[0]), p0 + 1, 0.0,
608 0.0);
609 Polylib::jacobd(nquad0, z0.data(), &(dLpm0[0]), p0 - 1, 0.0,
610 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,
614 0.0);
615 Polylib::jacobd(nquad1, z1.data(), &(dLpm1[0]), p1 - 1, 0.0,
616 0.0);
617
618 // Building the DG_c_Left
619 for (i = 0; i < nquad0; ++i)
620 {
621 m_dGL_xi1[n][i] = etap0 * dLpm0[i];
622 m_dGL_xi1[n][i] += dLpp0[i];
623 m_dGL_xi1[n][i] *= overeta0;
626 }
627
628 // Building the DG_c_Left
629 for (i = 0; i < nquad1; ++i)
630 {
631 m_dGL_xi2[n][i] = etap1 * dLpm1[i];
632 m_dGL_xi2[n][i] += dLpp1[i];
633 m_dGL_xi2[n][i] *= overeta1;
636 }
637
638 // Building the DG_c_Right
639 for (i = 0; i < nquad0; ++i)
640 {
641 m_dGR_xi1[n][i] = etap0 * dLpm0[i];
642 m_dGR_xi1[n][i] += dLpp0[i];
643 m_dGR_xi1[n][i] *= overeta0;
644 m_dGR_xi1[n][i] += dLp0[i];
646 }
647
648 // Building the DG_c_Right
649 for (i = 0; i < nquad1; ++i)
650 {
651 m_dGR_xi2[n][i] = etap1 * dLpm1[i];
652 m_dGR_xi2[n][i] += dLpp1[i];
653 m_dGR_xi2[n][i] *= overeta1;
654 m_dGR_xi2[n][i] += dLp1[i];
656 }
657 }
658 break;
659 }
660 case 3:
661 {
662 m_dGL_xi1 = Array<OneD, Array<OneD, NekDouble>>(nElements);
663 m_dGR_xi1 = Array<OneD, Array<OneD, NekDouble>>(nElements);
664 m_dGL_xi2 = Array<OneD, Array<OneD, NekDouble>>(nElements);
665 m_dGR_xi2 = Array<OneD, Array<OneD, NekDouble>>(nElements);
666 m_dGL_xi3 = Array<OneD, Array<OneD, NekDouble>>(nElements);
667 m_dGR_xi3 = Array<OneD, Array<OneD, NekDouble>>(nElements);
668
669 for (n = 0; n < nElements; ++n)
670 {
671 base = pFields[0]->GetExp(n)->GetBase();
672 nquad0 = base[0]->GetNumPoints();
673 nquad1 = base[1]->GetNumPoints();
674 nquad2 = base[2]->GetNumPoints();
675 nmodes0 = base[0]->GetNumModes();
676 nmodes1 = base[1]->GetNumModes();
677 nmodes2 = base[2]->GetNumModes();
678
679 Array<OneD, const NekDouble> z0;
680 Array<OneD, const NekDouble> w0;
681 Array<OneD, const NekDouble> z1;
682 Array<OneD, const NekDouble> w1;
683 Array<OneD, const NekDouble> z2;
684 Array<OneD, const NekDouble> w2;
685
686 base[0]->GetZW(z0, w0);
687 base[1]->GetZW(z1, w1);
688 base[1]->GetZW(z2, w2);
689
690 m_dGL_xi1[n] = Array<OneD, NekDouble>(nquad0);
691 m_dGR_xi1[n] = Array<OneD, NekDouble>(nquad0);
692 m_dGL_xi2[n] = Array<OneD, NekDouble>(nquad1);
693 m_dGR_xi2[n] = Array<OneD, NekDouble>(nquad1);
694 m_dGL_xi3[n] = Array<OneD, NekDouble>(nquad2);
695 m_dGR_xi3[n] = Array<OneD, NekDouble>(nquad2);
696
697 // Auxiliary vectors to build up the auxiliary
698 // derivatives of the Legendre polynomials
699 Array<OneD, NekDouble> dLp0(nquad0, 0.0);
700 Array<OneD, NekDouble> dLpp0(nquad0, 0.0);
701 Array<OneD, NekDouble> dLpm0(nquad0, 0.0);
702 Array<OneD, NekDouble> dLp1(nquad1, 0.0);
703 Array<OneD, NekDouble> dLpp1(nquad1, 0.0);
704 Array<OneD, NekDouble> dLpm1(nquad1, 0.0);
705 Array<OneD, NekDouble> dLp2(nquad2, 0.0);
706 Array<OneD, NekDouble> dLpp2(nquad2, 0.0);
707 Array<OneD, NekDouble> dLpm2(nquad2, 0.0);
708
709 // Degree of the correction functions
710 int p0 = nmodes0 - 1;
711 int p1 = nmodes1 - 1;
712 int p2 = nmodes2 - 1;
713
714 // Function sign to compute left correction function
715 NekDouble sign0 = pow(-1.0, p0);
716 NekDouble sign1 = pow(-1.0, p1);
717 NekDouble sign2 = pow(-1.0, p2);
718
719 // Factorial factor to build the scheme
720 NekDouble ap0 =
721 std::tgamma(2 * p0 + 1) /
722 (pow(2.0, p0) * std::tgamma(p0 + 1) * std::tgamma(p0 + 1));
723
724 // Factorial factor to build the scheme
725 NekDouble ap1 =
726 std::tgamma(2 * p1 + 1) /
727 (pow(2.0, p1) * std::tgamma(p1 + 1) * std::tgamma(p1 + 1));
728
729 // Factorial factor to build the scheme
730 NekDouble ap2 =
731 std::tgamma(2 * p2 + 1) /
732 (pow(2.0, p2) * std::tgamma(p2 + 1) * std::tgamma(p2 + 1));
733
734 // Scalar parameter which recovers the FR schemes
736 {
737 c0 = 0.0;
738 c1 = 0.0;
739 c2 = 0.0;
740 }
742 {
743 c0 = 2.0 * p0 /
744 ((2.0 * p0 + 1.0) * (p0 + 1.0) *
745 (ap0 * std::tgamma(p0 + 1)) *
746 (ap0 * std::tgamma(p0 + 1)));
747
748 c1 = 2.0 * p1 /
749 ((2.0 * p1 + 1.0) * (p1 + 1.0) *
750 (ap1 * std::tgamma(p1 + 1)) *
751 (ap1 * std::tgamma(p1 + 1)));
752
753 c2 = 2.0 * p2 /
754 ((2.0 * p2 + 1.0) * (p2 + 1.0) *
755 (ap2 * std::tgamma(p2 + 1)) *
756 (ap2 * std::tgamma(p2 + 1)));
757 }
759 {
760 c0 = 2.0 * (p0 + 1.0) /
761 ((2.0 * p0 + 1.0) * p0 * (ap0 * std::tgamma(p0 + 1)) *
762 (ap0 * std::tgamma(p0 + 1)));
763
764 c1 = 2.0 * (p1 + 1.0) /
765 ((2.0 * p1 + 1.0) * p1 * (ap1 * std::tgamma(p1 + 1)) *
766 (ap1 * std::tgamma(p1 + 1)));
767
768 c2 = 2.0 * (p2 + 1.0) /
769 ((2.0 * p2 + 1.0) * p2 * (ap2 * std::tgamma(p2 + 1)) *
770 (ap2 * std::tgamma(p2 + 1)));
771 }
773 {
774 c0 =
775 -2.0 / ((2.0 * p0 + 1.0) * (ap0 * std::tgamma(p0 + 1)) *
776 (ap0 * std::tgamma(p0 + 1)));
777
778 c1 =
779 -2.0 / ((2.0 * p1 + 1.0) * (ap1 * std::tgamma(p1 + 1)) *
780 (ap1 * std::tgamma(p1 + 1)));
781
782 c2 =
783 -2.0 / ((2.0 * p2 + 1.0) * (ap2 * std::tgamma(p2 + 1)) *
784 (ap2 * std::tgamma(p2 + 1)));
785 }
787 {
788 c0 = 10000000000000000.0;
789 c1 = 10000000000000000.0;
790 c2 = 10000000000000000.0;
791 }
792
793 NekDouble etap0 = 0.5 * c0 * (2.0 * p0 + 1.0) *
794 (ap0 * std::tgamma(p0 + 1)) *
795 (ap0 * std::tgamma(p0 + 1));
796
797 NekDouble etap1 = 0.5 * c1 * (2.0 * p1 + 1.0) *
798 (ap1 * std::tgamma(p1 + 1)) *
799 (ap1 * std::tgamma(p1 + 1));
800
801 NekDouble etap2 = 0.5 * c2 * (2.0 * p2 + 1.0) *
802 (ap2 * std::tgamma(p2 + 1)) *
803 (ap2 * std::tgamma(p2 + 1));
804
805 NekDouble overeta0 = 1.0 / (1.0 + etap0);
806 NekDouble overeta1 = 1.0 / (1.0 + etap1);
807 NekDouble overeta2 = 1.0 / (1.0 + etap2);
808
809 // Derivative of the Legendre polynomials
810 // dLp = derivative of the Legendre polynomial of order p
811 // dLpp = derivative of the Legendre polynomial of order p+1
812 // dLpm = derivative of the Legendre polynomial of order p-1
813 Polylib::jacobd(nquad0, z0.data(), &(dLp0[0]), p0, 0.0, 0.0);
814 Polylib::jacobd(nquad0, z0.data(), &(dLpp0[0]), p0 + 1, 0.0,
815 0.0);
816 Polylib::jacobd(nquad0, z0.data(), &(dLpm0[0]), p0 - 1, 0.0,
817 0.0);
818
819 Polylib::jacobd(nquad1, z1.data(), &(dLp1[0]), p1, 0.0, 0.0);
820 Polylib::jacobd(nquad1, z1.data(), &(dLpp1[0]), p1 + 1, 0.0,
821 0.0);
822 Polylib::jacobd(nquad1, z1.data(), &(dLpm1[0]), p1 - 1, 0.0,
823 0.0);
824
825 Polylib::jacobd(nquad2, z2.data(), &(dLp2[0]), p2, 0.0, 0.0);
826 Polylib::jacobd(nquad2, z2.data(), &(dLpp2[0]), p2 + 1, 0.0,
827 0.0);
828 Polylib::jacobd(nquad2, z2.data(), &(dLpm2[0]), p2 - 1, 0.0,
829 0.0);
830
831 // Building the DG_c_Left
832 for (i = 0; i < nquad0; ++i)
833 {
834 m_dGL_xi1[n][i] = etap0 * dLpm0[i];
835 m_dGL_xi1[n][i] += dLpp0[i];
836 m_dGL_xi1[n][i] *= overeta0;
839 }
840
841 // Building the DG_c_Left
842 for (i = 0; i < nquad1; ++i)
843 {
844 m_dGL_xi2[n][i] = etap1 * dLpm1[i];
845 m_dGL_xi2[n][i] += dLpp1[i];
846 m_dGL_xi2[n][i] *= overeta1;
849 }
850
851 // Building the DG_c_Left
852 for (i = 0; i < nquad2; ++i)
853 {
854 m_dGL_xi3[n][i] = etap2 * dLpm2[i];
855 m_dGL_xi3[n][i] += dLpp2[i];
856 m_dGL_xi3[n][i] *= overeta2;
859 }
860
861 // Building the DG_c_Right
862 for (i = 0; i < nquad0; ++i)
863 {
864 m_dGR_xi1[n][i] = etap0 * dLpm0[i];
865 m_dGR_xi1[n][i] += dLpp0[i];
866 m_dGR_xi1[n][i] *= overeta0;
867 m_dGR_xi1[n][i] += dLp0[i];
869 }
870
871 // Building the DG_c_Right
872 for (i = 0; i < nquad1; ++i)
873 {
874 m_dGR_xi2[n][i] = etap1 * dLpm1[i];
875 m_dGR_xi2[n][i] += dLpp1[i];
876 m_dGR_xi2[n][i] *= overeta1;
877 m_dGR_xi2[n][i] += dLp1[i];
879 }
880
881 // Building the DG_c_Right
882 for (i = 0; i < nquad2; ++i)
883 {
884 m_dGR_xi3[n][i] = etap2 * dLpm2[i];
885 m_dGR_xi3[n][i] += dLpp2[i];
886 m_dGR_xi3[n][i] *= overeta2;
887 m_dGR_xi3[n][i] += dLp2[i];
889 }
890 }
891 break;
892 }
893 default:
894 {
896 break;
897 }
898 }
899}
Array< OneD, Array< OneD, NekDouble > > m_dGR_xi3
Definition: DiffusionLFRNS.h:96
Array< OneD, Array< OneD, NekDouble > > m_dGL_xi3
Definition: DiffusionLFRNS.h:95
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:1378
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().
◆ SetupMetrics()
|
private |
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 #GetTraceQFactors 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
-
pSession Pointer to session reader. pFields Pointer to fields.
- Todo:
- Add the metric terms for 3D Hexahedra.
Definition at line 208 of file DiffusionLFRNS.cpp.
211{
212 int i, n;
213 int nquad0, nquad1;
214 int phys_offset;
215 int nLocalSolutionPts;
216 int nElements = pFields[0]->GetExpSize();
217 int nDimensions = pFields[0]->GetCoordim(0);
218 int nSolutionPts = pFields[0]->GetTotPoints();
219 int nTracePts = pFields[0]->GetTrace()->GetTotPoints();
220
221 m_traceNormals = Array<OneD, Array<OneD, NekDouble>>(nDimensions);
222 for (i = 0; i < nDimensions; ++i)
223 {
224 m_traceNormals[i] = Array<OneD, NekDouble>(nTracePts, 0.0);
225 }
226 pFields[0]->GetTrace()->GetNormals(m_traceNormals);
227
228 Array<TwoD, const NekDouble> gmat;
229 Array<OneD, const NekDouble> jac;
230
231 m_jac = Array<OneD, NekDouble>(nSolutionPts);
232
233 Array<OneD, NekDouble> auxArray1;
234 Array<OneD, LibUtilities::BasisSharedPtr> base;
235 LibUtilities::PointsKeyVector ptsKeys;
236
237 switch (nDimensions)
238 {
239 case 1:
240 {
241 for (n = 0; n < nElements; ++n)
242 {
243 ptsKeys = pFields[0]->GetExp(n)->GetPointsKeys();
244 nLocalSolutionPts = pFields[0]->GetExp(n)->GetTotPoints();
245 phys_offset = pFields[0]->GetPhys_Offset(n);
246 jac = pFields[0]
247 ->GetExp(n)
248 ->as<LocalRegions::Expansion1D>()
249 ->GetGeom1D()
250 ->GetMetricInfo()
251 ->GetJac(ptsKeys);
252 for (i = 0; i < nLocalSolutionPts; ++i)
253 {
254 m_jac[i + phys_offset] = jac[0];
255 }
256 }
257 break;
258 }
259 case 2:
260 {
261 m_gmat = Array<OneD, Array<OneD, NekDouble>>(4);
262 m_gmat[0] = Array<OneD, NekDouble>(nSolutionPts);
263 m_gmat[1] = Array<OneD, NekDouble>(nSolutionPts);
264 m_gmat[2] = Array<OneD, NekDouble>(nSolutionPts);
265 m_gmat[3] = Array<OneD, NekDouble>(nSolutionPts);
266
267 m_Q2D_e0 = Array<OneD, Array<OneD, NekDouble>>(nElements);
268 m_Q2D_e1 = Array<OneD, Array<OneD, NekDouble>>(nElements);
269 m_Q2D_e2 = Array<OneD, Array<OneD, NekDouble>>(nElements);
270 m_Q2D_e3 = Array<OneD, Array<OneD, NekDouble>>(nElements);
271
272 for (n = 0; n < nElements; ++n)
273 {
274 base = pFields[0]->GetExp(n)->GetBase();
275 nquad0 = base[0]->GetNumPoints();
276 nquad1 = base[1]->GetNumPoints();
277
278 m_Q2D_e0[n] = Array<OneD, NekDouble>(nquad0);
279 m_Q2D_e1[n] = Array<OneD, NekDouble>(nquad1);
280 m_Q2D_e2[n] = Array<OneD, NekDouble>(nquad0);
281 m_Q2D_e3[n] = Array<OneD, NekDouble>(nquad1);
282
283 // Extract the Q factors at each edge point
284 pFields[0]->GetExp(n)->GetTraceQFactors(0, auxArray1 =
285 m_Q2D_e0[n]);
286 pFields[0]->GetExp(n)->GetTraceQFactors(1, auxArray1 =
287 m_Q2D_e1[n]);
288 pFields[0]->GetExp(n)->GetTraceQFactors(2, auxArray1 =
289 m_Q2D_e2[n]);
290 pFields[0]->GetExp(n)->GetTraceQFactors(3, auxArray1 =
291 m_Q2D_e3[n]);
292
293 ptsKeys = pFields[0]->GetExp(n)->GetPointsKeys();
294 nLocalSolutionPts = pFields[0]->GetExp(n)->GetTotPoints();
295 phys_offset = pFields[0]->GetPhys_Offset(n);
296
297 jac = pFields[0]
298 ->GetExp(n)
299 ->as<LocalRegions::Expansion2D>()
300 ->GetGeom2D()
301 ->GetMetricInfo()
302 ->GetJac(ptsKeys);
303 gmat = pFields[0]
304 ->GetExp(n)
305 ->as<LocalRegions::Expansion2D>()
306 ->GetGeom2D()
307 ->GetMetricInfo()
308 ->GetDerivFactors(ptsKeys);
309
310 if (pFields[0]
311 ->GetExp(n)
312 ->as<LocalRegions::Expansion2D>()
313 ->GetGeom2D()
314 ->GetMetricInfo()
315 ->GetGtype() == SpatialDomains::eDeformed)
316 {
317 for (i = 0; i < nLocalSolutionPts; ++i)
318 {
319 m_jac[i + phys_offset] = jac[i];
320 m_gmat[0][i + phys_offset] = gmat[0][i];
321 m_gmat[1][i + phys_offset] = gmat[1][i];
322 m_gmat[2][i + phys_offset] = gmat[2][i];
323 m_gmat[3][i + phys_offset] = gmat[3][i];
324 }
325 }
326 else
327 {
328 for (i = 0; i < nLocalSolutionPts; ++i)
329 {
330 m_jac[i + phys_offset] = jac[0];
331 m_gmat[0][i + phys_offset] = gmat[0][0];
332 m_gmat[1][i + phys_offset] = gmat[1][0];
333 m_gmat[2][i + phys_offset] = gmat[2][0];
334 m_gmat[3][i + phys_offset] = gmat[3][0];
335 }
336 }
337 }
338 break;
339 }
340 case 3:
341 {
343 break;
344 }
345 default:
346 {
348 break;
349 }
350 }
351}
Array< OneD, NekDouble > m_jac
Definition: DiffusionLFRNS.h:83
Array< OneD, Array< OneD, NekDouble > > m_gmat
Definition: DiffusionLFRNS.h:84
References ASSERTL0, Nektar::SpatialDomains::eDeformed, Nektar::LocalRegions::Expansion::GetMetricInfo(), m_gmat, m_jac, m_Q2D_e0, m_Q2D_e1, m_Q2D_e2, m_Q2D_e3, and m_traceNormals.
Referenced by v_InitObject().
◆ v_Diffuse()
|
overrideprotectedvirtual |
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 909 of file DiffusionLFRNS.cpp.
916{
917 int i, j, n;
918 int phys_offset;
919
920 Array<OneD, NekDouble> auxArray1, auxArray2;
921
922 Array<OneD, LibUtilities::BasisSharedPtr> Basis;
923 Basis = fields[0]->GetExp(0)->GetBase();
924
925 int nElements = fields[0]->GetExpSize();
926 int nDim = fields[0]->GetCoordim(0);
927 int nScalars = inarray.size();
928 int nSolutionPts = fields[0]->GetTotPoints();
929 int nCoeffs = fields[0]->GetNcoeffs();
930
931 Array<OneD, Array<OneD, NekDouble>> outarrayCoeff(nConvectiveFields);
932 for (i = 0; i < nConvectiveFields; ++i)
933 {
934 outarrayCoeff[i] = Array<OneD, NekDouble>(nCoeffs);
935 }
936
937 // Compute interface numerical fluxes for inarray in physical space
939
940 switch (nDim)
941 {
942 // 1D problems
943 case 1:
944 {
945 for (i = 0; i < nScalars; ++i)
946 {
947 // Computing the physical first-order discountinuous
948 // derivative
949 for (n = 0; n < nElements; n++)
950 {
951 phys_offset = fields[0]->GetPhys_Offset(n);
952
953 fields[i]->GetExp(n)->PhysDeriv(
954 0, auxArray1 = inarray[i] + phys_offset,
955 auxArray2 = m_DU1[i][0] + phys_offset);
956 }
957
958 // Computing the standard first-order correction
959 // derivative
961 m_DFC1[i][0]);
962
963 // Back to the physical space using global operations
965 &m_DFC1[i][0][0], 1);
967 &m_DFC1[i][0][0], 1);
968
969 // Computing total first order derivatives
971 1, &m_D1[i][0][0], 1);
972
974 1);
975 }
976
977 // Computing the viscous tensor
979
980 // Compute u from q_{\eta} and q_{\xi}
981 // Obtain numerical fluxes
983
984 for (i = 0; i < nConvectiveFields; ++i)
985 {
986 // Computing the physical second-order discountinuous
987 // derivative
988 for (n = 0; n < nElements; n++)
989 {
990 phys_offset = fields[0]->GetPhys_Offset(n);
991
992 fields[i]->GetExp(n)->PhysDeriv(
993 0, auxArray1 = m_viscTensor[0][i] + phys_offset,
994 auxArray2 = m_DD1[i][0] + phys_offset);
995 }
996
997 // Computing the standard second-order correction
998 // derivative
1001
1002 // Back to the physical space using global operations
1004 &m_DFC2[i][0][0], 1);
1006 &m_DFC2[i][0][0], 1);
1007
1008 // Adding the total divergence to outarray (RHS)
1010 1, &outarray[i][0], 1);
1011
1012 // Primitive Dealiasing 1D
1013 if (!(Basis[0]->Collocation()))
1014 {
1015 fields[i]->FwdTrans(outarray[i], outarrayCoeff[i]);
1016 fields[i]->BwdTrans(outarrayCoeff[i], outarray[i]);
1017 }
1018 }
1019 break;
1020 }
1021 // 2D problems
1022 case 2:
1023 {
1024 for (i = 0; i < nScalars; ++i)
1025 {
1026 for (j = 0; j < nDim; ++j)
1027 {
1028 // Temporary vectors
1029 Array<OneD, NekDouble> u1_hat(nSolutionPts, 0.0);
1030 Array<OneD, NekDouble> u2_hat(nSolutionPts, 0.0);
1031
1032 if (j == 0)
1033 {
1034 Vmath::Vmul(nSolutionPts, &inarray[i][0], 1,
1035 &m_gmat[0][0], 1, &u1_hat[0], 1);
1036
1038 &u1_hat[0], 1);
1039
1040 Vmath::Vmul(nSolutionPts, &inarray[i][0], 1,
1041 &m_gmat[1][0], 1, &u2_hat[0], 1);
1042
1044 &u2_hat[0], 1);
1045 }
1046 else if (j == 1)
1047 {
1048 Vmath::Vmul(nSolutionPts, &inarray[i][0], 1,
1049 &m_gmat[2][0], 1, &u1_hat[0], 1);
1050
1052 &u1_hat[0], 1);
1053
1054 Vmath::Vmul(nSolutionPts, &inarray[i][0], 1,
1055 &m_gmat[3][0], 1, &u2_hat[0], 1);
1056
1058 &u2_hat[0], 1);
1059 }
1060
1061 for (n = 0; n < nElements; n++)
1062 {
1063 phys_offset = fields[0]->GetPhys_Offset(n);
1064
1065 fields[i]->GetExp(n)->StdPhysDeriv(
1066 0, auxArray1 = u1_hat + phys_offset,
1067 auxArray2 = m_tmp1[i][j] + phys_offset);
1068
1069 fields[i]->GetExp(n)->StdPhysDeriv(
1070 1, auxArray1 = u2_hat + phys_offset,
1071 auxArray2 = m_tmp2[i][j] + phys_offset);
1072 }
1073
1076
1077 // Divide by the metric jacobian
1079 &m_DU1[i][j][0], 1);
1080
1081 // Computing the standard first-order correction
1082 // derivatives
1083 DerCFlux_2D(nConvectiveFields, j, fields, inarray[i],
1085 }
1086
1087 // Multiplying derivatives by B matrix to get auxiliary
1088 // variables
1089 for (j = 0; j < nSolutionPts; ++j)
1090 {
1091 // std(q1)
1094 m_DFC1[i][0][j] +
1097 m_DFC1[i][1][j];
1098
1099 // std(q2)
1102 m_DFC1[i][0][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)
1115
1116 // q2 = A21^(-1)*std(q1) + A22^(-1)*std(q2)
1119 }
1120
1121 // Computing the physical first-order derivatives
1122 for (j = 0; j < nDim; ++j)
1123 {
1126 }
1127 } // Close loop on nScalars
1128
1129 // For 3D Homogeneous 1D only take derivatives in 3rd direction
1131 {
1132 for (i = 0; i < nScalars; ++i)
1133 {
1135 }
1136 }
1137
1138 // Computing the viscous tensor
1140
1141 // Compute u from q_{\eta} and q_{\xi}
1142 // Obtain numerical fluxes
1144
1145 // Computing the standard second-order discontinuous
1146 // derivatives
1147 for (i = 0; i < nConvectiveFields; ++i)
1148 {
1149 // Temporary vectors
1150 Array<OneD, NekDouble> f_hat(nSolutionPts, 0.0);
1151 Array<OneD, NekDouble> g_hat(nSolutionPts, 0.0);
1152
1153 for (j = 0; j < nSolutionPts; j++)
1154 {
1157 m_jac[j];
1158
1161 m_jac[j];
1162 }
1163
1164 for (n = 0; n < nElements; n++)
1165 {
1166 phys_offset = fields[0]->GetPhys_Offset(n);
1167
1168 fields[0]->GetExp(n)->StdPhysDeriv(
1169 0, auxArray1 = f_hat + phys_offset,
1170 auxArray2 = m_DD1[i][0] + phys_offset);
1171
1172 fields[0]->GetExp(n)->StdPhysDeriv(
1173 1, auxArray1 = g_hat + phys_offset,
1174 auxArray2 = m_DD1[i][1] + phys_offset);
1175 }
1176
1177 // Divergence of the standard discontinuous flux
1179 1, &m_divFD[i][0], 1);
1180
1181 // Divergence of the standard correction flux
1182 if (Basis[0]->GetPointsType() ==
1184 Basis[1]->GetPointsType() ==
1186 {
1187 DivCFlux_2D_Gauss(nConvectiveFields, fields, f_hat, g_hat,
1189 }
1190 else
1191 {
1194 }
1195
1196 // Divergence of the standard final flux
1198 &outarray[i][0], 1);
1199
1200 // Dividing by the jacobian to get back into
1201 // physical space
1203 &outarray[i][0], 1);
1204
1205 // Primitive Dealiasing 2D
1206 if (!(Basis[0]->Collocation()))
1207 {
1208 fields[i]->FwdTrans(outarray[i], outarrayCoeff[i]);
1209 fields[i]->BwdTrans(outarrayCoeff[i], outarray[i]);
1210 }
1211 }
1212 break;
1213 }
1214 // 3D problems
1215 case 3:
1216 {
1218 break;
1219 }
1220 }
1221}
DiffusionFluxVecCBNS m_fluxVectorNS
Definition: Diffusion.h:347
Array< OneD, Array< OneD, NekDouble > > m_divFD
Definition: DiffusionLFRNS.h:118
void 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:1650
void 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:1227
Array< OneD, Array< OneD, NekDouble > > m_viscFlux
Definition: DiffusionLFRNS.h:116
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > m_tmp2
Definition: DiffusionLFRNS.h:122
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > m_DFC2
Definition: DiffusionLFRNS.h:117
Array< OneD, Array< OneD, NekDouble > > m_divFC
Definition: DiffusionLFRNS.h:119
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > m_DFC1
Definition: DiffusionLFRNS.h:111
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > m_viscTensor
Definition: DiffusionLFRNS.h:115
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > m_tmp1
Definition: DiffusionLFRNS.h:121
void 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:1940
void 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:1744
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > m_IF1
Definition: DiffusionLFRNS.h:109
void 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:2128
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > m_D1
Definition: DiffusionLFRNS.h:113
Array< OneD, Array< OneD, NekDouble > > m_homoDerivs
Definition: DiffusionLFRNS.h:124
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > m_BD1
Definition: DiffusionLFRNS.h:112
void 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:1497
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > m_DU1
Definition: DiffusionLFRNS.h:110
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > m_DD1
Definition: DiffusionLFRNS.h:114
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.hpp:126
References ASSERTL0, DerCFlux_1D(), DerCFlux_2D(), DivCFlux_2D(), DivCFlux_2D_Gauss(), 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, NumericalFluxO1(), NumericalFluxO2(), Vmath::Vadd(), Vmath::Vcopy(), Vmath::Vdiv(), and Vmath::Vmul().
◆ v_GetFluxTensor()
|
inlineoverrideprotectedvirtual |
Reimplemented from Nektar::SolverUtils::Diffusion.
Definition at line 73 of file DiffusionLFRNS.h.
References m_viscTensor.
◆ v_InitObject()
|
overrideprotectedvirtual |
Initiliase DiffusionLFRNS objects and store them before starting the time-stepping.
This routine calls the functions SetupMetrics and SetupCFunctions to initialise the objects needed by DiffusionLFRNS.
- Parameters
-
pSession Pointer to session reader. pFields Pointer to fields.
Implements Nektar::SolverUtils::Diffusion.
Definition at line 80 of file DiffusionLFRNS.cpp.
83{
84 m_session = pSession;
91 0.0257);
94 SetupMetrics(pSession, pFields);
95 SetupCFunctions(pSession, pFields);
96
97 // Initialising arrays
98 int i, j;
99 int nConvectiveFields = pFields.size();
100 int nScalars = nConvectiveFields - 1;
101 int nDim = pFields[0]->GetCoordim(0);
102 int nSolutionPts = pFields[0]->GetTotPoints();
103 int nTracePts = pFields[0]->GetTrace()->GetTotPoints();
104
105 m_spaceDim = nDim;
106 if (pSession->DefinesSolverInfo("HOMOGENEOUS"))
107 {
108 m_spaceDim = 3;
109 }
110
112
114
116 {
117 m_traceVel[i] = Array<OneD, NekDouble>(nTracePts, 0.0);
118 }
119
120 m_IF1 = Array<OneD, Array<OneD, Array<OneD, NekDouble>>>(nScalars);
121 m_DU1 = Array<OneD, Array<OneD, Array<OneD, NekDouble>>>(nScalars);
122 m_DFC1 = Array<OneD, Array<OneD, Array<OneD, NekDouble>>>(nScalars);
123 m_tmp1 = Array<OneD, Array<OneD, Array<OneD, NekDouble>>>(nScalars);
124 m_tmp2 = Array<OneD, Array<OneD, Array<OneD, NekDouble>>>(nScalars);
125 m_BD1 = Array<OneD, Array<OneD, Array<OneD, NekDouble>>>(nScalars);
126
127 m_DFC2 =
128 Array<OneD, Array<OneD, Array<OneD, NekDouble>>>(nConvectiveFields);
129 m_DD1 = Array<OneD, Array<OneD, Array<OneD, NekDouble>>>(nConvectiveFields);
130 m_viscFlux = Array<OneD, Array<OneD, NekDouble>>(nConvectiveFields);
131 m_divFD = Array<OneD, Array<OneD, NekDouble>>(nConvectiveFields);
132 m_divFC = Array<OneD, Array<OneD, NekDouble>>(nConvectiveFields);
133
136 for (i = 0; i < nScalars; ++i)
137 {
139 m_DU1[i] = Array<OneD, Array<OneD, NekDouble>>(nDim);
140 m_DFC1[i] = Array<OneD, Array<OneD, NekDouble>>(nDim);
141 m_tmp1[i] = Array<OneD, Array<OneD, NekDouble>>(nDim);
142 m_tmp2[i] = Array<OneD, Array<OneD, NekDouble>>(nDim);
143 m_BD1[i] = Array<OneD, Array<OneD, NekDouble>>(nDim);
144
145 for (j = 0; j < nDim; ++j)
146 {
147 m_IF1[i][j] = Array<OneD, NekDouble>(nTracePts, 0.0);
148 m_DU1[i][j] = Array<OneD, NekDouble>(nSolutionPts, 0.0);
149 m_DFC1[i][j] = Array<OneD, NekDouble>(nSolutionPts, 0.0);
150 m_tmp1[i][j] = Array<OneD, NekDouble>(nSolutionPts, 0.0);
151 m_tmp2[i][j] = Array<OneD, NekDouble>(nSolutionPts, 0.0);
152 m_BD1[i][j] = Array<OneD, NekDouble>(nSolutionPts, 0.0);
153 }
154 }
155
156 for (i = 0; i < nConvectiveFields; ++i)
157 {
158 m_DFC2[i] = Array<OneD, Array<OneD, NekDouble>>(nDim);
159 m_DD1[i] = Array<OneD, Array<OneD, NekDouble>>(nDim);
160 m_viscFlux[i] = Array<OneD, NekDouble>(nTracePts, 0.0);
161 m_divFD[i] = Array<OneD, NekDouble>(nSolutionPts, 0.0);
162 m_divFC[i] = Array<OneD, NekDouble>(nSolutionPts, 0.0);
163
164 for (j = 0; j < nDim; ++j)
165 {
166 m_DFC2[i][j] = Array<OneD, NekDouble>(nSolutionPts, 0.0);
167 m_DD1[i][j] = Array<OneD, NekDouble>(nSolutionPts, 0.0);
168 }
169 }
170
172 {
173 m_D1[i] = Array<OneD, Array<OneD, NekDouble>>(nScalars);
174
175 for (j = 0; j < nScalars; ++j)
176 {
177 m_D1[i][j] = Array<OneD, NekDouble>(nSolutionPts, 0.0);
178 }
179 }
180
182 {
183 m_viscTensor[i] = Array<OneD, Array<OneD, NekDouble>>(nScalars + 1);
184
185 for (j = 0; j < nScalars + 1; ++j)
186 {
187 m_viscTensor[i][j] = Array<OneD, NekDouble>(nSolutionPts, 0.0);
188 }
189 }
190}
void 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:208
NekDouble m_thermalConductivity
Definition: DiffusionLFRNS.h:105
void 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...
Definition: DiffusionLFRNS.cpp:369
LibUtilities::SessionReaderSharedPtr m_session
Definition: DiffusionLFRNS.h:81
NekDouble m_gasConstant
Definition: DiffusionLFRNS.h:101
std::string m_ViscosityType
Definition: DiffusionLFRNS.h:103
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, SetupCFunctions(), and SetupMetrics().
◆ v_SetHomoDerivs()
|
inlineoverrideprotectedvirtual |
Reimplemented from Nektar::SolverUtils::Diffusion.
Definition at line 68 of file DiffusionLFRNS.h.
References m_homoDerivs.
◆ WeakPenaltyO1()
|
private |
Imposes appropriate bcs for the 1st order derivatives.
Definition at line 1282 of file DiffusionLFRNS.cpp.
1286{
1287 int cnt;
1288 int i, j, e;
1289 int id1, id2;
1290
1291 int nBndEdgePts, nBndEdges, nBndRegions;
1292
1293 int nTracePts = fields[0]->GetTrace()->GetTotPoints();
1294 int nScalars = inarray.size();
1295
1296 Array<OneD, NekDouble> tmp1(nTracePts, 0.0);
1297 Array<OneD, NekDouble> tmp2(nTracePts, 0.0);
1298 Array<OneD, NekDouble> Tw(nTracePts, m_Twall);
1299
1300 Array<OneD, Array<OneD, NekDouble>> scalarVariables(nScalars);
1301 Array<OneD, Array<OneD, NekDouble>> uplus(nScalars);
1302
1303 // Extract internal values of the scalar variables for Neumann bcs
1304 for (i = 0; i < nScalars; ++i)
1305 {
1306 scalarVariables[i] = Array<OneD, NekDouble>(nTracePts, 0.0);
1307
1308 uplus[i] = Array<OneD, NekDouble>(nTracePts, 0.0);
1309 fields[i]->ExtractTracePhys(inarray[i], uplus[i]);
1310 }
1311
1312 // Compute boundary conditions for velocities
1313 for (i = 0; i < nScalars - 1; ++i)
1314 {
1315 // Note that cnt has to loop on nBndRegions and nBndEdges
1316 // and has to be reset to zero per each equation
1317 cnt = 0;
1318 nBndRegions = fields[i + 1]->GetBndCondExpansions().size();
1319 for (j = 0; j < nBndRegions; ++j)
1320 {
1321 if (fields[i + 1]
1322 ->GetBndConditions()[j]
1323 ->GetBoundaryConditionType() == SpatialDomains::ePeriodic)
1324 {
1325 continue;
1326 }
1327
1328 nBndEdges = fields[i + 1]->GetBndCondExpansions()[j]->GetExpSize();
1329 for (e = 0; e < nBndEdges; ++e)
1330 {
1331 nBndEdgePts = fields[i + 1]
1332 ->GetBndCondExpansions()[j]
1333 ->GetExp(e)
1334 ->GetTotPoints();
1335
1336 id1 =
1337 fields[i + 1]->GetBndCondExpansions()[j]->GetPhys_Offset(e);
1338
1339 id2 = fields[0]->GetTrace()->GetPhys_Offset(
1340 fields[0]->GetTraceMap()->GetBndCondIDToGlobalTraceID(
1341 cnt++));
1342
1343 // Reinforcing bcs for velocity in case of Wall bcs
1344 if (boost::iequals(
1345 fields[i]->GetBndConditions()[j]->GetUserDefined(),
1346 "WallViscous") ||
1347 boost::iequals(
1348 fields[i]->GetBndConditions()[j]->GetUserDefined(),
1349 "WallAdiabatic"))
1350 {
1351 Vmath::Zero(nBndEdgePts, &scalarVariables[i][id2], 1);
1352 }
1353
1354 // Imposing velocity bcs if not Wall
1355 else if (fields[i]
1356 ->GetBndConditions()[j]
1357 ->GetBoundaryConditionType() ==
1359 {
1360 Vmath::Vdiv(nBndEdgePts,
1361 &(fields[i + 1]
1362 ->GetBndCondExpansions()[j]
1363 ->UpdatePhys())[id1],
1364 1,
1365 &(fields[0]
1366 ->GetBndCondExpansions()[j]
1367 ->UpdatePhys())[id1],
1368 1, &scalarVariables[i][id2], 1);
1369 }
1370
1371 // For Dirichlet boundary condition: uflux = u_bcs
1372 if (fields[i]
1373 ->GetBndConditions()[j]
1374 ->GetBoundaryConditionType() ==
1376 {
1377 Vmath::Vcopy(nBndEdgePts, &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, &uplus[i][id2], 1,
1387 &penaltyfluxO1[i][id2], 1);
1388 }
1389
1390 // Building kinetic energy to be used for T bcs
1391 Vmath::Vmul(nBndEdgePts, &scalarVariables[i][id2], 1,
1392 &scalarVariables[i][id2], 1, &tmp1[id2], 1);
1393
1394 Vmath::Smul(nBndEdgePts, 0.5, &tmp1[id2], 1, &tmp1[id2], 1);
1395
1396 Vmath::Vadd(nBndEdgePts, &tmp2[id2], 1, &tmp1[id2], 1,
1397 &tmp2[id2], 1);
1398 }
1399 }
1400 }
1401
1402 // Compute boundary conditions for temperature
1403 cnt = 0;
1404 nBndRegions = fields[nScalars]->GetBndCondExpansions().size();
1405 for (j = 0; j < nBndRegions; ++j)
1406 {
1407 nBndEdges = fields[nScalars]->GetBndCondExpansions()[j]->GetExpSize();
1408
1409 if (fields[nScalars]
1410 ->GetBndConditions()[j]
1411 ->GetBoundaryConditionType() == SpatialDomains::ePeriodic)
1412 {
1413 continue;
1414 }
1415
1416 for (e = 0; e < nBndEdges; ++e)
1417 {
1418 nBndEdgePts = fields[nScalars]
1419 ->GetBndCondExpansions()[j]
1420 ->GetExp(e)
1421 ->GetTotPoints();
1422
1423 id1 =
1424 fields[nScalars]->GetBndCondExpansions()[j]->GetPhys_Offset(e);
1425
1426 id2 = fields[0]->GetTrace()->GetPhys_Offset(
1427 fields[0]->GetTraceMap()->GetBndCondIDToGlobalTraceID(cnt++));
1428
1429 // Imposing Temperature Twall at the wall
1430 if (boost::iequals(
1431 fields[i]->GetBndConditions()[j]->GetUserDefined(),
1432 "WallViscous"))
1433 {
1434 Vmath::Vcopy(nBndEdgePts, &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]
1440 ->GetBndConditions()[j]
1441 ->GetBoundaryConditionType() ==
1443 {
1444 // Divide E by rho
1445 Vmath::Vdiv(
1446 nBndEdgePts,
1447 &(fields[nScalars]
1448 ->GetBndCondExpansions()[j]
1449 ->GetPhys())[id1],
1450 1, &(fields[0]->GetBndCondExpansions()[j]->GetPhys())[id1],
1451 1, &scalarVariables[nScalars - 1][id2], 1);
1452
1453 // Subtract kinetic energy to E/rho
1454 Vmath::Vsub(nBndEdgePts, &scalarVariables[nScalars - 1][id2], 1,
1455 &tmp2[id2], 1, &scalarVariables[nScalars - 1][id2],
1456 1);
1457
1458 // Multiply by constant factor (gamma-1)/R
1460 &scalarVariables[nScalars - 1][id2], 1,
1461 &scalarVariables[nScalars - 1][id2], 1);
1462 }
1463
1464 // For Dirichlet boundary condition: uflux = u_bcs
1465 if (fields[nScalars]
1466 ->GetBndConditions()[j]
1467 ->GetBoundaryConditionType() ==
1468 SpatialDomains::eDirichlet &&
1469 !boost::iequals(
1470 fields[nScalars]->GetBndConditions()[j]->GetUserDefined(),
1471 "WallAdiabatic"))
1472 {
1473 Vmath::Vcopy(nBndEdgePts, &scalarVariables[nScalars - 1][id2],
1474 1, &penaltyfluxO1[nScalars - 1][id2], 1);
1475 }
1476
1477 // For Neumann boundary condition: uflux = u_+
1478 else if (((fields[nScalars]->GetBndConditions()[j])
1479 ->GetBoundaryConditionType() ==
1480 SpatialDomains::eNeumann) ||
1481 boost::iequals(fields[nScalars]
1482 ->GetBndConditions()[j]
1483 ->GetUserDefined(),
1484 "WallAdiabatic"))
1485 {
1486 Vmath::Vcopy(nBndEdgePts, &uplus[nScalars - 1][id2], 1,
1487 &penaltyfluxO1[nScalars - 1][id2], 1);
1488 }
1489 }
1490 }
1491}
References Nektar::SpatialDomains::eDirichlet, Nektar::SpatialDomains::eNeumann, Nektar::SpatialDomains::ePeriodic, m_gamma, m_gasConstant, m_Twall, Vmath::Smul(), Vmath::Vadd(), Vmath::Vcopy(), Vmath::Vdiv(), Vmath::Vmul(), Vmath::Vsub(), and Vmath::Zero().
Referenced by NumericalFluxO1().
◆ WeakPenaltyO2()
|
private |
Imposes appropriate bcs for the 2nd order derivatives.
Definition at line 1558 of file DiffusionLFRNS.cpp.
1562{
1563 int cnt = 0;
1564 int nBndEdges, nBndEdgePts;
1565 int i, e;
1566 int id2;
1567
1568 int nTracePts = fields[0]->GetTrace()->GetTotPoints();
1569 int nBndRegions = fields[var]->GetBndCondExpansions().size();
1570
1571 Array<OneD, NekDouble> uterm(nTracePts);
1572 Array<OneD, NekDouble> qtemp(nTracePts);
1573
1574 // Extract the physical values of the solution at the boundaries
1575 fields[var]->ExtractTracePhys(qfield, qtemp);
1576
1577 // Loop on the boundary regions to apply appropriate bcs
1578 for (i = 0; i < nBndRegions; ++i)
1579 {
1580 // Number of boundary regions related to region 'i'
1581 nBndEdges = fields[var]->GetBndCondExpansions()[i]->GetExpSize();
1582
1583 if (fields[var]->GetBndConditions()[i]->GetBoundaryConditionType() ==
1585 {
1586 continue;
1587 }
1588
1589 // Weakly impose bcs by modifying flux values
1590 for (e = 0; e < nBndEdges; ++e)
1591 {
1592 nBndEdgePts = fields[var]
1593 ->GetBndCondExpansions()[i]
1594 ->GetExp(e)
1595 ->GetTotPoints();
1596
1597 id2 = fields[0]->GetTrace()->GetPhys_Offset(
1598 fields[0]->GetTraceMap()->GetBndCondIDToGlobalTraceID(cnt++));
1599
1600 // In case of Dirichlet bcs:
1601 // uflux = gD
1602 if (fields[var]
1603 ->GetBndConditions()[i]
1604 ->GetBoundaryConditionType() ==
1605 SpatialDomains::eDirichlet &&
1606 !boost::iequals(
1607 fields[var]->GetBndConditions()[i]->GetUserDefined(),
1608 "WallAdiabatic"))
1609 {
1611 &qtemp[id2], 1, &penaltyflux[id2], 1);
1612 }
1613 // 3.4) In case of Neumann bcs:
1614 // uflux = u+
1615 else if ((fields[var]->GetBndConditions()[i])
1616 ->GetBoundaryConditionType() ==
1617 SpatialDomains::eNeumann)
1618 {
1620 }
1621 else if (boost::iequals(
1622 fields[var]->GetBndConditions()[i]->GetUserDefined(),
1623 "WallAdiabatic"))
1624 {
1626 {
1627 Vmath::Zero(nBndEdgePts, &penaltyflux[id2], 1);
1628 }
1629 else
1630 {
1631
1633 &qtemp[id2], 1, &penaltyflux[id2], 1);
1634 }
1635 }
1636 }
1637 }
1638}
References ASSERTL0, Nektar::SpatialDomains::eDirichlet, Nektar::SpatialDomains::eNeumann, Nektar::SpatialDomains::ePeriodic, m_spaceDim, m_traceNormals, Vmath::Vmul(), and Vmath::Zero().
Referenced by NumericalFluxO2().
Member Data Documentation
◆ m_BD1
|
private |
Definition at line 112 of file DiffusionLFRNS.h.
Referenced by v_Diffuse(), and v_InitObject().
◆ m_D1
|
private |
Definition at line 113 of file DiffusionLFRNS.h.
Referenced by v_Diffuse(), and v_InitObject().
◆ m_DD1
|
private |
Definition at line 114 of file DiffusionLFRNS.h.
Referenced by v_Diffuse(), and v_InitObject().
◆ m_DFC1
|
private |
Definition at line 111 of file DiffusionLFRNS.h.
Referenced by v_Diffuse(), and v_InitObject().
◆ m_DFC2
|
private |
Definition at line 117 of file DiffusionLFRNS.h.
Referenced by v_Diffuse(), and v_InitObject().
◆ m_dGL_xi1
Definition at line 91 of file DiffusionLFRNS.h.
Referenced by DerCFlux_1D(), DerCFlux_2D(), DivCFlux_2D(), DivCFlux_2D_Gauss(), and SetupCFunctions().
◆ m_dGL_xi2
Definition at line 93 of file DiffusionLFRNS.h.
Referenced by DerCFlux_2D(), DivCFlux_2D(), DivCFlux_2D_Gauss(), and SetupCFunctions().
◆ m_dGL_xi3
Definition at line 95 of file DiffusionLFRNS.h.
Referenced by SetupCFunctions().
◆ m_dGR_xi1
Definition at line 92 of file DiffusionLFRNS.h.
Referenced by DerCFlux_1D(), DerCFlux_2D(), DivCFlux_2D(), DivCFlux_2D_Gauss(), and SetupCFunctions().
◆ m_dGR_xi2
Definition at line 94 of file DiffusionLFRNS.h.
Referenced by DerCFlux_2D(), DivCFlux_2D(), DivCFlux_2D_Gauss(), and SetupCFunctions().
◆ m_dGR_xi3
Definition at line 96 of file DiffusionLFRNS.h.
Referenced by SetupCFunctions().
◆ m_diffDim
|
private |
Definition at line 127 of file DiffusionLFRNS.h.
Referenced by v_Diffuse(), and v_InitObject().
◆ m_diffType
|
protected |
Definition at line 55 of file DiffusionLFRNS.h.
Referenced by SetupCFunctions().
◆ m_divFC
Definition at line 119 of file DiffusionLFRNS.h.
Referenced by v_Diffuse(), and v_InitObject().
◆ m_divFD
Definition at line 118 of file DiffusionLFRNS.h.
Referenced by v_Diffuse(), and v_InitObject().
◆ m_DU1
|
private |
Definition at line 110 of file DiffusionLFRNS.h.
Referenced by v_Diffuse(), and v_InitObject().
◆ m_gamma
|
private |
Definition at line 100 of file DiffusionLFRNS.h.
Referenced by v_InitObject(), and WeakPenaltyO1().
◆ m_gasConstant
|
private |
Definition at line 101 of file DiffusionLFRNS.h.
Referenced by v_InitObject(), and WeakPenaltyO1().
◆ m_gmat
Definition at line 84 of file DiffusionLFRNS.h.
Referenced by SetupMetrics(), and v_Diffuse().
◆ m_homoDerivs
Definition at line 124 of file DiffusionLFRNS.h.
Referenced by v_Diffuse(), and v_SetHomoDerivs().
◆ m_IF1
|
private |
Definition at line 109 of file DiffusionLFRNS.h.
Referenced by v_Diffuse(), and v_InitObject().
◆ m_Ixm
|
private |
Definition at line 97 of file DiffusionLFRNS.h.
◆ m_Ixp
|
private |
Definition at line 98 of file DiffusionLFRNS.h.
◆ m_jac
Definition at line 83 of file DiffusionLFRNS.h.
Referenced by SetupMetrics(), and v_Diffuse().
◆ m_mu
|
private |
Definition at line 104 of file DiffusionLFRNS.h.
Referenced by v_InitObject().
◆ m_pInf
|
private |
Definition at line 107 of file DiffusionLFRNS.h.
Referenced by v_InitObject().
◆ m_Q2D_e0
Definition at line 86 of file DiffusionLFRNS.h.
Referenced by DivCFlux_2D(), DivCFlux_2D_Gauss(), and SetupMetrics().
◆ m_Q2D_e1
Definition at line 87 of file DiffusionLFRNS.h.
Referenced by DivCFlux_2D(), DivCFlux_2D_Gauss(), and SetupMetrics().
◆ m_Q2D_e2
Definition at line 88 of file DiffusionLFRNS.h.
Referenced by DivCFlux_2D(), DivCFlux_2D_Gauss(), and SetupMetrics().
◆ m_Q2D_e3
Definition at line 89 of file DiffusionLFRNS.h.
Referenced by DivCFlux_2D(), DivCFlux_2D_Gauss(), and SetupMetrics().
◆ m_rhoInf
|
private |
Definition at line 106 of file DiffusionLFRNS.h.
Referenced by v_InitObject().
◆ m_session
|
private |
Definition at line 81 of file DiffusionLFRNS.h.
Referenced by v_InitObject().
◆ m_spaceDim
|
private |
Definition at line 126 of file DiffusionLFRNS.h.
Referenced by v_InitObject(), and WeakPenaltyO2().
◆ m_thermalConductivity
|
private |
Definition at line 105 of file DiffusionLFRNS.h.
Referenced by v_InitObject().
◆ m_tmp1
|
private |
Definition at line 121 of file DiffusionLFRNS.h.
Referenced by v_Diffuse(), and v_InitObject().
◆ m_tmp2
|
private |
Definition at line 122 of file DiffusionLFRNS.h.
Referenced by v_Diffuse(), and v_InitObject().
◆ m_traceNormals
Definition at line 80 of file DiffusionLFRNS.h.
Referenced by DivCFlux_2D(), DivCFlux_2D_Gauss(), NumericalFluxO1(), NumericalFluxO2(), SetupMetrics(), and WeakPenaltyO2().
◆ m_traceVel
Definition at line 79 of file DiffusionLFRNS.h.
Referenced by NumericalFluxO1(), NumericalFluxO2(), and v_InitObject().
◆ m_Twall
|
private |
Definition at line 102 of file DiffusionLFRNS.h.
Referenced by v_InitObject(), and WeakPenaltyO1().
◆ m_viscFlux
Definition at line 116 of file DiffusionLFRNS.h.
Referenced by v_Diffuse(), and v_InitObject().
◆ m_ViscosityType
|
private |
Definition at line 103 of file DiffusionLFRNS.h.
Referenced by v_InitObject().
◆ m_viscTensor
|
private |
Definition at line 115 of file DiffusionLFRNS.h.
Referenced by v_Diffuse(), v_GetFluxTensor(), and v_InitObject().
◆ type
|
static |
Initial value:
= {
tKey RegisterCreatorFunction(tKey idKey, CreatorFunction classCreator, std::string pDesc="")
Register a class with the factory.
Definition: NekFactory.hpp:197
static DiffusionSharedPtr create(std::string diffType)
Definition: DiffusionLFRNS.h:45
DiffusionFactory & GetDiffusionFactory()
Definition: Diffusion.cpp:39
Definition at line 50 of file DiffusionLFRNS.h.
