Nektar++
Diffusion3DHomogeneous1D.cpp
Go to the documentation of this file.
1///////////////////////////////////////////////////////////////////////////////
2//
3// File: Diffusion3DHomogeneous1D.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: LDG diffusion 3DHomogeneous1D class.
32//
33///////////////////////////////////////////////////////////////////////////////
34
35#include <iomanip>
36#include <iostream>
37
39
40namespace Nektar::SolverUtils
41{
42std::string Diffusion3DHomogeneous1D::type[] = {
44 "LDG3DHomogeneous1D", Diffusion3DHomogeneous1D::create),
46 "LFRDG3DHomogeneous1D", Diffusion3DHomogeneous1D::create),
48 "LFRSD3DHomogeneous1D", Diffusion3DHomogeneous1D::create),
50 "LFRHU3DHomogeneous1D", Diffusion3DHomogeneous1D::create),
52 "LFRcmin3DHomogeneous1D", Diffusion3DHomogeneous1D::create),
54 "LFRcinf3DHomogeneous1D", Diffusion3DHomogeneous1D::create),
56 "LDGNS3DHomogeneous1D", Diffusion3DHomogeneous1D::create),
58 "LFRDGNS3DHomogeneous1D", Diffusion3DHomogeneous1D::create),
60 "LFRSDNS3DHomogeneous1D", Diffusion3DHomogeneous1D::create),
62 "LFRHUNS3DHomogeneous1D", Diffusion3DHomogeneous1D::create),
64 "LFRcminNS3DHomogeneous1D", Diffusion3DHomogeneous1D::create),
66 "LFRcinfNS3DHomogeneous1D", Diffusion3DHomogeneous1D::create)};
67
68/**
69 * @brief Diffusion3DHomogeneous1D uses the 2D WeakDG approach
70 * to compute the diffusion term looping on the planes in the z
71 * direction and adding the flux in z direction at the end.
72 */
74{
75 // Strip trailing string "3DHomogeneous1D" to determine 2D diffusion
76 // type, and create a diffusion object for the plane.
77 m_diffType = diffType.substr(0, diffType.length() - 15);
79}
80
81/**
82 * @brief Initiliase Diffusion3DHomogeneous1D objects and store
83 * them before starting the time-stepping.
84 *
85 * @param pSession Pointer to session reader.
86 * @param pFields Pointer to fields.
87 */
91{
92 int nConvectiveFields = pFields.size();
93
95 nConvectiveFields);
96
97 // Initialise the plane advection object.
98 for (int i = 0; i < nConvectiveFields; ++i)
99 {
100 pFields_plane0[i] = pFields[i]->GetPlane(0);
101 }
102 m_planeDiff->InitObject(pSession, pFields_plane0);
103
104 m_numPoints = pFields[0]->GetTotPoints();
105 m_numPlanes = pFields[0]->GetZIDs().size();
107 m_homoLen = pFields[0]->GetHomoLen();
108 m_trans = pFields[0]->GetTransposition();
109 m_planeCounter = 0;
110
111 if (m_diffType == "LDG")
112 {
113 // Set viscous flux for LDG
114 m_planeDiff->SetFluxVector(m_fluxVector);
115 }
116 else if (m_diffType == "LDGNS")
117 {
118 // Set viscous flux for LDGNS
119 m_planeDiff->SetFluxVectorNS(m_fluxVectorNS);
120 // Set penalty flux
121 m_planeDiff->SetFluxPenaltyNS(m_fluxPenaltyNS);
122 }
123 else if (m_diffType == "LFRDGNS" || m_diffType == "LFRHUNS" ||
124 m_diffType == "LFRSDNS")
125 {
126 // Set viscous flux for FR cases
127 m_planeDiff->SetFluxVectorNS(m_fluxVectorNS);
128 }
129
132
133 if (m_fluxVectorNS)
134 {
136 Array<OneD, Array<OneD, NekDouble>>(nConvectiveFields - 1);
137 }
138 else
139 {
141 }
144
145 for (int i = 0; i < m_numPlanes; ++i)
146 {
148 }
149
150 if (m_fluxVectorNS)
151 {
153 Array<OneD, Array<OneD, NekDouble>>(nConvectiveFields);
156
157 for (int i = 0; i < nConvectiveFields; ++i)
158 {
160 }
161
162 for (int i = 0; i < m_numPlanes; ++i)
163 {
165 Array<OneD, Array<OneD, NekDouble>>(nConvectiveFields);
166
167 for (int j = 0; j < nConvectiveFields; ++j)
168 {
171 }
172 }
173 }
174}
175
176/**
177 * @brief Calculate WeakDG Diffusion for the linear problems
178 * using an LDG interface flux and the the flux in the third direction.
179 */
181 const std::size_t nConvectiveFields,
183 const Array<OneD, Array<OneD, NekDouble>> &inarray,
185 [[maybe_unused]] const Array<OneD, Array<OneD, NekDouble>> &pFwd,
186 [[maybe_unused]] const Array<OneD, Array<OneD, NekDouble>> &pBwd)
187{
190 const int nPointsTot = fields[0]->GetNpoints();
192
193 if (m_fluxVectorNS)
194 {
195 viscHComp = Array<OneD, Array<OneD, NekDouble>>(nConvectiveFields);
196 for (int i = 0; i < nConvectiveFields - 1; ++i)
197 {
198 fields[0]->PhysDeriv(2, inarray[i], m_homoDerivStore[i]);
199 viscHComp[i] = Array<OneD, NekDouble>(m_numPoints);
200 }
201 }
202
203 for (int i = 0; i < m_numPlanes; ++i)
204 {
205 // Set up memory references for fields, inarray and outarray for
206 // this plane.
207 for (int j = 0; j < inarray.size(); ++j)
208 {
210 m_numPointsPlane, tmp2 = inarray[j] + m_planePos[i]);
211 }
212
213 for (int j = 0; j < nConvectiveFields; ++j)
214 {
215 m_fieldsPlane[j] = fields[j]->GetPlane(i);
217 m_numPointsPlane, tmp2 = outarray[j] + m_planePos[i]);
218 }
219
220 if (m_fluxVectorNS)
221 {
222 m_planeDiff->SetHomoDerivs(m_homoDerivPlane[i]);
223 }
224
225 if (m_diffType == "LDGNS")
226 {
227 // Store plane Fwd/Bwd traces
228 std::size_t nTracePts =
229 m_fieldsPlane[0]->GetTrace()->GetTotPoints();
230 std::size_t nScalar = m_inarrayPlane.size();
233 {
234 for (std::size_t k = 0; k < nScalar; ++k)
235 {
236 Fwd[k] = Array<OneD, NekDouble>(nTracePts, 0.0);
237 Bwd[k] = Array<OneD, NekDouble>(nTracePts, 0.0);
238 m_fieldsPlane[k]->GetFwdBwdTracePhys(m_inarrayPlane[k],
239 Fwd[k], Bwd[k]);
240 }
241 }
242
243 m_planeDiff->Diffuse(nConvectiveFields, m_fieldsPlane,
245 }
246 else
247 {
248 m_planeDiff->Diffuse(nConvectiveFields, m_fieldsPlane,
250 }
251
252 if (m_fluxVectorNS)
253 {
255 m_planeDiff->GetFluxTensor();
256
257 // Extract H (viscTensor[2])
258 for (int j = 0; j < nConvectiveFields - 1; ++j)
259 {
260 Vmath::Vcopy(m_numPointsPlane, viscTensor[2][j + 1], 1,
261 tmp2 = viscHComp[j] + m_planePos[i], 1);
262 }
263 }
264 }
265
266 if (m_fluxVectorNS)
267 {
268 for (int j = 0; j < nConvectiveFields - 1; ++j)
269 {
270 fields[j + 1]->PhysDeriv(2, viscHComp[j], tmp);
271 Vmath::Vadd(nPointsTot, outarray[j + 1], 1, tmp, 1, outarray[j + 1],
272 1);
273 }
274 }
275 else
276 {
277 for (int j = 0; j < nConvectiveFields; ++j)
278 {
279 fields[j]->HomogeneousFwdTrans(m_numPoints, inarray[j], tmp);
280
281 for (int i = 0; i < m_numPlanes; ++i)
282 {
283 beta = 2 * M_PI * m_trans->GetK(i) / m_homoLen;
284 beta *= beta;
285
287 &tmp[0] + i * m_numPointsPlane, 1,
288 &tmp[0] + i * m_numPointsPlane, 1);
289 }
290
291 fields[0]->HomogeneousBwdTrans(m_numPoints, tmp, tmp);
292
293 Vmath::Vsub(nPointsTot, outarray[j], 1, tmp, 1, outarray[j], 1);
294 }
295 }
296}
297} // namespace Nektar::SolverUtils
tKey RegisterCreatorFunction(tKey idKey, CreatorFunction classCreator, std::string pDesc="")
Register a class with the factory.
tBaseSharedPtr CreateInstance(tKey idKey, tParam... args)
Create an instance of the class referred to by idKey.
Diffusion3DHomogeneous1D(std::string diffType)
Diffusion3DHomogeneous1D uses the 2D WeakDG approach to compute the diffusion term looping on the pla...
LibUtilities::TranspositionSharedPtr m_trans
Array< OneD, MultiRegions::ExpListSharedPtr > m_fieldsPlane
static DiffusionSharedPtr create(std::string diffType)
void v_InitObject(LibUtilities::SessionReaderSharedPtr pSession, Array< OneD, MultiRegions::ExpListSharedPtr > pFields) override
Initiliase Diffusion3DHomogeneous1D objects and store them before starting the time-stepping.
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > m_homoDerivPlane
Array< OneD, Array< OneD, NekDouble > > m_outarrayPlane
Array< OneD, Array< OneD, NekDouble > > m_homoDerivStore
Array< OneD, Array< OneD, NekDouble > > m_inarrayPlane
void v_Diffuse(const size_t nConvective, 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=NullNekDoubleArrayOfArray, const Array< OneD, Array< OneD, NekDouble > > &pBwd=NullNekDoubleArrayOfArray) override
Calculate WeakDG Diffusion for the linear problems using an LDG interface flux and the the flux in th...
DiffusionFluxVecCBNS m_fluxVectorNS
Definition: Diffusion.h:361
DiffusionFluxPenaltyNS m_fluxPenaltyNS
Definition: Diffusion.h:362
DiffusionFluxVecCB m_fluxVector
Definition: Diffusion.h:360
std::shared_ptr< SessionReader > SessionReaderSharedPtr
@ beta
Gauss Radau pinned at x=-1,.
Definition: PointsType.h:59
DiffusionFactory & GetDiffusionFactory()
Definition: Diffusion.cpp:39
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 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
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