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 // License for the specific language governing rights and limitations under
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: FR advection 3DHomogeneous1D class.
33 //
34 ///////////////////////////////////////////////////////////////////////////////
35 
38 #include <iostream>
39 #include <iomanip>
40 
41 namespace Nektar
42 {
43  namespace SolverUtils
44  {
45  std::string Advection3DHomogeneous1D::type[] = {
47  "WeakDG3DHomogeneous1D", Advection3DHomogeneous1D::create),
49  "FRDG3DHomogeneous1D", Advection3DHomogeneous1D::create),
51  "FRDG3DHomogeneous1D", Advection3DHomogeneous1D::create),
53  "FRSD3DHomogeneous1D", Advection3DHomogeneous1D::create),
55  "FRHU3DHomogeneous1D", Advection3DHomogeneous1D::create),
57  "FRcmin3DHomogeneous1D", Advection3DHomogeneous1D::create),
59  "FRcinf3DHomogeneous1D", Advection3DHomogeneous1D::create)
60  };
61 
62  /**
63  * @brief AdvectionFR uses the Flux Reconstruction (FR) approach to
64  * compute the advection term. The implementation is only for segments,
65  * quadrilaterals and hexahedra at the moment.
66  *
67  * \todo Extension to triangles, tetrahedra and other shapes.
68  * (Long term objective)
69  */
71  : m_advType(advType)
72  {
73  // Strip trailing string "3DHomogeneous1D" to determine 2D advection
74  // type, and create an advection object for the plane.
75  string advName = advType.substr(0, advType.length()-15);
76  m_planeAdv = GetAdvectionFactory().CreateInstance(advName, advName);
77  }
78 
79  /**
80  * @brief Initiliase Advection3DHomogeneous1D objects and store them
81  * before starting the time-stepping.
82  *
83  * @param pSession Pointer to session reader.
84  * @param pFields Pointer to fields.
85  */
89  {
90  int nConvectiveFields = pFields.num_elements();
91 
93  nConvectiveFields);
94 
95  // Initialise the plane advection object.
96  for (int i = 0; i < nConvectiveFields; ++i)
97  {
98  pFields_plane0[i] = pFields[i]->GetPlane(0);
99  }
100  m_planeAdv->InitObject(pSession, pFields_plane0);
101 
102  m_numPoints = pFields[0]->GetTotPoints();
103  m_planes = pFields[0]->GetZIDs();
104  m_numPlanes = m_planes.num_elements();
106 
107  // Set Riemann solver and flux vector callback for this plane.
108  m_planeAdv->SetRiemannSolver(m_riemann);
109  m_planeAdv->SetFluxVector (
111  m_planeCounter = 0;
112 
113  // Override Riemann solver scalar and vector callbacks.
116  map<string, RSScalarFuncType> scalars = m_riemann->GetScalars();
117  map<string, RSVecFuncType> vectors = m_riemann->GetVectors();
118 
119  for (it1 = scalars.begin(); it1 != scalars.end(); ++it1)
120  {
121  boost::shared_ptr<HomoRSScalar> tmp = MemoryManager<HomoRSScalar>
122  ::AllocateSharedPtr(it1->second, m_numPlanes);
123  m_riemann->SetScalar(it1->first, &HomoRSScalar::Exec, tmp);
124  }
125 
126  for (it2 = vectors.begin(); it2 != vectors.end(); ++it2)
127  {
128  boost::shared_ptr<HomoRSVector> tmp = MemoryManager<HomoRSVector>
129  ::AllocateSharedPtr(it2->second, m_numPlanes, it2->first);
130  m_riemann->SetVector(it2->first, &HomoRSVector::Exec, tmp);
131  }
132 
134  nConvectiveFields);
135 
136  // Set up storage for flux vector.
137  for (int i = 0; i < nConvectiveFields; ++i)
138  {
140  for (int j = 0; j < 3; ++j)
141  {
143  }
144  }
145 
146  m_fluxVecPlane = Array<OneD, Array<OneD,
149  (nConvectiveFields);
151  (nConvectiveFields);
153  (nConvectiveFields);
156 
157  // Set up memory reference which links fluxVecPlane to fluxVecStore.
158  for (int i = 0; i < m_numPlanes; ++i)
159  {
160  m_planePos[i] = i * m_numPointsPlane;
161  m_fluxVecPlane[i] =
163  nConvectiveFields);
164 
165  for (int j = 0; j < nConvectiveFields; ++j)
166  {
167  m_fluxVecPlane[i][j] =
169  for (int k = 0; k < 3; ++k)
170  {
173  m_fluxVecStore[j][k] + m_planePos[i]);
174  }
175  }
176  }
177  }
178 
179  /**
180  * @brief Compute the advection operator for a given input @a inarray
181  * and put the result in @a outarray.
182  *
183  * @param nConvectiveFields Number of fields to advect.
184  * @param fields Pointer to fields.
185  * @param advVel Advection velocities.
186  * @param inarray Input which will be advected.
187  * @param outarray Computed advection.
188  */
190  const int nConvectiveFields,
192  const Array<OneD, Array<OneD, NekDouble> > &advVel,
193  const Array<OneD, Array<OneD, NekDouble> > &inarray,
194  Array<OneD, Array<OneD, NekDouble> > &outarray,
195  const NekDouble &time)
196  {
198  int nVel = advVel.num_elements();
199 
200  // Call solver's flux vector function to compute the flux vector on
201  // the entire domain.
202  m_fluxVector(inarray, m_fluxVecStore);
203 
204  // Loop over each plane.
205  for (int i = 0; i < m_numPlanes; ++i)
206  {
207  // Set up memory references for fields, inarray and outarray for
208  // this plane.
209  for (int j = 0; j < nConvectiveFields; ++j)
210  {
211  m_fieldsPlane [j] = fields[j]->GetPlane(i);
213  m_numPointsPlane, tmp2 = inarray [j] + m_planePos[i]);
215  m_numPointsPlane, tmp2 = outarray[j] + m_planePos[i]);
216  }
217 
218  for (int j = 0; j < nVel; ++j)
219  {
220  if (advVel[j].num_elements() != 0)
221  {
223  m_numPointsPlane, tmp2 = advVel[j] + m_planePos[i]);
224  }
225  }
226 
227  // Compute advection term for this plane.
228  m_planeAdv->Advect(nConvectiveFields, m_fieldsPlane,
230  m_outarrayPlane, time);
231  }
232 
233  // Calculate Fourier derivative and add to final result.
234  for (int i = 0; i < nConvectiveFields; ++i)
235  {
236  fields[0]->PhysDeriv(2, m_fluxVecStore[i][2], tmp);
237 
238  Vmath::Vadd(m_numPoints, outarray[i], 1, tmp, 1,
239  outarray[i], 1);
240  }
241  }
242 
244  const Array<OneD, Array<OneD, NekDouble> > &inarray,
245  Array<OneD, Array<OneD, Array<OneD, NekDouble> > > &outarray)
246  {
247  // Return section of flux vector for this plane.
248  outarray = m_fluxVecPlane[m_planeCounter];
249 
250  // Increment the plane counter.
252  }
253  }
254 }
tBaseSharedPtr CreateInstance(tKey idKey BOOST_PP_COMMA_IF(MAX_PARAM) BOOST_PP_ENUM_BINARY_PARAMS(MAX_PARAM, tParam, x))
Create an instance of the class referred to by idKey.
Definition: NekFactory.hpp:162
Array< OneD, Array< OneD, NekDouble > > m_advVelPlane
static boost::shared_ptr< DataType > AllocateSharedPtr()
Allocate a shared pointer from the memory pool.
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > m_fluxVecStore
const Array< OneD, const NekDouble > & Exec()
Array< OneD, Array< OneD, Array< OneD, Array< OneD, NekDouble > > > > m_fluxVecPlane
RiemannSolverSharedPtr m_riemann
Riemann solver for DG-type schemes.
Definition: Advection.h:134
boost::shared_ptr< SessionReader > SessionReaderSharedPtr
Definition: MeshPartition.h:50
virtual 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)
Compute the advection operator for a given input inarray and put the result in outarray.
void ModifiedFluxVector(const Array< OneD, Array< OneD, NekDouble > > &physfield, Array< OneD, Array< OneD, Array< OneD, NekDouble > > > &flux)
Advection3DHomogeneous1D(std::string advType)
AdvectionFR uses the Flux Reconstruction (FR) approach to compute the advection term. The implementation is only for segments, quadrilaterals and hexahedra at the moment.
AdvectionFactory & GetAdvectionFactory()
Gets the factory for initialising advection objects.
Definition: Advection.cpp:46
static AdvectionSharedPtr create(std::string advType)
double NekDouble
const Array< OneD, const Array< OneD, NekDouble > > & Exec()
StandardMatrixTag boost::call_traits< LhsDataType >::const_reference rhs typedef NekMatrix< LhsDataType, StandardMatrixTag >::iterator iterator
Array< OneD, Array< OneD, NekDouble > > m_outarrayPlane
Array< OneD, Array< OneD, NekDouble > > m_inarrayPlane
AdvectionFluxVecCB m_fluxVector
Callback function to the flux vector (set when advection is in conservative form).
Definition: Advection.h:132
Array< OneD, MultiRegions::ExpListSharedPtr > m_fieldsPlane
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:285
virtual void v_InitObject(LibUtilities::SessionReaderSharedPtr pSession, Array< OneD, MultiRegions::ExpListSharedPtr > pFields)
Initiliase Advection3DHomogeneous1D objects and store them before starting the time-stepping.
tKey RegisterCreatorFunction(tKey idKey, CreatorFunction classCreator, tDescription pDesc="")
Register a class with the factory.
Definition: NekFactory.hpp:215