Nektar++
Advection.cpp
Go to the documentation of this file.
1///////////////////////////////////////////////////////////////////////////////
2//
3// File: Advection.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: Abstract base class for advection.
32//
33///////////////////////////////////////////////////////////////////////////////
34
36
38{
39
40/**
41 * @returns The advection factory.
42 */
44{
45 static AdvectionFactory instance;
46 return instance;
47}
48
49template <typename DataType, typename TypeNekBlkMatSharedPtr>
51 const int nConvectiveFields, const int nSpaceDim,
53 const Array<OneD, TypeNekBlkMatSharedPtr> &TracePntJacCons,
55 const Array<OneD, TypeNekBlkMatSharedPtr> &TracePntJacGrad,
56 const Array<OneD, Array<OneD, DataType>> &TracePntJacGradSign)
57{
58 MultiRegions::ExpListSharedPtr tracelist = pFields[0]->GetTrace();
59 std::shared_ptr<LocalRegions::ExpansionVector> traceExp =
60 tracelist->GetExp();
61 int ntotTrac = (*traceExp).size();
62 int nTracePnts = tracelist->GetTotPoints();
63 int nTracPnt, nTracCoef;
64
65 DNekMatSharedPtr tmp2Add;
66 Array<OneD, DataType> MatData0;
69
70 Array<OneD, DNekMatSharedPtr> TraceJacFwd(ntotTrac);
71 Array<OneD, DNekMatSharedPtr> TraceJacBwd(ntotTrac);
72
73 for (int nelmt = 0; nelmt < ntotTrac; nelmt++)
74 {
75 nTracCoef = (*traceExp)[nelmt]->GetNcoeffs();
76 nTracPnt = (*traceExp)[nelmt]->GetTotPoints();
77 TraceJacFwd[nelmt] =
78 MemoryManager<DNekMat>::AllocateSharedPtr(nTracCoef, nTracPnt, 0.0);
79 TraceJacBwd[nelmt] =
80 MemoryManager<DNekMat>::AllocateSharedPtr(nTracCoef, nTracPnt, 0.0);
81 }
82
83 std::shared_ptr<LocalRegions::ExpansionVector> fieldExp =
84 pFields[0]->GetExp();
85 int nTotElmt = (*fieldExp).size();
86 int nElmtPnt, nElmtCoef;
87
88 Array<OneD, DNekMatSharedPtr> ElmtJacQuad(nTotElmt);
89 Array<OneD, DNekMatSharedPtr> ElmtJacCoef(nTotElmt);
90
91 Array<OneD, NekDouble> SymmMatData;
92
93 for (int nelmt = 0; nelmt < nTotElmt; nelmt++)
94 {
95 nElmtCoef = (*fieldExp)[nelmt]->GetNcoeffs();
96 nElmtPnt = (*fieldExp)[nelmt]->GetTotPoints();
97
98 ElmtJacQuad[nelmt] =
99 MemoryManager<DNekMat>::AllocateSharedPtr(nElmtCoef, nElmtPnt, 0.0);
101 nElmtCoef, nElmtCoef, 0.0);
102 }
103
104 bool TracePntJacGradflag = true;
105
106 Array<OneD, Array<OneD, DataType>> TraceJacConsSign(2);
107 for (int i = 0; i < 2; i++)
108 {
109 TraceJacConsSign[i] = Array<OneD, DataType>(nTracePnts, 1.0);
110 }
111
112 if (0 == TracePntJacGrad.size())
113 {
114 TracePntJacGradflag = false;
115 }
116
117 for (int m = 0; m < nConvectiveFields; m++)
118 {
119 for (int n = 0; n < nConvectiveFields; n++)
120 {
121 // ElmtJacCons to set 0
122 for (int nelmt = 0; nelmt < nTotElmt; nelmt++)
123 {
124 (*ElmtJacCoef[nelmt]) = 0.0;
125 (*ElmtJacQuad[nelmt]) = 0.0;
126 }
127
128 if (TracePntJacGradflag)
129 {
130 // only BJac is stored here because BJac = - FJac
131 for (int ndir = 0; ndir < nSpaceDim; ndir++)
132 {
133 // ElmtJacGrad to set 0
134 for (int nelmt = 0; nelmt < nTotElmt; nelmt++)
135 {
136 (*ElmtJacQuad[nelmt]) = 0.0;
137 }
138 int ngrad = n * nSpaceDim + ndir;
139
140 CalcJacobTraceInteg(pFields, m, ngrad, TracePntJacGrad,
141 TracePntJacGradSign, TraceJacFwd,
142 TraceJacBwd);
143 pFields[0]->AddTraceJacToElmtJac(TraceJacFwd, TraceJacBwd,
144 ElmtJacQuad);
145 // cost can be lower down by doing 3 directions together
146 pFields[0]->AddRightIPTPhysDerivBase(ndir, ElmtJacQuad,
147 ElmtJacCoef);
148 }
149 for (int nelmt = 0; nelmt < nTotElmt; nelmt++)
150 {
151 nElmtCoef = (*fieldExp)[nelmt]->GetNcoeffs();
152 if (SymmMatData.size() < nElmtCoef)
153 {
154 SymmMatData =
155 Array<OneD, NekDouble>(nElmtCoef * nElmtCoef);
156 }
157 MatData1 = ElmtJacCoef[nelmt]->GetPtr();
158 for (int i = 0; i < nElmtCoef; i++)
159 {
160 Vmath::Vcopy(nElmtCoef, &MatData1[i * nElmtCoef], 1,
161 &SymmMatData[i], nElmtCoef);
162 }
163 Vmath::Vadd(nElmtCoef * nElmtCoef, SymmMatData, 1, MatData1,
164 1, MatData1, 1);
165 }
166 }
167
168 CalcJacobTraceInteg(pFields, m, n, TracePntJacCons,
169 TraceJacConsSign, TraceJacFwd, TraceJacBwd);
170
171 pFields[0]->AddTraceJacToElmtJac(TraceJacFwd, TraceJacBwd,
172 ElmtJacQuad);
173 // Memory can be save to explore reusing ElmtJacQuad as output
174 pFields[0]->AddRightIPTBaseMatrix(ElmtJacQuad, ElmtJacCoef);
175
176 // add ElmtJacCons to gmtxarray[m][n]
177 for (int nelmt = 0; nelmt < nTotElmt; nelmt++)
178 {
179 tmp2Add = ElmtJacCoef[nelmt];
180 MatData0 = gmtxarray[m][n]->GetBlock(nelmt, nelmt)->GetPtr();
181 MatData1 = tmp2Add->GetPtr();
182 for (int i = 0; i < MatData0.size(); i++)
183 {
184 MatData0[i] += DataType(MatData1[i]);
185 }
186 }
187 }
188 }
189}
190
192 const int nConvectiveFields, const int nSpaceDim,
194 const Array<OneD, DNekBlkMatSharedPtr> &TracePntJacCons,
196 const Array<OneD, DNekBlkMatSharedPtr> &TracePntJacGrad,
197 const Array<OneD, Array<OneD, NekDouble>> &TracePntJacGradSign);
199 const int nConvectiveFields, const int nSpaceDim,
201 const Array<OneD, SNekBlkMatSharedPtr> &TracePntJacCons,
203 const Array<OneD, SNekBlkMatSharedPtr> &TracePntJacGrad,
204 const Array<OneD, Array<OneD, NekSingle>> &TracePntJacGradSign);
205
206/**
207 *
208 */
209template <typename DataType, typename TypeNekBlkMatSharedPtr>
211 const Array<OneD, MultiRegions::ExpListSharedPtr> &pFields, const int m,
212 const int n, const Array<OneD, const TypeNekBlkMatSharedPtr> &PntJac,
213 const Array<OneD, const Array<OneD, DataType>> &PntJacSign,
216{
217 MultiRegions::ExpListSharedPtr tracelist = pFields[0]->GetTrace();
218 std::shared_ptr<LocalRegions::ExpansionVector> traceExp =
219 tracelist->GetExp();
220 int ntotTrac = (*traceExp).size();
221 int nTracPnt, noffset, pntoffset;
222
223 Array<OneD, int> tracepnts(ntotTrac);
224 Array<OneD, Array<OneD, NekDouble>> JacFwd(ntotTrac);
225 Array<OneD, Array<OneD, NekDouble>> JacBwd(ntotTrac);
226
227 for (int nelmt = 0; nelmt < ntotTrac; nelmt++)
228 {
229 nTracPnt = (*traceExp)[nelmt]->GetTotPoints();
230 tracepnts[nelmt] = nTracPnt;
231
232 JacFwd[nelmt] = Array<OneD, NekDouble>(nTracPnt, 0.0);
233 JacBwd[nelmt] = Array<OneD, NekDouble>(nTracPnt, 0.0);
234 }
235
236 DataType ftmp, btmp;
237 for (int nelmt = 0; nelmt < ntotTrac; nelmt++)
238 {
239 nTracPnt = tracepnts[nelmt];
240 noffset = tracelist->GetPhys_Offset(nelmt);
241 for (int npnt = 0; npnt < nTracPnt; npnt++)
242 {
243 pntoffset = noffset + npnt;
244 ftmp = (*PntJac[0]->GetBlock(pntoffset, pntoffset))(m, n);
245 JacFwd[nelmt][npnt] = NekDouble(PntJacSign[0][pntoffset] * ftmp);
246
247 btmp = (*PntJac[1]->GetBlock(pntoffset, pntoffset))(m, n);
248 JacBwd[nelmt][npnt] = NekDouble(PntJacSign[1][pntoffset] * btmp);
249 }
250 }
251 tracelist->GetDiagMatIpwrtBase(JacFwd, TraceJacFwd);
252 tracelist->GetDiagMatIpwrtBase(JacBwd, TraceJacBwd);
253}
254
255/**
256 * This function should be overridden in derived classes to initialise the
257 * specific advection data members. However, this base class function should
258 * be called as the first statement of the overridden function to ensure the
259 * base class is correctly initialised in order.
260 *
261 * @param pSession Session information.
262 * @param pFields Expansion lists for scalar fields.
263 */
267{
268 m_spaceDim = pFields[0]->GetCoordim(0);
269
270 if (pSession->DefinesSolverInfo("HOMOGENEOUS"))
271 {
272 std::string HomoStr = pSession->GetSolverInfo("HOMOGENEOUS");
273 if (HomoStr == "HOMOGENEOUS1D" || HomoStr == "Homogeneous1D" ||
274 HomoStr == "1D" || HomoStr == "Homo1D" ||
275 HomoStr == "HOMOGENEOUS2D" || HomoStr == "Homogeneous2D" ||
276 HomoStr == "2D" || HomoStr == "Homo2D")
277 {
278 m_spaceDim++;
279 }
280 else
281 {
283 "Only 1D homogeneous dimension supported.");
284 }
285 }
286}
287
288/**
289 *
290 */
292 [[maybe_unused]] const Array<OneD, Array<OneD, NekDouble>> &inarray,
293 [[maybe_unused]] const Array<OneD, MultiRegions::ExpListSharedPtr> &fields)
294{
296 "A baseflow is not appropriate for this advection type.");
297}
298
299} // namespace Nektar::SolverUtils
#define NEKERROR(type, msg)
Assert Level 0 – Fundamental assert which is used whether in FULLDEBUG, DEBUG or OPT compilation mode...
Definition: ErrorUtil.hpp:202
#define SOLVER_UTILS_EXPORT
size_type size() const
Returns the array's size.
Provides a generic Factory class.
static std::shared_ptr< DataType > AllocateSharedPtr(const Args &...args)
Allocate a shared pointer from the memory pool.
int m_spaceDim
Storage for space dimension. Used for homogeneous extension.
Definition: Advection.h:179
virtual SOLVER_UTILS_EXPORT void v_InitObject(LibUtilities::SessionReaderSharedPtr pSession, Array< OneD, MultiRegions::ExpListSharedPtr > pFields)
Initialises the advection object.
Definition: Advection.cpp:264
void CalcJacobTraceInteg(const Array< OneD, MultiRegions::ExpListSharedPtr > &pFields, const int m, const int n, const Array< OneD, const TypeNekBlkMatSharedPtr > &PntJac, const Array< OneD, const Array< OneD, DataType > > &PntJacSign, Array< OneD, DNekMatSharedPtr > &TraceJacFwd, Array< OneD, DNekMatSharedPtr > &TraceJacBwd)
Definition: Advection.cpp:210
virtual SOLVER_UTILS_EXPORT void v_SetBaseFlow(const Array< OneD, Array< OneD, NekDouble > > &inarray, const Array< OneD, MultiRegions::ExpListSharedPtr > &fields)
Overrides the base flow used during linearised advection.
Definition: Advection.cpp:291
SOLVER_UTILS_EXPORT void AddTraceJacToMat(const int nConvectiveFields, const int nSpaceDim, const Array< OneD, MultiRegions::ExpListSharedPtr > &pFields, const Array< OneD, TypeNekBlkMatSharedPtr > &TracePntJacCons, Array< OneD, Array< OneD, TypeNekBlkMatSharedPtr > > &gmtxarray, const Array< OneD, TypeNekBlkMatSharedPtr > &TracePntJacGrad, const Array< OneD, Array< OneD, DataType > > &TracePntJacGradSign)
std::shared_ptr< SessionReader > SessionReaderSharedPtr
std::shared_ptr< ExpList > ExpListSharedPtr
Shared pointer to an ExpList object.
AdvectionFactory & GetAdvectionFactory()
Gets the factory for initialising advection objects.
Definition: Advection.cpp:43
std::shared_ptr< DNekMat > DNekMatSharedPtr
Definition: NekTypeDefs.hpp:75
double NekDouble
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