Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
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 using namespace std;
43 
44 namespace Nektar
45 {
46 string APE::className = GetEquationSystemFactory().RegisterCreatorFunction(
47  "APE", APE::create,
48  "APE1/APE4 (Acoustic Perturbation Equations)");
49 
50 
51 APE::APE(
53  : UnsteadySystem(pSession)
54 {
55 }
56 
57 
58 /**
59  * @brief Initialization object for the APE class.
60  */
62 {
64 
65  // TODO: We have a bug somewhere in the 1D boundary conditions. Therefore 1D
66  // problems are currently disabled. This should get fixed in the future.
67  ASSERTL0(m_spacedim > 1, "1D problems currently not supported by the APE class.");
68 
70  "Only Projection=DisContinuous supported by the APE class.");
71 
72  // Load isentropic coefficient, Ratio of specific heats
73  m_session->LoadParameter("Gamma", m_gamma, 1.4);
74 
75  // Define Baseflow fields
77  m_basefield_names.push_back("p0");
78  m_basefield_names.push_back("rho0");
79  m_basefield_names.push_back("u0");
80  m_basefield_names.push_back("v0");
81  m_basefield_names.push_back("w0");
82 
83  // Resize the advection velocities vector to dimension of the problem
84  m_basefield_names.resize(m_spacedim + 2);
85 
86  // Initialize the sourceterm
88 
89  // Do not forwards transform initial condition
90  m_homoInitialFwd = false;
91 
92  // Define the normal velocity fields
93  if (m_fields[0]->GetTrace())
94  {
96  for (int i = 0; i < m_spacedim + 2; i++)
97  {
99  }
100  }
101 
102  // Set up locations of velocity and base velocity vectors.
105  for (int i = 0; i < m_spacedim; ++i)
106  {
107  // u', v', w'
108  m_vecLocs[0][i] = 1 + i;
109  }
110 
111  string riemName;
112  m_session->LoadSolverInfo("UpwindType", riemName, "APEUpwind");
114  riemName);
115  m_riemannSolver->SetVector("N", &APE::GetNormals, this);
116  m_riemannSolver->SetVector("basefield", &APE::GetBasefield, this);
117  m_riemannSolver->SetAuxVec("vecLocs", &APE::GetVecLocs, this);
118  m_riemannSolver->SetParam("Gamma", &APE::GetGamma, this);
119 
120  // Set up advection operator
121  string advName;
122  m_session->LoadSolverInfo("AdvectionType", advName, "WeakDG");
124  .CreateInstance(advName, advName);
125  m_advection->SetFluxVector(&APE::GetFluxVector, this);
126  m_advection->SetRiemannSolver(m_riemannSolver);
127  m_advection->InitObject(m_session, m_fields);
128 
130  {
133  }
134  else
135  {
136  ASSERTL0(false, "Implicit APE not set up.");
137  }
138 }
139 
140 
141 /**
142  * @brief Destructor for APE class.
143  */
145 {
146 
147 }
148 
149 
150 /**
151  * @brief Return the flux vector for the APE equations.
152  *
153  * @param physfield Fields.
154  * @param flux Resulting flux. flux[eq][dir][pt]
155  */
157  const Array<OneD, Array<OneD, NekDouble> > &physfield,
159 {
160  UpdateBasefield();
161 
162  int nq = physfield[0].num_elements();
163  Array<OneD, NekDouble> tmp1(nq);
164  Array<OneD, NekDouble> tmp2(nq);
165 
166  ASSERTL1(flux[0].num_elements() == m_spacedim,
167  "Dimension of flux array and velocity array do not match");
168 
169  // F_{adv,p',j} = \rho_0 u'_j + p' \bar{u}_j / c^2
170  for (int j = 0; j < m_spacedim; ++j)
171  {
172  Vmath::Zero(nq, flux[0][j], 1);
173 
174  // construct rho_0 u'_j term
175  Vmath::Vmul(nq, m_basefield[1], 1, physfield[j + 1], 1, flux[0][j], 1);
176 
177  // construct p' \bar{u}_j / c^2 term
178  // c^2
179  Vmath::Vdiv(nq, m_basefield[0], 1, m_basefield[1], 1, tmp1, 1);
180  Vmath::Smul(nq, m_gamma, tmp1, 1, tmp1, 1);
181 
182  // p' \bar{u}_j / c^2 term
183  Vmath::Vmul(nq, physfield[0], 1, m_basefield[j + 2], 1, tmp2, 1);
184  Vmath::Vdiv(nq, tmp2, 1, tmp1, 1, tmp2, 1);
185 
186  // \rho_0 u'_j + p' \bar{u}_j / c^2
187  Vmath::Vadd(nq, flux[0][j], 1, tmp2, 1, flux[0][j], 1);
188  }
189 
190  for (int i = 1; i < flux.num_elements(); ++i)
191  {
192  ASSERTL1(flux[i].num_elements() == m_spacedim,
193  "Dimension of flux array and velocity array do not match");
194 
195  // F_{adv,u'_i,j} = (p'/ \bar{rho} + \bar{u}_k u'_k) \delta_{ij}
196  for (int j = 0; j < m_spacedim; ++j)
197  {
198  Vmath::Zero(nq, flux[i][j], 1);
199 
200  if (i - 1 == j)
201  {
202  // contruct p'/ \bar{rho} term
203  Vmath::Vdiv(nq, physfield[0], 1, m_basefield[1], 1, flux[i][j], 1);
204 
205  // construct \bar{u}_k u'_k term
206  Vmath::Zero(nq, tmp1, 1);
207  for (int k = 0; k < m_spacedim; ++k)
208  {
209  Vmath::Vvtvp(nq, physfield[k + 1], 1, m_basefield[k + 2], 1, tmp1, 1, tmp1, 1);
210  }
211 
212  // add terms
213  Vmath::Vadd(nq, flux[i][j], 1, tmp1, 1, flux[i][j], 1);
214  }
215  }
216  }
217 }
218 
219 
220 /**
221  * @brief Compute the right-hand side.
222  */
223 void APE::DoOdeRhs(const Array<OneD, const Array<OneD, NekDouble> >&inarray,
224  Array<OneD, Array<OneD, NekDouble> >&outarray,
225  const NekDouble time)
226 {
227  int nVariables = inarray.num_elements();
228  int nq = GetTotPoints();
229  Array<OneD, NekDouble> tmp1(nq);
230 
231  // WeakDG does not use advVel, so we only provide a dummy array
233  m_advection->Advect(nVariables, m_fields, advVel, inarray, outarray, time);
234 
235  for (int i = 0; i < nVariables; ++i)
236  {
237  if (i == 0)
238  {
239  // c^2 = gamma*p0/rho0
240  Vmath::Vdiv(nq, m_basefield[0], 1, m_basefield[1], 1, tmp1, 1);
241  Vmath::Smul(nq, m_gamma, tmp1, 1, tmp1, 1);
242  Vmath::Vmul(nq, tmp1, 1, outarray[i], 1, outarray[i], 1);
243  }
244 
245  Vmath::Neg(nq, outarray[i], 1);
246  }
247 
248  AddSource(outarray);
249 }
250 
251 
252 /**
253  * @brief Compute the projection and call the method for imposing the
254  * boundary conditions in case of discontinuous projection.
255  */
257  Array<OneD, Array<OneD, NekDouble> >&outarray,
258  const NekDouble time)
259 {
260  int nvariables = inarray.num_elements();
261  int nq = m_fields[0]->GetNpoints();
262 
263  // deep copy
264  for (int i = 0; i < nvariables; ++i)
265  {
266  Vmath::Vcopy(nq, inarray[i], 1, outarray[i], 1);
267  }
268 
269  SetBoundaryConditions(outarray, time);
270 }
271 
272 
273 /**
274  * @brief Apply the Boundary Conditions to the APE equations.
275  */
277  NekDouble time)
278 {
279  std::string varName;
280  int nvariables = m_fields.num_elements();
281  int cnt = 0;
282  int nTracePts = GetTraceTotPoints();
283 
284  // Extract trace for boundaries. Needs to be done on all processors to avoid
285  // deadlock.
286  Array<OneD, Array<OneD, NekDouble> > Fwd(nvariables);
287  for (int i = 0; i < nvariables; ++i)
288  {
289  Fwd[i] = Array<OneD, NekDouble>(nTracePts);
290  m_fields[i]->ExtractTracePhys(inarray[i], Fwd[i]);
291  }
292 
293  // loop over Boundary Regions
294  for(int n = 0; n < m_fields[0]->GetBndConditions().num_elements(); ++n)
295  {
296  // Wall Boundary Condition
297  if (boost::iequals(m_fields[0]->GetBndConditions()[n]->GetUserDefined(),"Wall"))
298  {
299  WallBC(n, cnt, Fwd, inarray);
300  }
301 
302  // Time Dependent Boundary Condition (specified in meshfile)
303  if (m_fields[0]->GetBndConditions()[n]->IsTimeDependent())
304  {
305  for (int i = 0; i < nvariables; ++i)
306  {
307  varName = m_session->GetVariable(i);
308  m_fields[i]->EvaluateBoundaryConditions(time, varName);
309  }
310  }
311  cnt +=m_fields[0]->GetBndCondExpansions()[n]->GetExpSize();
312  }
313 }
314 
315 
316 /**
317  * @brief Wall boundary conditions for the APE equations.
318  */
319 void APE::WallBC(int bcRegion, int cnt,
321  Array<OneD, Array<OneD, NekDouble> > &physarray)
322 {
323  int nVariables = physarray.num_elements();
324 
325  const Array<OneD, const int> &traceBndMap = m_fields[0]->GetTraceBndMap();
326 
327  // Adjust the physical values of the trace to take
328  // user defined boundaries into account
329  int id1, id2, nBCEdgePts;
330  int eMax = m_fields[0]->GetBndCondExpansions()[bcRegion]->GetExpSize();
331 
332  for (int e = 0; e < eMax; ++e)
333  {
334  nBCEdgePts = m_fields[0]->GetBndCondExpansions()[bcRegion]->
335  GetExp(e)->GetTotPoints();
336  id1 = m_fields[0]->GetBndCondExpansions()[bcRegion]->GetPhys_Offset(e);
337  id2 = m_fields[0]->GetTrace()->GetPhys_Offset(traceBndMap[cnt+e]);
338 
339  // For 2D/3D, define: v* = v - 2(v.n)n
340  Array<OneD, NekDouble> tmp(nBCEdgePts, 0.0);
341 
342  // Calculate (v.n)
343  for (int i = 0; i < m_spacedim; ++i)
344  {
345  Vmath::Vvtvp(nBCEdgePts,
346  &Fwd[1+i][id2], 1,
347  &m_traceNormals[i][id2], 1,
348  &tmp[0], 1,
349  &tmp[0], 1);
350  }
351 
352  // Calculate 2.0(v.n)
353  Vmath::Smul(nBCEdgePts, -2.0, &tmp[0], 1, &tmp[0], 1);
354 
355  // Calculate v* = v - 2.0(v.n)n
356  for (int i = 0; i < m_spacedim; ++i)
357  {
358  Vmath::Vvtvp(nBCEdgePts,
359  &tmp[0], 1,
360  &m_traceNormals[i][id2], 1,
361  &Fwd[1+i][id2], 1,
362  &Fwd[1+i][id2], 1);
363  }
364 
365  // Copy boundary adjusted values into the boundary expansion
366  for (int i = 0; i < nVariables; ++i)
367  {
368  Vmath::Vcopy(nBCEdgePts,
369  &Fwd[i][id2], 1,
370  &(m_fields[i]->GetBndCondExpansions()[bcRegion]->UpdatePhys())[id1], 1);
371  }
372  }
373 }
374 
375 
376 /**
377  * @brief sourceterm for p' equation obtained from GetSource
378  */
380 {
382  Vmath::Vadd(GetTotPoints(), m_sourceTerms, 1, outarray[0], 1, outarray[0], 1);
383 }
384 
385 
387  std::vector<Array<OneD, NekDouble> > &fieldcoeffs,
388  std::vector<std::string> &variables)
389 {
390  UpdateBasefield();
391 
392  const int nCoeffs = m_fields[0]->GetNcoeffs();
393 
394  for (int i = 0; i < m_spacedim + 2; i++)
395  {
396  variables.push_back(m_basefield_names[i]);
397 
398  Array<OneD, NekDouble> tmpFwd(nCoeffs);
399  m_fields[0]->FwdTrans(m_basefield[i], tmpFwd);
400  fieldcoeffs.push_back(tmpFwd);
401  }
402 }
403 
404 
405 /**
406  * @brief Get the normal vectors.
407  */
409 {
410  return m_traceNormals;
411 }
412 
413 
414 /**
415  * @brief Get the locations of the components of the directed fields within the fields array.
416  */
418 {
419  return m_vecLocs;
420 }
421 
422 
423 /**
424  * @brief Get the baseflow field.
425  */
427 {
428  for (int i = 0; i < m_spacedim + 2; i++)
429  {
430  m_fields[0]->ExtractTracePhys(m_basefield[i], m_traceBasefield[i]);
431  }
432  return m_traceBasefield;
433 }
434 
435 
436 /**
437  * @brief Get the heat capacity ratio.
438  */
440 {
441  return m_gamma;
442 }
443 
444 
446 {
447  static NekDouble last_update = -1.0;
448 
449  if (m_time > last_update)
450  {
451  last_update = m_time;
453  }
454 }
455 
457 {
458  static NekDouble last_update = -1.0;
459 
460  if (m_time > last_update)
461  {
463 
464  EvaluateFunction("S", m_sourceTerms, "Source", m_time);
465 
466  m_fields[0]->IProductWRTBase(m_sourceTerms, sourceC);
467  m_fields[0]->MultiplyByElmtInvMass(sourceC, sourceC);
468  m_fields[0]->BwdTrans(sourceC, m_sourceTerms);
469 
470  last_update = m_time;
471  }
472 }
473 
474 
475 } //end of namespace
476 
void WallBC(int bcRegion, int cnt, Array< OneD, Array< OneD, NekDouble > > &Fwd, Array< OneD, Array< OneD, NekDouble > > &physarray)
Wall boundary conditions for the APE equations.
Definition: APE.cpp:319
const Array< OneD, const Array< OneD, NekDouble > > & GetNormals()
Get the normal vectors.
Definition: APE.cpp:408
#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:379
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
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.
STL namespace.
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:156
boost::shared_ptr< SessionReader > SessionReaderSharedPtr
Definition: MeshPartition.h:51
SOLVER_UTILS_EXPORT int GetTotPoints()
void UpdateBasefield()
Definition: APE.cpp:445
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:439
virtual void v_ExtraFldOutput(std::vector< Array< OneD, NekDouble > > &fieldcoeffs, std::vector< std::string > &variables)
Definition: APE.cpp:386
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:417
Array< OneD, Array< OneD, NekDouble > > m_traceBasefield
Definition: APE.h:74
void DefineOdeRhs(FuncPointerT func, ObjectPointerT obj)
virtual ~APE()
Destructor.
Definition: APE.cpp:144
Base class for unsteady solvers.
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:223
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
double NekDouble
virtual void v_InitObject()
Initialization object for the APE class.
Definition: APE.cpp:61
virtual SOLVER_UTILS_EXPORT void v_InitObject()
Init object for UnsteadySystem class.
void SetBoundaryConditions(Array< OneD, Array< OneD, NekDouble > > &physarray, NekDouble time)
Apply the Boundary Conditions to the APE equations.
Definition: APE.cpp:276
void UpdateSourceTerms()
Definition: APE.cpp:456
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:426
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:256
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:1047
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