#include <DiffusionLFR.h>
Inheritance diagram for Nektar::SolverUtils::DiffusionLFR:
Static Public Member Functions | |
| static DiffusionSharedPtr | create (std::string diffType) |
Static Public Attributes | |
| static std::string | type [] |
Protected Member Functions | |
| DiffusionLFR (std::string diffType) | |
| DiffusionLFR uses the Flux Reconstruction (FR) approach to compute the diffusion term. The implementation is only for segments, quadrilaterals and hexahedra at the moment. More... | |
| virtual void | v_InitObject (LibUtilities::SessionReaderSharedPtr pSession, Array< OneD, MultiRegions::ExpListSharedPtr > pFields) override |
| Initiliase DiffusionLFR objects and store them before starting the time-stepping. More... | |
| virtual void | v_Diffuse (const size_t 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, const Array< OneD, Array< OneD, NekDouble > > &pBwd) override |
| Calculate FR Diffusion for the linear problems using an LDG interface flux. More... | |
 Protected Member Functions inherited from Nektar::SolverUtils::Diffusion | |
| virtual SOLVER_UTILS_EXPORT void | v_InitObject (LibUtilities::SessionReaderSharedPtr pSession, Array< OneD, MultiRegions::ExpListSharedPtr > pFields)=0 |
| virtual SOLVER_UTILS_EXPORT void | v_Diffuse (const std::size_t 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, const Array< OneD, Array< OneD, NekDouble > > &pBwd)=0 |
| virtual SOLVER_UTILS_EXPORT void | v_DiffuseCoeffs (const std::size_t 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, const Array< OneD, Array< OneD, NekDouble > > &pBwd) |
| virtual SOLVER_UTILS_EXPORT void | v_DiffuseCalcDerivative (const Array< OneD, MultiRegions::ExpListSharedPtr > &fields, const Array< OneD, Array< OneD, NekDouble > > &inarray, TensorOfArray3D< NekDouble > &qfields, const Array< OneD, Array< OneD, NekDouble > > &pFwd, const Array< OneD, Array< OneD, NekDouble > > &pBwd) |
| Diffusion Flux, calculate the physical derivatives. More... | |
| virtual SOLVER_UTILS_EXPORT void | v_DiffuseVolumeFlux (const Array< OneD, MultiRegions::ExpListSharedPtr > &fields, const Array< OneD, Array< OneD, NekDouble > > &inarray, TensorOfArray3D< NekDouble > &qfields, TensorOfArray3D< NekDouble > &VolumeFlux, Array< OneD, int > &nonZeroIndex) |
| Diffusion Volume Flux. More... | |
| virtual SOLVER_UTILS_EXPORT void | v_DiffuseTraceFlux (const Array< OneD, MultiRegions::ExpListSharedPtr > &fields, const Array< OneD, Array< OneD, NekDouble > > &inarray, TensorOfArray3D< NekDouble > &Qfields, TensorOfArray3D< NekDouble > &VolumeFlux, Array< OneD, Array< OneD, NekDouble > > &TraceFlux, const Array< OneD, Array< OneD, NekDouble > > &pFwd, const Array< OneD, Array< OneD, NekDouble > > &pBwd, Array< OneD, int > &nonZeroIndex) |
| Diffusion term Trace Flux. More... | |
| virtual void | v_SetHomoDerivs (Array< OneD, Array< OneD, NekDouble > > &deriv) |
| virtual TensorOfArray3D< NekDouble > & | v_GetFluxTensor () |
| virtual SOLVER_UTILS_EXPORT const Array< OneD, const Array< OneD, NekDouble > > & | v_GetTraceNormal () |
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 | NumFluxforScalar (const Array< OneD, MultiRegions::ExpListSharedPtr > &fields, const Array< OneD, Array< OneD, NekDouble > > &ufield, Array< OneD, Array< OneD, Array< OneD, NekDouble > > > &uflux) |
| Builds the numerical flux for the 1st order derivatives. More... | |
| void | WeakPenaltyforScalar (const Array< OneD, MultiRegions::ExpListSharedPtr > &fields, const int var, const Array< OneD, const NekDouble > &ufield, Array< OneD, NekDouble > &penaltyflux) |
| Imposes appropriate bcs for the 1st order derivatives. More... | |
| void | NumFluxforVector (const Array< OneD, MultiRegions::ExpListSharedPtr > &fields, const Array< OneD, Array< OneD, NekDouble > > &ufield, Array< OneD, Array< OneD, Array< OneD, NekDouble > > > &qfield, Array< OneD, Array< OneD, NekDouble > > &qflux) |
| Build the numerical flux for the 2nd order derivatives. More... | |
| void | WeakPenaltyforVector (const Array< OneD, MultiRegions::ExpListSharedPtr > &fields, const int var, const int dir, const Array< OneD, const NekDouble > &qfield, Array< OneD, NekDouble > &penaltyflux, NekDouble C11) |
| Imposes appropriate bcs for the 2nd order derivatives. More... | |
| 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 44 of file DiffusionLFR.h.
Constructor & Destructor Documentation
◆ DiffusionLFR()
|
protected |
DiffusionLFR uses the Flux Reconstruction (FR) approach to compute the diffusion term. The implementation is only for segments, quadrilaterals and hexahedra at the moment.
- Todo:
- Extension to triangles, tetrahedra and other shapes. (Long term objective)
Definition at line 71 of file DiffusionLFR.cpp.
Referenced by create().
Member Function Documentation
◆ create()
|
inlinestatic |
Definition at line 47 of file DiffusionLFR.h.
48 {
50 }
DiffusionLFR(std::string diffType)
DiffusionLFR uses the Flux Reconstruction (FR) approach to compute the diffusion term....
Definition: DiffusionLFR.cpp:71
std::shared_ptr< Diffusion > DiffusionSharedPtr
A shared pointer to an EquationSystem object.
Definition: Diffusion.h:58
References DiffusionLFR().
◆ 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 1457 of file DiffusionLFR.cpp.
1462{
1463 boost::ignore_unused(nConvectiveFields);
1464
1465 int n;
1466 int nLocalSolutionPts, phys_offset, t_offset;
1467
1468 Array<OneD, NekDouble> auxArray1, auxArray2;
1469 Array<TwoD, const NekDouble> gmat;
1470 Array<OneD, const NekDouble> jac;
1471
1472 LibUtilities::PointsKeyVector ptsKeys;
1473 LibUtilities::BasisSharedPtr Basis;
1474 Basis = fields[0]->GetExp(0)->GetBasis(0);
1475
1476 int nElements = fields[0]->GetExpSize();
1477 int nSolutionPts = fields[0]->GetTotPoints();
1478
1479 vector<bool> leftAdjacentTraces =
1480 std::static_pointer_cast<MultiRegions::DisContField>(fields[0])
1481 ->GetLeftAdjacentTraces();
1482
1483 // Arrays to store the derivatives of the correction flux
1484 Array<OneD, NekDouble> DCL(nSolutionPts / nElements, 0.0);
1485 Array<OneD, NekDouble> DCR(nSolutionPts / nElements, 0.0);
1486
1487 // Arrays to store the intercell numerical flux jumps
1488 Array<OneD, NekDouble> JumpL(nElements);
1489 Array<OneD, NekDouble> JumpR(nElements);
1490
1491 Array<OneD, Array<OneD, LocalRegions::ExpansionSharedPtr>> &elmtToTrace =
1492 fields[0]->GetTraceMap()->GetElmtToTrace();
1493
1494 for (n = 0; n < nElements; ++n)
1495 {
1496 nLocalSolutionPts = fields[0]->GetExp(n)->GetTotPoints();
1497 phys_offset = fields[0]->GetPhys_Offset(n);
1498
1499 Array<OneD, NekDouble> tmparrayX1(nLocalSolutionPts, 0.0);
1500 Array<OneD, NekDouble> tmpFluxVertex(1, 0.0);
1501 Vmath::Vcopy(nLocalSolutionPts, &flux[phys_offset], 1, &tmparrayX1[0],
1502 1);
1503
1504 fields[0]->GetExp(n)->GetTracePhysVals(0, elmtToTrace[n][0], tmparrayX1,
1505 tmpFluxVertex);
1506
1507 t_offset = fields[0]->GetTrace()->GetPhys_Offset(
1508 elmtToTrace[n][0]->GetElmtId());
1509
1510 if (leftAdjacentTraces[2 * n])
1511 {
1512 JumpL[n] = -iFlux[t_offset] - tmpFluxVertex[0];
1513 }
1514 else
1515 {
1516 JumpL[n] = iFlux[t_offset] - tmpFluxVertex[0];
1517 }
1518
1519 t_offset = fields[0]->GetTrace()->GetPhys_Offset(
1520 elmtToTrace[n][1]->GetElmtId());
1521
1522 fields[0]->GetExp(n)->GetTracePhysVals(1, elmtToTrace[n][1], tmparrayX1,
1523 tmpFluxVertex);
1524 if (leftAdjacentTraces[2 * n + 1])
1525 {
1526 JumpR[n] = iFlux[t_offset] - tmpFluxVertex[0];
1527 }
1528 else
1529 {
1530 JumpR[n] = -iFlux[t_offset] - tmpFluxVertex[0];
1531 }
1532 }
1533
1534 for (n = 0; n < nElements; ++n)
1535 {
1536 ptsKeys = fields[0]->GetExp(n)->GetPointsKeys();
1537 nLocalSolutionPts = fields[0]->GetExp(n)->GetTotPoints();
1538 phys_offset = fields[0]->GetPhys_Offset(n);
1539 jac = fields[0]
1540 ->GetExp(n)
1541 ->as<LocalRegions::Expansion1D>()
1542 ->GetGeom1D()
1543 ->GetMetricInfo()
1544 ->GetJac(ptsKeys);
1545
1546 JumpL[n] = JumpL[n] * jac[0];
1547 JumpR[n] = JumpR[n] * jac[0];
1548
1549 // Left jump multiplied by left derivative of C function
1551 1);
1552
1553 // Right jump multiplied by right derivative of C function
1555 1);
1556
1557 // Assembling derivative of the correction flux
1558 Vmath::Vadd(nLocalSolutionPts, &DCL[0], 1, &DCR[0], 1,
1559 &derCFlux[phys_offset], 1);
1560 }
1561}
Array< OneD, Array< OneD, NekDouble > > m_dGR_xi1
Definition: DiffusionLFR.h:84
Array< OneD, Array< OneD, NekDouble > > m_dGL_xi1
Definition: DiffusionLFR.h:83
std::shared_ptr< Basis > BasisSharedPtr
Definition: FoundationsFwd.hpp:71
std::vector< PointsKey > PointsKeyVector
Definition: Points.h:236
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:354
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.cpp:245
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.cpp:1191
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 1578 of file DiffusionLFR.cpp.
1583{
1584 boost::ignore_unused(nConvectiveFields);
1585
1586 int n, e, i, j, cnt;
1587
1588 Array<OneD, const NekDouble> jac;
1589
1590 int nElements = fields[0]->GetExpSize();
1591 int trace_offset, phys_offset;
1592 int nLocalSolutionPts;
1593 int nquad0, nquad1;
1594 int nEdgePts;
1595
1596 LibUtilities::PointsKeyVector ptsKeys;
1597 Array<OneD, NekDouble> auxArray1, auxArray2;
1598 Array<OneD, LibUtilities::BasisSharedPtr> base;
1599 Array<OneD, Array<OneD, LocalRegions::ExpansionSharedPtr>> &elmtToTrace =
1600 fields[0]->GetTraceMap()->GetElmtToTrace();
1601
1602 // Loop on the elements
1603 for (n = 0; n < nElements; ++n)
1604 {
1605 // Offset of the element on the global vector
1606 phys_offset = fields[0]->GetPhys_Offset(n);
1607 nLocalSolutionPts = fields[0]->GetExp(n)->GetTotPoints();
1608 ptsKeys = fields[0]->GetExp(n)->GetPointsKeys();
1609
1610 jac = fields[0]
1611 ->GetExp(n)
1612 ->as<LocalRegions::Expansion2D>()
1613 ->GetGeom2D()
1614 ->GetMetricInfo()
1615 ->GetJac(ptsKeys);
1616
1617 base = fields[0]->GetExp(n)->GetBase();
1618 nquad0 = base[0]->GetNumPoints();
1619 nquad1 = base[1]->GetNumPoints();
1620
1621 Array<OneD, NekDouble> divCFluxE0(nLocalSolutionPts, 0.0);
1622 Array<OneD, NekDouble> divCFluxE1(nLocalSolutionPts, 0.0);
1623 Array<OneD, NekDouble> divCFluxE2(nLocalSolutionPts, 0.0);
1624 Array<OneD, NekDouble> divCFluxE3(nLocalSolutionPts, 0.0);
1625
1626 // Loop on the edges
1627 for (e = 0; e < fields[0]->GetExp(n)->GetNtraces(); ++e)
1628 {
1629 // Number of edge points of edge 'e'
1630 nEdgePts = fields[0]->GetExp(n)->GetTraceNumPoints(e);
1631
1632 // Array for storing volumetric fluxes on each edge
1633 Array<OneD, NekDouble> tmparray(nEdgePts, 0.0);
1634
1635 // Offset of the trace space correspondent to edge 'e'
1636 trace_offset = fields[0]->GetTrace()->GetPhys_Offset(
1637 elmtToTrace[n][e]->GetElmtId());
1638
1639 // Extract the edge values of the volumetric fluxes
1640 // on edge 'e' and order them accordingly to the order
1641 // of the trace space
1642 fields[0]->GetExp(n)->GetTracePhysVals(
1643 e, elmtToTrace[n][e], flux + phys_offset, auxArray1 = tmparray);
1644
1645 // Compute the fluxJumps per each edge 'e' and each
1646 // flux point
1647 Array<OneD, NekDouble> fluxJumps(nEdgePts, 0.0);
1648 Vmath::Vsub(nEdgePts, &iFlux[trace_offset], 1, &tmparray[0], 1,
1649 &fluxJumps[0], 1);
1650
1651 // Check the ordering of the fluxJumps and reverse
1652 // it in case of backward definition of edge 'e'
1653 if (fields[0]->GetExp(n)->GetTraceOrient(e) ==
1654 StdRegions::eBackwards)
1655 {
1656 Vmath::Reverse(nEdgePts, &fluxJumps[0], 1, &fluxJumps[0], 1);
1657 }
1658
1659 // Deformed elements
1660 if (fields[0]
1661 ->GetExp(n)
1662 ->as<LocalRegions::Expansion2D>()
1663 ->GetGeom2D()
1664 ->GetMetricInfo()
1665 ->GetGtype() == SpatialDomains::eDeformed)
1666 {
1667 // Extract the Jacobians along edge 'e'
1668 Array<OneD, NekDouble> jacEdge(nEdgePts, 0.0);
1669
1670 fields[0]->GetExp(n)->GetTracePhysVals(
1671 e, elmtToTrace[n][e], jac, auxArray1 = jacEdge);
1672
1673 // Check the ordering of the fluxJumps and reverse
1674 // it in case of backward definition of edge 'e'
1675 if (fields[0]->GetExp(n)->GetTraceOrient(e) ==
1676 StdRegions::eBackwards)
1677 {
1678 Vmath::Reverse(nEdgePts, &jacEdge[0], 1, &jacEdge[0], 1);
1679 }
1680
1681 // Multiply the fluxJumps by the edge 'e' Jacobians
1682 // to bring the fluxJumps into the standard space
1683 for (j = 0; j < nEdgePts; j++)
1684 {
1685 fluxJumps[j] = fluxJumps[j] * jacEdge[j];
1686 }
1687 }
1688 // Non-deformed elements
1689 else
1690 {
1691 // Multiply the fluxJumps by the edge 'e' Jacobians
1692 // to bring the fluxJumps into the standard space
1693 Vmath::Smul(nEdgePts, jac[0], fluxJumps, 1, fluxJumps, 1);
1694 }
1695
1696 // Multiply jumps by derivatives of the correction functions
1697 // All the quntities at this level should be defined into
1698 // the standard space
1699 switch (e)
1700 {
1701 case 0:
1702 for (i = 0; i < nquad0; ++i)
1703 {
1704 for (j = 0; j < nquad1; ++j)
1705 {
1706 cnt = i + j * nquad0;
1707 divCFluxE0[cnt] = fluxJumps[i] * m_dGL_xi2[n][j];
1708 }
1709 }
1710 break;
1711 case 1:
1712 for (i = 0; i < nquad1; ++i)
1713 {
1714 for (j = 0; j < nquad0; ++j)
1715 {
1716 cnt = (nquad0)*i + j;
1717 divCFluxE1[cnt] = fluxJumps[i] * m_dGR_xi1[n][j];
1718 }
1719 }
1720 break;
1721 case 2:
1722 for (i = 0; i < nquad0; ++i)
1723 {
1724 for (j = 0; j < nquad1; ++j)
1725 {
1726 cnt = j * nquad0 + i;
1727 divCFluxE2[cnt] = fluxJumps[i] * m_dGR_xi2[n][j];
1728 }
1729 }
1730 break;
1731 case 3:
1732 for (i = 0; i < nquad1; ++i)
1733 {
1734 for (j = 0; j < nquad0; ++j)
1735 {
1736 cnt = j + i * nquad0;
1737 divCFluxE3[cnt] = fluxJumps[i] * m_dGL_xi1[n][j];
1738 }
1739 }
1740 break;
1741
1742 default:
1744 break;
1745 }
1746 }
1747
1748 // Summing the component of the flux for calculating the
1749 // derivatives per each direction
1750 if (direction == 0)
1751 {
1752 Vmath::Vadd(nLocalSolutionPts, &divCFluxE1[0], 1, &divCFluxE3[0], 1,
1753 &derCFlux[phys_offset], 1);
1754 }
1755 else if (direction == 1)
1756 {
1757 Vmath::Vadd(nLocalSolutionPts, &divCFluxE0[0], 1, &divCFluxE2[0], 1,
1758 &derCFlux[phys_offset], 1);
1759 }
1760 }
1761}
Array< OneD, Array< OneD, NekDouble > > m_dGL_xi2
Definition: DiffusionLFR.h:85
Array< OneD, Array< OneD, NekDouble > > m_dGR_xi2
Definition: DiffusionLFR.h:86
void Reverse(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.cpp:1222
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:414
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 1776 of file DiffusionLFR.cpp.
1783{
1784 boost::ignore_unused(nConvectiveFields);
1785
1786 int n, e, i, j, cnt;
1787
1788 int nElements = fields[0]->GetExpSize();
1789
1790 int nLocalSolutionPts;
1791 int nEdgePts;
1792 int trace_offset;
1793 int phys_offset;
1794 int nquad0;
1795 int nquad1;
1796
1797 Array<OneD, NekDouble> auxArray1, auxArray2;
1798 Array<OneD, LibUtilities::BasisSharedPtr> base;
1799
1800 Array<OneD, Array<OneD, LocalRegions::ExpansionSharedPtr>> &elmtToTrace =
1801 fields[0]->GetTraceMap()->GetElmtToTrace();
1802
1803 // Loop on the elements
1804 for (n = 0; n < nElements; ++n)
1805 {
1806 // Offset of the element on the global vector
1807 phys_offset = fields[0]->GetPhys_Offset(n);
1808 nLocalSolutionPts = fields[0]->GetExp(n)->GetTotPoints();
1809
1810 base = fields[0]->GetExp(n)->GetBase();
1811 nquad0 = base[0]->GetNumPoints();
1812 nquad1 = base[1]->GetNumPoints();
1813
1814 Array<OneD, NekDouble> divCFluxE0(nLocalSolutionPts, 0.0);
1815 Array<OneD, NekDouble> divCFluxE1(nLocalSolutionPts, 0.0);
1816 Array<OneD, NekDouble> divCFluxE2(nLocalSolutionPts, 0.0);
1817 Array<OneD, NekDouble> divCFluxE3(nLocalSolutionPts, 0.0);
1818
1819 // Loop on the edges
1820 for (e = 0; e < fields[0]->GetExp(n)->GetNtraces(); ++e)
1821 {
1822 // Number of edge points of edge e
1823 nEdgePts = fields[0]->GetExp(n)->GetTraceNumPoints(e);
1824
1825 Array<OneD, NekDouble> tmparrayX1(nEdgePts, 0.0);
1826 Array<OneD, NekDouble> tmparrayX2(nEdgePts, 0.0);
1827 Array<OneD, NekDouble> fluxN(nEdgePts, 0.0);
1828 Array<OneD, NekDouble> fluxT(nEdgePts, 0.0);
1829 Array<OneD, NekDouble> fluxJumps(nEdgePts, 0.0);
1830
1831 // Offset of the trace space correspondent to edge e
1832 trace_offset = fields[0]->GetTrace()->GetPhys_Offset(
1833 elmtToTrace[n][e]->GetElmtId());
1834
1835 // Get the normals of edge e
1836 const Array<OneD, const Array<OneD, NekDouble>> &normals =
1837 fields[0]->GetExp(n)->GetTraceNormal(e);
1838
1839 // Extract the edge values of flux-x on edge e and order
1840 // them accordingly to the order of the trace space
1841 fields[0]->GetExp(n)->GetTracePhysVals(e, elmtToTrace[n][e],
1842 fluxX1 + phys_offset,
1843 auxArray1 = tmparrayX1);
1844
1845 // Extract the edge values of flux-y on edge e and order
1846 // them accordingly to the order of the trace space
1847 fields[0]->GetExp(n)->GetTracePhysVals(e, elmtToTrace[n][e],
1848 fluxX2 + phys_offset,
1849 auxArray1 = tmparrayX2);
1850
1851 // Multiply the edge components of the flux by the normal
1852 for (i = 0; i < nEdgePts; ++i)
1853 {
1854 fluxN[i] = tmparrayX1[i] * m_traceNormals[0][trace_offset + i] +
1855 tmparrayX2[i] * m_traceNormals[1][trace_offset + i];
1856 }
1857
1858 // Subtract to the Riemann flux the discontinuous flux
1859 Vmath::Vsub(nEdgePts, &numericalFlux[trace_offset], 1, &fluxN[0], 1,
1860 &fluxJumps[0], 1);
1861
1862 // Check the ordering of the jump vectors
1863 if (fields[0]->GetExp(n)->GetTraceOrient(e) ==
1864 StdRegions::eBackwards)
1865 {
1866 Vmath::Reverse(nEdgePts, auxArray1 = fluxJumps, 1,
1867 auxArray2 = fluxJumps, 1);
1868 }
1869
1870 for (i = 0; i < nEdgePts; ++i)
1871 {
1873 m_traceNormals[1][trace_offset + i] != normals[1][i])
1874 {
1875 fluxJumps[i] = -fluxJumps[i];
1876 }
1877 }
1878
1879 // Multiply jumps by derivatives of the correction functions
1880 switch (e)
1881 {
1882 case 0:
1883 for (i = 0; i < nquad0; ++i)
1884 {
1885 // Multiply fluxJumps by Q factors
1886 fluxJumps[i] = -(m_Q2D_e0[n][i]) * fluxJumps[i];
1887
1888 for (j = 0; j < nquad1; ++j)
1889 {
1890 cnt = i + j * nquad0;
1891 divCFluxE0[cnt] = fluxJumps[i] * m_dGL_xi2[n][j];
1892 }
1893 }
1894 break;
1895 case 1:
1896 for (i = 0; i < nquad1; ++i)
1897 {
1898 // Multiply fluxJumps by Q factors
1899 fluxJumps[i] = (m_Q2D_e1[n][i]) * fluxJumps[i];
1900
1901 for (j = 0; j < nquad0; ++j)
1902 {
1903 cnt = (nquad0)*i + j;
1904 divCFluxE1[cnt] = fluxJumps[i] * m_dGR_xi1[n][j];
1905 }
1906 }
1907 break;
1908 case 2:
1909 for (i = 0; i < nquad0; ++i)
1910 {
1911 // Multiply fluxJumps by Q factors
1912 fluxJumps[i] = (m_Q2D_e2[n][i]) * fluxJumps[i];
1913
1914 for (j = 0; j < nquad1; ++j)
1915 {
1916 cnt = j * nquad0 + i;
1917 divCFluxE2[cnt] = fluxJumps[i] * m_dGR_xi2[n][j];
1918 }
1919 }
1920 break;
1921 case 3:
1922 for (i = 0; i < nquad1; ++i)
1923 {
1924 // Multiply fluxJumps by Q factors
1925 fluxJumps[i] = -(m_Q2D_e3[n][i]) * fluxJumps[i];
1926 for (j = 0; j < nquad0; ++j)
1927 {
1928 cnt = j + i * nquad0;
1929 divCFluxE3[cnt] = fluxJumps[i] * m_dGL_xi1[n][j];
1930 }
1931 }
1932 break;
1933
1934 default:
1936 break;
1937 }
1938 }
1939
1940 // Sum each edge contribution
1941 Vmath::Vadd(nLocalSolutionPts, &divCFluxE0[0], 1, &divCFluxE1[0], 1,
1942 &divCFlux[phys_offset], 1);
1943
1944 Vmath::Vadd(nLocalSolutionPts, &divCFlux[phys_offset], 1,
1945 &divCFluxE2[0], 1, &divCFlux[phys_offset], 1);
1946
1947 Vmath::Vadd(nLocalSolutionPts, &divCFlux[phys_offset], 1,
1948 &divCFluxE3[0], 1, &divCFlux[phys_offset], 1);
1949 }
1950}
Array< OneD, Array< OneD, NekDouble > > m_traceNormals
Definition: DiffusionLFR.h:72
Array< OneD, Array< OneD, NekDouble > > m_Q2D_e0
Definition: DiffusionLFR.h:78
Array< OneD, Array< OneD, NekDouble > > m_Q2D_e3
Definition: DiffusionLFR.h:81
Array< OneD, Array< OneD, NekDouble > > m_Q2D_e2
Definition: DiffusionLFR.h:80
Array< OneD, Array< OneD, NekDouble > > m_Q2D_e1
Definition: DiffusionLFR.h:79
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 1966 of file DiffusionLFR.cpp.
1973{
1974 boost::ignore_unused(nConvectiveFields);
1975
1976 int n, e, i, j, cnt;
1977
1978 int nElements = fields[0]->GetExpSize();
1979 int nLocalSolutionPts;
1980 int nEdgePts;
1981 int trace_offset;
1982 int phys_offset;
1983 int nquad0;
1984 int nquad1;
1985
1986 Array<OneD, NekDouble> auxArray1, auxArray2;
1987 Array<OneD, LibUtilities::BasisSharedPtr> base;
1988
1989 Array<OneD, Array<OneD, LocalRegions::ExpansionSharedPtr>> &elmtToTrace =
1990 fields[0]->GetTraceMap()->GetElmtToTrace();
1991
1992 // Loop on the elements
1993 for (n = 0; n < nElements; ++n)
1994 {
1995 // Offset of the element on the global vector
1996 phys_offset = fields[0]->GetPhys_Offset(n);
1997 nLocalSolutionPts = fields[0]->GetExp(n)->GetTotPoints();
1998
1999 base = fields[0]->GetExp(n)->GetBase();
2000 nquad0 = base[0]->GetNumPoints();
2001 nquad1 = base[1]->GetNumPoints();
2002
2003 Array<OneD, NekDouble> divCFluxE0(nLocalSolutionPts, 0.0);
2004 Array<OneD, NekDouble> divCFluxE1(nLocalSolutionPts, 0.0);
2005 Array<OneD, NekDouble> divCFluxE2(nLocalSolutionPts, 0.0);
2006 Array<OneD, NekDouble> divCFluxE3(nLocalSolutionPts, 0.0);
2007
2008 // Loop on the edges
2009 for (e = 0; e < fields[0]->GetExp(n)->GetNtraces(); ++e)
2010 {
2011 // Number of edge points of edge e
2012 nEdgePts = fields[0]->GetExp(n)->GetTraceNumPoints(e);
2013
2014 Array<OneD, NekDouble> fluxN(nEdgePts, 0.0);
2015 Array<OneD, NekDouble> fluxT(nEdgePts, 0.0);
2016 Array<OneD, NekDouble> fluxN_R(nEdgePts, 0.0);
2017 Array<OneD, NekDouble> fluxN_D(nEdgePts, 0.0);
2018 Array<OneD, NekDouble> fluxJumps(nEdgePts, 0.0);
2019
2020 // Offset of the trace space correspondent to edge e
2021 trace_offset = fields[0]->GetTrace()->GetPhys_Offset(
2022 elmtToTrace[n][e]->GetElmtId());
2023
2024 // Get the normals of edge e
2025 const Array<OneD, const Array<OneD, NekDouble>> &normals =
2026 fields[0]->GetExp(n)->GetTraceNormal(e);
2027
2028 // Extract the trasformed normal flux at each edge
2029 switch (e)
2030 {
2031 case 0:
2032 // Extract the edge values of transformed flux-y on
2033 // edge e and order them accordingly to the order of
2034 // the trace space
2035 fields[0]->GetExp(n)->GetTracePhysVals(e, elmtToTrace[n][e],
2036 fluxX2 + phys_offset,
2037 auxArray1 = fluxN_D);
2038
2039 Vmath::Neg(nEdgePts, fluxN_D, 1);
2040
2041 // Extract the physical Rieamann flux at each edge
2042 Vmath::Vcopy(nEdgePts, &numericalFlux[trace_offset], 1,
2043 &fluxN[0], 1);
2044
2045 // Check the ordering of vectors
2046 if (fields[0]->GetExp(n)->GetTraceOrient(e) ==
2047 StdRegions::eBackwards)
2048 {
2049 Vmath::Reverse(nEdgePts, auxArray1 = fluxN, 1,
2050 auxArray2 = fluxN, 1);
2051
2052 Vmath::Reverse(nEdgePts, auxArray1 = fluxN_D, 1,
2053 auxArray2 = fluxN_D, 1);
2054 }
2055
2056 // Transform Riemann Fluxes in the standard element
2057 for (i = 0; i < nquad0; ++i)
2058 {
2059 // Multiply Riemann Flux by Q factors
2060 fluxN_R[i] = (m_Q2D_e0[n][i]) * fluxN[i];
2061 }
2062
2063 for (i = 0; i < nEdgePts; ++i)
2064 {
2066 normals[0][i] ||
2067 m_traceNormals[1][trace_offset + i] !=
2068 normals[1][i])
2069 {
2070 fluxN_R[i] = -fluxN_R[i];
2071 }
2072 }
2073
2074 // Subtract to the Riemann flux the discontinuous
2075 // flux
2076 Vmath::Vsub(nEdgePts, &fluxN_R[0], 1, &fluxN_D[0], 1,
2077 &fluxJumps[0], 1);
2078
2079 // Multiplicate the flux jumps for the correction
2080 // function
2081 for (i = 0; i < nquad0; ++i)
2082 {
2083 for (j = 0; j < nquad1; ++j)
2084 {
2085 cnt = i + j * nquad0;
2086 divCFluxE0[cnt] = -fluxJumps[i] * m_dGL_xi2[n][j];
2087 }
2088 }
2089 break;
2090 case 1:
2091 // Extract the edge values of transformed flux-x on
2092 // edge e and order them accordingly to the order of
2093 // the trace space
2094 fields[0]->GetExp(n)->GetTracePhysVals(e, elmtToTrace[n][e],
2095 fluxX1 + phys_offset,
2096 auxArray1 = fluxN_D);
2097
2098 // Extract the physical Rieamann flux at each edge
2099 Vmath::Vcopy(nEdgePts, &numericalFlux[trace_offset], 1,
2100 &fluxN[0], 1);
2101
2102 // Check the ordering of vectors
2103 if (fields[0]->GetExp(n)->GetTraceOrient(e) ==
2104 StdRegions::eBackwards)
2105 {
2106 Vmath::Reverse(nEdgePts, auxArray1 = fluxN, 1,
2107 auxArray2 = fluxN, 1);
2108
2109 Vmath::Reverse(nEdgePts, auxArray1 = fluxN_D, 1,
2110 auxArray2 = fluxN_D, 1);
2111 }
2112
2113 // Transform Riemann Fluxes in the standard element
2114 for (i = 0; i < nquad1; ++i)
2115 {
2116 // Multiply Riemann Flux by Q factors
2117 fluxN_R[i] = (m_Q2D_e1[n][i]) * fluxN[i];
2118 }
2119
2120 for (i = 0; i < nEdgePts; ++i)
2121 {
2123 normals[0][i] ||
2124 m_traceNormals[1][trace_offset + i] !=
2125 normals[1][i])
2126 {
2127 fluxN_R[i] = -fluxN_R[i];
2128 }
2129 }
2130
2131 // Subtract to the Riemann flux the discontinuous
2132 // flux
2133 Vmath::Vsub(nEdgePts, &fluxN_R[0], 1, &fluxN_D[0], 1,
2134 &fluxJumps[0], 1);
2135
2136 // Multiplicate the flux jumps for the correction
2137 // function
2138 for (i = 0; i < nquad1; ++i)
2139 {
2140 for (j = 0; j < nquad0; ++j)
2141 {
2142 cnt = (nquad0)*i + j;
2143 divCFluxE1[cnt] = fluxJumps[i] * m_dGR_xi1[n][j];
2144 }
2145 }
2146 break;
2147 case 2:
2148
2149 // Extract the edge values of transformed flux-y on
2150 // edge e and order them accordingly to the order of
2151 // the trace space
2152
2153 fields[0]->GetExp(n)->GetTracePhysVals(e, elmtToTrace[n][e],
2154 fluxX2 + phys_offset,
2155 auxArray1 = fluxN_D);
2156
2157 // Extract the physical Rieamann flux at each edge
2158 Vmath::Vcopy(nEdgePts, &numericalFlux[trace_offset], 1,
2159 &fluxN[0], 1);
2160
2161 // Check the ordering of vectors
2162 if (fields[0]->GetExp(n)->GetTraceOrient(e) ==
2163 StdRegions::eBackwards)
2164 {
2165 Vmath::Reverse(nEdgePts, auxArray1 = fluxN, 1,
2166 auxArray2 = fluxN, 1);
2167
2168 Vmath::Reverse(nEdgePts, auxArray1 = fluxN_D, 1,
2169 auxArray2 = fluxN_D, 1);
2170 }
2171
2172 // Transform Riemann Fluxes in the standard element
2173 for (i = 0; i < nquad0; ++i)
2174 {
2175 // Multiply Riemann Flux by Q factors
2176 fluxN_R[i] = (m_Q2D_e2[n][i]) * fluxN[i];
2177 }
2178
2179 for (i = 0; i < nEdgePts; ++i)
2180 {
2182 normals[0][i] ||
2183 m_traceNormals[1][trace_offset + i] !=
2184 normals[1][i])
2185 {
2186 fluxN_R[i] = -fluxN_R[i];
2187 }
2188 }
2189
2190 // Subtract to the Riemann flux the discontinuous
2191 // flux
2192
2193 Vmath::Vsub(nEdgePts, &fluxN_R[0], 1, &fluxN_D[0], 1,
2194 &fluxJumps[0], 1);
2195
2196 // Multiplicate the flux jumps for the correction
2197 // function
2198 for (i = 0; i < nquad0; ++i)
2199 {
2200 for (j = 0; j < nquad1; ++j)
2201 {
2202 cnt = j * nquad0 + i;
2203 divCFluxE2[cnt] = fluxJumps[i] * m_dGR_xi2[n][j];
2204 }
2205 }
2206 break;
2207 case 3:
2208 // Extract the edge values of transformed flux-x on
2209 // edge e and order them accordingly to the order of
2210 // the trace space
2211
2212 fields[0]->GetExp(n)->GetTracePhysVals(e, elmtToTrace[n][e],
2213 fluxX1 + phys_offset,
2214 auxArray1 = fluxN_D);
2215 Vmath::Neg(nEdgePts, fluxN_D, 1);
2216
2217 // Extract the physical Rieamann flux at each edge
2218 Vmath::Vcopy(nEdgePts, &numericalFlux[trace_offset], 1,
2219 &fluxN[0], 1);
2220
2221 // Check the ordering of vectors
2222 if (fields[0]->GetExp(n)->GetTraceOrient(e) ==
2223 StdRegions::eBackwards)
2224 {
2225 Vmath::Reverse(nEdgePts, auxArray1 = fluxN, 1,
2226 auxArray2 = fluxN, 1);
2227
2228 Vmath::Reverse(nEdgePts, auxArray1 = fluxN_D, 1,
2229 auxArray2 = fluxN_D, 1);
2230 }
2231
2232 // Transform Riemann Fluxes in the standard element
2233 for (i = 0; i < nquad1; ++i)
2234 {
2235 // Multiply Riemann Flux by Q factors
2236 fluxN_R[i] = (m_Q2D_e3[n][i]) * fluxN[i];
2237 }
2238
2239 for (i = 0; i < nEdgePts; ++i)
2240 {
2242 normals[0][i] ||
2243 m_traceNormals[1][trace_offset + i] !=
2244 normals[1][i])
2245 {
2246 fluxN_R[i] = -fluxN_R[i];
2247 }
2248 }
2249
2250 // Subtract to the Riemann flux the discontinuous
2251 // flux
2252
2253 Vmath::Vsub(nEdgePts, &fluxN_R[0], 1, &fluxN_D[0], 1,
2254 &fluxJumps[0], 1);
2255
2256 // Multiplicate the flux jumps for the correction
2257 // function
2258 for (i = 0; i < nquad1; ++i)
2259 {
2260 for (j = 0; j < nquad0; ++j)
2261 {
2262 cnt = j + i * nquad0;
2263 divCFluxE3[cnt] = -fluxJumps[i] * m_dGL_xi1[n][j];
2264 }
2265 }
2266 break;
2267 default:
2269 break;
2270 }
2271 }
2272
2273 // Sum each edge contribution
2274 Vmath::Vadd(nLocalSolutionPts, &divCFluxE0[0], 1, &divCFluxE1[0], 1,
2275 &divCFlux[phys_offset], 1);
2276
2277 Vmath::Vadd(nLocalSolutionPts, &divCFlux[phys_offset], 1,
2278 &divCFluxE2[0], 1, &divCFlux[phys_offset], 1);
2279
2280 Vmath::Vadd(nLocalSolutionPts, &divCFlux[phys_offset], 1,
2281 &divCFluxE3[0], 1, &divCFlux[phys_offset], 1);
2282 }
2283}
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().
◆ NumFluxforScalar()
|
private |
Builds the numerical flux for the 1st order derivatives.
Definition at line 1177 of file DiffusionLFR.cpp.
1181{
1182 int i, j;
1183 int nTracePts = fields[0]->GetTrace()->GetTotPoints();
1184 int nvariables = fields.size();
1185 int nDim = fields[0]->GetCoordim(0);
1186
1187 Array<OneD, NekDouble> Fwd(nTracePts);
1188 Array<OneD, NekDouble> Bwd(nTracePts);
1189 Array<OneD, NekDouble> Vn(nTracePts, 0.0);
1190 Array<OneD, NekDouble> fluxtemp(nTracePts, 0.0);
1191
1192 // Get the normal velocity Vn
1193 for (i = 0; i < nDim; ++i)
1194 {
1196 }
1197
1198 // Get the sign of (v \cdot n), v = an arbitrary vector
1199 // Evaluate upwind flux:
1200 // uflux = \hat{u} \phi \cdot u = u^{(+,-)} n
1201 for (j = 0; j < nDim; ++j)
1202 {
1203 for (i = 0; i < nvariables; ++i)
1204 {
1205 // Compute Fwd and Bwd value of ufield of i direction
1206 fields[i]->GetFwdBwdTracePhys(ufield[i], Fwd, Bwd);
1207
1208 // if Vn >= 0, flux = uFwd, i.e.,
1209 // edge::eForward, if V*n>=0 <=> V*n_F>=0, pick uflux = uFwd
1210 // edge::eBackward, if V*n>=0 <=> V*n_B<0, pick uflux = uFwd
1211
1212 // else if Vn < 0, flux = uBwd, i.e.,
1213 // edge::eForward, if V*n<0 <=> V*n_F<0, pick uflux = uBwd
1214 // edge::eBackward, if V*n<0 <=> V*n_B>=0, pick uflux = uBwd
1215
1216 fields[i]->GetTrace()->Upwind(Vn, Fwd, Bwd, fluxtemp);
1217
1218 // Imposing weak boundary condition with flux
1219 // if Vn >= 0, uflux = uBwd at Neumann, i.e.,
1220 // edge::eForward, if V*n>=0 <=> V*n_F>=0, pick uflux = uBwd
1221 // edge::eBackward, if V*n>=0 <=> V*n_B<0, pick uflux = uBwd
1222
1223 // if Vn >= 0, uflux = uFwd at Neumann, i.e.,
1224 // edge::eForward, if V*n<0 <=> V*n_F<0, pick uflux = uFwd
1225 // edge::eBackward, if V*n<0 <=> V*n_B>=0, pick uflux = uFwd
1226
1227 if (fields[0]->GetBndCondExpansions().size())
1228 {
1229 WeakPenaltyforScalar(fields, i, ufield[i], fluxtemp);
1230 }
1231
1232 // if Vn >= 0, flux = uFwd*(tan_{\xi}^- \cdot \vec{n}),
1233 // i.e,
1234 // edge::eForward, uFwd \(\tan_{\xi}^Fwd \cdot \vec{n})
1235 // edge::eBackward, uFwd \(\tan_{\xi}^Bwd \cdot \vec{n})
1236
1237 // else if Vn < 0, flux = uBwd*(tan_{\xi}^- \cdot \vec{n}),
1238 // i.e,
1239 // edge::eForward, uBwd \(\tan_{\xi}^Fwd \cdot \vec{n})
1240 // edge::eBackward, uBwd \(\tan_{\xi}^Bwd \cdot \vec{n})
1241
1242 Vmath::Vcopy(nTracePts, fluxtemp, 1, uflux[i][j], 1);
1243 }
1244 }
1245}
void WeakPenaltyforScalar(const Array< OneD, MultiRegions::ExpListSharedPtr > &fields, const int var, const Array< OneD, const NekDouble > &ufield, Array< OneD, NekDouble > &penaltyflux)
Imposes appropriate bcs for the 1st order derivatives.
Definition: DiffusionLFR.cpp:1251
void Svtvp(int n, const T alpha, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
svtvp (scalar times vector plus vector): z = alpha*x + y
Definition: Vmath.cpp:617
References m_traceNormals, Vmath::Svtvp(), Vmath::Vcopy(), and WeakPenaltyforScalar().
Referenced by v_Diffuse().
◆ NumFluxforVector()
|
private |
Build the numerical flux for the 2nd order derivatives.
Definition at line 1314 of file DiffusionLFR.cpp.
1319{
1320 boost::ignore_unused(ufield);
1321
1322 int i, j;
1323 int nTracePts = fields[0]->GetTrace()->GetTotPoints();
1324 int nvariables = fields.size();
1325 int nDim = fields[0]->GetCoordim(0);
1326
1327 NekDouble C11 = 0.0;
1328 Array<OneD, NekDouble> Fwd(nTracePts);
1329 Array<OneD, NekDouble> Bwd(nTracePts);
1330 Array<OneD, NekDouble> Vn(nTracePts, 0.0);
1331
1332 Array<OneD, NekDouble> qFwd(nTracePts);
1333 Array<OneD, NekDouble> qBwd(nTracePts);
1334 Array<OneD, NekDouble> qfluxtemp(nTracePts, 0.0);
1335 Array<OneD, NekDouble> uterm(nTracePts);
1336
1337 // Get the normal velocity Vn
1338 for (i = 0; i < nDim; ++i)
1339 {
1341 }
1342
1343 // Evaulate upwind flux:
1344 // qflux = \hat{q} \cdot u = q \cdot n - C_(11)*(u^+ - u^-)
1345 for (i = 0; i < nvariables; ++i)
1346 {
1347 qflux[i] = Array<OneD, NekDouble>(nTracePts, 0.0);
1348 for (j = 0; j < nDim; ++j)
1349 {
1350 // Compute Fwd and Bwd value of ufield of jth direction
1351 fields[i]->GetFwdBwdTracePhys(qfield[i][j], qFwd, qBwd);
1352
1353 // if Vn >= 0, flux = uFwd, i.e.,
1354 // edge::eForward, if V*n>=0 <=> V*n_F>=0, pick
1355 // qflux = qBwd = q+
1356 // edge::eBackward, if V*n>=0 <=> V*n_B<0, pick
1357 // qflux = qBwd = q-
1358
1359 // else if Vn < 0, flux = uBwd, i.e.,
1360 // edge::eForward, if V*n<0 <=> V*n_F<0, pick
1361 // qflux = qFwd = q-
1362 // edge::eBackward, if V*n<0 <=> V*n_B>=0, pick
1363 // qflux = qFwd = q+
1364
1365 fields[i]->GetTrace()->Upwind(Vn, qBwd, qFwd, qfluxtemp);
1366
1368 qfluxtemp, 1);
1369
1370 // Imposing weak boundary condition with flux
1371 if (fields[0]->GetBndCondExpansions().size())
1372 {
1373 WeakPenaltyforVector(fields, i, j, qfield[i][j], qfluxtemp,
1374 C11);
1375 }
1376
1377 // q_hat \cdot n = (q_xi \cdot n_xi) or (q_eta \cdot n_eta)
1378 // n_xi = n_x * tan_xi_x + n_y * tan_xi_y + n_z * tan_xi_z
1379 // n_xi = n_x * tan_eta_x + n_y * tan_eta_y + n_z*tan_eta_z
1380 Vmath::Vadd(nTracePts, qfluxtemp, 1, qflux[i], 1, qflux[i], 1);
1381 }
1382 }
1383}
void WeakPenaltyforVector(const Array< OneD, MultiRegions::ExpListSharedPtr > &fields, const int var, const int dir, const Array< OneD, const NekDouble > &qfield, Array< OneD, NekDouble > &penaltyflux, NekDouble C11)
Imposes appropriate bcs for the 2nd order derivatives.
Definition: DiffusionLFR.cpp:1389
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:207
References m_traceNormals, Vmath::Svtvp(), Vmath::Vadd(), Vmath::Vmul(), and WeakPenaltyforVector().
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 into the Functions #DivCFlux_1D, DivCFlux_2D, #DivCFlux_3D to compute the the divergence of the correction flux for 1D, 2D or 3D problems.
- Parameters
-
pSession Pointer to session reader. pFields Pointer to fields.
Definition at line 330 of file DiffusionLFR.cpp.
333{
334 boost::ignore_unused(pSession);
335
336 int i, n;
337 NekDouble c0 = 0.0;
338 NekDouble c1 = 0.0;
339 NekDouble c2 = 0.0;
340 int nquad0, nquad1, nquad2;
341 int nmodes0, nmodes1, nmodes2;
342 Array<OneD, LibUtilities::BasisSharedPtr> base;
343
344 int nElements = pFields[0]->GetExpSize();
345 int nDim = pFields[0]->GetCoordim(0);
346
347 switch (nDim)
348 {
349 case 1:
350 {
351 m_dGL_xi1 = Array<OneD, Array<OneD, NekDouble>>(nElements);
352 m_dGR_xi1 = Array<OneD, Array<OneD, NekDouble>>(nElements);
353
354 for (n = 0; n < nElements; ++n)
355 {
356 base = pFields[0]->GetExp(n)->GetBase();
357 nquad0 = base[0]->GetNumPoints();
358 nmodes0 = base[0]->GetNumModes();
359 Array<OneD, const NekDouble> z0;
360 Array<OneD, const NekDouble> w0;
361
362 base[0]->GetZW(z0, w0);
363
364 m_dGL_xi1[n] = Array<OneD, NekDouble>(nquad0);
365 m_dGR_xi1[n] = Array<OneD, NekDouble>(nquad0);
366
367 // Auxiliary vectors to build up the auxiliary
368 // derivatives of the Legendre polynomials
369 Array<OneD, NekDouble> dLp0(nquad0, 0.0);
370 Array<OneD, NekDouble> dLpp0(nquad0, 0.0);
371 Array<OneD, NekDouble> dLpm0(nquad0, 0.0);
372
373 // Degree of the correction functions
374 int p0 = nmodes0 - 1;
375
376 // Function sign to compute left correction function
377 NekDouble sign0 = pow(-1.0, p0);
378
379 // Factorial factor to build the scheme
380 NekDouble ap0 = boost::math::tgamma(2 * p0 + 1) /
381 (pow(2.0, p0) * boost::math::tgamma(p0 + 1) *
382 boost::math::tgamma(p0 + 1));
383
384 // Scalar parameter which recovers the FR schemes
386 {
387 c0 = 0.0;
388 }
390 {
391 c0 = 2.0 * p0 /
392 ((2.0 * p0 + 1.0) * (p0 + 1.0) *
393 (ap0 * boost::math::tgamma(p0 + 1)) *
394 (ap0 * boost::math::tgamma(p0 + 1)));
395 }
397 {
398 c0 = 2.0 * (p0 + 1.0) /
399 ((2.0 * p0 + 1.0) * p0 *
400 (ap0 * boost::math::tgamma(p0 + 1)) *
401 (ap0 * boost::math::tgamma(p0 + 1)));
402 }
404 {
405 c0 = -2.0 / ((2.0 * p0 + 1.0) *
406 (ap0 * boost::math::tgamma(p0 + 1)) *
407 (ap0 * boost::math::tgamma(p0 + 1)));
408 }
410 {
411 c0 = 10000000000000000.0;
412 }
413
414 NekDouble etap0 = 0.5 * c0 * (2.0 * p0 + 1.0) *
415 (ap0 * boost::math::tgamma(p0 + 1)) *
416 (ap0 * boost::math::tgamma(p0 + 1));
417
418 NekDouble overeta0 = 1.0 / (1.0 + etap0);
419
420 // Derivative of the Legendre polynomials
421 // dLp = derivative of the Legendre polynomial of order p
422 // dLpp = derivative of the Legendre polynomial of order p+1
423 // dLpm = derivative of the Legendre polynomial of order p-1
424 Polylib::jacobd(nquad0, z0.data(), &(dLp0[0]), p0, 0.0, 0.0);
425 Polylib::jacobd(nquad0, z0.data(), &(dLpp0[0]), p0 + 1, 0.0,
426 0.0);
427 Polylib::jacobd(nquad0, z0.data(), &(dLpm0[0]), p0 - 1, 0.0,
428 0.0);
429
430 // Building the DG_c_Left
431 for (i = 0; i < nquad0; ++i)
432 {
433 m_dGL_xi1[n][i] = etap0 * dLpm0[i];
434 m_dGL_xi1[n][i] += dLpp0[i];
435 m_dGL_xi1[n][i] *= overeta0;
438 }
439
440 // Building the DG_c_Right
441 for (i = 0; i < nquad0; ++i)
442 {
443 m_dGR_xi1[n][i] = etap0 * dLpm0[i];
444 m_dGR_xi1[n][i] += dLpp0[i];
445 m_dGR_xi1[n][i] *= overeta0;
446 m_dGR_xi1[n][i] += dLp0[i];
448 }
449 }
450 break;
451 }
452 case 2:
453 {
454 LibUtilities::BasisSharedPtr BasisFR_Left0;
455 LibUtilities::BasisSharedPtr BasisFR_Right0;
456 LibUtilities::BasisSharedPtr BasisFR_Left1;
457 LibUtilities::BasisSharedPtr BasisFR_Right1;
458
459 m_dGL_xi1 = Array<OneD, Array<OneD, NekDouble>>(nElements);
460 m_dGR_xi1 = Array<OneD, Array<OneD, NekDouble>>(nElements);
461 m_dGL_xi2 = Array<OneD, Array<OneD, NekDouble>>(nElements);
462 m_dGR_xi2 = Array<OneD, Array<OneD, NekDouble>>(nElements);
463
464 for (n = 0; n < nElements; ++n)
465 {
466 base = pFields[0]->GetExp(n)->GetBase();
467 nquad0 = base[0]->GetNumPoints();
468 nquad1 = base[1]->GetNumPoints();
469 nmodes0 = base[0]->GetNumModes();
470 nmodes1 = base[1]->GetNumModes();
471
472 Array<OneD, const NekDouble> z0;
473 Array<OneD, const NekDouble> w0;
474 Array<OneD, const NekDouble> z1;
475 Array<OneD, const NekDouble> w1;
476
477 base[0]->GetZW(z0, w0);
478 base[1]->GetZW(z1, w1);
479
480 m_dGL_xi1[n] = Array<OneD, NekDouble>(nquad0);
481 m_dGR_xi1[n] = Array<OneD, NekDouble>(nquad0);
482 m_dGL_xi2[n] = Array<OneD, NekDouble>(nquad1);
483 m_dGR_xi2[n] = Array<OneD, NekDouble>(nquad1);
484
485 // Auxiliary vectors to build up the auxiliary
486 // derivatives of the Legendre polynomials
487 Array<OneD, NekDouble> dLp0(nquad0, 0.0);
488 Array<OneD, NekDouble> dLpp0(nquad0, 0.0);
489 Array<OneD, NekDouble> dLpm0(nquad0, 0.0);
490 Array<OneD, NekDouble> dLp1(nquad1, 0.0);
491 Array<OneD, NekDouble> dLpp1(nquad1, 0.0);
492 Array<OneD, NekDouble> dLpm1(nquad1, 0.0);
493
494 // Degree of the correction functions
495 int p0 = nmodes0 - 1;
496 int p1 = nmodes1 - 1;
497
498 // Function sign to compute left correction function
499 NekDouble sign0 = pow(-1.0, p0);
500 NekDouble sign1 = pow(-1.0, p1);
501
502 // Factorial factor to build the scheme
503 NekDouble ap0 = boost::math::tgamma(2 * p0 + 1) /
504 (pow(2.0, p0) * boost::math::tgamma(p0 + 1) *
505 boost::math::tgamma(p0 + 1));
506
507 NekDouble ap1 = boost::math::tgamma(2 * p1 + 1) /
508 (pow(2.0, p1) * boost::math::tgamma(p1 + 1) *
509 boost::math::tgamma(p1 + 1));
510
511 // Scalar parameter which recovers the FR schemes
513 {
514 c0 = 0.0;
515 c1 = 0.0;
516 }
518 {
519 c0 = 2.0 * p0 /
520 ((2.0 * p0 + 1.0) * (p0 + 1.0) *
521 (ap0 * boost::math::tgamma(p0 + 1)) *
522 (ap0 * boost::math::tgamma(p0 + 1)));
523
524 c1 = 2.0 * p1 /
525 ((2.0 * p1 + 1.0) * (p1 + 1.0) *
526 (ap1 * boost::math::tgamma(p1 + 1)) *
527 (ap1 * boost::math::tgamma(p1 + 1)));
528 }
530 {
531 c0 = 2.0 * (p0 + 1.0) /
532 ((2.0 * p0 + 1.0) * p0 *
533 (ap0 * boost::math::tgamma(p0 + 1)) *
534 (ap0 * boost::math::tgamma(p0 + 1)));
535
536 c1 = 2.0 * (p1 + 1.0) /
537 ((2.0 * p1 + 1.0) * p1 *
538 (ap1 * boost::math::tgamma(p1 + 1)) *
539 (ap1 * boost::math::tgamma(p1 + 1)));
540 }
542 {
543 c0 = -2.0 / ((2.0 * p0 + 1.0) *
544 (ap0 * boost::math::tgamma(p0 + 1)) *
545 (ap0 * boost::math::tgamma(p0 + 1)));
546
547 c1 = -2.0 / ((2.0 * p1 + 1.0) *
548 (ap1 * boost::math::tgamma(p1 + 1)) *
549 (ap1 * boost::math::tgamma(p1 + 1)));
550 }
552 {
553 c0 = 10000000000000000.0;
554 c1 = 10000000000000000.0;
555 }
556
557 NekDouble etap0 = 0.5 * c0 * (2.0 * p0 + 1.0) *
558 (ap0 * boost::math::tgamma(p0 + 1)) *
559 (ap0 * boost::math::tgamma(p0 + 1));
560
561 NekDouble etap1 = 0.5 * c1 * (2.0 * p1 + 1.0) *
562 (ap1 * boost::math::tgamma(p1 + 1)) *
563 (ap1 * boost::math::tgamma(p1 + 1));
564
565 NekDouble overeta0 = 1.0 / (1.0 + etap0);
566 NekDouble overeta1 = 1.0 / (1.0 + etap1);
567
568 // Derivative of the Legendre polynomials
569 // dLp = derivative of the Legendre polynomial of order p
570 // dLpp = derivative of the Legendre polynomial of order p+1
571 // dLpm = derivative of the Legendre polynomial of order p-1
572 Polylib::jacobd(nquad0, z0.data(), &(dLp0[0]), p0, 0.0, 0.0);
573 Polylib::jacobd(nquad0, z0.data(), &(dLpp0[0]), p0 + 1, 0.0,
574 0.0);
575 Polylib::jacobd(nquad0, z0.data(), &(dLpm0[0]), p0 - 1, 0.0,
576 0.0);
577
578 Polylib::jacobd(nquad1, z1.data(), &(dLp1[0]), p1, 0.0, 0.0);
579 Polylib::jacobd(nquad1, z1.data(), &(dLpp1[0]), p1 + 1, 0.0,
580 0.0);
581 Polylib::jacobd(nquad1, z1.data(), &(dLpm1[0]), p1 - 1, 0.0,
582 0.0);
583
584 // Building the DG_c_Left
585 for (i = 0; i < nquad0; ++i)
586 {
587 m_dGL_xi1[n][i] = etap0 * dLpm0[i];
588 m_dGL_xi1[n][i] += dLpp0[i];
589 m_dGL_xi1[n][i] *= overeta0;
592 }
593
594 // Building the DG_c_Left
595 for (i = 0; i < nquad1; ++i)
596 {
597 m_dGL_xi2[n][i] = etap1 * dLpm1[i];
598 m_dGL_xi2[n][i] += dLpp1[i];
599 m_dGL_xi2[n][i] *= overeta1;
602 }
603
604 // Building the DG_c_Right
605 for (i = 0; i < nquad0; ++i)
606 {
607 m_dGR_xi1[n][i] = etap0 * dLpm0[i];
608 m_dGR_xi1[n][i] += dLpp0[i];
609 m_dGR_xi1[n][i] *= overeta0;
610 m_dGR_xi1[n][i] += dLp0[i];
612 }
613
614 // Building the DG_c_Right
615 for (i = 0; i < nquad1; ++i)
616 {
617 m_dGR_xi2[n][i] = etap1 * dLpm1[i];
618 m_dGR_xi2[n][i] += dLpp1[i];
619 m_dGR_xi2[n][i] *= overeta1;
620 m_dGR_xi2[n][i] += dLp1[i];
622 }
623 }
624 break;
625 }
626 case 3:
627 {
628 m_dGL_xi1 = Array<OneD, Array<OneD, NekDouble>>(nElements);
629 m_dGR_xi1 = Array<OneD, Array<OneD, NekDouble>>(nElements);
630 m_dGL_xi2 = Array<OneD, Array<OneD, NekDouble>>(nElements);
631 m_dGR_xi2 = Array<OneD, Array<OneD, NekDouble>>(nElements);
632 m_dGL_xi3 = Array<OneD, Array<OneD, NekDouble>>(nElements);
633 m_dGR_xi3 = Array<OneD, Array<OneD, NekDouble>>(nElements);
634
635 for (n = 0; n < nElements; ++n)
636 {
637 base = pFields[0]->GetExp(n)->GetBase();
638 nquad0 = base[0]->GetNumPoints();
639 nquad1 = base[1]->GetNumPoints();
640 nquad2 = base[2]->GetNumPoints();
641 nmodes0 = base[0]->GetNumModes();
642 nmodes1 = base[1]->GetNumModes();
643 nmodes2 = base[2]->GetNumModes();
644
645 Array<OneD, const NekDouble> z0;
646 Array<OneD, const NekDouble> w0;
647 Array<OneD, const NekDouble> z1;
648 Array<OneD, const NekDouble> w1;
649 Array<OneD, const NekDouble> z2;
650 Array<OneD, const NekDouble> w2;
651
652 base[0]->GetZW(z0, w0);
653 base[1]->GetZW(z1, w1);
654 base[1]->GetZW(z2, w2);
655
656 m_dGL_xi1[n] = Array<OneD, NekDouble>(nquad0);
657 m_dGR_xi1[n] = Array<OneD, NekDouble>(nquad0);
658 m_dGL_xi2[n] = Array<OneD, NekDouble>(nquad1);
659 m_dGR_xi2[n] = Array<OneD, NekDouble>(nquad1);
660 m_dGL_xi3[n] = Array<OneD, NekDouble>(nquad2);
661 m_dGR_xi3[n] = Array<OneD, NekDouble>(nquad2);
662
663 // Auxiliary vectors to build up the auxiliary
664 // derivatives of the Legendre polynomials
665 Array<OneD, NekDouble> dLp0(nquad0, 0.0);
666 Array<OneD, NekDouble> dLpp0(nquad0, 0.0);
667 Array<OneD, NekDouble> dLpm0(nquad0, 0.0);
668 Array<OneD, NekDouble> dLp1(nquad1, 0.0);
669 Array<OneD, NekDouble> dLpp1(nquad1, 0.0);
670 Array<OneD, NekDouble> dLpm1(nquad1, 0.0);
671 Array<OneD, NekDouble> dLp2(nquad2, 0.0);
672 Array<OneD, NekDouble> dLpp2(nquad2, 0.0);
673 Array<OneD, NekDouble> dLpm2(nquad2, 0.0);
674
675 // Degree of the correction functions
676 int p0 = nmodes0 - 1;
677 int p1 = nmodes1 - 1;
678 int p2 = nmodes2 - 1;
679
680 // Function sign to compute left correction function
681 NekDouble sign0 = pow(-1.0, p0);
682 NekDouble sign1 = pow(-1.0, p1);
683 NekDouble sign2 = pow(-1.0, p2);
684
685 // Factorial factor to build the scheme
686 NekDouble ap0 = boost::math::tgamma(2 * p0 + 1) /
687 (pow(2.0, p0) * boost::math::tgamma(p0 + 1) *
688 boost::math::tgamma(p0 + 1));
689
690 // Factorial factor to build the scheme
691 NekDouble ap1 = boost::math::tgamma(2 * p1 + 1) /
692 (pow(2.0, p1) * boost::math::tgamma(p1 + 1) *
693 boost::math::tgamma(p1 + 1));
694
695 // Factorial factor to build the scheme
696 NekDouble ap2 = boost::math::tgamma(2 * p2 + 1) /
697 (pow(2.0, p2) * boost::math::tgamma(p2 + 1) *
698 boost::math::tgamma(p2 + 1));
699
700 // Scalar parameter which recovers the FR schemes
702 {
703 c0 = 0.0;
704 c1 = 0.0;
705 c2 = 0.0;
706 }
708 {
709 c0 = 2.0 * p0 /
710 ((2.0 * p0 + 1.0) * (p0 + 1.0) *
711 (ap0 * boost::math::tgamma(p0 + 1)) *
712 (ap0 * boost::math::tgamma(p0 + 1)));
713
714 c1 = 2.0 * p1 /
715 ((2.0 * p1 + 1.0) * (p1 + 1.0) *
716 (ap1 * boost::math::tgamma(p1 + 1)) *
717 (ap1 * boost::math::tgamma(p1 + 1)));
718
719 c2 = 2.0 * p2 /
720 ((2.0 * p2 + 1.0) * (p2 + 1.0) *
721 (ap2 * boost::math::tgamma(p2 + 1)) *
722 (ap2 * boost::math::tgamma(p2 + 1)));
723 }
725 {
726 c0 = 2.0 * (p0 + 1.0) /
727 ((2.0 * p0 + 1.0) * p0 *
728 (ap0 * boost::math::tgamma(p0 + 1)) *
729 (ap0 * boost::math::tgamma(p0 + 1)));
730
731 c1 = 2.0 * (p1 + 1.0) /
732 ((2.0 * p1 + 1.0) * p1 *
733 (ap1 * boost::math::tgamma(p1 + 1)) *
734 (ap1 * boost::math::tgamma(p1 + 1)));
735
736 c2 = 2.0 * (p2 + 1.0) /
737 ((2.0 * p2 + 1.0) * p2 *
738 (ap2 * boost::math::tgamma(p2 + 1)) *
739 (ap2 * boost::math::tgamma(p2 + 1)));
740 }
742 {
743 c0 = -2.0 / ((2.0 * p0 + 1.0) *
744 (ap0 * boost::math::tgamma(p0 + 1)) *
745 (ap0 * boost::math::tgamma(p0 + 1)));
746
747 c1 = -2.0 / ((2.0 * p1 + 1.0) *
748 (ap1 * boost::math::tgamma(p1 + 1)) *
749 (ap1 * boost::math::tgamma(p1 + 1)));
750
751 c2 = -2.0 / ((2.0 * p2 + 1.0) *
752 (ap2 * boost::math::tgamma(p2 + 1)) *
753 (ap2 * boost::math::tgamma(p2 + 1)));
754 }
756 {
757 c0 = 10000000000000000.0;
758 c1 = 10000000000000000.0;
759 c2 = 10000000000000000.0;
760 }
761
762 NekDouble etap0 = 0.5 * c0 * (2.0 * p0 + 1.0) *
763 (ap0 * boost::math::tgamma(p0 + 1)) *
764 (ap0 * boost::math::tgamma(p0 + 1));
765
766 NekDouble etap1 = 0.5 * c1 * (2.0 * p1 + 1.0) *
767 (ap1 * boost::math::tgamma(p1 + 1)) *
768 (ap1 * boost::math::tgamma(p1 + 1));
769
770 NekDouble etap2 = 0.5 * c2 * (2.0 * p2 + 1.0) *
771 (ap2 * boost::math::tgamma(p2 + 1)) *
772 (ap2 * boost::math::tgamma(p2 + 1));
773
774 NekDouble overeta0 = 1.0 / (1.0 + etap0);
775 NekDouble overeta1 = 1.0 / (1.0 + etap1);
776 NekDouble overeta2 = 1.0 / (1.0 + etap2);
777
778 // Derivative of the Legendre polynomials
779 // dLp = derivative of the Legendre polynomial of order p
780 // dLpp = derivative of the Legendre polynomial of order p+1
781 // dLpm = derivative of the Legendre polynomial of order p-1
782 Polylib::jacobd(nquad0, z0.data(), &(dLp0[0]), p0, 0.0, 0.0);
783 Polylib::jacobd(nquad0, z0.data(), &(dLpp0[0]), p0 + 1, 0.0,
784 0.0);
785 Polylib::jacobd(nquad0, z0.data(), &(dLpm0[0]), p0 - 1, 0.0,
786 0.0);
787
788 Polylib::jacobd(nquad1, z1.data(), &(dLp1[0]), p1, 0.0, 0.0);
789 Polylib::jacobd(nquad1, z1.data(), &(dLpp1[0]), p1 + 1, 0.0,
790 0.0);
791 Polylib::jacobd(nquad1, z1.data(), &(dLpm1[0]), p1 - 1, 0.0,
792 0.0);
793
794 Polylib::jacobd(nquad2, z2.data(), &(dLp2[0]), p2, 0.0, 0.0);
795 Polylib::jacobd(nquad2, z2.data(), &(dLpp2[0]), p2 + 1, 0.0,
796 0.0);
797 Polylib::jacobd(nquad2, z2.data(), &(dLpm2[0]), p2 - 1, 0.0,
798 0.0);
799
800 // Building the DG_c_Left
801 for (i = 0; i < nquad0; ++i)
802 {
803 m_dGL_xi1[n][i] = etap0 * dLpm0[i];
804 m_dGL_xi1[n][i] += dLpp0[i];
805 m_dGL_xi1[n][i] *= overeta0;
808 }
809
810 // Building the DG_c_Left
811 for (i = 0; i < nquad1; ++i)
812 {
813 m_dGL_xi2[n][i] = etap1 * dLpm1[i];
814 m_dGL_xi2[n][i] += dLpp1[i];
815 m_dGL_xi2[n][i] *= overeta1;
818 }
819
820 // Building the DG_c_Left
821 for (i = 0; i < nquad2; ++i)
822 {
823 m_dGL_xi3[n][i] = etap2 * dLpm2[i];
824 m_dGL_xi3[n][i] += dLpp2[i];
825 m_dGL_xi3[n][i] *= overeta2;
828 }
829
830 // Building the DG_c_Right
831 for (i = 0; i < nquad0; ++i)
832 {
833 m_dGR_xi1[n][i] = etap0 * dLpm0[i];
834 m_dGR_xi1[n][i] += dLpp0[i];
835 m_dGR_xi1[n][i] *= overeta0;
836 m_dGR_xi1[n][i] += dLp0[i];
838 }
839
840 // Building the DG_c_Right
841 for (i = 0; i < nquad1; ++i)
842 {
843 m_dGR_xi2[n][i] = etap1 * dLpm1[i];
844 m_dGR_xi2[n][i] += dLpp1[i];
845 m_dGR_xi2[n][i] *= overeta1;
846 m_dGR_xi2[n][i] += dLp1[i];
848 }
849
850 // Building the DG_c_Right
851 for (i = 0; i < nquad2; ++i)
852 {
853 m_dGR_xi3[n][i] = etap2 * dLpm2[i];
854 m_dGR_xi3[n][i] += dLpp2[i];
855 m_dGR_xi3[n][i] *= overeta2;
856 m_dGR_xi3[n][i] += dLp2[i];
858 }
859 }
860 break;
861 }
862 default:
863 {
865 break;
866 }
867 }
868}
Array< OneD, Array< OneD, NekDouble > > m_dGR_xi3
Definition: DiffusionLFR.h:88
Array< OneD, Array< OneD, NekDouble > > m_dGL_xi3
Definition: DiffusionLFR.h:87
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:1318
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 166 of file DiffusionLFR.cpp.
169{
170 boost::ignore_unused(pSession);
171
172 int i, n;
173 int nquad0, nquad1;
174 int phys_offset;
175 int nLocalSolutionPts;
176 int nElements = pFields[0]->GetExpSize();
177 int nDimensions = pFields[0]->GetCoordim(0);
178 int nSolutionPts = pFields[0]->GetTotPoints();
179 int nTracePts = pFields[0]->GetTrace()->GetTotPoints();
180
181 m_traceNormals = Array<OneD, Array<OneD, NekDouble>>(nDimensions);
182 for (i = 0; i < nDimensions; ++i)
183 {
184 m_traceNormals[i] = Array<OneD, NekDouble>(nTracePts);
185 }
186 pFields[0]->GetTrace()->GetNormals(m_traceNormals);
187
188 Array<TwoD, const NekDouble> gmat;
189 Array<OneD, const NekDouble> jac;
190
191 m_jac = Array<OneD, NekDouble>(nSolutionPts);
192
193 Array<OneD, NekDouble> auxArray1;
194 Array<OneD, LibUtilities::BasisSharedPtr> base;
195 LibUtilities::PointsKeyVector ptsKeys;
196
197 switch (nDimensions)
198 {
199 case 1:
200 {
201 for (n = 0; n < nElements; ++n)
202 {
203 ptsKeys = pFields[0]->GetExp(n)->GetPointsKeys();
204 nLocalSolutionPts = pFields[0]->GetExp(n)->GetTotPoints();
205 phys_offset = pFields[0]->GetPhys_Offset(n);
206 jac = pFields[0]
207 ->GetExp(n)
208 ->as<LocalRegions::Expansion1D>()
209 ->GetGeom1D()
210 ->GetMetricInfo()
211 ->GetJac(ptsKeys);
212 for (i = 0; i < nLocalSolutionPts; ++i)
213 {
214 m_jac[i + phys_offset] = jac[0];
215 }
216 }
217 break;
218 }
219 case 2:
220 {
221 m_gmat = Array<OneD, Array<OneD, NekDouble>>(4);
222 m_gmat[0] = Array<OneD, NekDouble>(nSolutionPts);
223 m_gmat[1] = Array<OneD, NekDouble>(nSolutionPts);
224 m_gmat[2] = Array<OneD, NekDouble>(nSolutionPts);
225 m_gmat[3] = Array<OneD, NekDouble>(nSolutionPts);
226
227 m_Q2D_e0 = Array<OneD, Array<OneD, NekDouble>>(nElements);
228 m_Q2D_e1 = Array<OneD, Array<OneD, NekDouble>>(nElements);
229 m_Q2D_e2 = Array<OneD, Array<OneD, NekDouble>>(nElements);
230 m_Q2D_e3 = Array<OneD, Array<OneD, NekDouble>>(nElements);
231
232 for (n = 0; n < nElements; ++n)
233 {
234 base = pFields[0]->GetExp(n)->GetBase();
235 nquad0 = base[0]->GetNumPoints();
236 nquad1 = base[1]->GetNumPoints();
237
238 m_Q2D_e0[n] = Array<OneD, NekDouble>(nquad0);
239 m_Q2D_e1[n] = Array<OneD, NekDouble>(nquad1);
240 m_Q2D_e2[n] = Array<OneD, NekDouble>(nquad0);
241 m_Q2D_e3[n] = Array<OneD, NekDouble>(nquad1);
242
243 // Extract the Q factors at each edge point
244 pFields[0]->GetExp(n)->GetTraceQFactors(0, auxArray1 =
245 m_Q2D_e0[n]);
246 pFields[0]->GetExp(n)->GetTraceQFactors(1, auxArray1 =
247 m_Q2D_e1[n]);
248 pFields[0]->GetExp(n)->GetTraceQFactors(2, auxArray1 =
249 m_Q2D_e2[n]);
250 pFields[0]->GetExp(n)->GetTraceQFactors(3, auxArray1 =
251 m_Q2D_e3[n]);
252
253 ptsKeys = pFields[0]->GetExp(n)->GetPointsKeys();
254 nLocalSolutionPts = pFields[0]->GetExp(n)->GetTotPoints();
255 phys_offset = pFields[0]->GetPhys_Offset(n);
256
257 jac = pFields[0]
258 ->GetExp(n)
259 ->as<LocalRegions::Expansion2D>()
260 ->GetGeom2D()
261 ->GetMetricInfo()
262 ->GetJac(ptsKeys);
263 gmat = pFields[0]
264 ->GetExp(n)
265 ->as<LocalRegions::Expansion2D>()
266 ->GetGeom2D()
267 ->GetMetricInfo()
268 ->GetDerivFactors(ptsKeys);
269
270 if (pFields[0]
271 ->GetExp(n)
272 ->as<LocalRegions::Expansion2D>()
273 ->GetGeom2D()
274 ->GetMetricInfo()
275 ->GetGtype() == SpatialDomains::eDeformed)
276 {
277 for (i = 0; i < nLocalSolutionPts; ++i)
278 {
279 m_jac[i + phys_offset] = jac[i];
280 m_gmat[0][i + phys_offset] = gmat[0][i];
281 m_gmat[1][i + phys_offset] = gmat[1][i];
282 m_gmat[2][i + phys_offset] = gmat[2][i];
283 m_gmat[3][i + phys_offset] = gmat[3][i];
284 }
285 }
286 else
287 {
288 for (i = 0; i < nLocalSolutionPts; ++i)
289 {
290 m_jac[i + phys_offset] = jac[0];
291 m_gmat[0][i + phys_offset] = gmat[0][0];
292 m_gmat[1][i + phys_offset] = gmat[1][0];
293 m_gmat[2][i + phys_offset] = gmat[2][0];
294 m_gmat[3][i + phys_offset] = gmat[3][0];
295 }
296 }
297 }
298 break;
299 }
300 case 3:
301 {
303 break;
304 }
305 default:
306 {
308 break;
309 }
310 }
311}
Array< OneD, Array< OneD, NekDouble > > m_gmat
Definition: DiffusionLFR.h:76
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 linear problems using an LDG interface flux.
Definition at line 875 of file DiffusionLFR.cpp.
882{
883 boost::ignore_unused(pFwd, pBwd);
884
885 int i, j, n;
886 int phys_offset;
887 Array<OneD, NekDouble> auxArray1, auxArray2;
888
889 Array<OneD, LibUtilities::BasisSharedPtr> Basis;
890 Basis = fields[0]->GetExp(0)->GetBase();
891
892 int nElements = fields[0]->GetExpSize();
893 int nDim = fields[0]->GetCoordim(0);
894 int nSolutionPts = fields[0]->GetTotPoints();
895 int nCoeffs = fields[0]->GetNcoeffs();
896
897 Array<OneD, Array<OneD, NekDouble>> outarrayCoeff(nConvectiveFields);
898 for (i = 0; i < nConvectiveFields; ++i)
899 {
900 outarrayCoeff[i] = Array<OneD, NekDouble>(nCoeffs);
901 }
902
903 // Compute interface numerical fluxes for inarray in physical space
905
906 switch (nDim)
907 {
908 // 1D problems
909 case 1:
910 {
911 for (i = 0; i < nConvectiveFields; ++i)
912 {
913 // Computing the physical first-order discountinuous
914 // derivative
915 for (n = 0; n < nElements; n++)
916 {
917 phys_offset = fields[0]->GetPhys_Offset(n);
918
919 fields[i]->GetExp(n)->PhysDeriv(
920 0, auxArray1 = inarray[i] + phys_offset,
921 auxArray2 = m_DU1[i][0] + phys_offset);
922 }
923
924 // Computing the standard first-order correction
925 // derivative
927 m_DFC1[i][0]);
928
929 // Back to the physical space using global operations
931 &m_DFC1[i][0][0], 1);
933 &m_DFC1[i][0][0], 1);
934
935 // Computing total first order derivatives
937 1, &m_D1[i][0][0], 1);
938
940 1);
941 }
942
943 // Computing interface numerical fluxes for m_D1
944 // in physical space
946
947 for (i = 0; i < nConvectiveFields; ++i)
948 {
949 // Computing the physical second-order discountinuous
950 // derivative
951 for (n = 0; n < nElements; n++)
952 {
953 phys_offset = fields[0]->GetPhys_Offset(n);
954
955 fields[i]->GetExp(n)->PhysDeriv(
956 0, auxArray1 = m_D1[i][0] + phys_offset,
957 auxArray2 = m_DD1[i][0] + phys_offset);
958 }
959
960 // Computing the standard second-order correction
961 // derivative
963 m_DFC2[i][0]);
964
965 // Back to the physical space using global operations
967 &m_DFC2[i][0][0], 1);
969 &m_DFC2[i][0][0], 1);
970
971 // Adding the total divergence to outarray (RHS)
973 1, &outarray[i][0], 1);
974
975 // Primitive Dealiasing 1D
976 if (!(Basis[0]->Collocation()))
977 {
978 fields[i]->FwdTrans(outarray[i], outarrayCoeff[i]);
979 fields[i]->BwdTrans(outarrayCoeff[i], outarray[i]);
980 }
981 }
982 break;
983 }
984 // 2D problems
985 case 2:
986 {
987 for (i = 0; i < nConvectiveFields; ++i)
988 {
989 for (j = 0; j < nDim; ++j)
990 {
991 // Temporary vectors
992 Array<OneD, NekDouble> u1_hat(nSolutionPts, 0.0);
993 Array<OneD, NekDouble> u2_hat(nSolutionPts, 0.0);
994
995 if (j == 0)
996 {
997 Vmath::Vmul(nSolutionPts, &inarray[i][0], 1,
998 &m_gmat[0][0], 1, &u1_hat[0], 1);
999
1001 &u1_hat[0], 1);
1002
1003 Vmath::Vmul(nSolutionPts, &inarray[i][0], 1,
1004 &m_gmat[1][0], 1, &u2_hat[0], 1);
1005
1007 &u2_hat[0], 1);
1008 }
1009 else if (j == 1)
1010 {
1011 Vmath::Vmul(nSolutionPts, &inarray[i][0], 1,
1012 &m_gmat[2][0], 1, &u1_hat[0], 1);
1013
1015 &u1_hat[0], 1);
1016
1017 Vmath::Vmul(nSolutionPts, &inarray[i][0], 1,
1018 &m_gmat[3][0], 1, &u2_hat[0], 1);
1019
1021 &u2_hat[0], 1);
1022 }
1023
1024 for (n = 0; n < nElements; n++)
1025 {
1026 phys_offset = fields[0]->GetPhys_Offset(n);
1027
1028 fields[i]->GetExp(n)->StdPhysDeriv(
1029 0, auxArray1 = u1_hat + phys_offset,
1030 auxArray2 = m_tmp1[i][j] + phys_offset);
1031
1032 fields[i]->GetExp(n)->StdPhysDeriv(
1033 1, auxArray1 = u2_hat + phys_offset,
1034 auxArray2 = m_tmp2[i][j] + phys_offset);
1035 }
1036
1039
1040 // Divide by the metric jacobian
1042 &m_DU1[i][j][0], 1);
1043
1044 // Computing the standard first-order correction
1045 // derivatives
1046 DerCFlux_2D(nConvectiveFields, j, fields, inarray[i],
1048 }
1049
1050 // Multiplying derivatives by B matrix to get auxiliary
1051 // variables
1052 for (j = 0; j < nSolutionPts; ++j)
1053 {
1054 // std(q1)
1057 m_DFC1[i][0][j] +
1060 m_DFC1[i][1][j];
1061
1062 // std(q2)
1065 m_DFC1[i][0][j] +
1068 m_DFC1[i][1][j];
1069 }
1070
1071 // Multiplying derivatives by A^(-1) to get back
1072 // into the physical space
1073 for (j = 0; j < nSolutionPts; j++)
1074 {
1075 // q1 = A11^(-1)*std(q1) + A12^(-1)*std(q2)
1078
1079 // q2 = A21^(-1)*std(q1) + A22^(-1)*std(q2)
1082 }
1083
1084 // Computing the physical first-order derivatives
1085 for (j = 0; j < nDim; ++j)
1086 {
1089 }
1090 }
1091
1092 // Computing interface numerical fluxes for m_D1
1093 // in physical space
1095
1096 // Computing the standard second-order discontinuous
1097 // derivatives
1098 for (i = 0; i < nConvectiveFields; ++i)
1099 {
1100 Array<OneD, NekDouble> f_hat(nSolutionPts, 0.0);
1101 Array<OneD, NekDouble> g_hat(nSolutionPts, 0.0);
1102
1103 for (j = 0; j < nSolutionPts; j++)
1104 {
1107 m_jac[j];
1108
1111 m_jac[j];
1112 }
1113
1114 for (n = 0; n < nElements; n++)
1115 {
1116 phys_offset = fields[0]->GetPhys_Offset(n);
1117
1118 fields[0]->GetExp(n)->StdPhysDeriv(
1119 0, auxArray1 = f_hat + phys_offset,
1120 auxArray2 = m_DD1[i][0] + phys_offset);
1121
1122 fields[0]->GetExp(n)->StdPhysDeriv(
1123 1, auxArray1 = g_hat + phys_offset,
1124 auxArray2 = m_DD1[i][1] + phys_offset);
1125 }
1126
1127 // Divergence of the standard discontinuous flux
1129 1, &m_divFD[i][0], 1);
1130
1131 // Divergence of the standard correction flux
1132 if (Basis[0]->GetPointsType() ==
1134 Basis[1]->GetPointsType() ==
1136 {
1137 DivCFlux_2D_Gauss(nConvectiveFields, fields, f_hat, g_hat,
1139 }
1140 else
1141 {
1144 }
1145
1146 // Divergence of the standard final flux
1148 &outarray[i][0], 1);
1149
1150 // Dividing by the jacobian to get back into
1151 // physical space
1153 &outarray[i][0], 1);
1154
1155 // Primitive Dealiasing 2D
1156 if (!(Basis[0]->Collocation()))
1157 {
1158 fields[i]->FwdTrans(outarray[i], outarrayCoeff[i]);
1159 fields[i]->BwdTrans(outarrayCoeff[i], outarray[i]);
1160 }
1161 } // Close loop on nConvectiveFields
1162 break;
1163 }
1164 // 3D problems
1165 case 3:
1166 {
1168 break;
1169 }
1170 }
1171}
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > m_tmp2
Definition: DiffusionLFR.h:104
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > m_IF1
Definition: DiffusionLFR.h:92
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: DiffusionLFR.cpp:1457
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > m_DFC1
Definition: DiffusionLFR.h:94
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > m_DD1
Definition: DiffusionLFR.h:97
Array< OneD, Array< OneD, NekDouble > > m_divFC
Definition: DiffusionLFR.h:101
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: DiffusionLFR.cpp:1578
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: DiffusionLFR.cpp:1776
void NumFluxforVector(const Array< OneD, MultiRegions::ExpListSharedPtr > &fields, const Array< OneD, Array< OneD, NekDouble > > &ufield, Array< OneD, Array< OneD, Array< OneD, NekDouble > > > &qfield, Array< OneD, Array< OneD, NekDouble > > &qflux)
Build the numerical flux for the 2nd order derivatives.
Definition: DiffusionLFR.cpp:1314
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > m_DU1
Definition: DiffusionLFR.h:93
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > m_DFC2
Definition: DiffusionLFR.h:99
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > m_tmp1
Definition: DiffusionLFR.h:103
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > m_D1
Definition: DiffusionLFR.h:96
void NumFluxforScalar(const Array< OneD, MultiRegions::ExpListSharedPtr > &fields, const Array< OneD, Array< OneD, NekDouble > > &ufield, Array< OneD, Array< OneD, Array< OneD, NekDouble > > > &uflux)
Builds the numerical flux for the 1st order derivatives.
Definition: DiffusionLFR.cpp:1177
Array< OneD, Array< OneD, NekDouble > > m_IF2
Definition: DiffusionLFR.h:98
Array< OneD, Array< OneD, NekDouble > > m_divFD
Definition: DiffusionLFR.h:100
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > m_BD1
Definition: DiffusionLFR.h:95
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: DiffusionLFR.cpp:1966
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:280
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_divFC, m_divFD, m_DU1, m_gmat, m_IF1, m_IF2, m_jac, m_tmp1, m_tmp2, NumFluxforScalar(), NumFluxforVector(), Vmath::Vadd(), Vmath::Vcopy(), Vmath::Vdiv(), and Vmath::Vmul().
◆ v_InitObject()
|
overrideprotectedvirtual |
Initiliase DiffusionLFR objects and store them before starting the time-stepping.
This routine calls the functions SetupMetrics and SetupCFunctions to initialise the objects needed by DiffusionLFR.
- Parameters
-
pSession Pointer to session reader. pFields Pointer to fields.
Implements Nektar::SolverUtils::Diffusion.
Definition at line 85 of file DiffusionLFR.cpp.
88{
89 m_session = pSession;
90 SetupMetrics(pSession, pFields);
91 SetupCFunctions(pSession, pFields);
92
93 // Initialising arrays
94 int i, j;
95 int nConvectiveFields = pFields.size();
96 int nDim = pFields[0]->GetCoordim(0);
97 int nSolutionPts = pFields[0]->GetTotPoints();
98 int nTracePts = pFields[0]->GetTrace()->GetTotPoints();
99
100 m_IF1 = Array<OneD, Array<OneD, Array<OneD, NekDouble>>>(nConvectiveFields);
101 m_DU1 = Array<OneD, Array<OneD, Array<OneD, NekDouble>>>(nConvectiveFields);
102 m_DFC1 =
103 Array<OneD, Array<OneD, Array<OneD, NekDouble>>>(nConvectiveFields);
104 m_BD1 = Array<OneD, Array<OneD, Array<OneD, NekDouble>>>(nConvectiveFields);
105 m_D1 = Array<OneD, Array<OneD, Array<OneD, NekDouble>>>(nConvectiveFields);
106 m_DD1 = Array<OneD, Array<OneD, Array<OneD, NekDouble>>>(nConvectiveFields);
107 m_IF2 = Array<OneD, Array<OneD, NekDouble>>(nConvectiveFields);
108 m_DFC2 =
109 Array<OneD, Array<OneD, Array<OneD, NekDouble>>>(nConvectiveFields);
110 m_divFD = Array<OneD, Array<OneD, NekDouble>>(nConvectiveFields);
111 m_divFC = Array<OneD, Array<OneD, NekDouble>>(nConvectiveFields);
112
113 m_tmp1 =
114 Array<OneD, Array<OneD, Array<OneD, NekDouble>>>(nConvectiveFields);
115 m_tmp2 =
116 Array<OneD, Array<OneD, Array<OneD, NekDouble>>>(nConvectiveFields);
117
118 for (i = 0; i < nConvectiveFields; ++i)
119 {
120 m_IF1[i] = Array<OneD, Array<OneD, NekDouble>>(nDim);
121 m_DU1[i] = Array<OneD, Array<OneD, NekDouble>>(nDim);
122 m_DFC1[i] = Array<OneD, Array<OneD, NekDouble>>(nDim);
123 m_BD1[i] = Array<OneD, Array<OneD, NekDouble>>(nDim);
124 m_D1[i] = Array<OneD, Array<OneD, NekDouble>>(nDim);
125 m_DD1[i] = Array<OneD, Array<OneD, NekDouble>>(nDim);
126 m_IF2[i] = Array<OneD, NekDouble>(nTracePts, 0.0);
127 m_DFC2[i] = Array<OneD, Array<OneD, NekDouble>>(nDim);
128 m_divFD[i] = Array<OneD, NekDouble>(nSolutionPts, 0.0);
129 m_divFC[i] = Array<OneD, NekDouble>(nSolutionPts, 0.0);
130
131 m_tmp1[i] = Array<OneD, Array<OneD, NekDouble>>(nDim);
132 m_tmp2[i] = Array<OneD, Array<OneD, NekDouble>>(nDim);
133
134 for (j = 0; j < nDim; ++j)
135 {
136 m_IF1[i][j] = Array<OneD, NekDouble>(nTracePts, 0.0);
137 m_DU1[i][j] = Array<OneD, NekDouble>(nSolutionPts, 0.0);
138 m_DFC1[i][j] = Array<OneD, NekDouble>(nSolutionPts, 0.0);
139 m_BD1[i][j] = Array<OneD, NekDouble>(nSolutionPts, 0.0);
140 m_D1[i][j] = Array<OneD, NekDouble>(nSolutionPts, 0.0);
141 m_DD1[i][j] = Array<OneD, NekDouble>(nSolutionPts, 0.0);
142 m_DFC2[i][j] = Array<OneD, NekDouble>(nSolutionPts, 0.0);
143
144 m_tmp1[i][j] = Array<OneD, NekDouble>(nSolutionPts, 0.0);
145 m_tmp2[i][j] = Array<OneD, NekDouble>(nSolutionPts, 0.0);
146 }
147 }
148}
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: DiffusionLFR.cpp:166
LibUtilities::SessionReaderSharedPtr m_session
Definition: DiffusionLFR.h:73
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: DiffusionLFR.cpp:330
References m_BD1, m_D1, m_DD1, m_DFC1, m_DFC2, m_divFC, m_divFD, m_DU1, m_IF1, m_IF2, m_session, m_tmp1, m_tmp2, SetupCFunctions(), and SetupMetrics().
◆ WeakPenaltyforScalar()
|
private |
Imposes appropriate bcs for the 1st order derivatives.
Definition at line 1251 of file DiffusionLFR.cpp.
1255{
1256 // Variable initialisation
1257 int i, e, id1, id2;
1258 int nBndEdgePts, nBndEdges;
1259 int cnt = 0;
1260 int nBndRegions = fields[var]->GetBndCondExpansions().size();
1261 int nTracePts = fields[0]->GetTrace()->GetTotPoints();
1262
1263 // Extract physical values of the fields at the trace space
1264 Array<OneD, NekDouble> uplus(nTracePts);
1265 fields[var]->ExtractTracePhys(ufield, uplus);
1266
1267 // Impose boundary conditions
1268 for (i = 0; i < nBndRegions; ++i)
1269 {
1270 // Number of boundary edges related to bcRegion 'i'
1271 nBndEdges = fields[var]->GetBndCondExpansions()[i]->GetExpSize();
1272
1273 // Weakly impose boundary conditions by modifying flux values
1274 for (e = 0; e < nBndEdges; ++e)
1275 {
1276 // Number of points on the expansion
1277 nBndEdgePts = fields[var]
1278 ->GetBndCondExpansions()[i]
1279 ->GetExp(e)
1280 ->GetTotPoints();
1281
1282 // Offset of the boundary expansion
1283 id1 = fields[var]->GetBndCondExpansions()[i]->GetPhys_Offset(e);
1284
1285 // Offset of the trace space related to boundary expansion
1286 id2 = fields[0]->GetTrace()->GetPhys_Offset(
1287 fields[0]->GetTraceMap()->GetBndCondIDToGlobalTraceID(cnt++));
1288
1289 // Dirichlet bcs ==> uflux = gD
1290 if (fields[var]
1291 ->GetBndConditions()[i]
1292 ->GetBoundaryConditionType() == SpatialDomains::eDirichlet)
1293 {
1294 Vmath::Vcopy(
1295 nBndEdgePts,
1296 &(fields[var]->GetBndCondExpansions()[i]->GetPhys())[id1],
1297 1, &penaltyflux[id2], 1);
1298 }
1299 // Neumann bcs ==> uflux = u+
1300 else if ((fields[var]->GetBndConditions()[i])
1301 ->GetBoundaryConditionType() ==
1302 SpatialDomains::eNeumann)
1303 {
1304 Vmath::Vcopy(nBndEdgePts, &uplus[id2], 1, &penaltyflux[id2], 1);
1305 }
1306 }
1307 }
1308}
References Nektar::SpatialDomains::eDirichlet, Nektar::SpatialDomains::eNeumann, and Vmath::Vcopy().
Referenced by NumFluxforScalar().
◆ WeakPenaltyforVector()
|
private |
Imposes appropriate bcs for the 2nd order derivatives.
Definition at line 1389 of file DiffusionLFR.cpp.
1393{
1394 boost::ignore_unused(C11);
1395
1396 int i, e, id1, id2;
1397 int nBndEdges, nBndEdgePts;
1398 int nBndRegions = fields[var]->GetBndCondExpansions().size();
1399 int nTracePts = fields[0]->GetTrace()->GetTotPoints();
1400
1401 Array<OneD, NekDouble> uterm(nTracePts);
1402 Array<OneD, NekDouble> qtemp(nTracePts);
1403 int cnt = 0;
1404
1405 fields[var]->ExtractTracePhys(qfield, qtemp);
1406
1407 for (i = 0; i < nBndRegions; ++i)
1408 {
1409 nBndEdges = fields[var]->GetBndCondExpansions()[i]->GetExpSize();
1410
1411 // Weakly impose boundary conditions by modifying flux values
1412 for (e = 0; e < nBndEdges; ++e)
1413 {
1414 nBndEdgePts = fields[var]
1415 ->GetBndCondExpansions()[i]
1416 ->GetExp(e)
1417 ->GetTotPoints();
1418
1419 id1 = fields[var]->GetBndCondExpansions()[i]->GetPhys_Offset(e);
1420
1421 id2 = fields[0]->GetTrace()->GetPhys_Offset(
1422 fields[0]->GetTraceMap()->GetBndCondIDToGlobalTraceID(cnt++));
1423
1424 // For Dirichlet boundary condition:
1425 // qflux = q+ - C_11 (u+ - g_D) (nx, ny)
1426 if (fields[var]
1427 ->GetBndConditions()[i]
1428 ->GetBoundaryConditionType() == SpatialDomains::eDirichlet)
1429 {
1431 &qtemp[id2], 1, &penaltyflux[id2], 1);
1432 }
1433 // For Neumann boundary condition: qflux = g_N
1434 else if ((fields[var]->GetBndConditions()[i])
1435 ->GetBoundaryConditionType() ==
1436 SpatialDomains::eNeumann)
1437 {
1438 Vmath::Vmul(
1439 nBndEdgePts, &m_traceNormals[dir][id2], 1,
1440 &(fields[var]->GetBndCondExpansions()[i]->GetPhys())[id1],
1441 1, &penaltyflux[id2], 1);
1442 }
1443 }
1444 }
1445}
References Nektar::SpatialDomains::eDirichlet, Nektar::SpatialDomains::eNeumann, m_traceNormals, and Vmath::Vmul().
Referenced by NumFluxforVector().
Member Data Documentation
◆ m_BD1
Definition at line 95 of file DiffusionLFR.h.
Referenced by v_Diffuse(), and v_InitObject().
◆ m_D1
Definition at line 96 of file DiffusionLFR.h.
Referenced by v_Diffuse(), and v_InitObject().
◆ m_DD1
Definition at line 97 of file DiffusionLFR.h.
Referenced by v_Diffuse(), and v_InitObject().
◆ m_DFC1
|
private |
Definition at line 94 of file DiffusionLFR.h.
Referenced by v_Diffuse(), and v_InitObject().
◆ m_DFC2
|
private |
Definition at line 99 of file DiffusionLFR.h.
Referenced by v_Diffuse(), and v_InitObject().
◆ m_dGL_xi1
Definition at line 83 of file DiffusionLFR.h.
Referenced by DerCFlux_1D(), DerCFlux_2D(), DivCFlux_2D(), DivCFlux_2D_Gauss(), and SetupCFunctions().
◆ m_dGL_xi2
Definition at line 85 of file DiffusionLFR.h.
Referenced by DerCFlux_2D(), DivCFlux_2D(), DivCFlux_2D_Gauss(), and SetupCFunctions().
◆ m_dGL_xi3
Definition at line 87 of file DiffusionLFR.h.
Referenced by SetupCFunctions().
◆ m_dGR_xi1
Definition at line 84 of file DiffusionLFR.h.
Referenced by DerCFlux_1D(), DerCFlux_2D(), DivCFlux_2D(), DivCFlux_2D_Gauss(), and SetupCFunctions().
◆ m_dGR_xi2
Definition at line 86 of file DiffusionLFR.h.
Referenced by DerCFlux_2D(), DivCFlux_2D(), DivCFlux_2D_Gauss(), and SetupCFunctions().
◆ m_dGR_xi3
Definition at line 88 of file DiffusionLFR.h.
Referenced by SetupCFunctions().
◆ m_diffType
|
protected |
Definition at line 57 of file DiffusionLFR.h.
Referenced by SetupCFunctions().
◆ m_divFC
Definition at line 101 of file DiffusionLFR.h.
Referenced by v_Diffuse(), and v_InitObject().
◆ m_divFD
Definition at line 100 of file DiffusionLFR.h.
Referenced by v_Diffuse(), and v_InitObject().
◆ m_DU1
Definition at line 93 of file DiffusionLFR.h.
Referenced by v_Diffuse(), and v_InitObject().
◆ m_gmat
Definition at line 76 of file DiffusionLFR.h.
Referenced by SetupMetrics(), and v_Diffuse().
◆ m_IF1
Definition at line 92 of file DiffusionLFR.h.
Referenced by v_Diffuse(), and v_InitObject().
◆ m_IF2
Definition at line 98 of file DiffusionLFR.h.
Referenced by v_Diffuse(), and v_InitObject().
◆ m_Ixm
|
private |
Definition at line 89 of file DiffusionLFR.h.
◆ m_Ixp
|
private |
Definition at line 90 of file DiffusionLFR.h.
◆ m_jac
Definition at line 75 of file DiffusionLFR.h.
Referenced by SetupMetrics(), and v_Diffuse().
◆ m_Q2D_e0
Definition at line 78 of file DiffusionLFR.h.
Referenced by DivCFlux_2D(), DivCFlux_2D_Gauss(), and SetupMetrics().
◆ m_Q2D_e1
Definition at line 79 of file DiffusionLFR.h.
Referenced by DivCFlux_2D(), DivCFlux_2D_Gauss(), and SetupMetrics().
◆ m_Q2D_e2
Definition at line 80 of file DiffusionLFR.h.
Referenced by DivCFlux_2D(), DivCFlux_2D_Gauss(), and SetupMetrics().
◆ m_Q2D_e3
Definition at line 81 of file DiffusionLFR.h.
Referenced by DivCFlux_2D(), DivCFlux_2D_Gauss(), and SetupMetrics().
◆ m_session
|
private |
Definition at line 73 of file DiffusionLFR.h.
Referenced by v_InitObject().
◆ m_tmp1
|
private |
Definition at line 103 of file DiffusionLFR.h.
Referenced by v_Diffuse(), and v_InitObject().
◆ m_tmp2
|
private |
Definition at line 104 of file DiffusionLFR.h.
Referenced by v_Diffuse(), and v_InitObject().
◆ m_traceNormals
Definition at line 72 of file DiffusionLFR.h.
Referenced by DivCFlux_2D(), DivCFlux_2D_Gauss(), NumFluxforScalar(), NumFluxforVector(), SetupMetrics(), and WeakPenaltyforVector().
◆ type
|
static |
Initial value:
= {
"LFRDG"),
"LRFSD"),
"LFRHU"),
tKey RegisterCreatorFunction(tKey idKey, CreatorFunction classCreator, std::string pDesc="")
Register a class with the factory.
Definition: NekFactory.hpp:198
static DiffusionSharedPtr create(std::string diffType)
Definition: DiffusionLFR.h:47
DiffusionFactory & GetDiffusionFactory()
Definition: Diffusion.cpp:41
Definition at line 52 of file DiffusionLFR.h.
