Nektar++
Advection3DHomogeneous1D.cpp
Go to the documentation of this file.
1///////////////////////////////////////////////////////////////////////////////
2//
3// File: Advection3DHomogeneous1D.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: FR advection 3DHomogeneous1D class.
32//
33///////////////////////////////////////////////////////////////////////////////
34
37#include <iomanip>
38#include <iostream>
39
40namespace Nektar::SolverUtils
41{
42std::string Advection3DHomogeneous1D::type[] = {
44 "WeakDG3DHomogeneous1D", Advection3DHomogeneous1D::create),
46 "FRDG3DHomogeneous1D", Advection3DHomogeneous1D::create),
48 "FRDG3DHomogeneous1D", Advection3DHomogeneous1D::create),
50 "FRSD3DHomogeneous1D", Advection3DHomogeneous1D::create),
52 "FRHU3DHomogeneous1D", Advection3DHomogeneous1D::create),
54 "FRcmin3DHomogeneous1D", Advection3DHomogeneous1D::create),
56 "FRcinf3DHomogeneous1D", Advection3DHomogeneous1D::create)};
57
58/**
59 * @brief AdvectionFR uses the Flux Reconstruction (FR) approach to
60 * compute the advection term. The implementation is only for segments,
61 * quadrilaterals and hexahedra at the moment.
62 *
63 * \todo Extension to triangles, tetrahedra and other shapes.
64 * (Long term objective)
65 */
67 : m_advType(advType)
68{
69 // Strip trailing string "3DHomogeneous1D" to determine 2D advection
70 // type, and create an advection object for the plane.
71 std::string advName = advType.substr(0, advType.length() - 15);
72 m_planeAdv = GetAdvectionFactory().CreateInstance(advName, advName);
73}
74
75/**
76 * @brief Initiliase Advection3DHomogeneous1D objects and store them
77 * before starting the time-stepping.
78 *
79 * @param pSession Pointer to session reader.
80 * @param pFields Pointer to fields.
81 */
85{
86 int nConvectiveFields = pFields.size();
87
89 nConvectiveFields);
90
91 // Initialise the plane advection object.
92 for (int i = 0; i < nConvectiveFields; ++i)
93 {
94 pFields_plane0[i] = pFields[i]->GetPlane(0);
95 }
96 m_planeAdv->InitObject(pSession, pFields_plane0);
97
98 m_numPoints = pFields[0]->GetTotPoints();
99 m_planes = pFields[0]->GetZIDs();
100 m_numPlanes = m_planes.size();
102
103 // Set Riemann solver and flux vector callback for this plane.
104 m_planeAdv->SetRiemannSolver(m_riemann);
106 this);
107 m_planeCounter = 0;
108
109 // Override Riemann solver scalar and vector callbacks.
110 std::map<std::string, RSScalarFuncType> scalars = m_riemann->GetScalars();
111 std::map<std::string, RSVecFuncType> vectors = m_riemann->GetVectors();
112
113 for (auto &it1 : scalars)
114 {
115 std::shared_ptr<HomoRSScalar> tmp =
118 m_riemann->SetScalar(it1.first, &HomoRSScalar::Exec, tmp);
119 }
120
121 for (auto &it2 : vectors)
122 {
123 std::shared_ptr<HomoRSVector> tmp =
125 it2.second, m_numPlanes, it2.first);
126 m_riemann->SetVector(it2.first, &HomoRSVector::Exec, tmp);
127 }
128
131
132 // Set up storage for flux vector.
133 for (int i = 0; i < nConvectiveFields; ++i)
134 {
136 for (int j = 0; j < 3; ++j)
137 {
139 }
140 }
141
151
152 // Set up memory reference which links fluxVecPlane to fluxVecStore.
153 for (int i = 0; i < m_numPlanes; ++i)
154 {
156 m_fluxVecPlane[i] =
158
159 for (int j = 0; j < nConvectiveFields; ++j)
160 {
162 for (int k = 0; k < 3; ++k)
163 {
166 }
167 }
168 }
169}
170
171/**
172 * @brief Compute the advection operator for a given input @a inarray
173 * and put the result in @a outarray.
174 *
175 * @param nConvectiveFields Number of fields to advect.
176 * @param fields Pointer to fields.
177 * @param advVel Advection velocities.
178 * @param inarray Input which will be advected.
179 * @param outarray Computed advection.
180 */
182 const int nConvectiveFields,
184 const Array<OneD, Array<OneD, NekDouble>> &advVel,
185 const Array<OneD, Array<OneD, NekDouble>> &inarray,
186 Array<OneD, Array<OneD, NekDouble>> &outarray, const NekDouble &time,
187 [[maybe_unused]] const Array<OneD, Array<OneD, NekDouble>> &pFwd,
188 [[maybe_unused]] const Array<OneD, Array<OneD, NekDouble>> &pBwd)
189{
191 int nVel = advVel.size();
192
193 // Call solver's flux vector function to compute the flux vector on
194 // the entire domain.
196
197 // Loop over each plane.
198 for (int i = 0; i < m_numPlanes; ++i)
199 {
200 // Set up memory references for fields, inarray and outarray for
201 // this plane.
202 for (int j = 0; j < nConvectiveFields; ++j)
203 {
204 m_fieldsPlane[j] = fields[j]->GetPlane(i);
206 m_numPointsPlane, tmp2 = inarray[j] + m_planePos[i]);
208 m_numPointsPlane, tmp2 = outarray[j] + m_planePos[i]);
209 }
210
211 for (int j = 0; j < nVel; ++j)
212 {
213 if (advVel[j].size() != 0)
214 {
216 m_numPointsPlane, tmp2 = advVel[j] + m_planePos[i]);
217 }
218 }
219
220 // Compute advection term for this plane.
221 m_planeAdv->Advect(nConvectiveFields, m_fieldsPlane, m_advVelPlane,
223 }
224
225 // Calculate Fourier derivative and add to final result.
226 for (int i = 0; i < nConvectiveFields; ++i)
227 {
228 fields[0]->PhysDeriv(2, m_fluxVecStore[i][2], tmp);
229
230 Vmath::Vadd(m_numPoints, outarray[i], 1, tmp, 1, outarray[i], 1);
231 }
232}
233
234/**
235 *
236 */
238 [[maybe_unused]] const Array<OneD, Array<OneD, NekDouble>> &inarray,
240{
241 // Return section of flux vector for this plane.
242 outarray = m_fluxVecPlane[m_planeCounter];
243
244 // Increment the plane counter.
246}
247} // 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.
static std::shared_ptr< DataType > AllocateSharedPtr(const Args &...args)
Allocate a shared pointer from the memory pool.
Array< OneD, Array< OneD, Array< OneD, Array< OneD, NekDouble > > > > m_fluxVecPlane
void ModifiedFluxVector(const Array< OneD, Array< OneD, NekDouble > > &physfield, Array< OneD, Array< OneD, Array< OneD, NekDouble > > > &flux)
Array< OneD, Array< OneD, NekDouble > > m_outarrayPlane
void v_Advect(const int nConvField, const Array< OneD, MultiRegions::ExpListSharedPtr > &fields, const Array< OneD, Array< OneD, NekDouble > > &advVel, const Array< OneD, Array< OneD, NekDouble > > &inarray, Array< OneD, Array< OneD, NekDouble > > &outarray, const NekDouble &time, const Array< OneD, Array< OneD, NekDouble > > &pFwd=NullNekDoubleArrayOfArray, const Array< OneD, Array< OneD, NekDouble > > &pBwd=NullNekDoubleArrayOfArray) override
Compute the advection operator for a given input inarray and put the result in outarray.
Array< OneD, Array< OneD, NekDouble > > m_advVelPlane
Advection3DHomogeneous1D(std::string advType)
AdvectionFR uses the Flux Reconstruction (FR) approach to compute the advection term....
Array< OneD, Array< OneD, NekDouble > > m_inarrayPlane
static AdvectionSharedPtr create(std::string advType)
Array< OneD, MultiRegions::ExpListSharedPtr > m_fieldsPlane
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > m_fluxVecStore
void v_InitObject(LibUtilities::SessionReaderSharedPtr pSession, Array< OneD, MultiRegions::ExpListSharedPtr > pFields) override
Initiliase Advection3DHomogeneous1D objects and store them before starting the time-stepping.
AdvectionFluxVecCB m_fluxVector
Callback function to the flux vector (set when advection is in conservative form).
Definition: Advection.h:175
RiemannSolverSharedPtr m_riemann
Riemann solver for DG-type schemes.
Definition: Advection.h:177
const Array< OneD, const NekDouble > & Exec()
const Array< OneD, const Array< OneD, NekDouble > > & Exec()
std::shared_ptr< SessionReader > SessionReaderSharedPtr
AdvectionFactory & GetAdvectionFactory()
Gets the factory for initialising advection objects.
Definition: Advection.cpp:43
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