Nektar++
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
ForcingQuasi1D.cpp
Go to the documentation of this file.
1 ///////////////////////////////////////////////////////////////////////////////
2 //
3 // File: ForcingQuasi1D.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 quasi-1D nozzle flow.
32 //
33 ///////////////////////////////////////////////////////////////////////////////
34 
35 #include <boost/core/ignore_unused.hpp>
36 
38 
39 using namespace std;
40 
41 namespace Nektar
42 {
43 std::string ForcingQuasi1D::className =
45  "Quasi1D", ForcingQuasi1D::create, "Quasi-1D nozzle Forcing");
46 
47 ForcingQuasi1D::ForcingQuasi1D(
49  const std::weak_ptr<SolverUtils::EquationSystem> &pEquation)
50  : Forcing(pSession, pEquation)
51 {
52 }
53 
56  const unsigned int &pNumForcingFields, const TiXmlElement *pForce)
57 {
58  m_NumVariable = pNumForcingFields;
59  m_varConv =
61 
62  ASSERTL0(pFields[0]->GetGraph()->GetSpaceDimension() == 1,
63  "ForcingQuasi1D requires a 1D problem.");
64 
65  const TiXmlElement *funcNameElmt = pForce->FirstChildElement("AREAFCN");
66  if (!funcNameElmt)
67  {
68  ASSERTL0(funcNameElmt,
69  "Requires AREAFCN tag "
70  "specifying function name which prescribes nozzle area.");
71  }
72 
73  string funcName = funcNameElmt->GetText();
74  ASSERTL0(m_session->DefinesFunction(funcName),
75  "Function '" + funcName + "' not defined.");
76 
77  // Evaluate geometrical term -Ax/A for forcing
78  m_geomFactor = Array<OneD, NekDouble>(pFields[0]->GetTotPoints(), 0.0);
79  Array<OneD, NekDouble> tmp(pFields[0]->GetTotPoints(), 0.0);
80 
81  std::string sFieldStr = m_session->GetVariable(0);
82  ASSERTL0(m_session->DefinesFunction(funcName, sFieldStr),
83  "Variable '" + sFieldStr + "' not defined.");
84  GetFunction(pFields, m_session, funcName, true)
85  ->Evaluate(sFieldStr, m_geomFactor, 0.0);
86 
87  // Check if DADXFCN is defined
88  const TiXmlElement *dAFuncNameElmt = pForce->FirstChildElement("DADXFCN");
89  if (dAFuncNameElmt)
90  {
91  funcName = dAFuncNameElmt->GetText();
92  ASSERTL0(m_session->DefinesFunction(funcName),
93  "Function '" + funcName + "' not defined.");
94  ASSERTL0(m_session->DefinesFunction(funcName, sFieldStr),
95  "Variable '" + sFieldStr + "' not defined.");
96  GetFunction(pFields, m_session, funcName, true)
97  ->Evaluate(sFieldStr, tmp, 0.0);
98  }
99  else
100  {
101  // Numerically evaluate dA/dX
102  pFields[0]->PhysDeriv(MultiRegions::DirCartesianMap[0], m_geomFactor,
103  tmp);
104  }
105 
106  Vmath::Vdiv(pFields[0]->GetTotPoints(), tmp, 1, m_geomFactor, 1,
107  m_geomFactor, 1);
108  Vmath::Neg(pFields[0]->GetTotPoints(), m_geomFactor, 1);
109 
110  // Project m_geomFactor to solution space
111  Array<OneD, NekDouble> tmpCoeff(pFields[0]->GetNcoeffs(), 0.0);
112  pFields[0]->FwdTransLocalElmt(m_geomFactor, tmpCoeff);
113  pFields[0]->BwdTrans(tmpCoeff, m_geomFactor);
114 
116  for (int i = 0; i < m_NumVariable; ++i)
117  {
118  m_Forcing[i] = Array<OneD, NekDouble>(pFields[0]->GetTotPoints(), 0.0);
119  }
120 }
121 
124  const Array<OneD, Array<OneD, NekDouble>> &inarray,
125  Array<OneD, Array<OneD, NekDouble>> &outarray, const NekDouble &time)
126 {
127  boost::ignore_unused(time);
128 
129  int nPoints = pFields[0]->GetTotPoints();
130 
131  // Get (E+p)
132  Array<OneD, NekDouble> tmp(nPoints, 0.0);
133  m_varConv->GetPressure(inarray, tmp);
134  Vmath::Vadd(nPoints, tmp, 1, inarray[2], 1, tmp, 1);
135 
136  // F-rho = -Ax/A *rhou
137  Vmath::Vmul(nPoints, m_geomFactor, 1, inarray[1], 1, m_Forcing[0], 1);
138 
139  // F-rhou = -Ax/A *rhou*u
140  Vmath::Vmul(nPoints, m_geomFactor, 1, inarray[1], 1, m_Forcing[1], 1);
141  Vmath::Vmul(nPoints, m_Forcing[1], 1, inarray[1], 1, m_Forcing[1], 1);
142  Vmath::Vdiv(nPoints, m_Forcing[1], 1, inarray[0], 1, m_Forcing[1], 1);
143 
144  // F-E = -Ax/A *(E+p)*u
145  Vmath::Vmul(nPoints, m_geomFactor, 1, inarray[1], 1, m_Forcing[2], 1);
146  Vmath::Vmul(nPoints, m_Forcing[2], 1, tmp, 1, m_Forcing[2], 1);
147  Vmath::Vdiv(nPoints, m_Forcing[2], 1, inarray[0], 1, m_Forcing[2], 1);
148 
149  // Apply forcing
150  for (int i = 0; i < m_NumVariable; i++)
151  {
152  Vmath::Vadd(nPoints, outarray[i], 1, m_Forcing[i], 1, outarray[i], 1);
153  }
154 }
155 
156 } // namespace Nektar
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:215
virtual void v_InitObject(const Array< OneD, MultiRegions::ExpListSharedPtr > &pFields, const unsigned int &pNumForcingFields, const TiXmlElement *pForce) override
VariableConverterSharedPtr m_varConv
Array< OneD, NekDouble > m_geomFactor
virtual 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
tKey RegisterCreatorFunction(tKey idKey, CreatorFunction classCreator, std::string pDesc="")
Register a class with the factory.
Definition: NekFactory.hpp:198
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:73
int m_NumVariable
Number of variables.
Definition: Forcing.h:122
Array< OneD, Array< OneD, NekDouble > > m_Forcing
Evaluated forcing function.
Definition: Forcing.h:120
SOLVER_UTILS_EXPORT SessionFunctionSharedPtr GetFunction(const Array< OneD, MultiRegions::ExpListSharedPtr > &pFields, const LibUtilities::SessionReaderSharedPtr &pSession, std::string pName, bool pCache=false)
Get a SessionFunction by name.
Definition: Forcing.cpp:190
LibUtilities::SessionReaderSharedPtr m_session
Session reader.
Definition: Forcing.h:116
std::shared_ptr< SessionReader > SessionReaderSharedPtr
MultiRegions::Direction const DirCartesianMap[]
Definition: ExpList.h:91
ForcingFactory & GetForcingFactory()
Declaration of the forcing factory singleton.
Definition: Forcing.cpp:44
The above copyright notice and this permission notice shall be included.
Definition: CoupledSolver.h:2
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.cpp:209
void Neg(int n, T *x, const int incx)
Negate x = -x.
Definition: Vmath.cpp:518
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:359
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:284