Nektar++
ForcingAbsorption.cpp
Go to the documentation of this file.
1///////////////////////////////////////////////////////////////////////////////
2//
3// File: ForcingAbsorption.cpp
4//
5// For more information, please see: http://www.nektar.info
6//
7// The MIT License
8//
9// Copyright (c) 2016 Kilian Lackhove
10// Copyright (c) 2006 Division of Applied Mathematics, Brown University (USA),
11// Department of Aeronautics, Imperial College London (UK), and Scientific
12// Computing and Imaging Institute, University of Utah (USA).
13//
14// Permission is hereby granted, free of charge, to any person obtaining a
15// copy of this software and associated documentation files (the "Software"),
16// to deal in the Software without restriction, including without limitation
17// the rights to use, copy, modify, merge, publish, distribute, sublicense,
18// and/or sell copies of the Software, and to permit persons to whom the
19// Software is furnished to do so, subject to the following conditions:
20//
21// The above copyright notice and this permission notice shall be included
22// in all copies or substantial portions of the Software.
23//
24// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
25// OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
26// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
27// THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
28// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
29// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
30// DEALINGS IN THE SOFTWARE.
31//
32// Description: Absorption layer forcing
33//
34///////////////////////////////////////////////////////////////////////////////
35
37
39
40namespace Nektar::SolverUtils
41{
42
45 "Absorption", ForcingAbsorption::create, "Forcing Absorption");
46
49 const std::weak_ptr<EquationSystem> &pEquation)
50 : Forcing(pSession, pEquation), m_hasRefFlow(false), m_hasRefFlowTime(false)
51{
52}
53
56 const unsigned int &pNumForcingFields, const TiXmlElement *pForce)
57{
58 m_NumVariable = pNumForcingFields;
59 int npts = pFields[0]->GetTotPoints();
60
61 // Check for homogeneous expansion
62 // for later transformation of evaluated Absorption (SessionFunction)
63 m_homogeneous = pFields[0]->GetExpType() == MultiRegions::e3DH1D ||
64 pFields[0]->GetExpType() == MultiRegions::e3DH2D;
65
66 CalcAbsorption(pFields, pForce);
67
69 for (int i = 0; i < m_NumVariable; ++i)
70 {
71 m_Forcing[i] = Array<OneD, NekDouble>(npts, 0.0);
72 }
73
74 const TiXmlElement *funcNameElmt = pForce->FirstChildElement("REFFLOW");
75 if (funcNameElmt)
76 {
77 std::string funcName = funcNameElmt->GetText();
78 ASSERTL0(m_session->DefinesFunction(funcName),
79 "Function '" + funcName + "' not defined.");
81 for (int i = 0; i < m_NumVariable; ++i)
82 {
83 std::string s_FieldStr = m_session->GetVariable(i);
84 ASSERTL0(m_session->DefinesFunction(funcName, s_FieldStr),
85 "Variable '" + s_FieldStr + "' not defined.");
86 m_Refflow[i] = Array<OneD, NekDouble>(npts, 0.0);
87 GetFunction(pFields, m_session, funcName)
88 ->Evaluate(s_FieldStr, m_Refflow[i]);
89 }
90 m_hasRefFlow = true;
91 }
92
93 funcNameElmt = pForce->FirstChildElement("REFFLOWTIME");
94 if (funcNameElmt)
95 {
96 m_funcNameTime = funcNameElmt->GetText();
97 m_hasRefFlowTime = true;
98 m_hasRefFlow = true;
100 for (int i = 0; i < m_NumVariable; ++i)
101 {
102 std::string s_FieldStr = m_session->GetVariable(i);
103 ASSERTL0(m_session->DefinesFunction(m_funcNameTime, s_FieldStr),
104 "Variable '" + s_FieldStr + "' not defined.");
105 m_Refflow[i] = Array<OneD, NekDouble>(npts, 0.0);
106 }
107 }
108}
109
112 &pFields,
113 const TiXmlElement *pForce)
114{
115 const TiXmlElement *funcNameElmt = pForce->FirstChildElement("COEFF");
116 ASSERTL0(funcNameElmt,
117 "Requires COEFF tag, specifying function "
118 "name which prescribes absorption layer coefficient.");
119 std::string funcName = funcNameElmt->GetText();
120 ASSERTL0(m_session->DefinesFunction(funcName),
121 "Function '" + funcName + "' not defined.");
122
123 int npts = pFields[0]->GetTotPoints();
124
126 for (int i = 0; i < m_NumVariable; ++i)
127 {
128 m_Absorption[i] = Array<OneD, NekDouble>(npts, 0.0);
129 }
130
131 funcNameElmt = pForce->FirstChildElement("BOUNDARYREGIONS");
132 if (funcNameElmt)
133 {
134 ASSERTL0(
135 ParseUtils::GenerateVector(funcNameElmt->GetText(), m_bRegions),
136 "Unable to process list of BOUNDARYREGIONS in Absorption "
137 "Forcing: " +
138 std::string(funcNameElmt->GetText()));
139
140 // alter m_bRegions so that it contains the boundaryRegions of this rank
141 std::vector<unsigned int> localBRegions;
143 pFields[0]->GetGraph());
145 bcs.GetBoundaryRegions();
146 SpatialDomains::BoundaryRegionCollection::iterator it1;
147 int n = 0;
148 for (it1 = regions.begin(); it1 != regions.end(); ++it1)
149 {
150 if (std::find(m_bRegions.begin(), m_bRegions.end(), it1->first) !=
151 m_bRegions.end())
152 {
153 localBRegions.push_back(n);
154 }
155 n++;
156 }
157 m_bRegions = localBRegions;
158
159 if (m_bRegions.size() == 0)
160 {
161 return;
162 }
163
164 std::vector<Array<OneD, const NekDouble>> points;
165
167 for (int i = 0; i < 3; i++)
168 {
169 x[i] = Array<OneD, NekDouble>(npts, 0.0);
170 }
171 pFields[0]->GetCoords(x[0], x[1], x[2]);
172 for (int i = 0; i < 3; i++)
173 {
174 points.push_back(x[i]);
175 }
176
177 Array<OneD, NekDouble> t(npts, 0.0);
178 points.push_back(t);
179
180 Array<OneD, NekDouble> r(npts, 0.0);
181 std::vector<unsigned int>::iterator it;
182 std::vector<BPointPair> inPoints;
184 for (it = m_bRegions.begin(); it != m_bRegions.end(); ++it)
185 {
186 int bpts = pFields[0]->GetBndCondExpansions()[*it]->GetNpoints();
187 for (int i = 0; i < 3; i++)
188 {
189 b[i] = Array<OneD, NekDouble>(bpts, 0.0);
190 }
191 pFields[0]->GetBndCondExpansions()[*it]->GetCoords(b[0], b[1],
192 b[2]);
193 for (int i = 0;
194 i < pFields[0]->GetBndCondExpansions()[*it]->GetNpoints(); ++i)
195 {
196 inPoints.push_back(
197 BPointPair(BPoint(b[0][i], b[1][i], b[2][i]), i));
198 }
199 }
201 m_rtree->insert(inPoints.begin(), inPoints.end());
202
203 for (int i = 0; i < npts; ++i)
204 {
205 std::vector<BPointPair> result;
206 BPoint sPoint(x[0][i], x[1][i], x[2][i]);
207 m_rtree->query(bgi::nearest(sPoint, 1), std::back_inserter(result));
208 r[i] = bg::distance(sPoint, result[0].first);
209 }
210 points.push_back(r);
211
212 std::string s_FieldStr;
213 for (int i = 0; i < m_NumVariable; ++i)
214 {
215 s_FieldStr = m_session->GetVariable(i);
216 ASSERTL0(m_session->DefinesFunction(funcName, s_FieldStr),
217 "Variable '" + s_FieldStr + "' not defined.");
218
220 m_session->GetFunction(funcName, s_FieldStr);
221 ASSERTL0(ffunc->GetVlist() == "x y z t r",
222 "EVARS of " + funcName + " must be 'r'");
223
224 ffunc->Evaluate(points, m_Absorption[i]);
225 }
226 }
227 else
228 {
229 for (int i = 0; i < m_NumVariable; ++i)
230 {
231 std::string s_FieldStr = m_session->GetVariable(i);
232 GetFunction(pFields, m_session, funcName)
233 ->Evaluate(s_FieldStr, m_Absorption[i]);
234 }
235 }
236
237 // If homogeneous expansion is used, transform the forcing term to
238 // be in the Fourier space
239 if (m_homogeneous)
240 {
241 for (int i = 0; i < m_NumVariable; ++i)
242 {
243 pFields[i]->HomogeneousFwdTrans(pFields[i]->GetTotPoints(),
245 }
246 }
247}
248
251 const Array<OneD, Array<OneD, NekDouble>> &inarray,
252 Array<OneD, Array<OneD, NekDouble>> &outarray, const NekDouble &time)
253{
254 int nq = m_Forcing[0].size();
255 CalculateForcing(fields, inarray, time);
256 for (int i = 0; i < m_NumVariable; ++i)
257 {
258 Vmath::Vadd(nq, m_Forcing[i], 1, outarray[i], 1, outarray[i], 1);
259 }
260}
261
264 const Array<OneD, Array<OneD, NekDouble>> &inarray,
265 Array<OneD, Array<OneD, NekDouble>> &outarray, const NekDouble &time)
266{
267 int ncoeff = outarray[m_NumVariable - 1].size();
268 Array<OneD, NekDouble> tmp(ncoeff, 0.0);
269 CalculateForcing(fields, inarray, time);
270
271 for (int i = 0; i < m_NumVariable; ++i)
272 {
273 fields[i]->FwdTrans(m_Forcing[i], tmp);
274 Vmath::Vadd(ncoeff, tmp, 1, outarray[i], 1, outarray[i], 1);
275 }
276}
277
279 [[maybe_unused]] const Array<OneD, MultiRegions::ExpListSharedPtr> &fields,
280 const Array<OneD, Array<OneD, NekDouble>> &inarray, const NekDouble &time)
281{
282 int nq = m_Forcing[0].size();
283
284 std::string s_FieldStr;
285 Array<OneD, NekDouble> TimeScale(1);
287
288 if (m_hasRefFlow)
289 {
290 for (int i = 0; i < m_NumVariable; i++)
291 {
292 RefflowScaled[i] = Array<OneD, NekDouble>(nq);
294 {
295 s_FieldStr = m_session->GetVariable(i);
296
297 std::string s_FieldStr = m_session->GetVariable(i);
299 ->Evaluate(s_FieldStr, m_Refflow[i], time);
300 Vmath::Vcopy(nq, m_Refflow[i], 1, RefflowScaled[i], 1);
301 }
302 else
303 {
304 Vmath::Vcopy(nq, m_Refflow[i], 1, RefflowScaled[i], 1);
305 }
306
307 Vmath::Vsub(nq, inarray[i], 1, RefflowScaled[i], 1, m_Forcing[i],
308 1);
309 Vmath::Vmul(nq, m_Absorption[i], 1, m_Forcing[i], 1, m_Forcing[i],
310 1);
311 }
312 }
313 else
314 {
315 for (int i = 0; i < m_NumVariable; i++)
316 {
317 Vmath::Vmul(nq, m_Absorption[i], 1, inarray[i], 1, m_Forcing[i], 1);
318 }
319 }
320}
321
322} // namespace Nektar::SolverUtils
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:208
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.
static bool GenerateVector(const std::string &str, std::vector< T > &out)
Takes a comma-separated string and converts it to entries in a vector.
Definition: ParseUtils.cpp:130
bg::model::point< NekDouble, 3, bg::cs::cartesian > BPoint
SOLVER_UTILS_EXPORT void v_InitObject(const Array< OneD, MultiRegions::ExpListSharedPtr > &pFields, const unsigned int &pNumForcingFields, const TiXmlElement *pForce) override
SOLVER_UTILS_EXPORT void v_ApplyCoeff(const Array< OneD, MultiRegions::ExpListSharedPtr > &fields, const Array< OneD, Array< OneD, NekDouble > > &inarray, Array< OneD, Array< OneD, NekDouble > > &outarray, const NekDouble &time) override
ForcingAbsorption(const LibUtilities::SessionReaderSharedPtr &pSession, const std::weak_ptr< EquationSystem > &pEquation)
static SOLVER_UTILS_EXPORT ForcingSharedPtr create(const LibUtilities::SessionReaderSharedPtr &pSession, const std::weak_ptr< EquationSystem > &pEquation, const Array< OneD, MultiRegions::ExpListSharedPtr > &pFields, const unsigned int &pNumForcingFields, const TiXmlElement *pForce)
Creates an instance of this class.
SOLVER_UTILS_EXPORT 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
Array< OneD, Array< OneD, NekDouble > > m_Refflow
std::pair< BPoint, unsigned int > BPointPair
void CalcAbsorption(const Array< OneD, MultiRegions::ExpListSharedPtr > &pFields, const TiXmlElement *pForce)
void CalculateForcing(const Array< OneD, MultiRegions::ExpListSharedPtr > &fields, const Array< OneD, Array< OneD, NekDouble > > &inarray, const NekDouble &time)
static std::string className
Name of the class.
std::vector< unsigned int > m_bRegions
Array< OneD, Array< OneD, NekDouble > > m_Absorption
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
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:150
LibUtilities::SessionReaderSharedPtr m_session
Session reader.
Definition: Forcing.h:123
const BoundaryRegionCollection & GetBoundaryRegions(void) const
Definition: Conditions.h:235
std::shared_ptr< SessionReader > SessionReaderSharedPtr
std::shared_ptr< Equation > EquationSharedPtr
Definition: Equation.h:125
ForcingFactory & GetForcingFactory()
Declaration of the forcing factory singleton.
Definition: Forcing.cpp:41
std::map< int, BoundaryRegionShPtr > BoundaryRegionCollection
Definition: Conditions.h:211
InputIterator find(InputIterator first, InputIterator last, InputIterator startingpoint, const EqualityComparable &value)
Definition: StdRegions.hpp:475
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 Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.hpp:825
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