Nektar++
APE.cpp
Go to the documentation of this file.
1 ///////////////////////////////////////////////////////////////////////////////
2 //
3 // File APE.cpp
4 //
5 // For more information, please see: http://www.nektar.info
6 //
7 // The MIT License
8 //
9 // Copyright (c) 2015 Kilian Lackhove
10 // Copyright (c) 2006 Division of Applied Mathematics, Brown University (USA),
11 // Department of Aeronautics, Imperial College London (UK), and Scientific
12 // Computing and Imaging Institute, University of Utah (USA).
13 //
14 // License for the specific language governing rights and limitations under
15 // Permission is hereby granted, free of charge, to any person obtaining a
16 // copy of this software and associated documentation files (the "Software"),
17 // to deal in the Software without restriction, including without limitation
18 // the rights to use, copy, modify, merge, publish, distribute, sublicense,
19 // and/or sell copies of the Software, and to permit persons to whom the
20 // Software is furnished to do so, subject to the following conditions:
21 //
22 // The above copyright notice and this permission notice shall be included
23 // in all copies or substantial portions of the Software.
24 //
25 // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
26 // OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
27 // FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
28 // THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
29 // LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
30 // FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
31 // DEALINGS IN THE SOFTWARE.
32 //
33 // Description: APE1/APE4 (Acoustic Perturbation Equations)
34 //
35 ///////////////////////////////////////////////////////////////////////////////
36 
37 #include <iostream>
38 
41 
42 namespace Nektar
43 {
45  "APE", APE::create,
46  "APE1/APE4 (Acoustic Perturbation Equations)");
47 
48 
51  : UnsteadySystem(pSession)
52 {
53 }
54 
55 
56 /**
57  * @brief Initialization object for the APE class.
58  */
60 {
61  UnsteadySystem::v_InitObject();
62 
63  // TODO: We have a bug somewhere in the 1D boundary conditions. Therefore 1D
64  // problems are currently disabled. This should get fixed in the future.
65  ASSERTL0(m_spacedim > 1, "1D problems currently not supported by the APE class.");
66 
68  "Only Projection=DisContinuous supported by the APE class.");
69 
70  // Load isentropic coefficient, Ratio of specific heats
71  m_session->LoadParameter("Gamma", m_gamma, 1.4);
72 
73  // Define Baseflow fields
75  m_basefield_names.push_back("p0");
76  m_basefield_names.push_back("rho0");
77  m_basefield_names.push_back("u0");
78  m_basefield_names.push_back("v0");
79  m_basefield_names.push_back("w0");
80 
81  // Resize the advection velocities vector to dimension of the problem
82  m_basefield_names.resize(m_spacedim + 2);
83 
84  // Initialize the sourceterm
86 
87  // Do not forwards transform initial condition
88  m_homoInitialFwd = false;
89 
90  // Define the normal velocity fields
91  if (m_fields[0]->GetTrace())
92  {
94  for (int i = 0; i < m_spacedim + 2; i++)
95  {
97  }
98  }
99 
100  // Set up locations of velocity and base velocity vectors.
103  for (int i = 0; i < m_spacedim; ++i)
104  {
105  // u', v', w'
106  m_vecLocs[0][i] = 1 + i;
107  }
108 
109  string riemName;
110  m_session->LoadSolverInfo("UpwindType", riemName, "APEUpwind");
112  riemName);
113  m_riemannSolver->SetVector("N", &APE::GetNormals, this);
114  m_riemannSolver->SetVector("basefield", &APE::GetBasefield, this);
115  m_riemannSolver->SetAuxVec("vecLocs", &APE::GetVecLocs, this);
116  m_riemannSolver->SetParam("Gamma", &APE::GetGamma, this);
117 
118  // Set up advection operator
119  string advName;
120  m_session->LoadSolverInfo("AdvectionType", advName, "WeakDG");
122  .CreateInstance(advName, advName);
123  m_advection->SetFluxVector(&APE::GetFluxVector, this);
124  m_advection->SetRiemannSolver(m_riemannSolver);
125  m_advection->InitObject(m_session, m_fields);
126 
128  {
131  }
132  else
133  {
134  ASSERTL0(false, "Implicit APE not set up.");
135  }
136 }
137 
138 
139 /**
140  * @brief Destructor for APE class.
141  */
143 {
144 
145 }
146 
147 
148 /**
149  *
150  */
152 {
154 }
155 
156 
157 /**
158  * @brief Return the flux vector for the APE equations.
159  *
160  * @param physfield Fields.
161  * @param flux Resulting flux. flux[eq][dir][pt]
162  */
164  const Array<OneD, Array<OneD, NekDouble> > &physfield,
166 {
167  UpdateBasefield();
168 
169  int nq = physfield[0].num_elements();
170  Array<OneD, NekDouble> tmp1(nq);
171  Array<OneD, NekDouble> tmp2(nq);
172 
173  ASSERTL1(flux[0].num_elements() == m_spacedim,
174  "Dimension of flux array and velocity array do not match");
175 
176  // F_{adv,p',j} = \rho_0 u'_j + p' \bar{u}_j / c^2
177  for (int j = 0; j < m_spacedim; ++j)
178  {
179  Vmath::Zero(nq, flux[0][j], 1);
180 
181  // construct rho_0 u'_j term
182  Vmath::Vmul(nq, m_basefield[1], 1, physfield[j + 1], 1, flux[0][j], 1);
183 
184  // construct p' \bar{u}_j / c^2 term
185  // c^2
186  Vmath::Vdiv(nq, m_basefield[0], 1, m_basefield[1], 1, tmp1, 1);
187  Vmath::Smul(nq, m_gamma, tmp1, 1, tmp1, 1);
188 
189  // p' \bar{u}_j / c^2 term
190  Vmath::Vmul(nq, physfield[0], 1, m_basefield[j + 2], 1, tmp2, 1);
191  Vmath::Vdiv(nq, tmp2, 1, tmp1, 1, tmp2, 1);
192 
193  // \rho_0 u'_j + p' \bar{u}_j / c^2
194  Vmath::Vadd(nq, flux[0][j], 1, tmp2, 1, flux[0][j], 1);
195  }
196 
197  for (int i = 1; i < flux.num_elements(); ++i)
198  {
199  ASSERTL1(flux[i].num_elements() == m_spacedim,
200  "Dimension of flux array and velocity array do not match");
201 
202  // F_{adv,u'_i,j} = (p'/ \bar{rho} + \bar{u}_k u'_k) \delta_{ij}
203  for (int j = 0; j < m_spacedim; ++j)
204  {
205  Vmath::Zero(nq, flux[i][j], 1);
206 
207  if (i - 1 == j)
208  {
209  // contruct p'/ \bar{rho} term
210  Vmath::Vdiv(nq, physfield[0], 1, m_basefield[1], 1, flux[i][j], 1);
211 
212  // construct \bar{u}_k u'_k term
213  Vmath::Zero(nq, tmp1, 1);
214  for (int k = 0; k < m_spacedim; ++k)
215  {
216  Vmath::Vvtvp(nq, physfield[k + 1], 1, m_basefield[k + 2], 1, tmp1, 1, tmp1, 1);
217  }
218 
219  // add terms
220  Vmath::Vadd(nq, flux[i][j], 1, tmp1, 1, flux[i][j], 1);
221  }
222  }
223  }
224 }
225 
226 
227 /**
228  * @brief Compute the right-hand side.
229  */
230 void APE::DoOdeRhs(const Array<OneD, const Array<OneD, NekDouble> >&inarray,
231  Array<OneD, Array<OneD, NekDouble> >&outarray,
232  const NekDouble time)
233 {
234  int nVariables = inarray.num_elements();
235  int nq = GetTotPoints();
236  Array<OneD, NekDouble> tmp1(nq);
237 
238  // WeakDG does not use advVel, so we only provide a dummy array
240  m_advection->Advect(nVariables, m_fields, advVel, inarray, outarray, time);
241 
242  for (int i = 0; i < nVariables; ++i)
243  {
244  if (i == 0)
245  {
246  // c^2 = gamma*p0/rho0
247  Vmath::Vdiv(nq, m_basefield[0], 1, m_basefield[1], 1, tmp1, 1);
248  Vmath::Smul(nq, m_gamma, tmp1, 1, tmp1, 1);
249  Vmath::Vmul(nq, tmp1, 1, outarray[i], 1, outarray[i], 1);
250  }
251 
252  Vmath::Neg(nq, outarray[i], 1);
253  }
254 
255  AddSource(outarray);
256 }
257 
258 
259 /**
260  * @brief Compute the projection and call the method for imposing the
261  * boundary conditions in case of discontinuous projection.
262  */
264  Array<OneD, Array<OneD, NekDouble> >&outarray,
265  const NekDouble time)
266 {
267  int nvariables = inarray.num_elements();
268  int nq = m_fields[0]->GetNpoints();
269 
270  // deep copy
271  for (int i = 0; i < nvariables; ++i)
272  {
273  Vmath::Vcopy(nq, inarray[i], 1, outarray[i], 1);
274  }
275 
276  SetBoundaryConditions(outarray, time);
277 }
278 
279 
280 /**
281  * @brief Apply the Boundary Conditions to the APE equations.
282  */
284  NekDouble time)
285 {
286  std::string varName;
287  int nvariables = m_fields.num_elements();
288  int cnt = 0;
289 
290  // loop over Boundary Regions
291  for(int n = 0; n < m_fields[0]->GetBndConditions().num_elements(); ++n)
292  {
293  // Wall Boundary Condition
294  if (boost::iequals(m_fields[0]->GetBndConditions()[n]->GetUserDefined(),"Wall"))
295  {
296  WallBC(n, cnt, inarray);
297  }
298 
299  // Time Dependent Boundary Condition (specified in meshfile)
300  if (m_fields[0]->GetBndConditions()[n]->IsTimeDependent())
301  {
302  for (int i = 0; i < nvariables; ++i)
303  {
304  varName = m_session->GetVariable(i);
305  m_fields[i]->EvaluateBoundaryConditions(time, varName);
306  }
307  }
308  cnt +=m_fields[0]->GetBndCondExpansions()[n]->GetExpSize();
309  }
310 }
311 
312 
313 /**
314  * @brief Wall boundary conditions for the APE equations.
315  */
316 void APE::WallBC(int bcRegion, int cnt,
317  Array<OneD, Array<OneD, NekDouble> > &physarray)
318 {
319  int nTracePts = GetTraceTotPoints();
320  int nVariables = physarray.num_elements();
321 
322  const Array<OneD, const int> &traceBndMap = m_fields[0]->GetTraceBndMap();
323 
324  // Get physical values of the forward trace
325  Array<OneD, Array<OneD, NekDouble> > Fwd(nVariables);
326  for (int i = 0; i < nVariables; ++i)
327  {
328  Fwd[i] = Array<OneD, NekDouble>(nTracePts);
329  m_fields[i]->ExtractTracePhys(physarray[i], Fwd[i]);
330  }
331 
332  // Adjust the physical values of the trace to take
333  // user defined boundaries into account
334  int id1, id2, nBCEdgePts;
335  int eMax = m_fields[0]->GetBndCondExpansions()[bcRegion]->GetExpSize();
336 
337  for (int e = 0; e < eMax; ++e)
338  {
339  nBCEdgePts = m_fields[0]->GetBndCondExpansions()[bcRegion]->
340  GetExp(e)->GetTotPoints();
341  id1 = m_fields[0]->GetBndCondExpansions()[bcRegion]->GetPhys_Offset(e);
342  id2 = m_fields[0]->GetTrace()->GetPhys_Offset(traceBndMap[cnt+e]);
343 
344  // For 2D/3D, define: v* = v - 2(v.n)n
345  Array<OneD, NekDouble> tmp(nBCEdgePts, 0.0);
346 
347  // Calculate (v.n)
348  for (int i = 0; i < m_spacedim; ++i)
349  {
350  Vmath::Vvtvp(nBCEdgePts,
351  &Fwd[1+i][id2], 1,
352  &m_traceNormals[i][id2], 1,
353  &tmp[0], 1,
354  &tmp[0], 1);
355  }
356 
357  // Calculate 2.0(v.n)
358  Vmath::Smul(nBCEdgePts, -2.0, &tmp[0], 1, &tmp[0], 1);
359 
360  // Calculate v* = v - 2.0(v.n)n
361  for (int i = 0; i < m_spacedim; ++i)
362  {
363  Vmath::Vvtvp(nBCEdgePts,
364  &tmp[0], 1,
365  &m_traceNormals[i][id2], 1,
366  &Fwd[1+i][id2], 1,
367  &Fwd[1+i][id2], 1);
368  }
369 
370  // Copy boundary adjusted values into the boundary expansion
371  for (int i = 0; i < nVariables; ++i)
372  {
373  Vmath::Vcopy(nBCEdgePts,
374  &Fwd[i][id2], 1,
375  &(m_fields[i]->GetBndCondExpansions()[bcRegion]->UpdatePhys())[id1], 1);
376  }
377  }
378 }
379 
380 
381 /**
382  * @brief sourceterm for p' equation obtained from GetSource
383  */
385 {
387  Vmath::Vadd(GetTotPoints(), m_sourceTerms, 1, outarray[0], 1, outarray[0], 1);
388 }
389 
390 
392  std::vector<Array<OneD, NekDouble> > &fieldcoeffs,
393  std::vector<std::string> &variables)
394 {
395  UpdateBasefield();
396 
397  const int nCoeffs = m_fields[0]->GetNcoeffs();
398 
399  for (int i = 0; i < m_spacedim + 2; i++)
400  {
401  variables.push_back(m_basefield_names[i]);
402 
403  Array<OneD, NekDouble> tmpFwd(nCoeffs);
404  m_fields[0]->FwdTrans(m_basefield[i], tmpFwd);
405  fieldcoeffs.push_back(tmpFwd);
406  }
407 }
408 
409 
410 /**
411  * @brief Get the normal vectors.
412  */
414 {
415  return m_traceNormals;
416 }
417 
418 
419 /**
420  * @brief Get the locations of the components of the directed fields within the fields array.
421  */
423 {
424  return m_vecLocs;
425 }
426 
427 
428 /**
429  * @brief Get the baseflow field.
430  */
432 {
433  for (int i = 0; i < m_spacedim + 2; i++)
434  {
435  m_fields[0]->ExtractTracePhys(m_basefield[i], m_traceBasefield[i]);
436  }
437  return m_traceBasefield;
438 }
439 
440 
441 /**
442  * @brief Get the heat capacity ratio.
443  */
445 {
446  return m_gamma;
447 }
448 
449 
451 {
452  static NekDouble last_update = -1.0;
453 
454  if (m_time > last_update)
455  {
457  last_update = m_time;
458  }
459 }
460 
462 {
463  static NekDouble last_update = -1.0;
464 
465  if (m_time > last_update)
466  {
468 
469  EvaluateFunction("S", m_sourceTerms, "Source", m_time);
470 
471  m_fields[0]->IProductWRTBase(m_sourceTerms, sourceC);
472  m_fields[0]->MultiplyByElmtInvMass(sourceC, sourceC);
473  m_fields[0]->BwdTrans(sourceC, m_sourceTerms);
474 
475  last_update = m_time;
476  }
477 }
478 
479 
480 } //end of namespace
481 
const Array< OneD, const Array< OneD, NekDouble > > & GetNormals()
Get the normal vectors.
Definition: APE.cpp:413
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:161
SolverUtils::AdvectionSharedPtr m_advection
Definition: APE.h:72
void AddSource(Array< OneD, Array< OneD, NekDouble > > &outarray)
sourceterm for p' equation obtained from GetSource
Definition: APE.cpp:384
bool m_homoInitialFwd
Flag to determine if simulation should start in homogeneous forward transformed state.
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
NekDouble m_time
Current time of simulation.
std::vector< std::string > m_basefield_names
Definition: APE.h:80
LibUtilities::TimeIntegrationSchemeOperators m_ode
The time integration scheme operators to use.
Array< OneD, Array< OneD, NekDouble > > m_basefield
Definition: APE.h:78
Array< OneD, Array< OneD, NekDouble > > m_vecLocs
Definition: APE.h:75
static EquationSystemSharedPtr create(const LibUtilities::SessionReaderSharedPtr &pSession)
Creates an instance of this class.
Definition: APE.h:56
void Vvtvp(int n, const T *w, const int incw, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
vvtvp (vector times vector plus vector): z = w*x + y
Definition: Vmath.cpp:428
enum MultiRegions::ProjectionType m_projectionType
Type of projection; e.g continuous or discontinuous.
void Vdiv(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Multiply vector z = x/y.
Definition: Vmath.cpp:227
void GetFluxVector(const Array< OneD, Array< OneD, NekDouble > > &physfield, Array< OneD, Array< OneD, Array< OneD, NekDouble > > > &flux)
Return the flux vector for the APE equations.
Definition: APE.cpp:163
boost::shared_ptr< SessionReader > SessionReaderSharedPtr
Definition: MeshPartition.h:50
SOLVER_UTILS_EXPORT int GetTotPoints()
void UpdateBasefield()
Definition: APE.cpp:450
Array< OneD, NekDouble > m_sourceTerms
Definition: APE.h:79
Array< OneD, Array< OneD, NekDouble > > m_traceNormals
Array holding trace normals for DG simulations in the forwards direction.
void DefineProjection(FuncPointerT func, ObjectPointerT obj)
NekDouble GetGamma()
Get the heat capacity ratio.
Definition: APE.cpp:444
virtual void v_ExtraFldOutput(std::vector< Array< OneD, NekDouble > > &fieldcoeffs, std::vector< std::string > &variables)
Definition: APE.cpp:391
void Smul(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Scalar multiply y = alpha*y.
Definition: Vmath.cpp:199
bool m_explicitAdvection
Indicates if explicit or implicit treatment of advection is used.
const Array< OneD, const Array< OneD, NekDouble > > & GetVecLocs()
Get the locations of the components of the directed fields within the fields array.
Definition: APE.cpp:422
Array< OneD, Array< OneD, NekDouble > > m_traceBasefield
Definition: APE.h:74
void DefineOdeRhs(FuncPointerT func, ObjectPointerT obj)
APE(const LibUtilities::SessionReaderSharedPtr &pSession)
Initialises UnsteadySystem class members.
Definition: APE.cpp:49
virtual ~APE()
Destructor.
Definition: APE.cpp:142
Base class for unsteady solvers.
virtual void v_DoInitialise()
Sets up initial conditions.
Definition: APE.cpp:151
void DoOdeRhs(const Array< OneD, const Array< OneD, NekDouble > > &inarray, Array< OneD, Array< OneD, NekDouble > > &outarray, const NekDouble time)
Compute the right-hand side.
Definition: APE.cpp:230
RiemannSolverFactory & GetRiemannSolverFactory()
int m_spacedim
Spatial dimension (>= expansion dim).
AdvectionFactory & GetAdvectionFactory()
Gets the factory for initialising advection objects.
Definition: Advection.cpp:46
void Neg(int n, T *x, const int incx)
Negate x = -x.
Definition: Vmath.cpp:382
static std::string className
Name of class.
Definition: APE.h:64
double NekDouble
virtual void v_InitObject()
Initialization object for the APE class.
Definition: APE.cpp:59
void WallBC(int bcRegion, int cnt, Array< OneD, Array< OneD, NekDouble > > &physarray)
Wall boundary conditions for the APE equations.
Definition: APE.cpp:316
void SetBoundaryConditions(Array< OneD, Array< OneD, NekDouble > > &physarray, NekDouble time)
Apply the Boundary Conditions to the APE equations.
Definition: APE.cpp:283
void UpdateSourceTerms()
Definition: APE.cpp:461
SOLVER_UTILS_EXPORT void EvaluateFunction(Array< OneD, Array< OneD, NekDouble > > &pArray, std::string pFunctionName, const NekDouble pTime=0.0, const int domain=0)
Evaluates a function as specified in the session file.
NekDouble m_gamma
Isentropic coefficient, Ratio of specific heats (APE)
Definition: APE.h:77
const Array< OneD, const Array< OneD, NekDouble > > & GetBasefield()
Get the baseflow field.
Definition: APE.cpp:431
EquationSystemFactory & GetEquationSystemFactory()
SOLVER_UTILS_EXPORT int GetTraceTotPoints()
SolverUtils::RiemannSolverSharedPtr m_riemannSolver
Definition: APE.h:73
Array< OneD, MultiRegions::ExpListSharedPtr > m_fields
Array holding all dependent variables.
LibUtilities::SessionReaderSharedPtr m_session
The session reader.
SOLVER_UTILS_EXPORT int GetTraceNpoints()
SOLVER_UTILS_EXPORT int GetNcoeffs()
void DoOdeProjection(const Array< OneD, const Array< OneD, NekDouble > > &inarray, Array< OneD, Array< OneD, NekDouble > > &outarray, const NekDouble time)
Compute the projection and call the method for imposing the boundary conditions in case of discontinu...
Definition: APE.cpp:263
SOLVER_UTILS_EXPORT void SetInitialConditions(NekDouble initialtime=0.0, bool dumpInitialConditions=true, const int domain=0)
Initialise the data in the dependent fields.
void Zero(int n, T *x, const int incx)
Zero vector.
Definition: Vmath.cpp:359
#define ASSERTL1(condition, msg)
Assert Level 1 – Debugging which is used whether in FULLDEBUG or DEBUG compilation mode...
Definition: ErrorUtil.hpp:191
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.cpp:1038
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
void Vmul(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Multiply vector z = x*y.
Definition: Vmath.cpp:169
tKey RegisterCreatorFunction(tKey idKey, CreatorFunction classCreator, tDescription pDesc="")
Register a class with the factory.
Definition: NekFactory.hpp:215