Nektar++
ForcingAxiSymmetric.cpp
Go to the documentation of this file.
1///////////////////////////////////////////////////////////////////////////////
2//
3// File: ForcingAxiSymmetric.cpp
4//
5// For more information, please see: http://www.nektar.info
6//
7// The MIT License
8//
9// Copyright (c) 2006 Division of Applied Mathematics, Brown University (USA),
10// Department of Aeronautics, Imperial College London (UK), and Scientific
11// Computing and Imaging Institute, University of Utah (USA).
12//
13// Permission is hereby granted, free of charge, to any person obtaining a
14// copy of this software and associated documentation files (the "Software"),
15// to deal in the Software without restriction, including without limitation
16// the rights to use, copy, modify, merge, publish, distribute, sublicense,
17// and/or sell copies of the Software, and to permit persons to whom the
18// Software is furnished to do so, subject to the following conditions:
19//
20// The above copyright notice and this permission notice shall be included
21// in all copies or substantial portions of the Software.
22//
23// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
24// OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
25// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
26// THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
27// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
28// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
29// DEALINGS IN THE SOFTWARE.
30//
31// Description: Forcing for axi-symmetric flow.
32//
33///////////////////////////////////////////////////////////////////////////////
34
36
37namespace Nektar
38{
39
42 "AxiSymmetric", ForcingAxiSymmetric::create,
43 "Forcing for axi-symmetric flow (around x=0)");
44
47 const std::weak_ptr<SolverUtils::EquationSystem> &pEquation)
48 : Forcing(pSession, pEquation)
49{
50}
51
54 const unsigned int &pNumForcingFields,
55 [[maybe_unused]] const TiXmlElement *pForce)
56{
57 int spacedim = pFields[0]->GetGraph()->GetSpaceDimension();
58 int nPoints = pFields[0]->GetTotPoints();
59
60 m_NumVariable = pNumForcingFields;
62 spacedim);
63
64 // Get coordinates
66 for (int i = 0; i < 3; i++)
67 {
68 coords[i] = Array<OneD, NekDouble>(nPoints);
69 }
70 pFields[0]->GetCoords(coords[0], coords[1], coords[2]);
71
72 // Calculate fac = -1/r if r!=0, fac = 0 if r == 0
74 for (int i = 0; i < nPoints; ++i)
75 {
76 if (coords[0][i] < NekConstants::kNekZeroTol)
77 {
78 m_geomFactor[i] = 0;
79 }
80 else
81 {
82 m_geomFactor[i] = -1.0 / coords[0][i];
83 }
84 }
85
86 // Project m_geomFactor to solution space
88 for (int i = 0; i < m_NumVariable; ++i)
89 {
90 m_Forcing[i] = Array<OneD, NekDouble>(pFields[0]->GetTotPoints(), 0.0);
91 }
92}
93
96 const Array<OneD, Array<OneD, NekDouble>> &inarray,
98 [[maybe_unused]] const NekDouble &time)
99{
100 int nPoints = pFields[0]->GetTotPoints();
101
102 // Get (E+p)
103 Array<OneD, NekDouble> tmp(nPoints, 0.0);
104 m_varConv->GetPressure(inarray, tmp);
105 Vmath::Vadd(nPoints, tmp, 1, inarray[m_NumVariable - 1], 1, tmp, 1);
106
107 // F-rho = -1/r *rhou
108 Vmath::Vmul(nPoints, m_geomFactor, 1, inarray[1], 1, m_Forcing[0], 1);
109
110 // F-rhou_r = -1/r *rhou_r * u_r and F-rhou_y = -1/r *rhou_y * u_r
111 for (int i = 1; i < 3; ++i)
112 {
113 Vmath::Vmul(nPoints, inarray[1], 1, inarray[i], 1, m_Forcing[i], 1);
114 Vmath::Vdiv(nPoints, m_Forcing[i], 1, inarray[0], 1, m_Forcing[i], 1);
115 Vmath::Vmul(nPoints, m_Forcing[i], 1, m_geomFactor, 1, m_Forcing[i], 1);
116 }
117
118 // F-E = -1/r *(E+p)*u
119 Vmath::Vmul(nPoints, inarray[1], 1, tmp, 1, m_Forcing[m_NumVariable - 1],
120 1);
121 Vmath::Vdiv(nPoints, m_Forcing[m_NumVariable - 1], 1, inarray[0], 1,
122 m_Forcing[m_NumVariable - 1], 1);
123 Vmath::Vmul(nPoints, m_Forcing[m_NumVariable - 1], 1, m_geomFactor, 1,
124 m_Forcing[m_NumVariable - 1], 1);
125
126 // Swirl
127 if (m_NumVariable == 5)
128 {
129 // F-rhou_r -= (-1/r) * rho * u_theta * u_theta
130 Vmath::Vmul(nPoints, inarray[3], 1, inarray[3], 1, tmp, 1);
131 Vmath::Vdiv(nPoints, tmp, 1, inarray[0], 1, tmp, 1);
132 Vmath::Vmul(nPoints, tmp, 1, m_geomFactor, 1, tmp, 1);
133 Vmath::Vsub(nPoints, m_Forcing[1], 1, tmp, 1, m_Forcing[1], 1);
134
135 // F-rhou_theta = 2 * (-1/r *rhou_theta * u_r)
136 Vmath::Vmul(nPoints, inarray[1], 1, inarray[3], 1, m_Forcing[3], 1);
137 Vmath::Vdiv(nPoints, m_Forcing[3], 1, inarray[0], 1, m_Forcing[3], 1);
138 Vmath::Vmul(nPoints, m_Forcing[3], 1, m_geomFactor, 1, m_Forcing[3], 1);
139 Vmath::Smul(nPoints, 2.0, m_Forcing[3], 1, m_Forcing[3], 1);
140 }
141
142 // Apply forcing
143 for (int i = 0; i < m_NumVariable; i++)
144 {
145 Vmath::Vadd(nPoints, outarray[i], 1, m_Forcing[i], 1, outarray[i], 1);
146 }
147}
148
149} // namespace Nektar
void v_Apply(const Array< OneD, MultiRegions::ExpListSharedPtr > &fields, const Array< OneD, Array< OneD, NekDouble > > &inarray, Array< OneD, Array< OneD, NekDouble > > &outarray, const NekDouble &time) override
VariableConverterSharedPtr m_varConv
static std::string className
Name of the class.
Array< OneD, NekDouble > m_geomFactor
void v_InitObject(const Array< OneD, MultiRegions::ExpListSharedPtr > &pFields, const unsigned int &pNumForcingFields, const TiXmlElement *pForce) override
ForcingAxiSymmetric(const LibUtilities::SessionReaderSharedPtr &pSession, const std::weak_ptr< SolverUtils::EquationSystem > &pEquation)
static SolverUtils::ForcingSharedPtr create(const LibUtilities::SessionReaderSharedPtr &pSession, const std::weak_ptr< SolverUtils::EquationSystem > &pEquation, const Array< OneD, MultiRegions::ExpListSharedPtr > &pFields, const unsigned int &pNumForcingFields, const TiXmlElement *pForce)
Creates an instance of this class.
tKey RegisterCreatorFunction(tKey idKey, CreatorFunction classCreator, std::string pDesc="")
Register a class with the factory.
static std::shared_ptr< DataType > AllocateSharedPtr(const Args &...args)
Allocate a shared pointer from the memory pool.
Defines a forcing term to be explicitly applied.
Definition: Forcing.h:71
int m_NumVariable
Number of variables.
Definition: Forcing.h:129
Array< OneD, Array< OneD, NekDouble > > m_Forcing
Evaluated forcing function.
Definition: Forcing.h:127
LibUtilities::SessionReaderSharedPtr m_session
Session reader.
Definition: Forcing.h:123
std::shared_ptr< SessionReader > SessionReaderSharedPtr
static const NekDouble kNekZeroTol
ForcingFactory & GetForcingFactory()
Declaration of the forcing factory singleton.
Definition: Forcing.cpp:41
double NekDouble
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
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 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
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