Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
MappingExtrapolate.cpp
Go to the documentation of this file.
1 ///////////////////////////////////////////////////////////////////////////////
2 //
3 // File: StandardExtrapolate.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: Abstract base class for StandardExtrapolate.
33 //
34 ///////////////////////////////////////////////////////////////////////////////
35 
38 
39 namespace Nektar
40 {
41  /**
42  * Registers the class with the Factory.
43  */
45  "Mapping",
47  "Mapping");
48 
53  const Array<OneD, int> pVel,
54  const SolverUtils::AdvectionSharedPtr advObject)
55  : StandardExtrapolate(pSession,pFields,pPressure,pVel,advObject)
56  {
58 
59  // Load solve parameters related to the mapping
60  // Flags determining if pressure/viscous terms should be treated implicitly
61  m_session->MatchSolverInfo("MappingImplicitPressure","True",m_implicitPressure,false);
62  m_session->MatchSolverInfo("MappingImplicitViscous","True",m_implicitViscous,false);
63 
64  // Relaxation parameter for pressure system
65  m_session->LoadParameter("MappingPressureRelaxation",m_pressureRelaxation,1.0);
66  }
67 
69  {
70  }
71 
72  /**
73  *
74  */
76  const Array<OneD, NekDouble> &pressure)
77  {
78  if(m_HBCdata.num_elements()>0)
79  {
80  int cnt, n;
81  int physTot = m_fields[0]->GetTotPoints();
82  int nvel = m_fields.num_elements()-1;
83 
85  // Remove previous correction
86  for(cnt = n = 0; n < m_PBndConds.num_elements(); ++n)
87  {
88  if(m_PBndConds[n]->GetUserDefined() == "H")
89  {
90  int nq = m_PBndExp[n]->GetNcoeffs();
91  Vmath::Vsub(nq, &(m_PBndExp[n]->GetCoeffs()[0]), 1,
92  &(m_bcCorrection[cnt]), 1,
93  &(m_PBndExp[n]->UpdateCoeffs()[0]), 1);
94  cnt += nq;
95  }
96  }
97 
98  // Calculate new correction
99  Array<OneD, NekDouble> Jac(physTot, 0.0);
100  m_mapping->GetJacobian(Jac);
101 
102  Array<OneD, Array<OneD, NekDouble> > correction(nvel);
106  for (int i=0; i<nvel; i++)
107  {
108  wk[i] = Array<OneD, NekDouble> (physTot, 0.0);
109  gradP[i] = Array<OneD, NekDouble> (physTot, 0.0);
110  correction[i] = Array<OneD, NekDouble> (physTot, 0.0);
111  }
112 
113  // Calculate G(p)
114  for(int i = 0; i < nvel; ++i)
115  {
116  m_fields[0]->PhysDeriv(MultiRegions::DirCartesianMap[i], pressure, gradP[i]);
117  if(m_fields[0]->GetWaveSpace())
118  {
119  m_fields[0]->HomogeneousBwdTrans(gradP[i], wk[i]);
120  }
121  else
122  {
123  Vmath::Vcopy(physTot, gradP[i], 1, wk[i], 1);
124  }
125  }
126  m_mapping->RaiseIndex(wk, correction); // G(p)
127 
128  // alpha*J*(G(p))
129  if (!m_mapping->HasConstantJacobian())
130  {
131  for(int i = 0; i < nvel; ++i)
132  {
133  Vmath::Vmul(physTot, correction[i], 1, Jac, 1, correction[i], 1);
134  }
135  }
136  for(int i = 0; i < nvel; ++i)
137  {
138  Vmath::Smul(physTot, m_pressureRelaxation, correction[i], 1, correction[i], 1);
139  }
140 
141  if(m_pressure->GetWaveSpace())
142  {
143  for(int i = 0; i < nvel; ++i)
144  {
145  m_pressure->HomogeneousFwdTrans(correction[i], correction[i]);
146  }
147  }
148  // p_i - alpha*J*div(G(p))
149  for (int i = 0; i < nvel; ++i)
150  {
151  Vmath::Vsub(physTot, gradP[i], 1, correction[i], 1, correction[i], 1);
152  }
153 
154  // Get value at boundary and calculate Inner product
159  for(int i = 0; i < m_bnd_dim; i++)
160  {
161  BndValues[i] = Array<OneD, NekDouble> (m_pressureBCsMaxPts,0.0);
162  }
163  for(int j = 0 ; j < m_HBCdata.num_elements() ; j++)
164  {
165  /// Casting the boundary expansion to the specific case
166  Pbc = boost::dynamic_pointer_cast<StdRegions::StdExpansion>
167  (m_PBndExp[m_HBCdata[j].m_bndryID]
168  ->GetExp(m_HBCdata[j].m_bndElmtID));
169 
170  /// Picking up the element where the HOPBc is located
171  elmt = m_pressure->GetExp(m_HBCdata[j].m_globalElmtID);
172 
173  /// Assigning
174  for(int i = 0; i < m_bnd_dim; i++)
175  {
176  correctionElmt[i] = correction[i] + m_HBCdata[j].m_physOffset;
177  }
178  Vals = m_bcCorrection + m_HBCdata[j].m_coeffOffset;
179  // Getting values on the edge and filling the correction
180  switch(m_pressure->GetExpType())
181  {
182  case MultiRegions::e2D:
184  {
185  elmt->GetEdgePhysVals(m_HBCdata[j].m_elmtTraceID, Pbc,
186  correctionElmt[0], BndValues[0]);
187  elmt->GetEdgePhysVals(m_HBCdata[j].m_elmtTraceID, Pbc,
188  correctionElmt[1], BndValues[1]);
189 
190  // InnerProduct
191  Pbc->NormVectorIProductWRTBase(BndValues[0], BndValues[1],
192  Vals);
193  }
194  break;
195 
196  case MultiRegions::e3D:
197  {
198  elmt->GetFacePhysVals(m_HBCdata[j].m_elmtTraceID, Pbc,
199  correctionElmt[0], BndValues[0]);
200  elmt->GetFacePhysVals(m_HBCdata[j].m_elmtTraceID, Pbc,
201  correctionElmt[1], BndValues[1]);
202  elmt->GetFacePhysVals(m_HBCdata[j].m_elmtTraceID, Pbc,
203  correctionElmt[2], BndValues[2]);
204  Pbc->NormVectorIProductWRTBase(BndValues[0], BndValues[1],
205  BndValues[2], Vals);
206  }
207  break;
208  default:
209  ASSERTL0(0,"Dimension not supported");
210  break;
211  }
212  }
213 
214  // Apply new correction
215  for(cnt = n = 0; n < m_PBndConds.num_elements(); ++n)
216  {
217  if(m_PBndConds[n]->GetUserDefined() == "H")
218  {
219  int nq = m_PBndExp[n]->GetNcoeffs();
220  Vmath::Vadd(nq, &(m_PBndExp[n]->GetCoeffs()[0]), 1,
221  &(m_bcCorrection[cnt]), 1,
222  &(m_PBndExp[n]->UpdateCoeffs()[0]), 1);
223  cnt += nq;
224  }
225  }
226  }
227  }
228 
230  const Array<OneD, const Array<OneD, NekDouble> > &fields,
231  const Array<OneD, const Array<OneD, NekDouble> > &N,
232  NekDouble kinvis)
233  {
234  if (m_mapping->HasConstantJacobian() && !m_implicitViscous)
235  {
236  Extrapolate::v_CalcNeumannPressureBCs( fields, N, kinvis);
237  }
238  else
239  {
240  int physTot = m_fields[0]->GetTotPoints();
241  int nvel = m_fields.num_elements()-1;
242 
246 
249  // Get transformation Jacobian
250  Array<OneD, NekDouble> Jac(physTot,0.0);
251  m_mapping->GetJacobian(Jac);
252  // Declare variables
255  Array<OneD, Array<OneD, NekDouble> > Q_field(nvel);
256  Array<OneD, Array<OneD, NekDouble> > fields_new(nvel);
258  // Temporary variables
259  Array<OneD, NekDouble> tmp(physTot,0.0);
260  Array<OneD, NekDouble> tmp2(physTot,0.0);
261  for(int i = 0; i < m_bnd_dim; i++)
262  {
263  BndValues[i] = Array<OneD, NekDouble> (m_pressureBCsMaxPts,0.0);
265  N_new[i] = Array<OneD, NekDouble> (physTot,0.0);
266  }
267  for(int i = 0; i < nvel; i++)
268  {
269  Q_field[i] = Array<OneD, NekDouble> (physTot,0.0);
270  fields_new[i] = Array<OneD, NekDouble> (physTot,0.0);
271  }
272 
273  // Multiply convective terms by Jacobian
274  for(int i = 0; i < m_bnd_dim; i++)
275  {
276  if (m_fields[0]->GetWaveSpace())
277  {
278  m_fields[0]->HomogeneousBwdTrans(N[i],N_new[i]);
279  }
280  else
281  {
282  Vmath::Vcopy(physTot, N[i], 1, N_new[i], 1);
283  }
284  Vmath::Vmul(physTot, Jac, 1, N_new[i], 1, N_new[i], 1);
285  if (m_fields[0]->GetWaveSpace())
286  {
287  m_fields[0]->HomogeneousFwdTrans(N_new[i],N_new[i]);
288  }
289  }
290 
291  // Get velocity in physical space
292  for(int i = 0; i < nvel; i++)
293  {
294  if (m_fields[0]->GetWaveSpace())
295  {
296  m_fields[0]->HomogeneousBwdTrans(fields[i],fields_new[i]);
297  }
298  else
299  {
300  Vmath::Vcopy(physTot, fields[i], 1, fields_new[i], 1);
301  }
302  }
303 
304  // Calculate appropriate form of the CurlCurl operator
305  m_mapping->CurlCurlField(fields_new, Q_field, m_implicitViscous);
306 
307  // If viscous terms are treated explicitly,
308  // add grad(U/J \dot grad J) to CurlCurl
309  if ( !m_implicitViscous)
310  {
311  m_mapping->DotGradJacobian(fields_new, tmp);
312  Vmath::Vdiv(physTot, tmp, 1, Jac, 1, tmp, 1);
313 
314  bool wavespace = m_fields[0]->GetWaveSpace();
315  m_fields[0]->SetWaveSpace(false);
316  for(int i = 0; i < m_bnd_dim; i++)
317  {
318  m_fields[0]->PhysDeriv(MultiRegions::DirCartesianMap[i],
319  tmp, tmp2);
320  Vmath::Vadd(physTot, Q_field[i], 1, tmp2, 1, Q_field[i], 1);
321  }
322  m_fields[0]->SetWaveSpace(wavespace);
323  }
324 
325  // Multiply by Jacobian and convert to wavespace (if necessary)
326  for(int i = 0; i < m_bnd_dim; i++)
327  {
328  Vmath::Vmul(physTot, Jac, 1, fields_new[i], 1, fields_new[i], 1);
329  Vmath::Vmul(physTot, Jac, 1, Q_field[i] , 1, Q_field[i] , 1);
330  if (m_fields[0]->GetWaveSpace())
331  {
332  m_fields[0]->HomogeneousFwdTrans(fields_new[i],fields_new[i]);
333  m_fields[0]->HomogeneousFwdTrans(Q_field[i],Q_field[i]);
334  }
335  }
336 
337  for(int j = 0 ; j < m_HBCdata.num_elements() ; j++)
338  {
339  /// Casting the boundary expansion to the specific case
340  Pbc = boost::dynamic_pointer_cast<StdRegions::StdExpansion>
341  (m_PBndExp[m_HBCdata[j].m_bndryID]
342  ->GetExp(m_HBCdata[j].m_bndElmtID));
343 
344  /// Picking up the element where the HOPBc is located
345  elmt = m_pressure->GetExp(m_HBCdata[j].m_globalElmtID);
346 
347  /// Assigning
348  for(int i = 0; i < m_bnd_dim; i++)
349  {
350  Velocity[i] = fields_new[i] + m_HBCdata[j].m_physOffset;
351  Advection[i] = N_new[i] + m_HBCdata[j].m_physOffset;
352  Q[i] = Q_field[i] + m_HBCdata[j].m_physOffset;
353  }
354 
355  // Mounting advection component into the high-order condition
356  for(int i = 0; i < m_bnd_dim; i++)
357  {
358  MountHOPBCs(m_HBCdata[j].m_ptsInElmt,kinvis,Q[i],Advection[i]);
359  }
360 
361  Pvals = m_pressureHBCs[0] + m_HBCdata[j].m_coeffOffset;
362 
363  // Getting values on the edge and filling the pressure boundary
364  // expansion and the acceleration term. Multiplication by the
365  // normal is required
366  switch(m_pressure->GetExpType())
367  {
368  case MultiRegions::e2D:
370  {
371  elmt->GetEdgePhysVals(m_HBCdata[j].m_elmtTraceID, Pbc,
372  Q[0], BndValues[0]);
373  elmt->GetEdgePhysVals(m_HBCdata[j].m_elmtTraceID, Pbc,
374  Q[1], BndValues[1]);
375 
376  // InnerProduct
377  Pbc->NormVectorIProductWRTBase(BndValues[0], BndValues[1],
378  Pvals);
379  }
380  break;
382  {
383  if(m_HBCdata[j].m_elmtTraceID == 0)
384  {
385  (m_PBndExp[m_HBCdata[j].m_bndryID]->UpdateCoeffs()
386  + m_PBndExp[m_HBCdata[j].m_bndryID]
387  ->GetCoeff_Offset(
388  m_HBCdata[j].m_bndElmtID))[0]
389  = -1.0*Q[0][0];
390  }
391  else if (m_HBCdata[j].m_elmtTraceID == 1)
392  {
393  (m_PBndExp[m_HBCdata[j].m_bndryID]->UpdateCoeffs()
394  + m_PBndExp[m_HBCdata[j].m_bndryID]
395  ->GetCoeff_Offset(
396  m_HBCdata[j].m_bndElmtID))[0]
397  = Q[0][m_HBCdata[j].m_ptsInElmt-1];
398  }
399  else
400  {
401  ASSERTL0(false,
402  "In the 3D homogeneous 2D approach BCs edge "
403  "ID can be just 0 or 1 ");
404  }
405  }
406  break;
407  case MultiRegions::e3D:
408  {
409  elmt->GetFacePhysVals(m_HBCdata[j].m_elmtTraceID, Pbc,
410  Q[0], BndValues[0]);
411  elmt->GetFacePhysVals(m_HBCdata[j].m_elmtTraceID, Pbc,
412  Q[1], BndValues[1]);
413  elmt->GetFacePhysVals(m_HBCdata[j].m_elmtTraceID, Pbc,
414  Q[2], BndValues[2]);
415  Pbc->NormVectorIProductWRTBase(BndValues[0], BndValues[1],
416  BndValues[2], Pvals);
417  }
418  break;
419  default:
420  ASSERTL0(0,"Dimension not supported");
421  break;
422  }
423  }
424  }
425  // If pressure terms are treated implicitly, we need to multiply
426  // by the relaxation parameter, and zero the correction term
427  if (m_implicitPressure)
428  {
430  m_pressureHBCs[0], 1,
431  m_pressureHBCs[0], 1);
432  }
433  m_bcCorrection = Array<OneD, NekDouble> (m_pressureHBCs[0].num_elements(), 0.0);
434  }
435 
436 }
437 
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:188
MultiRegions::ExpListSharedPtr m_pressure
Pointer to field holding pressure field.
Definition: Extrapolate.h:215
LibUtilities::SessionReaderSharedPtr m_session
Definition: Extrapolate.h:208
MappingExtrapolate(const LibUtilities::SessionReaderSharedPtr pSession, Array< OneD, MultiRegions::ExpListSharedPtr > pFields, MultiRegions::ExpListSharedPtr pPressure, const Array< OneD, int > pVel, const SolverUtils::AdvectionSharedPtr advObject)
int m_pressureBCsMaxPts
Maximum points used in pressure BC evaluation.
Definition: Extrapolate.h:241
Array< OneD, MultiRegions::ExpListSharedPtr > m_PBndExp
pressure boundary conditions expansion container
Definition: Extrapolate.h:235
ExtrapolateFactory & GetExtrapolateFactory()
Definition: Extrapolate.cpp:50
static ExtrapolateSharedPtr create(const LibUtilities::SessionReaderSharedPtr &pSession, Array< OneD, MultiRegions::ExpListSharedPtr > &pFields, MultiRegions::ExpListSharedPtr &pPressure, const Array< OneD, int > &pVel, const SolverUtils::AdvectionSharedPtr &advObject)
Creates an instance of this class.
int m_pressureBCsElmtMaxPts
Maximum points used in Element adjacent to pressure BC evaluation.
Definition: Extrapolate.h:244
boost::shared_ptr< Advection > AdvectionSharedPtr
A shared pointer to an Advection object.
Definition: Advection.h:158
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
boost::shared_ptr< SessionReader > SessionReaderSharedPtr
Definition: MeshPartition.h:51
Array< OneD, MultiRegions::ExpListSharedPtr > m_fields
Definition: Extrapolate.h:212
Array< OneD, NekDouble > m_bcCorrection
GlobalMapping::MappingSharedPtr m_mapping
int m_bnd_dim
bounday dimensionality
Definition: Extrapolate.h:229
void MountHOPBCs(int HBCdata, NekDouble kinvis, Array< OneD, NekDouble > &Q, Array< OneD, const NekDouble > &Advection)
Definition: Extrapolate.h:376
The base class for all shapes.
Definition: StdExpansion.h:69
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
boost::shared_ptr< ExpList > ExpListSharedPtr
Shared pointer to an ExpList object.
double NekDouble
virtual void v_CalcNeumannPressureBCs(const Array< OneD, const Array< OneD, NekDouble > > &fields, const Array< OneD, const Array< OneD, NekDouble > > &N, NekDouble kinvis)
void Vsub(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Subtract vector z = x-y.
Definition: Vmath.cpp:329
MultiRegions::Direction const DirCartesianMap[]
Definition: ExpList.h:86
Array< OneD, Array< OneD, NekDouble > > m_pressureHBCs
Storage for current and previous levels of high order pressure boundary conditions.
Definition: Extrapolate.h:271
Array< OneD, HBCInfo > m_HBCdata
data structure to old all the information regarding High order pressure BCs
Definition: Extrapolate.h:277
virtual void v_CalcNeumannPressureBCs(const Array< OneD, const Array< OneD, NekDouble > > &fields, const Array< OneD, const Array< OneD, NekDouble > > &N, NekDouble kinvis)
static std::string className
Name of class.
boost::shared_ptr< StdExpansion > StdExpansionSharedPtr
Array< OneD, const SpatialDomains::BoundaryConditionShPtr > m_PBndConds
pressure boundary conditions container
Definition: Extrapolate.h:232
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.cpp:1047
virtual void v_CorrectPressureBCs(const Array< OneD, NekDouble > &pressure)
static GLOBAL_MAPPING_EXPORT MappingSharedPtr Load(const LibUtilities::SessionReaderSharedPtr &pSession, const Array< OneD, MultiRegions::ExpListSharedPtr > &pFields)
Return a pointer to the mapping, creating it on first call.
Definition: Mapping.cpp:266
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
Defines a callback function which evaluates the flux vector.
Definition: Advection.h:69
tKey RegisterCreatorFunction(tKey idKey, CreatorFunction classCreator, tDescription pDesc="")
Register a class with the factory.
Definition: NekFactory.hpp:215