#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... | |
| void | v_InitObject (LibUtilities::SessionReaderSharedPtr pSession, Array< OneD, MultiRegions::ExpListSharedPtr > pFields) override |
| Initiliase DiffusionLFR objects and store them before starting the time-stepping. More... | |
| 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 42 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 67 of file DiffusionLFR.cpp.
Referenced by create().
Member Function Documentation
◆ create()
|
inlinestatic |
Definition at line 45 of file DiffusionLFR.h.
46 {
48 }
DiffusionLFR(std::string diffType)
DiffusionLFR uses the Flux Reconstruction (FR) approach to compute the diffusion term....
Definition: DiffusionLFR.cpp:67
std::shared_ptr< Diffusion > DiffusionSharedPtr
A shared pointer to an EquationSystem object.
Definition: Diffusion.h:54
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 1437 of file DiffusionLFR.cpp.
1442{
1443 int n;
1444 int nLocalSolutionPts, phys_offset, t_offset;
1445
1446 Array<OneD, NekDouble> auxArray1, auxArray2;
1447 Array<TwoD, const NekDouble> gmat;
1448 Array<OneD, const NekDouble> jac;
1449
1450 LibUtilities::PointsKeyVector ptsKeys;
1451 LibUtilities::BasisSharedPtr Basis;
1452 Basis = fields[0]->GetExp(0)->GetBasis(0);
1453
1454 int nElements = fields[0]->GetExpSize();
1455 int nSolutionPts = fields[0]->GetTotPoints();
1456
1457 vector<bool> leftAdjacentTraces =
1458 std::static_pointer_cast<MultiRegions::DisContField>(fields[0])
1459 ->GetLeftAdjacentTraces();
1460
1461 // Arrays to store the derivatives of the correction flux
1462 Array<OneD, NekDouble> DCL(nSolutionPts / nElements, 0.0);
1463 Array<OneD, NekDouble> DCR(nSolutionPts / nElements, 0.0);
1464
1465 // Arrays to store the intercell numerical flux jumps
1466 Array<OneD, NekDouble> JumpL(nElements);
1467 Array<OneD, NekDouble> JumpR(nElements);
1468
1469 Array<OneD, Array<OneD, LocalRegions::ExpansionSharedPtr>> &elmtToTrace =
1470 fields[0]->GetTraceMap()->GetElmtToTrace();
1471
1472 for (n = 0; n < nElements; ++n)
1473 {
1474 nLocalSolutionPts = fields[0]->GetExp(n)->GetTotPoints();
1475 phys_offset = fields[0]->GetPhys_Offset(n);
1476
1477 Array<OneD, NekDouble> tmparrayX1(nLocalSolutionPts, 0.0);
1478 Array<OneD, NekDouble> tmpFluxVertex(1, 0.0);
1479 Vmath::Vcopy(nLocalSolutionPts, &flux[phys_offset], 1, &tmparrayX1[0],
1480 1);
1481
1482 fields[0]->GetExp(n)->GetTracePhysVals(0, elmtToTrace[n][0], tmparrayX1,
1483 tmpFluxVertex);
1484
1485 t_offset = fields[0]->GetTrace()->GetPhys_Offset(
1486 elmtToTrace[n][0]->GetElmtId());
1487
1488 if (leftAdjacentTraces[2 * n])
1489 {
1490 JumpL[n] = -iFlux[t_offset] - tmpFluxVertex[0];
1491 }
1492 else
1493 {
1494 JumpL[n] = iFlux[t_offset] - tmpFluxVertex[0];
1495 }
1496
1497 t_offset = fields[0]->GetTrace()->GetPhys_Offset(
1498 elmtToTrace[n][1]->GetElmtId());
1499
1500 fields[0]->GetExp(n)->GetTracePhysVals(1, elmtToTrace[n][1], tmparrayX1,
1501 tmpFluxVertex);
1502 if (leftAdjacentTraces[2 * n + 1])
1503 {
1504 JumpR[n] = iFlux[t_offset] - tmpFluxVertex[0];
1505 }
1506 else
1507 {
1508 JumpR[n] = -iFlux[t_offset] - tmpFluxVertex[0];
1509 }
1510 }
1511
1512 for (n = 0; n < nElements; ++n)
1513 {
1514 ptsKeys = fields[0]->GetExp(n)->GetPointsKeys();
1515 nLocalSolutionPts = fields[0]->GetExp(n)->GetTotPoints();
1516 phys_offset = fields[0]->GetPhys_Offset(n);
1517 jac = fields[0]
1518 ->GetExp(n)
1519 ->as<LocalRegions::Expansion1D>()
1520 ->GetGeom1D()
1521 ->GetMetricInfo()
1522 ->GetJac(ptsKeys);
1523
1524 JumpL[n] = JumpL[n] * jac[0];
1525 JumpR[n] = JumpR[n] * jac[0];
1526
1527 // Left jump multiplied by left derivative of C function
1529 1);
1530
1531 // Right jump multiplied by right derivative of C function
1533 1);
1534
1535 // Assembling derivative of the correction flux
1536 Vmath::Vadd(nLocalSolutionPts, &DCL[0], 1, &DCR[0], 1,
1537 &derCFlux[phys_offset], 1);
1538 }
1539}
Array< OneD, Array< OneD, NekDouble > > m_dGR_xi1
Definition: DiffusionLFR.h:81
Array< OneD, Array< OneD, NekDouble > > m_dGL_xi1
Definition: DiffusionLFR.h:80
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 1556 of file DiffusionLFR.cpp.
1561{
1562 int n, e, i, j, cnt;
1563
1564 Array<OneD, const NekDouble> jac;
1565
1566 int nElements = fields[0]->GetExpSize();
1567 int trace_offset, phys_offset;
1568 int nLocalSolutionPts;
1569 int nquad0, nquad1;
1570 int nEdgePts;
1571
1572 LibUtilities::PointsKeyVector ptsKeys;
1573 Array<OneD, NekDouble> auxArray1, auxArray2;
1574 Array<OneD, LibUtilities::BasisSharedPtr> base;
1575 Array<OneD, Array<OneD, LocalRegions::ExpansionSharedPtr>> &elmtToTrace =
1576 fields[0]->GetTraceMap()->GetElmtToTrace();
1577
1578 // Loop on the elements
1579 for (n = 0; n < nElements; ++n)
1580 {
1581 // Offset of the element on the global vector
1582 phys_offset = fields[0]->GetPhys_Offset(n);
1583 nLocalSolutionPts = fields[0]->GetExp(n)->GetTotPoints();
1584 ptsKeys = fields[0]->GetExp(n)->GetPointsKeys();
1585
1586 jac = fields[0]
1587 ->GetExp(n)
1588 ->as<LocalRegions::Expansion2D>()
1589 ->GetGeom2D()
1590 ->GetMetricInfo()
1591 ->GetJac(ptsKeys);
1592
1593 base = fields[0]->GetExp(n)->GetBase();
1594 nquad0 = base[0]->GetNumPoints();
1595 nquad1 = base[1]->GetNumPoints();
1596
1597 Array<OneD, NekDouble> divCFluxE0(nLocalSolutionPts, 0.0);
1598 Array<OneD, NekDouble> divCFluxE1(nLocalSolutionPts, 0.0);
1599 Array<OneD, NekDouble> divCFluxE2(nLocalSolutionPts, 0.0);
1600 Array<OneD, NekDouble> divCFluxE3(nLocalSolutionPts, 0.0);
1601
1602 // Loop on the edges
1603 for (e = 0; e < fields[0]->GetExp(n)->GetNtraces(); ++e)
1604 {
1605 // Number of edge points of edge 'e'
1606 nEdgePts = fields[0]->GetExp(n)->GetTraceNumPoints(e);
1607
1608 // Array for storing volumetric fluxes on each edge
1609 Array<OneD, NekDouble> tmparray(nEdgePts, 0.0);
1610
1611 // Offset of the trace space correspondent to edge 'e'
1612 trace_offset = fields[0]->GetTrace()->GetPhys_Offset(
1613 elmtToTrace[n][e]->GetElmtId());
1614
1615 // Extract the edge values of the volumetric fluxes
1616 // on edge 'e' and order them accordingly to the order
1617 // of the trace space
1618 fields[0]->GetExp(n)->GetTracePhysVals(
1619 e, elmtToTrace[n][e], flux + phys_offset, auxArray1 = tmparray);
1620
1621 // Compute the fluxJumps per each edge 'e' and each
1622 // flux point
1623 Array<OneD, NekDouble> fluxJumps(nEdgePts, 0.0);
1624 Vmath::Vsub(nEdgePts, &iFlux[trace_offset], 1, &tmparray[0], 1,
1625 &fluxJumps[0], 1);
1626
1627 // Check the ordering of the fluxJumps and reverse
1628 // it in case of backward definition of edge 'e'
1629 if (fields[0]->GetExp(n)->GetTraceOrient(e) ==
1630 StdRegions::eBackwards)
1631 {
1632 Vmath::Reverse(nEdgePts, &fluxJumps[0], 1, &fluxJumps[0], 1);
1633 }
1634
1635 // Deformed elements
1636 if (fields[0]
1637 ->GetExp(n)
1638 ->as<LocalRegions::Expansion2D>()
1639 ->GetGeom2D()
1640 ->GetMetricInfo()
1641 ->GetGtype() == SpatialDomains::eDeformed)
1642 {
1643 // Extract the Jacobians along edge 'e'
1644 Array<OneD, NekDouble> jacEdge(nEdgePts, 0.0);
1645
1646 fields[0]->GetExp(n)->GetTracePhysVals(
1647 e, elmtToTrace[n][e], jac, auxArray1 = jacEdge);
1648
1649 // Check the ordering of the fluxJumps and reverse
1650 // it in case of backward definition of edge 'e'
1651 if (fields[0]->GetExp(n)->GetTraceOrient(e) ==
1652 StdRegions::eBackwards)
1653 {
1654 Vmath::Reverse(nEdgePts, &jacEdge[0], 1, &jacEdge[0], 1);
1655 }
1656
1657 // Multiply the fluxJumps by the edge 'e' Jacobians
1658 // to bring the fluxJumps into the standard space
1659 for (j = 0; j < nEdgePts; j++)
1660 {
1661 fluxJumps[j] = fluxJumps[j] * jacEdge[j];
1662 }
1663 }
1664 // Non-deformed elements
1665 else
1666 {
1667 // Multiply the fluxJumps by the edge 'e' Jacobians
1668 // to bring the fluxJumps into the standard space
1669 Vmath::Smul(nEdgePts, jac[0], fluxJumps, 1, fluxJumps, 1);
1670 }
1671
1672 // Multiply jumps by derivatives of the correction functions
1673 // All the quntities at this level should be defined into
1674 // the standard space
1675 switch (e)
1676 {
1677 case 0:
1678 for (i = 0; i < nquad0; ++i)
1679 {
1680 for (j = 0; j < nquad1; ++j)
1681 {
1682 cnt = i + j * nquad0;
1683 divCFluxE0[cnt] = fluxJumps[i] * m_dGL_xi2[n][j];
1684 }
1685 }
1686 break;
1687 case 1:
1688 for (i = 0; i < nquad1; ++i)
1689 {
1690 for (j = 0; j < nquad0; ++j)
1691 {
1692 cnt = (nquad0)*i + j;
1693 divCFluxE1[cnt] = fluxJumps[i] * m_dGR_xi1[n][j];
1694 }
1695 }
1696 break;
1697 case 2:
1698 for (i = 0; i < nquad0; ++i)
1699 {
1700 for (j = 0; j < nquad1; ++j)
1701 {
1702 cnt = j * nquad0 + i;
1703 divCFluxE2[cnt] = fluxJumps[i] * m_dGR_xi2[n][j];
1704 }
1705 }
1706 break;
1707 case 3:
1708 for (i = 0; i < nquad1; ++i)
1709 {
1710 for (j = 0; j < nquad0; ++j)
1711 {
1712 cnt = j + i * nquad0;
1713 divCFluxE3[cnt] = fluxJumps[i] * m_dGL_xi1[n][j];
1714 }
1715 }
1716 break;
1717
1718 default:
1720 break;
1721 }
1722 }
1723
1724 // Summing the component of the flux for calculating the
1725 // derivatives per each direction
1726 if (direction == 0)
1727 {
1728 Vmath::Vadd(nLocalSolutionPts, &divCFluxE1[0], 1, &divCFluxE3[0], 1,
1729 &derCFlux[phys_offset], 1);
1730 }
1731 else if (direction == 1)
1732 {
1733 Vmath::Vadd(nLocalSolutionPts, &divCFluxE0[0], 1, &divCFluxE2[0], 1,
1734 &derCFlux[phys_offset], 1);
1735 }
1736 }
1737}
Array< OneD, Array< OneD, NekDouble > > m_dGL_xi2
Definition: DiffusionLFR.h:82
Array< OneD, Array< OneD, NekDouble > > m_dGR_xi2
Definition: DiffusionLFR.h:83
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 1752 of file DiffusionLFR.cpp.
1759{
1760 int n, e, i, j, cnt;
1761
1762 int nElements = fields[0]->GetExpSize();
1763
1764 int nLocalSolutionPts;
1765 int nEdgePts;
1766 int trace_offset;
1767 int phys_offset;
1768 int nquad0;
1769 int nquad1;
1770
1771 Array<OneD, NekDouble> auxArray1, auxArray2;
1772 Array<OneD, LibUtilities::BasisSharedPtr> base;
1773
1774 Array<OneD, Array<OneD, LocalRegions::ExpansionSharedPtr>> &elmtToTrace =
1775 fields[0]->GetTraceMap()->GetElmtToTrace();
1776
1777 // Loop on the elements
1778 for (n = 0; n < nElements; ++n)
1779 {
1780 // Offset of the element on the global vector
1781 phys_offset = fields[0]->GetPhys_Offset(n);
1782 nLocalSolutionPts = fields[0]->GetExp(n)->GetTotPoints();
1783
1784 base = fields[0]->GetExp(n)->GetBase();
1785 nquad0 = base[0]->GetNumPoints();
1786 nquad1 = base[1]->GetNumPoints();
1787
1788 Array<OneD, NekDouble> divCFluxE0(nLocalSolutionPts, 0.0);
1789 Array<OneD, NekDouble> divCFluxE1(nLocalSolutionPts, 0.0);
1790 Array<OneD, NekDouble> divCFluxE2(nLocalSolutionPts, 0.0);
1791 Array<OneD, NekDouble> divCFluxE3(nLocalSolutionPts, 0.0);
1792
1793 // Loop on the edges
1794 for (e = 0; e < fields[0]->GetExp(n)->GetNtraces(); ++e)
1795 {
1796 // Number of edge points of edge e
1797 nEdgePts = fields[0]->GetExp(n)->GetTraceNumPoints(e);
1798
1799 Array<OneD, NekDouble> tmparrayX1(nEdgePts, 0.0);
1800 Array<OneD, NekDouble> tmparrayX2(nEdgePts, 0.0);
1801 Array<OneD, NekDouble> fluxN(nEdgePts, 0.0);
1802 Array<OneD, NekDouble> fluxT(nEdgePts, 0.0);
1803 Array<OneD, NekDouble> fluxJumps(nEdgePts, 0.0);
1804
1805 // Offset of the trace space correspondent to edge e
1806 trace_offset = fields[0]->GetTrace()->GetPhys_Offset(
1807 elmtToTrace[n][e]->GetElmtId());
1808
1809 // Get the normals of edge e
1810 const Array<OneD, const Array<OneD, NekDouble>> &normals =
1811 fields[0]->GetExp(n)->GetTraceNormal(e);
1812
1813 // Extract the edge values of flux-x on edge e and order
1814 // them accordingly to the order of the trace space
1815 fields[0]->GetExp(n)->GetTracePhysVals(e, elmtToTrace[n][e],
1816 fluxX1 + phys_offset,
1817 auxArray1 = tmparrayX1);
1818
1819 // Extract the edge values of flux-y on edge e and order
1820 // them accordingly to the order of the trace space
1821 fields[0]->GetExp(n)->GetTracePhysVals(e, elmtToTrace[n][e],
1822 fluxX2 + phys_offset,
1823 auxArray1 = tmparrayX2);
1824
1825 // Multiply the edge components of the flux by the normal
1826 for (i = 0; i < nEdgePts; ++i)
1827 {
1828 fluxN[i] = tmparrayX1[i] * m_traceNormals[0][trace_offset + i] +
1829 tmparrayX2[i] * m_traceNormals[1][trace_offset + i];
1830 }
1831
1832 // Subtract to the Riemann flux the discontinuous flux
1833 Vmath::Vsub(nEdgePts, &numericalFlux[trace_offset], 1, &fluxN[0], 1,
1834 &fluxJumps[0], 1);
1835
1836 // Check the ordering of the jump vectors
1837 if (fields[0]->GetExp(n)->GetTraceOrient(e) ==
1838 StdRegions::eBackwards)
1839 {
1840 Vmath::Reverse(nEdgePts, auxArray1 = fluxJumps, 1,
1841 auxArray2 = fluxJumps, 1);
1842 }
1843
1844 for (i = 0; i < nEdgePts; ++i)
1845 {
1847 m_traceNormals[1][trace_offset + i] != normals[1][i])
1848 {
1849 fluxJumps[i] = -fluxJumps[i];
1850 }
1851 }
1852
1853 // Multiply jumps by derivatives of the correction functions
1854 switch (e)
1855 {
1856 case 0:
1857 for (i = 0; i < nquad0; ++i)
1858 {
1859 // Multiply fluxJumps by Q factors
1860 fluxJumps[i] = -(m_Q2D_e0[n][i]) * fluxJumps[i];
1861
1862 for (j = 0; j < nquad1; ++j)
1863 {
1864 cnt = i + j * nquad0;
1865 divCFluxE0[cnt] = fluxJumps[i] * m_dGL_xi2[n][j];
1866 }
1867 }
1868 break;
1869 case 1:
1870 for (i = 0; i < nquad1; ++i)
1871 {
1872 // Multiply fluxJumps by Q factors
1873 fluxJumps[i] = (m_Q2D_e1[n][i]) * fluxJumps[i];
1874
1875 for (j = 0; j < nquad0; ++j)
1876 {
1877 cnt = (nquad0)*i + j;
1878 divCFluxE1[cnt] = fluxJumps[i] * m_dGR_xi1[n][j];
1879 }
1880 }
1881 break;
1882 case 2:
1883 for (i = 0; i < nquad0; ++i)
1884 {
1885 // Multiply fluxJumps by Q factors
1886 fluxJumps[i] = (m_Q2D_e2[n][i]) * fluxJumps[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 // Multiply fluxJumps by Q factors
1899 fluxJumps[i] = -(m_Q2D_e3[n][i]) * fluxJumps[i];
1900 for (j = 0; j < nquad0; ++j)
1901 {
1902 cnt = j + i * nquad0;
1903 divCFluxE3[cnt] = fluxJumps[i] * m_dGL_xi1[n][j];
1904 }
1905 }
1906 break;
1907
1908 default:
1910 break;
1911 }
1912 }
1913
1914 // Sum each edge contribution
1915 Vmath::Vadd(nLocalSolutionPts, &divCFluxE0[0], 1, &divCFluxE1[0], 1,
1916 &divCFlux[phys_offset], 1);
1917
1918 Vmath::Vadd(nLocalSolutionPts, &divCFlux[phys_offset], 1,
1919 &divCFluxE2[0], 1, &divCFlux[phys_offset], 1);
1920
1921 Vmath::Vadd(nLocalSolutionPts, &divCFlux[phys_offset], 1,
1922 &divCFluxE3[0], 1, &divCFlux[phys_offset], 1);
1923 }
1924}
Array< OneD, Array< OneD, NekDouble > > m_traceNormals
Definition: DiffusionLFR.h:69
Array< OneD, Array< OneD, NekDouble > > m_Q2D_e0
Definition: DiffusionLFR.h:75
Array< OneD, Array< OneD, NekDouble > > m_Q2D_e3
Definition: DiffusionLFR.h:78
Array< OneD, Array< OneD, NekDouble > > m_Q2D_e2
Definition: DiffusionLFR.h:77
Array< OneD, Array< OneD, NekDouble > > m_Q2D_e1
Definition: DiffusionLFR.h:76
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 1940 of file DiffusionLFR.cpp.
1947{
1948 int n, e, i, j, cnt;
1949
1950 int nElements = fields[0]->GetExpSize();
1951 int nLocalSolutionPts;
1952 int nEdgePts;
1953 int trace_offset;
1954 int phys_offset;
1955 int nquad0;
1956 int nquad1;
1957
1958 Array<OneD, NekDouble> auxArray1, auxArray2;
1959 Array<OneD, LibUtilities::BasisSharedPtr> base;
1960
1961 Array<OneD, Array<OneD, LocalRegions::ExpansionSharedPtr>> &elmtToTrace =
1962 fields[0]->GetTraceMap()->GetElmtToTrace();
1963
1964 // Loop on the elements
1965 for (n = 0; n < nElements; ++n)
1966 {
1967 // Offset of the element on the global vector
1968 phys_offset = fields[0]->GetPhys_Offset(n);
1969 nLocalSolutionPts = fields[0]->GetExp(n)->GetTotPoints();
1970
1971 base = fields[0]->GetExp(n)->GetBase();
1972 nquad0 = base[0]->GetNumPoints();
1973 nquad1 = base[1]->GetNumPoints();
1974
1975 Array<OneD, NekDouble> divCFluxE0(nLocalSolutionPts, 0.0);
1976 Array<OneD, NekDouble> divCFluxE1(nLocalSolutionPts, 0.0);
1977 Array<OneD, NekDouble> divCFluxE2(nLocalSolutionPts, 0.0);
1978 Array<OneD, NekDouble> divCFluxE3(nLocalSolutionPts, 0.0);
1979
1980 // Loop on the edges
1981 for (e = 0; e < fields[0]->GetExp(n)->GetNtraces(); ++e)
1982 {
1983 // Number of edge points of edge e
1984 nEdgePts = fields[0]->GetExp(n)->GetTraceNumPoints(e);
1985
1986 Array<OneD, NekDouble> fluxN(nEdgePts, 0.0);
1987 Array<OneD, NekDouble> fluxT(nEdgePts, 0.0);
1988 Array<OneD, NekDouble> fluxN_R(nEdgePts, 0.0);
1989 Array<OneD, NekDouble> fluxN_D(nEdgePts, 0.0);
1990 Array<OneD, NekDouble> fluxJumps(nEdgePts, 0.0);
1991
1992 // Offset of the trace space correspondent to edge e
1993 trace_offset = fields[0]->GetTrace()->GetPhys_Offset(
1994 elmtToTrace[n][e]->GetElmtId());
1995
1996 // Get the normals of edge e
1997 const Array<OneD, const Array<OneD, NekDouble>> &normals =
1998 fields[0]->GetExp(n)->GetTraceNormal(e);
1999
2000 // Extract the trasformed normal flux at each edge
2001 switch (e)
2002 {
2003 case 0:
2004 // Extract the edge values of transformed flux-y on
2005 // edge e and order them accordingly to the order of
2006 // the trace space
2007 fields[0]->GetExp(n)->GetTracePhysVals(e, elmtToTrace[n][e],
2008 fluxX2 + phys_offset,
2009 auxArray1 = fluxN_D);
2010
2011 Vmath::Neg(nEdgePts, fluxN_D, 1);
2012
2013 // Extract the physical Rieamann flux at each edge
2014 Vmath::Vcopy(nEdgePts, &numericalFlux[trace_offset], 1,
2015 &fluxN[0], 1);
2016
2017 // Check the ordering of vectors
2018 if (fields[0]->GetExp(n)->GetTraceOrient(e) ==
2019 StdRegions::eBackwards)
2020 {
2021 Vmath::Reverse(nEdgePts, auxArray1 = fluxN, 1,
2022 auxArray2 = fluxN, 1);
2023
2024 Vmath::Reverse(nEdgePts, auxArray1 = fluxN_D, 1,
2025 auxArray2 = fluxN_D, 1);
2026 }
2027
2028 // Transform Riemann Fluxes in the standard element
2029 for (i = 0; i < nquad0; ++i)
2030 {
2031 // Multiply Riemann Flux by Q factors
2032 fluxN_R[i] = (m_Q2D_e0[n][i]) * fluxN[i];
2033 }
2034
2035 for (i = 0; i < nEdgePts; ++i)
2036 {
2038 normals[0][i] ||
2039 m_traceNormals[1][trace_offset + i] !=
2040 normals[1][i])
2041 {
2042 fluxN_R[i] = -fluxN_R[i];
2043 }
2044 }
2045
2046 // Subtract to the Riemann flux the discontinuous
2047 // flux
2048 Vmath::Vsub(nEdgePts, &fluxN_R[0], 1, &fluxN_D[0], 1,
2049 &fluxJumps[0], 1);
2050
2051 // Multiplicate the flux jumps for the correction
2052 // function
2053 for (i = 0; i < nquad0; ++i)
2054 {
2055 for (j = 0; j < nquad1; ++j)
2056 {
2057 cnt = i + j * nquad0;
2058 divCFluxE0[cnt] = -fluxJumps[i] * m_dGL_xi2[n][j];
2059 }
2060 }
2061 break;
2062 case 1:
2063 // Extract the edge values of transformed flux-x on
2064 // edge e and order them accordingly to the order of
2065 // the trace space
2066 fields[0]->GetExp(n)->GetTracePhysVals(e, elmtToTrace[n][e],
2067 fluxX1 + phys_offset,
2068 auxArray1 = fluxN_D);
2069
2070 // Extract the physical Rieamann flux at each edge
2071 Vmath::Vcopy(nEdgePts, &numericalFlux[trace_offset], 1,
2072 &fluxN[0], 1);
2073
2074 // Check the ordering of vectors
2075 if (fields[0]->GetExp(n)->GetTraceOrient(e) ==
2076 StdRegions::eBackwards)
2077 {
2078 Vmath::Reverse(nEdgePts, auxArray1 = fluxN, 1,
2079 auxArray2 = fluxN, 1);
2080
2081 Vmath::Reverse(nEdgePts, auxArray1 = fluxN_D, 1,
2082 auxArray2 = fluxN_D, 1);
2083 }
2084
2085 // Transform Riemann Fluxes in the standard element
2086 for (i = 0; i < nquad1; ++i)
2087 {
2088 // Multiply Riemann Flux by Q factors
2089 fluxN_R[i] = (m_Q2D_e1[n][i]) * fluxN[i];
2090 }
2091
2092 for (i = 0; i < nEdgePts; ++i)
2093 {
2095 normals[0][i] ||
2096 m_traceNormals[1][trace_offset + i] !=
2097 normals[1][i])
2098 {
2099 fluxN_R[i] = -fluxN_R[i];
2100 }
2101 }
2102
2103 // Subtract to the Riemann flux the discontinuous
2104 // flux
2105 Vmath::Vsub(nEdgePts, &fluxN_R[0], 1, &fluxN_D[0], 1,
2106 &fluxJumps[0], 1);
2107
2108 // Multiplicate the flux jumps for the correction
2109 // function
2110 for (i = 0; i < nquad1; ++i)
2111 {
2112 for (j = 0; j < nquad0; ++j)
2113 {
2114 cnt = (nquad0)*i + j;
2115 divCFluxE1[cnt] = fluxJumps[i] * m_dGR_xi1[n][j];
2116 }
2117 }
2118 break;
2119 case 2:
2120
2121 // Extract the edge values of transformed flux-y on
2122 // edge e and order them accordingly to the order of
2123 // the trace space
2124
2125 fields[0]->GetExp(n)->GetTracePhysVals(e, elmtToTrace[n][e],
2126 fluxX2 + phys_offset,
2127 auxArray1 = fluxN_D);
2128
2129 // Extract the physical Rieamann flux at each edge
2130 Vmath::Vcopy(nEdgePts, &numericalFlux[trace_offset], 1,
2131 &fluxN[0], 1);
2132
2133 // Check the ordering of vectors
2134 if (fields[0]->GetExp(n)->GetTraceOrient(e) ==
2135 StdRegions::eBackwards)
2136 {
2137 Vmath::Reverse(nEdgePts, auxArray1 = fluxN, 1,
2138 auxArray2 = fluxN, 1);
2139
2140 Vmath::Reverse(nEdgePts, auxArray1 = fluxN_D, 1,
2141 auxArray2 = fluxN_D, 1);
2142 }
2143
2144 // Transform Riemann Fluxes in the standard element
2145 for (i = 0; i < nquad0; ++i)
2146 {
2147 // Multiply Riemann Flux by Q factors
2148 fluxN_R[i] = (m_Q2D_e2[n][i]) * fluxN[i];
2149 }
2150
2151 for (i = 0; i < nEdgePts; ++i)
2152 {
2154 normals[0][i] ||
2155 m_traceNormals[1][trace_offset + i] !=
2156 normals[1][i])
2157 {
2158 fluxN_R[i] = -fluxN_R[i];
2159 }
2160 }
2161
2162 // Subtract to the Riemann flux the discontinuous
2163 // flux
2164
2165 Vmath::Vsub(nEdgePts, &fluxN_R[0], 1, &fluxN_D[0], 1,
2166 &fluxJumps[0], 1);
2167
2168 // Multiplicate the flux jumps for the correction
2169 // function
2170 for (i = 0; i < nquad0; ++i)
2171 {
2172 for (j = 0; j < nquad1; ++j)
2173 {
2174 cnt = j * nquad0 + i;
2175 divCFluxE2[cnt] = fluxJumps[i] * m_dGR_xi2[n][j];
2176 }
2177 }
2178 break;
2179 case 3:
2180 // Extract the edge values of transformed flux-x on
2181 // edge e and order them accordingly to the order of
2182 // the trace space
2183
2184 fields[0]->GetExp(n)->GetTracePhysVals(e, elmtToTrace[n][e],
2185 fluxX1 + phys_offset,
2186 auxArray1 = fluxN_D);
2187 Vmath::Neg(nEdgePts, fluxN_D, 1);
2188
2189 // Extract the physical Rieamann flux at each edge
2190 Vmath::Vcopy(nEdgePts, &numericalFlux[trace_offset], 1,
2191 &fluxN[0], 1);
2192
2193 // Check the ordering of vectors
2194 if (fields[0]->GetExp(n)->GetTraceOrient(e) ==
2195 StdRegions::eBackwards)
2196 {
2197 Vmath::Reverse(nEdgePts, auxArray1 = fluxN, 1,
2198 auxArray2 = fluxN, 1);
2199
2200 Vmath::Reverse(nEdgePts, auxArray1 = fluxN_D, 1,
2201 auxArray2 = fluxN_D, 1);
2202 }
2203
2204 // Transform Riemann Fluxes in the standard element
2205 for (i = 0; i < nquad1; ++i)
2206 {
2207 // Multiply Riemann Flux by Q factors
2208 fluxN_R[i] = (m_Q2D_e3[n][i]) * fluxN[i];
2209 }
2210
2211 for (i = 0; i < nEdgePts; ++i)
2212 {
2214 normals[0][i] ||
2215 m_traceNormals[1][trace_offset + i] !=
2216 normals[1][i])
2217 {
2218 fluxN_R[i] = -fluxN_R[i];
2219 }
2220 }
2221
2222 // Subtract to the Riemann flux the discontinuous
2223 // flux
2224
2225 Vmath::Vsub(nEdgePts, &fluxN_R[0], 1, &fluxN_D[0], 1,
2226 &fluxJumps[0], 1);
2227
2228 // Multiplicate the flux jumps for the correction
2229 // function
2230 for (i = 0; i < nquad1; ++i)
2231 {
2232 for (j = 0; j < nquad0; ++j)
2233 {
2234 cnt = j + i * nquad0;
2235 divCFluxE3[cnt] = -fluxJumps[i] * m_dGL_xi1[n][j];
2236 }
2237 }
2238 break;
2239 default:
2241 break;
2242 }
2243 }
2244
2245 // Sum each edge contribution
2246 Vmath::Vadd(nLocalSolutionPts, &divCFluxE0[0], 1, &divCFluxE1[0], 1,
2247 &divCFlux[phys_offset], 1);
2248
2249 Vmath::Vadd(nLocalSolutionPts, &divCFlux[phys_offset], 1,
2250 &divCFluxE2[0], 1, &divCFlux[phys_offset], 1);
2251
2252 Vmath::Vadd(nLocalSolutionPts, &divCFlux[phys_offset], 1,
2253 &divCFluxE3[0], 1, &divCFlux[phys_offset], 1);
2254 }
2255}
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 1161 of file DiffusionLFR.cpp.
1165{
1166 int i, j;
1167 int nTracePts = fields[0]->GetTrace()->GetTotPoints();
1168 int nvariables = fields.size();
1169 int nDim = fields[0]->GetCoordim(0);
1170
1171 Array<OneD, NekDouble> Fwd(nTracePts);
1172 Array<OneD, NekDouble> Bwd(nTracePts);
1173 Array<OneD, NekDouble> Vn(nTracePts, 0.0);
1174 Array<OneD, NekDouble> fluxtemp(nTracePts, 0.0);
1175
1176 // Get the normal velocity Vn
1177 for (i = 0; i < nDim; ++i)
1178 {
1180 }
1181
1182 // Get the sign of (v \cdot n), v = an arbitrary vector
1183 // Evaluate upwind flux:
1184 // uflux = \hat{u} \phi \cdot u = u^{(+,-)} n
1185 for (j = 0; j < nDim; ++j)
1186 {
1187 for (i = 0; i < nvariables; ++i)
1188 {
1189 // Compute Fwd and Bwd value of ufield of i direction
1190 fields[i]->GetFwdBwdTracePhys(ufield[i], Fwd, Bwd);
1191
1192 // if Vn >= 0, flux = uFwd, i.e.,
1193 // edge::eForward, if V*n>=0 <=> V*n_F>=0, pick uflux = uFwd
1194 // edge::eBackward, if V*n>=0 <=> V*n_B<0, pick uflux = uFwd
1195
1196 // else if Vn < 0, flux = uBwd, i.e.,
1197 // edge::eForward, if V*n<0 <=> V*n_F<0, pick uflux = uBwd
1198 // edge::eBackward, if V*n<0 <=> V*n_B>=0, pick uflux = uBwd
1199
1200 fields[i]->GetTrace()->Upwind(Vn, Fwd, Bwd, fluxtemp);
1201
1202 // Imposing weak boundary condition with flux
1203 // if Vn >= 0, uflux = uBwd at Neumann, i.e.,
1204 // edge::eForward, if V*n>=0 <=> V*n_F>=0, pick uflux = uBwd
1205 // edge::eBackward, if V*n>=0 <=> V*n_B<0, pick uflux = uBwd
1206
1207 // if Vn >= 0, uflux = uFwd at Neumann, i.e.,
1208 // edge::eForward, if V*n<0 <=> V*n_F<0, pick uflux = uFwd
1209 // edge::eBackward, if V*n<0 <=> V*n_B>=0, pick uflux = uFwd
1210
1211 if (fields[0]->GetBndCondExpansions().size())
1212 {
1213 WeakPenaltyforScalar(fields, i, ufield[i], fluxtemp);
1214 }
1215
1216 // if Vn >= 0, flux = uFwd*(tan_{\xi}^- \cdot \vec{n}),
1217 // i.e,
1218 // edge::eForward, uFwd \(\tan_{\xi}^Fwd \cdot \vec{n})
1219 // edge::eBackward, uFwd \(\tan_{\xi}^Bwd \cdot \vec{n})
1220
1221 // else if Vn < 0, flux = uBwd*(tan_{\xi}^- \cdot \vec{n}),
1222 // i.e,
1223 // edge::eForward, uBwd \(\tan_{\xi}^Fwd \cdot \vec{n})
1224 // edge::eBackward, uBwd \(\tan_{\xi}^Bwd \cdot \vec{n})
1225
1226 Vmath::Vcopy(nTracePts, fluxtemp, 1, uflux[i][j], 1);
1227 }
1228 }
1229}
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:1235
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.hpp:396
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 1298 of file DiffusionLFR.cpp.
1303{
1304 int i, j;
1305 int nTracePts = fields[0]->GetTrace()->GetTotPoints();
1306 int nvariables = fields.size();
1307 int nDim = fields[0]->GetCoordim(0);
1308
1309 NekDouble C11 = 0.0;
1310 Array<OneD, NekDouble> Fwd(nTracePts);
1311 Array<OneD, NekDouble> Bwd(nTracePts);
1312 Array<OneD, NekDouble> Vn(nTracePts, 0.0);
1313
1314 Array<OneD, NekDouble> qFwd(nTracePts);
1315 Array<OneD, NekDouble> qBwd(nTracePts);
1316 Array<OneD, NekDouble> qfluxtemp(nTracePts, 0.0);
1317 Array<OneD, NekDouble> uterm(nTracePts);
1318
1319 // Get the normal velocity Vn
1320 for (i = 0; i < nDim; ++i)
1321 {
1323 }
1324
1325 // Evaulate upwind flux:
1326 // qflux = \hat{q} \cdot u = q \cdot n - C_(11)*(u^+ - u^-)
1327 for (i = 0; i < nvariables; ++i)
1328 {
1329 qflux[i] = Array<OneD, NekDouble>(nTracePts, 0.0);
1330 for (j = 0; j < nDim; ++j)
1331 {
1332 // Compute Fwd and Bwd value of ufield of jth direction
1333 fields[i]->GetFwdBwdTracePhys(qfield[i][j], qFwd, qBwd);
1334
1335 // if Vn >= 0, flux = uFwd, i.e.,
1336 // edge::eForward, if V*n>=0 <=> V*n_F>=0, pick
1337 // qflux = qBwd = q+
1338 // edge::eBackward, if V*n>=0 <=> V*n_B<0, pick
1339 // qflux = qBwd = q-
1340
1341 // else if Vn < 0, flux = uBwd, i.e.,
1342 // edge::eForward, if V*n<0 <=> V*n_F<0, pick
1343 // qflux = qFwd = q-
1344 // edge::eBackward, if V*n<0 <=> V*n_B>=0, pick
1345 // qflux = qFwd = q+
1346
1347 fields[i]->GetTrace()->Upwind(Vn, qBwd, qFwd, qfluxtemp);
1348
1350 qfluxtemp, 1);
1351
1352 // Imposing weak boundary condition with flux
1353 if (fields[0]->GetBndCondExpansions().size())
1354 {
1355 WeakPenaltyforVector(fields, i, j, qfield[i][j], qfluxtemp,
1356 C11);
1357 }
1358
1359 // q_hat \cdot n = (q_xi \cdot n_xi) or (q_eta \cdot n_eta)
1360 // n_xi = n_x * tan_xi_x + n_y * tan_xi_y + n_z * tan_xi_z
1361 // n_xi = n_x * tan_eta_x + n_y * tan_eta_y + n_z*tan_eta_z
1362 Vmath::Vadd(nTracePts, qfluxtemp, 1, qflux[i], 1, qflux[i], 1);
1363 }
1364 }
1365}
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:1371
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, 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 324 of file DiffusionLFR.cpp.
327{
328 int i, n;
329 NekDouble c0 = 0.0;
330 NekDouble c1 = 0.0;
331 NekDouble c2 = 0.0;
332 int nquad0, nquad1, nquad2;
333 int nmodes0, nmodes1, nmodes2;
334 Array<OneD, LibUtilities::BasisSharedPtr> base;
335
336 int nElements = pFields[0]->GetExpSize();
337 int nDim = pFields[0]->GetCoordim(0);
338
339 switch (nDim)
340 {
341 case 1:
342 {
343 m_dGL_xi1 = Array<OneD, Array<OneD, NekDouble>>(nElements);
344 m_dGR_xi1 = Array<OneD, Array<OneD, NekDouble>>(nElements);
345
346 for (n = 0; n < nElements; ++n)
347 {
348 base = pFields[0]->GetExp(n)->GetBase();
349 nquad0 = base[0]->GetNumPoints();
350 nmodes0 = base[0]->GetNumModes();
351 Array<OneD, const NekDouble> z0;
352 Array<OneD, const NekDouble> w0;
353
354 base[0]->GetZW(z0, w0);
355
356 m_dGL_xi1[n] = Array<OneD, NekDouble>(nquad0);
357 m_dGR_xi1[n] = Array<OneD, NekDouble>(nquad0);
358
359 // Auxiliary vectors to build up the auxiliary
360 // derivatives of the Legendre polynomials
361 Array<OneD, NekDouble> dLp0(nquad0, 0.0);
362 Array<OneD, NekDouble> dLpp0(nquad0, 0.0);
363 Array<OneD, NekDouble> dLpm0(nquad0, 0.0);
364
365 // Degree of the correction functions
366 int p0 = nmodes0 - 1;
367
368 // Function sign to compute left correction function
369 NekDouble sign0 = pow(-1.0, p0);
370
371 // Factorial factor to build the scheme
372 NekDouble ap0 =
373 std::tgamma(2 * p0 + 1) /
374 (pow(2.0, p0) * std::tgamma(p0 + 1) * std::tgamma(p0 + 1));
375
376 // Scalar parameter which recovers the FR schemes
378 {
379 c0 = 0.0;
380 }
382 {
383 c0 = 2.0 * p0 /
384 ((2.0 * p0 + 1.0) * (p0 + 1.0) *
385 (ap0 * std::tgamma(p0 + 1)) *
386 (ap0 * std::tgamma(p0 + 1)));
387 }
389 {
390 c0 = 2.0 * (p0 + 1.0) /
391 ((2.0 * p0 + 1.0) * p0 * (ap0 * std::tgamma(p0 + 1)) *
392 (ap0 * std::tgamma(p0 + 1)));
393 }
395 {
396 c0 =
397 -2.0 / ((2.0 * p0 + 1.0) * (ap0 * std::tgamma(p0 + 1)) *
398 (ap0 * std::tgamma(p0 + 1)));
399 }
401 {
402 c0 = 10000000000000000.0;
403 }
404
405 NekDouble etap0 = 0.5 * c0 * (2.0 * p0 + 1.0) *
406 (ap0 * std::tgamma(p0 + 1)) *
407 (ap0 * std::tgamma(p0 + 1));
408
409 NekDouble overeta0 = 1.0 / (1.0 + etap0);
410
411 // Derivative of the Legendre polynomials
412 // dLp = derivative of the Legendre polynomial of order p
413 // dLpp = derivative of the Legendre polynomial of order p+1
414 // dLpm = derivative of the Legendre polynomial of order p-1
415 Polylib::jacobd(nquad0, z0.data(), &(dLp0[0]), p0, 0.0, 0.0);
416 Polylib::jacobd(nquad0, z0.data(), &(dLpp0[0]), p0 + 1, 0.0,
417 0.0);
418 Polylib::jacobd(nquad0, z0.data(), &(dLpm0[0]), p0 - 1, 0.0,
419 0.0);
420
421 // Building the DG_c_Left
422 for (i = 0; i < nquad0; ++i)
423 {
424 m_dGL_xi1[n][i] = etap0 * dLpm0[i];
425 m_dGL_xi1[n][i] += dLpp0[i];
426 m_dGL_xi1[n][i] *= overeta0;
429 }
430
431 // Building the DG_c_Right
432 for (i = 0; i < nquad0; ++i)
433 {
434 m_dGR_xi1[n][i] = etap0 * dLpm0[i];
435 m_dGR_xi1[n][i] += dLpp0[i];
436 m_dGR_xi1[n][i] *= overeta0;
437 m_dGR_xi1[n][i] += dLp0[i];
439 }
440 }
441 break;
442 }
443 case 2:
444 {
445 LibUtilities::BasisSharedPtr BasisFR_Left0;
446 LibUtilities::BasisSharedPtr BasisFR_Right0;
447 LibUtilities::BasisSharedPtr BasisFR_Left1;
448 LibUtilities::BasisSharedPtr BasisFR_Right1;
449
450 m_dGL_xi1 = Array<OneD, Array<OneD, NekDouble>>(nElements);
451 m_dGR_xi1 = Array<OneD, Array<OneD, NekDouble>>(nElements);
452 m_dGL_xi2 = Array<OneD, Array<OneD, NekDouble>>(nElements);
453 m_dGR_xi2 = Array<OneD, Array<OneD, NekDouble>>(nElements);
454
455 for (n = 0; n < nElements; ++n)
456 {
457 base = pFields[0]->GetExp(n)->GetBase();
458 nquad0 = base[0]->GetNumPoints();
459 nquad1 = base[1]->GetNumPoints();
460 nmodes0 = base[0]->GetNumModes();
461 nmodes1 = base[1]->GetNumModes();
462
463 Array<OneD, const NekDouble> z0;
464 Array<OneD, const NekDouble> w0;
465 Array<OneD, const NekDouble> z1;
466 Array<OneD, const NekDouble> w1;
467
468 base[0]->GetZW(z0, w0);
469 base[1]->GetZW(z1, w1);
470
471 m_dGL_xi1[n] = Array<OneD, NekDouble>(nquad0);
472 m_dGR_xi1[n] = Array<OneD, NekDouble>(nquad0);
473 m_dGL_xi2[n] = Array<OneD, NekDouble>(nquad1);
474 m_dGR_xi2[n] = Array<OneD, NekDouble>(nquad1);
475
476 // Auxiliary vectors to build up the auxiliary
477 // derivatives of the Legendre polynomials
478 Array<OneD, NekDouble> dLp0(nquad0, 0.0);
479 Array<OneD, NekDouble> dLpp0(nquad0, 0.0);
480 Array<OneD, NekDouble> dLpm0(nquad0, 0.0);
481 Array<OneD, NekDouble> dLp1(nquad1, 0.0);
482 Array<OneD, NekDouble> dLpp1(nquad1, 0.0);
483 Array<OneD, NekDouble> dLpm1(nquad1, 0.0);
484
485 // Degree of the correction functions
486 int p0 = nmodes0 - 1;
487 int p1 = nmodes1 - 1;
488
489 // Function sign to compute left correction function
490 NekDouble sign0 = pow(-1.0, p0);
491 NekDouble sign1 = pow(-1.0, p1);
492
493 // Factorial factor to build the scheme
494 NekDouble ap0 =
495 std::tgamma(2 * p0 + 1) /
496 (pow(2.0, p0) * std::tgamma(p0 + 1) * std::tgamma(p0 + 1));
497
498 NekDouble ap1 =
499 std::tgamma(2 * p1 + 1) /
500 (pow(2.0, p1) * std::tgamma(p1 + 1) * std::tgamma(p1 + 1));
501
502 // Scalar parameter which recovers the FR schemes
504 {
505 c0 = 0.0;
506 c1 = 0.0;
507 }
509 {
510 c0 = 2.0 * p0 /
511 ((2.0 * p0 + 1.0) * (p0 + 1.0) *
512 (ap0 * std::tgamma(p0 + 1)) *
513 (ap0 * std::tgamma(p0 + 1)));
514
515 c1 = 2.0 * p1 /
516 ((2.0 * p1 + 1.0) * (p1 + 1.0) *
517 (ap1 * std::tgamma(p1 + 1)) *
518 (ap1 * std::tgamma(p1 + 1)));
519 }
521 {
522 c0 = 2.0 * (p0 + 1.0) /
523 ((2.0 * p0 + 1.0) * p0 * (ap0 * std::tgamma(p0 + 1)) *
524 (ap0 * std::tgamma(p0 + 1)));
525
526 c1 = 2.0 * (p1 + 1.0) /
527 ((2.0 * p1 + 1.0) * p1 * (ap1 * std::tgamma(p1 + 1)) *
528 (ap1 * std::tgamma(p1 + 1)));
529 }
531 {
532 c0 =
533 -2.0 / ((2.0 * p0 + 1.0) * (ap0 * std::tgamma(p0 + 1)) *
534 (ap0 * std::tgamma(p0 + 1)));
535
536 c1 =
537 -2.0 / ((2.0 * p1 + 1.0) * (ap1 * std::tgamma(p1 + 1)) *
538 (ap1 * std::tgamma(p1 + 1)));
539 }
541 {
542 c0 = 10000000000000000.0;
543 c1 = 10000000000000000.0;
544 }
545
546 NekDouble etap0 = 0.5 * c0 * (2.0 * p0 + 1.0) *
547 (ap0 * std::tgamma(p0 + 1)) *
548 (ap0 * std::tgamma(p0 + 1));
549
550 NekDouble etap1 = 0.5 * c1 * (2.0 * p1 + 1.0) *
551 (ap1 * std::tgamma(p1 + 1)) *
552 (ap1 * std::tgamma(p1 + 1));
553
554 NekDouble overeta0 = 1.0 / (1.0 + etap0);
555 NekDouble overeta1 = 1.0 / (1.0 + etap1);
556
557 // Derivative of the Legendre polynomials
558 // dLp = derivative of the Legendre polynomial of order p
559 // dLpp = derivative of the Legendre polynomial of order p+1
560 // dLpm = derivative of the Legendre polynomial of order p-1
561 Polylib::jacobd(nquad0, z0.data(), &(dLp0[0]), p0, 0.0, 0.0);
562 Polylib::jacobd(nquad0, z0.data(), &(dLpp0[0]), p0 + 1, 0.0,
563 0.0);
564 Polylib::jacobd(nquad0, z0.data(), &(dLpm0[0]), p0 - 1, 0.0,
565 0.0);
566
567 Polylib::jacobd(nquad1, z1.data(), &(dLp1[0]), p1, 0.0, 0.0);
568 Polylib::jacobd(nquad1, z1.data(), &(dLpp1[0]), p1 + 1, 0.0,
569 0.0);
570 Polylib::jacobd(nquad1, z1.data(), &(dLpm1[0]), p1 - 1, 0.0,
571 0.0);
572
573 // Building the DG_c_Left
574 for (i = 0; i < nquad0; ++i)
575 {
576 m_dGL_xi1[n][i] = etap0 * dLpm0[i];
577 m_dGL_xi1[n][i] += dLpp0[i];
578 m_dGL_xi1[n][i] *= overeta0;
581 }
582
583 // Building the DG_c_Left
584 for (i = 0; i < nquad1; ++i)
585 {
586 m_dGL_xi2[n][i] = etap1 * dLpm1[i];
587 m_dGL_xi2[n][i] += dLpp1[i];
588 m_dGL_xi2[n][i] *= overeta1;
591 }
592
593 // Building the DG_c_Right
594 for (i = 0; i < nquad0; ++i)
595 {
596 m_dGR_xi1[n][i] = etap0 * dLpm0[i];
597 m_dGR_xi1[n][i] += dLpp0[i];
598 m_dGR_xi1[n][i] *= overeta0;
599 m_dGR_xi1[n][i] += dLp0[i];
601 }
602
603 // Building the DG_c_Right
604 for (i = 0; i < nquad1; ++i)
605 {
606 m_dGR_xi2[n][i] = etap1 * dLpm1[i];
607 m_dGR_xi2[n][i] += dLpp1[i];
608 m_dGR_xi2[n][i] *= overeta1;
609 m_dGR_xi2[n][i] += dLp1[i];
611 }
612 }
613 break;
614 }
615 case 3:
616 {
617 m_dGL_xi1 = Array<OneD, Array<OneD, NekDouble>>(nElements);
618 m_dGR_xi1 = Array<OneD, Array<OneD, NekDouble>>(nElements);
619 m_dGL_xi2 = Array<OneD, Array<OneD, NekDouble>>(nElements);
620 m_dGR_xi2 = Array<OneD, Array<OneD, NekDouble>>(nElements);
621 m_dGL_xi3 = Array<OneD, Array<OneD, NekDouble>>(nElements);
622 m_dGR_xi3 = Array<OneD, Array<OneD, NekDouble>>(nElements);
623
624 for (n = 0; n < nElements; ++n)
625 {
626 base = pFields[0]->GetExp(n)->GetBase();
627 nquad0 = base[0]->GetNumPoints();
628 nquad1 = base[1]->GetNumPoints();
629 nquad2 = base[2]->GetNumPoints();
630 nmodes0 = base[0]->GetNumModes();
631 nmodes1 = base[1]->GetNumModes();
632 nmodes2 = base[2]->GetNumModes();
633
634 Array<OneD, const NekDouble> z0;
635 Array<OneD, const NekDouble> w0;
636 Array<OneD, const NekDouble> z1;
637 Array<OneD, const NekDouble> w1;
638 Array<OneD, const NekDouble> z2;
639 Array<OneD, const NekDouble> w2;
640
641 base[0]->GetZW(z0, w0);
642 base[1]->GetZW(z1, w1);
643 base[1]->GetZW(z2, w2);
644
645 m_dGL_xi1[n] = Array<OneD, NekDouble>(nquad0);
646 m_dGR_xi1[n] = Array<OneD, NekDouble>(nquad0);
647 m_dGL_xi2[n] = Array<OneD, NekDouble>(nquad1);
648 m_dGR_xi2[n] = Array<OneD, NekDouble>(nquad1);
649 m_dGL_xi3[n] = Array<OneD, NekDouble>(nquad2);
650 m_dGR_xi3[n] = Array<OneD, NekDouble>(nquad2);
651
652 // Auxiliary vectors to build up the auxiliary
653 // derivatives of the Legendre polynomials
654 Array<OneD, NekDouble> dLp0(nquad0, 0.0);
655 Array<OneD, NekDouble> dLpp0(nquad0, 0.0);
656 Array<OneD, NekDouble> dLpm0(nquad0, 0.0);
657 Array<OneD, NekDouble> dLp1(nquad1, 0.0);
658 Array<OneD, NekDouble> dLpp1(nquad1, 0.0);
659 Array<OneD, NekDouble> dLpm1(nquad1, 0.0);
660 Array<OneD, NekDouble> dLp2(nquad2, 0.0);
661 Array<OneD, NekDouble> dLpp2(nquad2, 0.0);
662 Array<OneD, NekDouble> dLpm2(nquad2, 0.0);
663
664 // Degree of the correction functions
665 int p0 = nmodes0 - 1;
666 int p1 = nmodes1 - 1;
667 int p2 = nmodes2 - 1;
668
669 // Function sign to compute left correction function
670 NekDouble sign0 = pow(-1.0, p0);
671 NekDouble sign1 = pow(-1.0, p1);
672 NekDouble sign2 = pow(-1.0, p2);
673
674 // Factorial factor to build the scheme
675 NekDouble ap0 =
676 std::tgamma(2 * p0 + 1) /
677 (pow(2.0, p0) * std::tgamma(p0 + 1) * std::tgamma(p0 + 1));
678
679 // Factorial factor to build the scheme
680 NekDouble ap1 =
681 std::tgamma(2 * p1 + 1) /
682 (pow(2.0, p1) * std::tgamma(p1 + 1) * std::tgamma(p1 + 1));
683
684 // Factorial factor to build the scheme
685 NekDouble ap2 =
686 std::tgamma(2 * p2 + 1) /
687 (pow(2.0, p2) * std::tgamma(p2 + 1) * std::tgamma(p2 + 1));
688
689 // Scalar parameter which recovers the FR schemes
691 {
692 c0 = 0.0;
693 c1 = 0.0;
694 c2 = 0.0;
695 }
697 {
698 c0 = 2.0 * p0 /
699 ((2.0 * p0 + 1.0) * (p0 + 1.0) *
700 (ap0 * std::tgamma(p0 + 1)) *
701 (ap0 * std::tgamma(p0 + 1)));
702
703 c1 = 2.0 * p1 /
704 ((2.0 * p1 + 1.0) * (p1 + 1.0) *
705 (ap1 * std::tgamma(p1 + 1)) *
706 (ap1 * std::tgamma(p1 + 1)));
707
708 c2 = 2.0 * p2 /
709 ((2.0 * p2 + 1.0) * (p2 + 1.0) *
710 (ap2 * std::tgamma(p2 + 1)) *
711 (ap2 * std::tgamma(p2 + 1)));
712 }
714 {
715 c0 = 2.0 * (p0 + 1.0) /
716 ((2.0 * p0 + 1.0) * p0 * (ap0 * std::tgamma(p0 + 1)) *
717 (ap0 * std::tgamma(p0 + 1)));
718
719 c1 = 2.0 * (p1 + 1.0) /
720 ((2.0 * p1 + 1.0) * p1 * (ap1 * std::tgamma(p1 + 1)) *
721 (ap1 * std::tgamma(p1 + 1)));
722
723 c2 = 2.0 * (p2 + 1.0) /
724 ((2.0 * p2 + 1.0) * p2 * (ap2 * std::tgamma(p2 + 1)) *
725 (ap2 * std::tgamma(p2 + 1)));
726 }
728 {
729 c0 =
730 -2.0 / ((2.0 * p0 + 1.0) * (ap0 * std::tgamma(p0 + 1)) *
731 (ap0 * std::tgamma(p0 + 1)));
732
733 c1 =
734 -2.0 / ((2.0 * p1 + 1.0) * (ap1 * std::tgamma(p1 + 1)) *
735 (ap1 * std::tgamma(p1 + 1)));
736
737 c2 =
738 -2.0 / ((2.0 * p2 + 1.0) * (ap2 * std::tgamma(p2 + 1)) *
739 (ap2 * std::tgamma(p2 + 1)));
740 }
742 {
743 c0 = 10000000000000000.0;
744 c1 = 10000000000000000.0;
745 c2 = 10000000000000000.0;
746 }
747
748 NekDouble etap0 = 0.5 * c0 * (2.0 * p0 + 1.0) *
749 (ap0 * std::tgamma(p0 + 1)) *
750 (ap0 * std::tgamma(p0 + 1));
751
752 NekDouble etap1 = 0.5 * c1 * (2.0 * p1 + 1.0) *
753 (ap1 * std::tgamma(p1 + 1)) *
754 (ap1 * std::tgamma(p1 + 1));
755
756 NekDouble etap2 = 0.5 * c2 * (2.0 * p2 + 1.0) *
757 (ap2 * std::tgamma(p2 + 1)) *
758 (ap2 * std::tgamma(p2 + 1));
759
760 NekDouble overeta0 = 1.0 / (1.0 + etap0);
761 NekDouble overeta1 = 1.0 / (1.0 + etap1);
762 NekDouble overeta2 = 1.0 / (1.0 + etap2);
763
764 // Derivative of the Legendre polynomials
765 // dLp = derivative of the Legendre polynomial of order p
766 // dLpp = derivative of the Legendre polynomial of order p+1
767 // dLpm = derivative of the Legendre polynomial of order p-1
768 Polylib::jacobd(nquad0, z0.data(), &(dLp0[0]), p0, 0.0, 0.0);
769 Polylib::jacobd(nquad0, z0.data(), &(dLpp0[0]), p0 + 1, 0.0,
770 0.0);
771 Polylib::jacobd(nquad0, z0.data(), &(dLpm0[0]), p0 - 1, 0.0,
772 0.0);
773
774 Polylib::jacobd(nquad1, z1.data(), &(dLp1[0]), p1, 0.0, 0.0);
775 Polylib::jacobd(nquad1, z1.data(), &(dLpp1[0]), p1 + 1, 0.0,
776 0.0);
777 Polylib::jacobd(nquad1, z1.data(), &(dLpm1[0]), p1 - 1, 0.0,
778 0.0);
779
780 Polylib::jacobd(nquad2, z2.data(), &(dLp2[0]), p2, 0.0, 0.0);
781 Polylib::jacobd(nquad2, z2.data(), &(dLpp2[0]), p2 + 1, 0.0,
782 0.0);
783 Polylib::jacobd(nquad2, z2.data(), &(dLpm2[0]), p2 - 1, 0.0,
784 0.0);
785
786 // Building the DG_c_Left
787 for (i = 0; i < nquad0; ++i)
788 {
789 m_dGL_xi1[n][i] = etap0 * dLpm0[i];
790 m_dGL_xi1[n][i] += dLpp0[i];
791 m_dGL_xi1[n][i] *= overeta0;
794 }
795
796 // Building the DG_c_Left
797 for (i = 0; i < nquad1; ++i)
798 {
799 m_dGL_xi2[n][i] = etap1 * dLpm1[i];
800 m_dGL_xi2[n][i] += dLpp1[i];
801 m_dGL_xi2[n][i] *= overeta1;
804 }
805
806 // Building the DG_c_Left
807 for (i = 0; i < nquad2; ++i)
808 {
809 m_dGL_xi3[n][i] = etap2 * dLpm2[i];
810 m_dGL_xi3[n][i] += dLpp2[i];
811 m_dGL_xi3[n][i] *= overeta2;
814 }
815
816 // Building the DG_c_Right
817 for (i = 0; i < nquad0; ++i)
818 {
819 m_dGR_xi1[n][i] = etap0 * dLpm0[i];
820 m_dGR_xi1[n][i] += dLpp0[i];
821 m_dGR_xi1[n][i] *= overeta0;
822 m_dGR_xi1[n][i] += dLp0[i];
824 }
825
826 // Building the DG_c_Right
827 for (i = 0; i < nquad1; ++i)
828 {
829 m_dGR_xi2[n][i] = etap1 * dLpm1[i];
830 m_dGR_xi2[n][i] += dLpp1[i];
831 m_dGR_xi2[n][i] *= overeta1;
832 m_dGR_xi2[n][i] += dLp1[i];
834 }
835
836 // Building the DG_c_Right
837 for (i = 0; i < nquad2; ++i)
838 {
839 m_dGR_xi3[n][i] = etap2 * dLpm2[i];
840 m_dGR_xi3[n][i] += dLpp2[i];
841 m_dGR_xi3[n][i] *= overeta2;
842 m_dGR_xi3[n][i] += dLp2[i];
844 }
845 }
846 break;
847 }
848 default:
849 {
851 break;
852 }
853 }
854}
Array< OneD, Array< OneD, NekDouble > > m_dGR_xi3
Definition: DiffusionLFR.h:85
Array< OneD, Array< OneD, NekDouble > > m_dGL_xi3
Definition: DiffusionLFR.h:84
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 162 of file DiffusionLFR.cpp.
165{
166 int i, n;
167 int nquad0, nquad1;
168 int phys_offset;
169 int nLocalSolutionPts;
170 int nElements = pFields[0]->GetExpSize();
171 int nDimensions = pFields[0]->GetCoordim(0);
172 int nSolutionPts = pFields[0]->GetTotPoints();
173 int nTracePts = pFields[0]->GetTrace()->GetTotPoints();
174
175 m_traceNormals = Array<OneD, Array<OneD, NekDouble>>(nDimensions);
176 for (i = 0; i < nDimensions; ++i)
177 {
178 m_traceNormals[i] = Array<OneD, NekDouble>(nTracePts);
179 }
180 pFields[0]->GetTrace()->GetNormals(m_traceNormals);
181
182 Array<TwoD, const NekDouble> gmat;
183 Array<OneD, const NekDouble> jac;
184
185 m_jac = Array<OneD, NekDouble>(nSolutionPts);
186
187 Array<OneD, NekDouble> auxArray1;
188 Array<OneD, LibUtilities::BasisSharedPtr> base;
189 LibUtilities::PointsKeyVector ptsKeys;
190
191 switch (nDimensions)
192 {
193 case 1:
194 {
195 for (n = 0; n < nElements; ++n)
196 {
197 ptsKeys = pFields[0]->GetExp(n)->GetPointsKeys();
198 nLocalSolutionPts = pFields[0]->GetExp(n)->GetTotPoints();
199 phys_offset = pFields[0]->GetPhys_Offset(n);
200 jac = pFields[0]
201 ->GetExp(n)
202 ->as<LocalRegions::Expansion1D>()
203 ->GetGeom1D()
204 ->GetMetricInfo()
205 ->GetJac(ptsKeys);
206 for (i = 0; i < nLocalSolutionPts; ++i)
207 {
208 m_jac[i + phys_offset] = jac[0];
209 }
210 }
211 break;
212 }
213 case 2:
214 {
215 m_gmat = Array<OneD, Array<OneD, NekDouble>>(4);
216 m_gmat[0] = Array<OneD, NekDouble>(nSolutionPts);
217 m_gmat[1] = Array<OneD, NekDouble>(nSolutionPts);
218 m_gmat[2] = Array<OneD, NekDouble>(nSolutionPts);
219 m_gmat[3] = Array<OneD, NekDouble>(nSolutionPts);
220
221 m_Q2D_e0 = Array<OneD, Array<OneD, NekDouble>>(nElements);
222 m_Q2D_e1 = Array<OneD, Array<OneD, NekDouble>>(nElements);
223 m_Q2D_e2 = Array<OneD, Array<OneD, NekDouble>>(nElements);
224 m_Q2D_e3 = Array<OneD, Array<OneD, NekDouble>>(nElements);
225
226 for (n = 0; n < nElements; ++n)
227 {
228 base = pFields[0]->GetExp(n)->GetBase();
229 nquad0 = base[0]->GetNumPoints();
230 nquad1 = base[1]->GetNumPoints();
231
232 m_Q2D_e0[n] = Array<OneD, NekDouble>(nquad0);
233 m_Q2D_e1[n] = Array<OneD, NekDouble>(nquad1);
234 m_Q2D_e2[n] = Array<OneD, NekDouble>(nquad0);
235 m_Q2D_e3[n] = Array<OneD, NekDouble>(nquad1);
236
237 // Extract the Q factors at each edge point
238 pFields[0]->GetExp(n)->GetTraceQFactors(0, auxArray1 =
239 m_Q2D_e0[n]);
240 pFields[0]->GetExp(n)->GetTraceQFactors(1, auxArray1 =
241 m_Q2D_e1[n]);
242 pFields[0]->GetExp(n)->GetTraceQFactors(2, auxArray1 =
243 m_Q2D_e2[n]);
244 pFields[0]->GetExp(n)->GetTraceQFactors(3, auxArray1 =
245 m_Q2D_e3[n]);
246
247 ptsKeys = pFields[0]->GetExp(n)->GetPointsKeys();
248 nLocalSolutionPts = pFields[0]->GetExp(n)->GetTotPoints();
249 phys_offset = pFields[0]->GetPhys_Offset(n);
250
251 jac = pFields[0]
252 ->GetExp(n)
253 ->as<LocalRegions::Expansion2D>()
254 ->GetGeom2D()
255 ->GetMetricInfo()
256 ->GetJac(ptsKeys);
257 gmat = pFields[0]
258 ->GetExp(n)
259 ->as<LocalRegions::Expansion2D>()
260 ->GetGeom2D()
261 ->GetMetricInfo()
262 ->GetDerivFactors(ptsKeys);
263
264 if (pFields[0]
265 ->GetExp(n)
266 ->as<LocalRegions::Expansion2D>()
267 ->GetGeom2D()
268 ->GetMetricInfo()
269 ->GetGtype() == SpatialDomains::eDeformed)
270 {
271 for (i = 0; i < nLocalSolutionPts; ++i)
272 {
273 m_jac[i + phys_offset] = jac[i];
274 m_gmat[0][i + phys_offset] = gmat[0][i];
275 m_gmat[1][i + phys_offset] = gmat[1][i];
276 m_gmat[2][i + phys_offset] = gmat[2][i];
277 m_gmat[3][i + phys_offset] = gmat[3][i];
278 }
279 }
280 else
281 {
282 for (i = 0; i < nLocalSolutionPts; ++i)
283 {
284 m_jac[i + phys_offset] = jac[0];
285 m_gmat[0][i + phys_offset] = gmat[0][0];
286 m_gmat[1][i + phys_offset] = gmat[1][0];
287 m_gmat[2][i + phys_offset] = gmat[2][0];
288 m_gmat[3][i + phys_offset] = gmat[3][0];
289 }
290 }
291 }
292 break;
293 }
294 case 3:
295 {
297 break;
298 }
299 default:
300 {
302 break;
303 }
304 }
305}
Array< OneD, Array< OneD, NekDouble > > m_gmat
Definition: DiffusionLFR.h:73
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()
|
overrideprotected |
Calculate FR Diffusion for the linear problems using an LDG interface flux.
Definition at line 861 of file DiffusionLFR.cpp.
868{
869 int i, j, n;
870 int phys_offset;
871 Array<OneD, NekDouble> auxArray1, auxArray2;
872
873 Array<OneD, LibUtilities::BasisSharedPtr> Basis;
874 Basis = fields[0]->GetExp(0)->GetBase();
875
876 int nElements = fields[0]->GetExpSize();
877 int nDim = fields[0]->GetCoordim(0);
878 int nSolutionPts = fields[0]->GetTotPoints();
879 int nCoeffs = fields[0]->GetNcoeffs();
880
881 Array<OneD, Array<OneD, NekDouble>> outarrayCoeff(nConvectiveFields);
882 for (i = 0; i < nConvectiveFields; ++i)
883 {
884 outarrayCoeff[i] = Array<OneD, NekDouble>(nCoeffs);
885 }
886
887 // Compute interface numerical fluxes for inarray in physical space
889
890 switch (nDim)
891 {
892 // 1D problems
893 case 1:
894 {
895 for (i = 0; i < nConvectiveFields; ++i)
896 {
897 // Computing the physical first-order discountinuous
898 // derivative
899 for (n = 0; n < nElements; n++)
900 {
901 phys_offset = fields[0]->GetPhys_Offset(n);
902
903 fields[i]->GetExp(n)->PhysDeriv(
904 0, auxArray1 = inarray[i] + phys_offset,
905 auxArray2 = m_DU1[i][0] + phys_offset);
906 }
907
908 // Computing the standard first-order correction
909 // derivative
911 m_DFC1[i][0]);
912
913 // Back to the physical space using global operations
915 &m_DFC1[i][0][0], 1);
917 &m_DFC1[i][0][0], 1);
918
919 // Computing total first order derivatives
921 1, &m_D1[i][0][0], 1);
922
924 1);
925 }
926
927 // Computing interface numerical fluxes for m_D1
928 // in physical space
930
931 for (i = 0; i < nConvectiveFields; ++i)
932 {
933 // Computing the physical second-order discountinuous
934 // derivative
935 for (n = 0; n < nElements; n++)
936 {
937 phys_offset = fields[0]->GetPhys_Offset(n);
938
939 fields[i]->GetExp(n)->PhysDeriv(
940 0, auxArray1 = m_D1[i][0] + phys_offset,
941 auxArray2 = m_DD1[i][0] + phys_offset);
942 }
943
944 // Computing the standard second-order correction
945 // derivative
947 m_DFC2[i][0]);
948
949 // Back to the physical space using global operations
951 &m_DFC2[i][0][0], 1);
953 &m_DFC2[i][0][0], 1);
954
955 // Adding the total divergence to outarray (RHS)
957 1, &outarray[i][0], 1);
958
959 // Primitive Dealiasing 1D
960 if (!(Basis[0]->Collocation()))
961 {
962 fields[i]->FwdTrans(outarray[i], outarrayCoeff[i]);
963 fields[i]->BwdTrans(outarrayCoeff[i], outarray[i]);
964 }
965 }
966 break;
967 }
968 // 2D problems
969 case 2:
970 {
971 for (i = 0; i < nConvectiveFields; ++i)
972 {
973 for (j = 0; j < nDim; ++j)
974 {
975 // Temporary vectors
976 Array<OneD, NekDouble> u1_hat(nSolutionPts, 0.0);
977 Array<OneD, NekDouble> u2_hat(nSolutionPts, 0.0);
978
979 if (j == 0)
980 {
981 Vmath::Vmul(nSolutionPts, &inarray[i][0], 1,
982 &m_gmat[0][0], 1, &u1_hat[0], 1);
983
985 &u1_hat[0], 1);
986
987 Vmath::Vmul(nSolutionPts, &inarray[i][0], 1,
988 &m_gmat[1][0], 1, &u2_hat[0], 1);
989
991 &u2_hat[0], 1);
992 }
993 else if (j == 1)
994 {
995 Vmath::Vmul(nSolutionPts, &inarray[i][0], 1,
996 &m_gmat[2][0], 1, &u1_hat[0], 1);
997
999 &u1_hat[0], 1);
1000
1001 Vmath::Vmul(nSolutionPts, &inarray[i][0], 1,
1002 &m_gmat[3][0], 1, &u2_hat[0], 1);
1003
1005 &u2_hat[0], 1);
1006 }
1007
1008 for (n = 0; n < nElements; n++)
1009 {
1010 phys_offset = fields[0]->GetPhys_Offset(n);
1011
1012 fields[i]->GetExp(n)->StdPhysDeriv(
1013 0, auxArray1 = u1_hat + phys_offset,
1014 auxArray2 = m_tmp1[i][j] + phys_offset);
1015
1016 fields[i]->GetExp(n)->StdPhysDeriv(
1017 1, auxArray1 = u2_hat + phys_offset,
1018 auxArray2 = m_tmp2[i][j] + phys_offset);
1019 }
1020
1023
1024 // Divide by the metric jacobian
1026 &m_DU1[i][j][0], 1);
1027
1028 // Computing the standard first-order correction
1029 // derivatives
1030 DerCFlux_2D(nConvectiveFields, j, fields, inarray[i],
1032 }
1033
1034 // Multiplying derivatives by B matrix to get auxiliary
1035 // variables
1036 for (j = 0; j < nSolutionPts; ++j)
1037 {
1038 // std(q1)
1041 m_DFC1[i][0][j] +
1044 m_DFC1[i][1][j];
1045
1046 // std(q2)
1049 m_DFC1[i][0][j] +
1052 m_DFC1[i][1][j];
1053 }
1054
1055 // Multiplying derivatives by A^(-1) to get back
1056 // into the physical space
1057 for (j = 0; j < nSolutionPts; j++)
1058 {
1059 // q1 = A11^(-1)*std(q1) + A12^(-1)*std(q2)
1062
1063 // q2 = A21^(-1)*std(q1) + A22^(-1)*std(q2)
1066 }
1067
1068 // Computing the physical first-order derivatives
1069 for (j = 0; j < nDim; ++j)
1070 {
1073 }
1074 }
1075
1076 // Computing interface numerical fluxes for m_D1
1077 // in physical space
1079
1080 // Computing the standard second-order discontinuous
1081 // derivatives
1082 for (i = 0; i < nConvectiveFields; ++i)
1083 {
1084 Array<OneD, NekDouble> f_hat(nSolutionPts, 0.0);
1085 Array<OneD, NekDouble> g_hat(nSolutionPts, 0.0);
1086
1087 for (j = 0; j < nSolutionPts; j++)
1088 {
1091 m_jac[j];
1092
1095 m_jac[j];
1096 }
1097
1098 for (n = 0; n < nElements; n++)
1099 {
1100 phys_offset = fields[0]->GetPhys_Offset(n);
1101
1102 fields[0]->GetExp(n)->StdPhysDeriv(
1103 0, auxArray1 = f_hat + phys_offset,
1104 auxArray2 = m_DD1[i][0] + phys_offset);
1105
1106 fields[0]->GetExp(n)->StdPhysDeriv(
1107 1, auxArray1 = g_hat + phys_offset,
1108 auxArray2 = m_DD1[i][1] + phys_offset);
1109 }
1110
1111 // Divergence of the standard discontinuous flux
1113 1, &m_divFD[i][0], 1);
1114
1115 // Divergence of the standard correction flux
1116 if (Basis[0]->GetPointsType() ==
1118 Basis[1]->GetPointsType() ==
1120 {
1121 DivCFlux_2D_Gauss(nConvectiveFields, fields, f_hat, g_hat,
1123 }
1124 else
1125 {
1128 }
1129
1130 // Divergence of the standard final flux
1132 &outarray[i][0], 1);
1133
1134 // Dividing by the jacobian to get back into
1135 // physical space
1137 &outarray[i][0], 1);
1138
1139 // Primitive Dealiasing 2D
1140 if (!(Basis[0]->Collocation()))
1141 {
1142 fields[i]->FwdTrans(outarray[i], outarrayCoeff[i]);
1143 fields[i]->BwdTrans(outarrayCoeff[i], outarray[i]);
1144 }
1145 } // Close loop on nConvectiveFields
1146 break;
1147 }
1148 // 3D problems
1149 case 3:
1150 {
1152 break;
1153 }
1154 }
1155}
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > m_tmp2
Definition: DiffusionLFR.h:101
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > m_IF1
Definition: DiffusionLFR.h:89
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:1437
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > m_DFC1
Definition: DiffusionLFR.h:91
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > m_DD1
Definition: DiffusionLFR.h:94
Array< OneD, Array< OneD, NekDouble > > m_divFC
Definition: DiffusionLFR.h:98
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:1556
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:1752
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:1298
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > m_DU1
Definition: DiffusionLFR.h:90
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > m_DFC2
Definition: DiffusionLFR.h:96
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > m_tmp1
Definition: DiffusionLFR.h:100
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > m_D1
Definition: DiffusionLFR.h:93
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:1161
Array< OneD, Array< OneD, NekDouble > > m_IF2
Definition: DiffusionLFR.h:95
Array< OneD, Array< OneD, NekDouble > > m_divFD
Definition: DiffusionLFR.h:97
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > m_BD1
Definition: DiffusionLFR.h:92
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:1940
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_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 81 of file DiffusionLFR.cpp.
84{
85 m_session = pSession;
86 SetupMetrics(pSession, pFields);
87 SetupCFunctions(pSession, pFields);
88
89 // Initialising arrays
90 int i, j;
91 int nConvectiveFields = pFields.size();
92 int nDim = pFields[0]->GetCoordim(0);
93 int nSolutionPts = pFields[0]->GetTotPoints();
94 int nTracePts = pFields[0]->GetTrace()->GetTotPoints();
95
96 m_IF1 = Array<OneD, Array<OneD, Array<OneD, NekDouble>>>(nConvectiveFields);
97 m_DU1 = Array<OneD, Array<OneD, Array<OneD, NekDouble>>>(nConvectiveFields);
98 m_DFC1 =
99 Array<OneD, Array<OneD, Array<OneD, NekDouble>>>(nConvectiveFields);
100 m_BD1 = Array<OneD, Array<OneD, Array<OneD, NekDouble>>>(nConvectiveFields);
101 m_D1 = Array<OneD, Array<OneD, Array<OneD, NekDouble>>>(nConvectiveFields);
102 m_DD1 = Array<OneD, Array<OneD, Array<OneD, NekDouble>>>(nConvectiveFields);
103 m_IF2 = Array<OneD, Array<OneD, NekDouble>>(nConvectiveFields);
104 m_DFC2 =
105 Array<OneD, Array<OneD, Array<OneD, NekDouble>>>(nConvectiveFields);
106 m_divFD = Array<OneD, Array<OneD, NekDouble>>(nConvectiveFields);
107 m_divFC = Array<OneD, Array<OneD, NekDouble>>(nConvectiveFields);
108
109 m_tmp1 =
110 Array<OneD, Array<OneD, Array<OneD, NekDouble>>>(nConvectiveFields);
111 m_tmp2 =
112 Array<OneD, Array<OneD, Array<OneD, NekDouble>>>(nConvectiveFields);
113
114 for (i = 0; i < nConvectiveFields; ++i)
115 {
116 m_IF1[i] = Array<OneD, Array<OneD, NekDouble>>(nDim);
117 m_DU1[i] = Array<OneD, Array<OneD, NekDouble>>(nDim);
118 m_DFC1[i] = Array<OneD, Array<OneD, NekDouble>>(nDim);
119 m_BD1[i] = Array<OneD, Array<OneD, NekDouble>>(nDim);
120 m_D1[i] = Array<OneD, Array<OneD, NekDouble>>(nDim);
121 m_DD1[i] = Array<OneD, Array<OneD, NekDouble>>(nDim);
122 m_IF2[i] = Array<OneD, NekDouble>(nTracePts, 0.0);
123 m_DFC2[i] = Array<OneD, Array<OneD, NekDouble>>(nDim);
124 m_divFD[i] = Array<OneD, NekDouble>(nSolutionPts, 0.0);
125 m_divFC[i] = Array<OneD, NekDouble>(nSolutionPts, 0.0);
126
127 m_tmp1[i] = Array<OneD, Array<OneD, NekDouble>>(nDim);
128 m_tmp2[i] = Array<OneD, Array<OneD, NekDouble>>(nDim);
129
130 for (j = 0; j < nDim; ++j)
131 {
132 m_IF1[i][j] = Array<OneD, NekDouble>(nTracePts, 0.0);
133 m_DU1[i][j] = Array<OneD, NekDouble>(nSolutionPts, 0.0);
134 m_DFC1[i][j] = Array<OneD, NekDouble>(nSolutionPts, 0.0);
135 m_BD1[i][j] = Array<OneD, NekDouble>(nSolutionPts, 0.0);
136 m_D1[i][j] = Array<OneD, NekDouble>(nSolutionPts, 0.0);
137 m_DD1[i][j] = Array<OneD, NekDouble>(nSolutionPts, 0.0);
138 m_DFC2[i][j] = Array<OneD, NekDouble>(nSolutionPts, 0.0);
139
140 m_tmp1[i][j] = Array<OneD, NekDouble>(nSolutionPts, 0.0);
141 m_tmp2[i][j] = Array<OneD, NekDouble>(nSolutionPts, 0.0);
142 }
143 }
144}
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:162
LibUtilities::SessionReaderSharedPtr m_session
Definition: DiffusionLFR.h:70
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:324
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 1235 of file DiffusionLFR.cpp.
1239{
1240 // Variable initialisation
1241 int i, e, id1, id2;
1242 int nBndEdgePts, nBndEdges;
1243 int cnt = 0;
1244 int nBndRegions = fields[var]->GetBndCondExpansions().size();
1245 int nTracePts = fields[0]->GetTrace()->GetTotPoints();
1246
1247 // Extract physical values of the fields at the trace space
1248 Array<OneD, NekDouble> uplus(nTracePts);
1249 fields[var]->ExtractTracePhys(ufield, uplus);
1250
1251 // Impose boundary conditions
1252 for (i = 0; i < nBndRegions; ++i)
1253 {
1254 // Number of boundary edges related to bcRegion 'i'
1255 nBndEdges = fields[var]->GetBndCondExpansions()[i]->GetExpSize();
1256
1257 // Weakly impose boundary conditions by modifying flux values
1258 for (e = 0; e < nBndEdges; ++e)
1259 {
1260 // Number of points on the expansion
1261 nBndEdgePts = fields[var]
1262 ->GetBndCondExpansions()[i]
1263 ->GetExp(e)
1264 ->GetTotPoints();
1265
1266 // Offset of the boundary expansion
1267 id1 = fields[var]->GetBndCondExpansions()[i]->GetPhys_Offset(e);
1268
1269 // Offset of the trace space related to boundary expansion
1270 id2 = fields[0]->GetTrace()->GetPhys_Offset(
1271 fields[0]->GetTraceMap()->GetBndCondIDToGlobalTraceID(cnt++));
1272
1273 // Dirichlet bcs ==> uflux = gD
1274 if (fields[var]
1275 ->GetBndConditions()[i]
1276 ->GetBoundaryConditionType() == SpatialDomains::eDirichlet)
1277 {
1278 Vmath::Vcopy(
1279 nBndEdgePts,
1280 &(fields[var]->GetBndCondExpansions()[i]->GetPhys())[id1],
1281 1, &penaltyflux[id2], 1);
1282 }
1283 // Neumann bcs ==> uflux = u+
1284 else if ((fields[var]->GetBndConditions()[i])
1285 ->GetBoundaryConditionType() ==
1286 SpatialDomains::eNeumann)
1287 {
1288 Vmath::Vcopy(nBndEdgePts, &uplus[id2], 1, &penaltyflux[id2], 1);
1289 }
1290 }
1291 }
1292}
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 1371 of file DiffusionLFR.cpp.
1375{
1376 int i, e, id1, id2;
1377 int nBndEdges, nBndEdgePts;
1378 int nBndRegions = fields[var]->GetBndCondExpansions().size();
1379 int nTracePts = fields[0]->GetTrace()->GetTotPoints();
1380
1381 Array<OneD, NekDouble> uterm(nTracePts);
1382 Array<OneD, NekDouble> qtemp(nTracePts);
1383 int cnt = 0;
1384
1385 fields[var]->ExtractTracePhys(qfield, qtemp);
1386
1387 for (i = 0; i < nBndRegions; ++i)
1388 {
1389 nBndEdges = fields[var]->GetBndCondExpansions()[i]->GetExpSize();
1390
1391 // Weakly impose boundary conditions by modifying flux values
1392 for (e = 0; e < nBndEdges; ++e)
1393 {
1394 nBndEdgePts = fields[var]
1395 ->GetBndCondExpansions()[i]
1396 ->GetExp(e)
1397 ->GetTotPoints();
1398
1399 id1 = fields[var]->GetBndCondExpansions()[i]->GetPhys_Offset(e);
1400
1401 id2 = fields[0]->GetTrace()->GetPhys_Offset(
1402 fields[0]->GetTraceMap()->GetBndCondIDToGlobalTraceID(cnt++));
1403
1404 // For Dirichlet boundary condition:
1405 // qflux = q+ - C_11 (u+ - g_D) (nx, ny)
1406 if (fields[var]
1407 ->GetBndConditions()[i]
1408 ->GetBoundaryConditionType() == SpatialDomains::eDirichlet)
1409 {
1411 &qtemp[id2], 1, &penaltyflux[id2], 1);
1412 }
1413 // For Neumann boundary condition: qflux = g_N
1414 else if ((fields[var]->GetBndConditions()[i])
1415 ->GetBoundaryConditionType() ==
1416 SpatialDomains::eNeumann)
1417 {
1418 Vmath::Vmul(
1419 nBndEdgePts, &m_traceNormals[dir][id2], 1,
1420 &(fields[var]->GetBndCondExpansions()[i]->GetPhys())[id1],
1421 1, &penaltyflux[id2], 1);
1422 }
1423 }
1424 }
1425}
References Nektar::SpatialDomains::eDirichlet, Nektar::SpatialDomains::eNeumann, m_traceNormals, and Vmath::Vmul().
Referenced by NumFluxforVector().
Member Data Documentation
◆ m_BD1
Definition at line 92 of file DiffusionLFR.h.
Referenced by v_Diffuse(), and v_InitObject().
◆ m_D1
Definition at line 93 of file DiffusionLFR.h.
Referenced by v_Diffuse(), and v_InitObject().
◆ m_DD1
Definition at line 94 of file DiffusionLFR.h.
Referenced by v_Diffuse(), and v_InitObject().
◆ m_DFC1
|
private |
Definition at line 91 of file DiffusionLFR.h.
Referenced by v_Diffuse(), and v_InitObject().
◆ m_DFC2
|
private |
Definition at line 96 of file DiffusionLFR.h.
Referenced by v_Diffuse(), and v_InitObject().
◆ m_dGL_xi1
Definition at line 80 of file DiffusionLFR.h.
Referenced by DerCFlux_1D(), DerCFlux_2D(), DivCFlux_2D(), DivCFlux_2D_Gauss(), and SetupCFunctions().
◆ m_dGL_xi2
Definition at line 82 of file DiffusionLFR.h.
Referenced by DerCFlux_2D(), DivCFlux_2D(), DivCFlux_2D_Gauss(), and SetupCFunctions().
◆ m_dGL_xi3
Definition at line 84 of file DiffusionLFR.h.
Referenced by SetupCFunctions().
◆ m_dGR_xi1
Definition at line 81 of file DiffusionLFR.h.
Referenced by DerCFlux_1D(), DerCFlux_2D(), DivCFlux_2D(), DivCFlux_2D_Gauss(), and SetupCFunctions().
◆ m_dGR_xi2
Definition at line 83 of file DiffusionLFR.h.
Referenced by DerCFlux_2D(), DivCFlux_2D(), DivCFlux_2D_Gauss(), and SetupCFunctions().
◆ m_dGR_xi3
Definition at line 85 of file DiffusionLFR.h.
Referenced by SetupCFunctions().
◆ m_diffType
|
protected |
Definition at line 55 of file DiffusionLFR.h.
Referenced by SetupCFunctions().
◆ m_divFC
Definition at line 98 of file DiffusionLFR.h.
Referenced by v_Diffuse(), and v_InitObject().
◆ m_divFD
Definition at line 97 of file DiffusionLFR.h.
Referenced by v_Diffuse(), and v_InitObject().
◆ m_DU1
Definition at line 90 of file DiffusionLFR.h.
Referenced by v_Diffuse(), and v_InitObject().
◆ m_gmat
Definition at line 73 of file DiffusionLFR.h.
Referenced by SetupMetrics(), and v_Diffuse().
◆ m_IF1
Definition at line 89 of file DiffusionLFR.h.
Referenced by v_Diffuse(), and v_InitObject().
◆ m_IF2
Definition at line 95 of file DiffusionLFR.h.
Referenced by v_Diffuse(), and v_InitObject().
◆ m_Ixm
|
private |
Definition at line 86 of file DiffusionLFR.h.
◆ m_Ixp
|
private |
Definition at line 87 of file DiffusionLFR.h.
◆ m_jac
Definition at line 72 of file DiffusionLFR.h.
Referenced by SetupMetrics(), and v_Diffuse().
◆ m_Q2D_e0
Definition at line 75 of file DiffusionLFR.h.
Referenced by DivCFlux_2D(), DivCFlux_2D_Gauss(), and SetupMetrics().
◆ m_Q2D_e1
Definition at line 76 of file DiffusionLFR.h.
Referenced by DivCFlux_2D(), DivCFlux_2D_Gauss(), and SetupMetrics().
◆ m_Q2D_e2
Definition at line 77 of file DiffusionLFR.h.
Referenced by DivCFlux_2D(), DivCFlux_2D_Gauss(), and SetupMetrics().
◆ m_Q2D_e3
Definition at line 78 of file DiffusionLFR.h.
Referenced by DivCFlux_2D(), DivCFlux_2D_Gauss(), and SetupMetrics().
◆ m_session
|
private |
Definition at line 70 of file DiffusionLFR.h.
Referenced by v_InitObject().
◆ m_tmp1
|
private |
Definition at line 100 of file DiffusionLFR.h.
Referenced by v_Diffuse(), and v_InitObject().
◆ m_tmp2
|
private |
Definition at line 101 of file DiffusionLFR.h.
Referenced by v_Diffuse(), and v_InitObject().
◆ m_traceNormals
Definition at line 69 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:197
static DiffusionSharedPtr create(std::string diffType)
Definition: DiffusionLFR.h:45
DiffusionFactory & GetDiffusionFactory()
Definition: Diffusion.cpp:39
Definition at line 50 of file DiffusionLFR.h.
