Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
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 using namespace std;
42 
43 namespace Nektar
44 {
45  namespace SolverUtils
46  {
47  std::string Advection3DHomogeneous1D::type[] = {
49  "WeakDG3DHomogeneous1D", Advection3DHomogeneous1D::create),
51  "FRDG3DHomogeneous1D", Advection3DHomogeneous1D::create),
53  "FRDG3DHomogeneous1D", Advection3DHomogeneous1D::create),
55  "FRSD3DHomogeneous1D", Advection3DHomogeneous1D::create),
57  "FRHU3DHomogeneous1D", Advection3DHomogeneous1D::create),
59  "FRcmin3DHomogeneous1D", Advection3DHomogeneous1D::create),
61  "FRcinf3DHomogeneous1D", Advection3DHomogeneous1D::create)
62  };
63 
64  /**
65  * @brief AdvectionFR uses the Flux Reconstruction (FR) approach to
66  * compute the advection term. The implementation is only for segments,
67  * quadrilaterals and hexahedra at the moment.
68  *
69  * \todo Extension to triangles, tetrahedra and other shapes.
70  * (Long term objective)
71  */
72  Advection3DHomogeneous1D::Advection3DHomogeneous1D(std::string advType)
73  : m_advType(advType)
74  {
75  // Strip trailing string "3DHomogeneous1D" to determine 2D advection
76  // type, and create an advection object for the plane.
77  string advName = advType.substr(0, advType.length()-15);
78  m_planeAdv = GetAdvectionFactory().CreateInstance(advName, advName);
79  }
80 
81  /**
82  * @brief Initiliase Advection3DHomogeneous1D objects and store them
83  * 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.num_elements();
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_planeAdv->InitObject(pSession, pFields_plane0);
103 
104  m_numPoints = pFields[0]->GetTotPoints();
105  m_planes = pFields[0]->GetZIDs();
106  m_numPlanes = m_planes.num_elements();
108 
109  // Set Riemann solver and flux vector callback for this plane.
110  m_planeAdv->SetRiemannSolver(m_riemann);
111  m_planeAdv->SetFluxVector (
113  m_planeCounter = 0;
114 
115  // Override Riemann solver scalar and vector callbacks.
118  map<string, RSScalarFuncType> scalars = m_riemann->GetScalars();
119  map<string, RSVecFuncType> vectors = m_riemann->GetVectors();
120 
121  for (it1 = scalars.begin(); it1 != scalars.end(); ++it1)
122  {
123  boost::shared_ptr<HomoRSScalar> tmp = MemoryManager<HomoRSScalar>
124  ::AllocateSharedPtr(it1->second, m_numPlanes);
125  m_riemann->SetScalar(it1->first, &HomoRSScalar::Exec, tmp);
126  }
127 
128  for (it2 = vectors.begin(); it2 != vectors.end(); ++it2)
129  {
130  boost::shared_ptr<HomoRSVector> tmp = MemoryManager<HomoRSVector>
131  ::AllocateSharedPtr(it2->second, m_numPlanes, it2->first);
132  m_riemann->SetVector(it2->first, &HomoRSVector::Exec, tmp);
133  }
134 
136  nConvectiveFields);
137 
138  // Set up storage for flux vector.
139  for (int i = 0; i < nConvectiveFields; ++i)
140  {
142  for (int j = 0; j < 3; ++j)
143  {
145  }
146  }
147 
148  m_fluxVecPlane = Array<OneD, Array<OneD,
151  (nConvectiveFields);
153  (nConvectiveFields);
155  (nConvectiveFields);
158 
159  // Set up memory reference which links fluxVecPlane to fluxVecStore.
160  for (int i = 0; i < m_numPlanes; ++i)
161  {
162  m_planePos[i] = i * m_numPointsPlane;
163  m_fluxVecPlane[i] =
165  nConvectiveFields);
166 
167  for (int j = 0; j < nConvectiveFields; ++j)
168  {
169  m_fluxVecPlane[i][j] =
171  for (int k = 0; k < 3; ++k)
172  {
175  m_fluxVecStore[j][k] + m_planePos[i]);
176  }
177  }
178  }
179  }
180 
181  /**
182  * @brief Compute the advection operator for a given input @a inarray
183  * and put the result in @a outarray.
184  *
185  * @param nConvectiveFields Number of fields to advect.
186  * @param fields Pointer to fields.
187  * @param advVel Advection velocities.
188  * @param inarray Input which will be advected.
189  * @param outarray Computed advection.
190  */
192  const int nConvectiveFields,
194  const Array<OneD, Array<OneD, NekDouble> > &advVel,
195  const Array<OneD, Array<OneD, NekDouble> > &inarray,
196  Array<OneD, Array<OneD, NekDouble> > &outarray,
197  const NekDouble &time,
198  const Array<OneD, Array<OneD, NekDouble> > &pFwd,
199  const Array<OneD, Array<OneD, NekDouble> > &pBwd)
200  {
202  int nVel = advVel.num_elements();
203 
204  // Call solver's flux vector function to compute the flux vector on
205  // the entire domain.
206  m_fluxVector(inarray, m_fluxVecStore);
207 
208  // Loop over each plane.
209  for (int i = 0; i < m_numPlanes; ++i)
210  {
211  // Set up memory references for fields, inarray and outarray for
212  // this plane.
213  for (int j = 0; j < nConvectiveFields; ++j)
214  {
215  m_fieldsPlane [j] = fields[j]->GetPlane(i);
217  m_numPointsPlane, tmp2 = inarray [j] + m_planePos[i]);
219  m_numPointsPlane, tmp2 = outarray[j] + m_planePos[i]);
220  }
221 
222  for (int j = 0; j < nVel; ++j)
223  {
224  if (advVel[j].num_elements() != 0)
225  {
227  m_numPointsPlane, tmp2 = advVel[j] + m_planePos[i]);
228  }
229  }
230 
231  // Compute advection term for this plane.
232  m_planeAdv->Advect(nConvectiveFields, m_fieldsPlane,
234  m_outarrayPlane, time);
235  }
236 
237  // Calculate Fourier derivative and add to final result.
238  for (int i = 0; i < nConvectiveFields; ++i)
239  {
240  fields[0]->PhysDeriv(2, m_fluxVecStore[i][2], tmp);
241 
242  Vmath::Vadd(m_numPoints, outarray[i], 1, tmp, 1,
243  outarray[i], 1);
244  }
245  }
246 
248  const Array<OneD, Array<OneD, NekDouble> > &inarray,
249  Array<OneD, Array<OneD, Array<OneD, NekDouble> > > &outarray)
250  {
251  // Return section of flux vector for this plane.
252  outarray = m_fluxVecPlane[m_planeCounter];
253 
254  // Increment the plane counter.
256  }
257  }
258 }
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
STL namespace.
RiemannSolverSharedPtr m_riemann
Riemann solver for DG-type schemes.
Definition: Advection.h:138
boost::shared_ptr< SessionReader > SessionReaderSharedPtr
Definition: MeshPartition.h:51
void ModifiedFluxVector(const Array< OneD, Array< OneD, NekDouble > > &physfield, Array< OneD, Array< OneD, Array< OneD, NekDouble > > > &flux)
AdvectionFactory & GetAdvectionFactory()
Gets the factory for initialising advection objects.
Definition: Advection.cpp:46
double NekDouble
const Array< OneD, const Array< OneD, NekDouble > > & Exec()
StandardMatrixTag boost::call_traits< LhsDataType >::const_reference rhs typedef NekMatrix< LhsDataType, StandardMatrixTag >::iterator iterator
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, const Array< OneD, Array< OneD, NekDouble > > &pFwd=NullNekDoubleArrayofArray, const Array< OneD, Array< OneD, NekDouble > > &pBwd=NullNekDoubleArrayofArray)
Compute the advection operator for a given input inarray and put the result in outarray.
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:136
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:299
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