Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties 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  */
46  "Mapping", MappingExtrapolate::create, "Mapping");
47 
52  const Array<OneD, int> pVel,
53  const SolverUtils::AdvectionSharedPtr advObject)
54  : StandardExtrapolate(pSession, pFields, pPressure, pVel, advObject)
55 {
57 
58  // Load solve parameters related to the mapping
59  // Flags determining if pressure/viscous terms should be treated implicitly
60  m_session->MatchSolverInfo(
61  "MappingImplicitPressure", "True", m_implicitPressure, false);
62  m_session->MatchSolverInfo(
63  "MappingImplicitViscous", "True", m_implicitViscous, false);
64 
65  // Relaxation parameter for pressure system
66  m_session->LoadParameter(
67  "MappingPressureRelaxation", m_pressureRelaxation, 1.0);
68 }
69 
71 {
72 }
73 
74 /**
75  *
76  */
78  const Array<OneD, NekDouble> &pressure)
79 {
80  if (m_HBCnumber > 0)
81  {
82  int cnt, n;
83  int physTot = m_fields[0]->GetTotPoints();
84  int nvel = m_fields.num_elements() - 1;
85 
87  // Remove previous correction
88  for (cnt = n = 0; n < m_PBndConds.num_elements(); ++n)
89  {
90  if (m_PBndConds[n]->GetUserDefined() == "H")
91  {
92  int nq = m_PBndExp[n]->GetNcoeffs();
93  Vmath::Vsub(nq,
94  &(m_PBndExp[n]->GetCoeffs()[0]),
95  1,
96  &(m_bcCorrection[cnt]),
97  1,
98  &(m_PBndExp[n]->UpdateCoeffs()[0]),
99  1);
100  cnt += nq;
101  }
102  }
103 
104  // Calculate new correction
105  Array<OneD, NekDouble> Jac(physTot, 0.0);
106  m_mapping->GetJacobian(Jac);
107 
108  Array<OneD, Array<OneD, NekDouble> > correction(nvel);
112  for (int i = 0; i < nvel; i++)
113  {
114  wk[i] = Array<OneD, NekDouble>(physTot, 0.0);
115  gradP[i] = Array<OneD, NekDouble>(physTot, 0.0);
116  correction[i] = Array<OneD, NekDouble>(physTot, 0.0);
117  }
118 
119  // Calculate G(p)
120  for (int i = 0; i < nvel; ++i)
121  {
122  m_fields[0]->PhysDeriv(
123  MultiRegions::DirCartesianMap[i], pressure, gradP[i]);
124  if (m_fields[0]->GetWaveSpace())
125  {
126  m_fields[0]->HomogeneousBwdTrans(gradP[i], wk[i]);
127  }
128  else
129  {
130  Vmath::Vcopy(physTot, gradP[i], 1, wk[i], 1);
131  }
132  }
133  m_mapping->RaiseIndex(wk, correction); // G(p)
134 
135  // alpha*J*(G(p))
136  if (!m_mapping->HasConstantJacobian())
137  {
138  for (int i = 0; i < nvel; ++i)
139  {
140  Vmath::Vmul(
141  physTot, correction[i], 1, Jac, 1, correction[i], 1);
142  }
143  }
144  for (int i = 0; i < nvel; ++i)
145  {
146  Vmath::Smul(physTot,
148  correction[i],
149  1,
150  correction[i],
151  1);
152  }
153 
154  if (m_pressure->GetWaveSpace())
155  {
156  for (int i = 0; i < nvel; ++i)
157  {
158  m_pressure->HomogeneousFwdTrans(correction[i], correction[i]);
159  }
160  }
161  // p_i - alpha*J*div(G(p))
162  for (int i = 0; i < nvel; ++i)
163  {
164  Vmath::Vsub(
165  physTot, gradP[i], 1, correction[i], 1, correction[i], 1);
166  }
167 
168  // Get value at boundary and calculate Inner product
172  for (n = cnt = 0; n < m_PBndConds.num_elements(); ++n)
173  {
174  // High order boundary condition;
175  if (boost::iequals(m_PBndConds[n]->GetUserDefined(), "H"))
176  {
177  m_fields[0]->GetBndElmtExpansion(n, BndElmtExp);
178 
179  // Obtaining fields on BndElmtExp
180  for (int i = 0; i < m_bnd_dim; i++)
181  {
182  m_fields[0]->ExtractPhysToBndElmt(
183  n, correction[i], correctionElmt[i]);
184  }
185 
186  Vals = m_bcCorrection + cnt;
187  // Getting values on the edge and filling the correction
188  for (int i = 0; i < m_bnd_dim; i++)
189  {
190  m_fields[0]->ExtractElmtToBndPhys(
191  n, correctionElmt[i], BndValues[i]);
192  }
193  m_PBndExp[n]->NormVectorIProductWRTBase(BndValues, Vals);
194 
195  // Get offset for next terms
196  cnt += m_PBndExp[n]->GetNcoeffs();
197  }
198  }
199 
200  // Apply new correction
201  for (cnt = n = 0; n < m_PBndConds.num_elements(); ++n)
202  {
203  if (m_PBndConds[n]->GetUserDefined() == "H")
204  {
205  int nq = m_PBndExp[n]->GetNcoeffs();
206  Vmath::Vadd(nq,
207  &(m_PBndExp[n]->GetCoeffs()[0]),
208  1,
209  &(m_bcCorrection[cnt]),
210  1,
211  &(m_PBndExp[n]->UpdateCoeffs()[0]),
212  1);
213  cnt += nq;
214  }
215  }
216  }
217 }
218 
220  const Array<OneD, const Array<OneD, NekDouble> > &fields,
221  const Array<OneD, const Array<OneD, NekDouble> > &N,
222  NekDouble kinvis)
223 {
224  if (m_mapping->HasConstantJacobian() && !m_implicitViscous)
225  {
226  Extrapolate::v_CalcNeumannPressureBCs(fields, N, kinvis);
227  }
228  else
229  {
230  int physTot = m_fields[0]->GetTotPoints();
231  int nvel = m_fields.num_elements() - 1;
232  int i, n, cnt;
233 
236 
239  // Get transformation Jacobian
240  Array<OneD, NekDouble> Jac(physTot, 0.0);
241  m_mapping->GetJacobian(Jac);
242  // Declare variables
245  Array<OneD, Array<OneD, NekDouble> > Q_field(nvel);
246  Array<OneD, Array<OneD, NekDouble> > fields_new(nvel);
248  // Temporary variables
249  Array<OneD, NekDouble> tmp(physTot, 0.0);
250  Array<OneD, NekDouble> tmp2(physTot, 0.0);
251  for (i = 0; i < m_bnd_dim; i++)
252  {
253  N_new[i] = Array<OneD, NekDouble>(physTot, 0.0);
254  }
255  for (int i = 0; i < nvel; i++)
256  {
257  Q_field[i] = Array<OneD, NekDouble>(physTot, 0.0);
258  fields_new[i] = Array<OneD, NekDouble>(physTot, 0.0);
259  }
260 
261  // Multiply convective terms by Jacobian
262  for (i = 0; i < m_bnd_dim; i++)
263  {
264  if (m_fields[0]->GetWaveSpace())
265  {
266  m_fields[0]->HomogeneousBwdTrans(N[i], N_new[i]);
267  }
268  else
269  {
270  Vmath::Vcopy(physTot, N[i], 1, N_new[i], 1);
271  }
272  Vmath::Vmul(physTot, Jac, 1, N_new[i], 1, N_new[i], 1);
273  if (m_fields[0]->GetWaveSpace())
274  {
275  m_fields[0]->HomogeneousFwdTrans(N_new[i], N_new[i]);
276  }
277  }
278 
279  // Get velocity in physical space
280  for (i = 0; i < nvel; i++)
281  {
282  if (m_fields[0]->GetWaveSpace())
283  {
284  m_fields[0]->HomogeneousBwdTrans(fields[i], fields_new[i]);
285  }
286  else
287  {
288  Vmath::Vcopy(physTot, fields[i], 1, fields_new[i], 1);
289  }
290  }
291 
292  // Calculate appropriate form of the CurlCurl operator
293  m_mapping->CurlCurlField(fields_new, Q_field, m_implicitViscous);
294 
295  // If viscous terms are treated explicitly,
296  // add grad(U/J \dot grad J) to CurlCurl
297  if (!m_implicitViscous)
298  {
299  m_mapping->DotGradJacobian(fields_new, tmp);
300  Vmath::Vdiv(physTot, tmp, 1, Jac, 1, tmp, 1);
301 
302  bool wavespace = m_fields[0]->GetWaveSpace();
303  m_fields[0]->SetWaveSpace(false);
304  for (int i = 0; i < m_bnd_dim; i++)
305  {
306  m_fields[0]->PhysDeriv(
307  MultiRegions::DirCartesianMap[i], tmp, tmp2);
308  Vmath::Vadd(physTot, Q_field[i], 1, tmp2, 1, Q_field[i], 1);
309  }
310  m_fields[0]->SetWaveSpace(wavespace);
311  }
312 
313  // Multiply by Jacobian and convert to wavespace (if necessary)
314  for (i = 0; i < m_bnd_dim; i++)
315  {
316  Vmath::Vmul(physTot, Jac, 1, fields_new[i], 1, fields_new[i], 1);
317  Vmath::Vmul(physTot, Jac, 1, Q_field[i], 1, Q_field[i], 1);
318  if (m_fields[0]->GetWaveSpace())
319  {
320  m_fields[0]->HomogeneousFwdTrans(fields_new[i], fields_new[i]);
321  m_fields[0]->HomogeneousFwdTrans(Q_field[i], Q_field[i]);
322  }
323  }
324 
326  for (n = cnt = 0; n < m_PBndConds.num_elements(); ++n)
327  {
328  // High order boundary condition;
329  if (boost::iequals(m_PBndConds[n]->GetUserDefined(), "H"))
330  {
331  m_fields[0]->GetBndElmtExpansion(n, BndElmtExp);
332  int nq = BndElmtExp->GetTotPoints();
333 
334  // Obtaining fields on BndElmtExp
335  for (int i = 0; i < m_bnd_dim; i++)
336  {
337  m_fields[0]->ExtractPhysToBndElmt(
338  n, fields_new[i], Velocity[i]);
339  m_fields[0]->ExtractPhysToBndElmt(
340  n, N_new[i], Advection[i]);
341  m_fields[0]->ExtractPhysToBndElmt(n, Q_field[i], Q[i]);
342  }
343 
344  // Mounting advection component into the high-order condition
345  for (int i = 0; i < m_bnd_dim; i++)
346  {
347  MountHOPBCs(nq, kinvis, Q[i], Advection[i]);
348  }
349 
350  Pvals = (m_pressureHBCs[m_intSteps - 1]) + cnt;
351  Uvals = (m_iprodnormvel[m_intSteps]) + cnt;
352 
353  // Getting values on the edge and filling the pressure boundary
354  // expansion and the acceleration term. Multiplication by the
355  // normal is required
356  for (int i = 0; i < m_bnd_dim; i++)
357  {
358  m_fields[0]->ExtractElmtToBndPhys(n, Q[i], BndValues[i]);
359  }
360  m_PBndExp[n]->NormVectorIProductWRTBase(BndValues, Pvals);
361 
362  for (int i = 0; i < m_bnd_dim; i++)
363  {
364  m_fields[0]->ExtractElmtToBndPhys(
365  n, Velocity[i], BndValues[i]);
366  }
367  m_PBndExp[n]->NormVectorIProductWRTBase(BndValues, Uvals);
368 
369  // Get offset for next terms
370  cnt += m_PBndExp[n]->GetNcoeffs();
371  }
372  }
373  }
374  // If pressure terms are treated implicitly, we need to multiply
375  // by the relaxation parameter, and zero the correction term
376  if (m_implicitPressure)
377  {
381  1,
383  1);
384  }
386 }
387 }
MultiRegions::ExpListSharedPtr m_pressure
Pointer to field holding pressure field.
Definition: Extrapolate.h:222
LibUtilities::SessionReaderSharedPtr m_session
Definition: Extrapolate.h:211
MappingExtrapolate(const LibUtilities::SessionReaderSharedPtr pSession, Array< OneD, MultiRegions::ExpListSharedPtr > pFields, MultiRegions::ExpListSharedPtr pPressure, const Array< OneD, int > pVel, const SolverUtils::AdvectionSharedPtr advObject)
Array< OneD, MultiRegions::ExpListSharedPtr > m_PBndExp
pressure boundary conditions expansion container
Definition: Extrapolate.h:242
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.
boost::shared_ptr< Advection > AdvectionSharedPtr
A shared pointer to an Advection object.
Definition: Advection.h:165
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:241
boost::shared_ptr< SessionReader > SessionReaderSharedPtr
Definition: MeshPartition.h:51
Array< OneD, MultiRegions::ExpListSharedPtr > m_fields
Velocity fields.
Definition: Extrapolate.h:219
Array< OneD, NekDouble > m_bcCorrection
GlobalMapping::MappingSharedPtr m_mapping
int m_bnd_dim
bounday dimensionality
Definition: Extrapolate.h:236
void MountHOPBCs(int HBCdata, NekDouble kinvis, Array< OneD, NekDouble > &Q, Array< OneD, const NekDouble > &Advection)
Definition: Extrapolate.h:392
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:213
int m_intSteps
Maximum points used in pressure BC evaluation.
Definition: Extrapolate.h:254
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:343
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:259
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.
Array< OneD, Array< OneD, NekDouble > > m_iprodnormvel
Storage for current and previous levels of the inner product of normal velocity.
Definition: Extrapolate.h:262
Array< OneD, const SpatialDomains::BoundaryConditionShPtr > m_PBndConds
pressure boundary conditions container
Definition: Extrapolate.h:239
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.cpp:1061
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:265
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
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:183
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